<figure><img src="image_1.png" /></figure>

HAL
open science

# Kalman-based approaches for online estimation of bioreactor dynamics from fluorescent reporter measurements

Rand Asswad, Eugenio Cinquemani, Jean-Luc Gouzé

► To cite this version:

Rand Asswad, Eugenio Cinquemani, Jean-Luc Gouzé. Kalman-based approaches for online estimation of bioreactor dynamics from fluorescent reporter measurements. ECC 2024 - European Control Conference, Jun 2024, Stockholm, Sweden. pp.1-6. hal-04556535

HAL Id: hal-04556535

https://hal.science/hal-04556535

Submitted on 23 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Kalman-based approaches for online estimation of bioreactor dynamics from fluorescent reporter measurements*

Rand Asswad<fnref n="1" />,†, Eugenio Cinquemani<fnref n="1" />, Jean-Luc Gouzé<fnref n="2" />

**Abstract**— We address online estimation of microbial growth dynamics in bioreactors from measurements of a fluorescent reporter protein synthesized along with microbial growth. We consider an extended version of standard growth models that accounts for the dynamics of reporter synthesis. We develop state estimation from sampled, noisy measurements in the cases of known and unknown growth rate functions. Leveraging conservation laws and regularized estimation techniques, we reduce these nonlinear estimation problems to linear time-varying ones, and solve them via Kalman filtering. We establish convergence results in absence of noise and show performance on noisy data in simulation.

## I. INTRODUCTION

Estimation and control of cellular growth in bioreactors have been dedicated a significant amount of work over several decades [1]. Nowadays, these challenges receive renewed attention in relation with new biotechnological developments. Among the new research frontiers is the feedback control of microbial consortia [2]. Motivated by the potential to outperform single species in biosynthetic processes and other applications [3], control of microbial consortia requires in the first place to be able to discern single species dynamics by suitable real-time monitoring, e.g. the use of fluorescent reporter proteins expressed along with the growth of the different species.

In a recent paper [4], we have explored (local) observability of the prototypical microbial consortium model from [5] under different assumptions about the observed variables. The construction and analysis of observers for (noisy) time-sampled measurements (as encountered frequently, e.g. with the use of pipetting robots [6]) have not been explored in detail. Given the nonlinear nature of these models, the problem is nontrivial. General observer synthesis approaches exist [7], [8], however their viability depends on perfect knowledge of the system, which hardly applies to biological systems, and/or on idealized assumptions about the measurement process. Robust estimation approaches can tackle in part modelling uncertainty [9], but they are generally not conceived for sampled data. Nonlinear generalizations of Kalman filtering exist that cope explicitly with sampled noisy data [10], [11], but performance is usually not guaranteed even with a perfectly known model.

In this paper, we focus on the construction of state estimators for a single species growing in a bioreactor from indirect,

sampled fluorescent reporter measurements of the species abundance. We consider classical nonlinear growth models [1], [12]) modified to account for the dynamical relationship between fluorescence abundance and biomass synthesis. As explained above, besides its interest per se, the problem constitutes a key step toward estimation of dynamics of microbial communities. Much literature has been dedicated to the problem in the case of direct observations of biomass abundance [1], [13] or gaseous outflow [14]. None of these approaches applies to the augmented model we consider with sampled, noisy measurements.

We consider two cases, perfect knowledge of the functional form of the growth law, and lack of this knowledge. We apply hybrid (continuous-time dynamics, discrete-time measurements) Kalman filtering to linear, time-varying versions of the original nonlinear estimation problems, obtained by treating nonlinearities as unknown time-varying inputs. In the first case, we show that the unknown input can be pre-estimated thanks to the system conservation laws. In the second case, we use a Bayesian regularization approach to linearly estimate the unknown input along with the other system state variables. We demonstrate in simulation that our approaches perform well in presence of noise, and we prove that the estimator for a known growth law is an observer, i.e., it enjoys deterministic convergence in absence of noise.

The paper is organized as follows. In Sec. II, we describe bioreactor growth models and illustrate the properties of our model. In Sec. III we review concepts of observability and Kalman filtering estimation that are used in the sequel. In Sec. IV we introduce our estimation method for a perfectly known growth law, and prove its convergence properties. In Sec. V we present our method to address state estimation for unknown growth laws. In Sec. VI, based on noisy simulated data, we compare performance of our methods with the reference approach of Extended Kalman filtering. Conclusions and perspectives are drawn in Sec. VII.

## II. MODEL AND PROPERTIES

The model of a continuously stirred-tank bioreactor (CSTR) for microbial growth in constant volume without accounting for gaseous outflow is generally expressed as

$$ \dot{x} = Nr(x, t) + d(x_{\text{in}} - x) \quad (1) $$

where $x(t) \in \Omega = \mathbb{R}_+^n$ represents the reactants concentrations vector, $x_{\text{in}} \in \mathbb{R}_+^n$ is the reactor feed concentrations vector, $N \in \mathbb{R}^{n \times q}$ is the stoichiometric matrix, $r(x, t) \in \mathbb{R}_+^q$ is the reaction rates vector ($gL^{-1}h^{-1}$), and $d \in \mathbb{R}_+^*$ is the dilution rate ($h^{-1}$) [1]. Concentrations (in $gL^{-1}$) are relative

<footnote>*Work supported in part by project Ctrl-AB [ANR-20-CE45-0014]</footnote>

<footnote n="1">¹Université Grenoble Alpes, Inria, Grenoble, France</footnote>

<footnote n="2">²Université Côte d'Azur, Inria, INRAE, CNRS, Sorbonne Université, Macbes Team, Sophia Antipolis, France</footnote>

<footnote>†Corresponding author, rand.asswad@inria.fr</footnote>
# Kalman-based approaches for online estimation of bioreactor dynamics from fluorescent reporter measurements*

Rand Asswad<fnref n="1" />,†, Eugenio Cinquemani<fnref n="1" />, Jean-Luc Gouzé<fnref n="2" />

**Abstract**— We address online estimation of microbial growth dynamics in bioreactors from measurements of a fluorescent reporter protein synthesized along with microbial growth. We consider an extended version of standard growth models that accounts for the dynamics of reporter synthesis. We develop state estimation from sampled, noisy measurements in the cases of known and unknown growth rate functions. Leveraging conservation laws and regularized estimation techniques, we reduce these nonlinear estimation problems to linear time-varying ones, and solve them via Kalman filtering. We establish convergence results in absence of noise and show performance on noisy data in simulation.

## I. INTRODUCTION

Estimation and control of cellular growth in bioreactors have been dedicated a significant amount of work over several decades [1]. Nowadays, these challenges receive renewed attention in relation with new biotechnological developments. Among the new research frontiers is the feedback control of microbial consortia [2]. Motivated by the potential to outperform single species in biosynthetic processes and other applications [3], control of microbial consortia requires in the first place to be able to discern single species dynamics by suitable real-time monitoring, e.g. the use of fluorescent reporter proteins expressed along with the growth of the different species.

In a recent paper [4], we have explored (local) observability of the prototypical microbial consortium model from [5] under different assumptions about the observed variables. The construction and analysis of observers for (noisy) time-sampled measurements (as encountered frequently, e.g. with the use of pipetting robots [6]) have not been explored in detail. Given the nonlinear nature of these models, the problem is nontrivial. General observer synthesis approaches exist [7], [8], however their viability depends on perfect knowledge of the system, which hardly applies to biological systems, and/or on idealized assumptions about the measurement process. Robust estimation approaches can tackle in part modelling uncertainty [9], but they are generally not conceived for sampled data. Nonlinear generalizations of Kalman filtering exist that cope explicitly with sampled noisy data [10], [11], but performance is usually not guaranteed even with a perfectly known model.

In this paper, we focus on the construction of state estimators for a single species growing in a bioreactor from indirect,

sampled fluorescent reporter measurements of the species abundance. We consider classical nonlinear growth models [1], [12]) modified to account for the dynamical relationship between fluorescence abundance and biomass synthesis. As explained above, besides its interest per se, the problem constitutes a key step toward estimation of dynamics of microbial communities. Much literature has been dedicated to the problem in the case of direct observations of biomass abundance [1], [13] or gaseous outflow [14]. None of these approaches applies to the augmented model we consider with sampled, noisy measurements.

We consider two cases, perfect knowledge of the functional form of the growth law, and lack of this knowledge. We apply hybrid (continuous-time dynamics, discrete-time measurements) Kalman filtering to linear, time-varying versions of the original nonlinear estimation problems, obtained by treating nonlinearities as unknown time-varying inputs. In the first case, we show that the unknown input can be pre-estimated thanks to the system conservation laws. In the second case, we use a Bayesian regularization approach to linearly estimate the unknown input along with the other system state variables. We demonstrate in simulation that our approaches perform well in presence of noise, and we prove that the estimator for a known growth law is an observer, i.e., it enjoys deterministic convergence in absence of noise.

The paper is organized as follows. In Sec. II, we describe bioreactor growth models and illustrate the properties of our model. In Sec. III we review concepts of observability and Kalman filtering estimation that are used in the sequel. In Sec. IV we introduce our estimation method for a perfectly known growth law, and prove its convergence properties. In Sec. V we present our method to address state estimation for unknown growth laws. In Sec. VI, based on noisy simulated data, we compare performance of our methods with the reference approach of Extended Kalman filtering. Conclusions and perspectives are drawn in Sec. VII.

## II. MODEL AND PROPERTIES

The model of a continuously stirred-tank bioreactor (CSTR) for microbial growth in constant volume without accounting for gaseous outflow is generally expressed as

$$ \dot{x} = Nr(x, t) + d(x_{in} - x) \quad (1) $$

where $x(t) \in \Omega = \mathbb{R}_+^n$ represents the reactants concentrations vector, $x_{in} \in \mathbb{R}_+^n$ is the reactor feed concentrations vector, $N \in \mathbb{R}^{n \times q}$ is the stoichiometric matrix, $r(x, t) \in \mathbb{R}_+^q$ is the reaction rates vector ($gL^{-1}h^{-1}$), and $d \in \mathbb{R}_+^*$ is the dilution rate ($h^{-1}$) [1]. Concentrations (in $gL^{-1}$) are relative

<footnote>*Work supported in part by project Ctrl-AB [ANR-20-CE45-0014]</footnote>

<footnote n="1">¹Université Grenoble Alpes, Inria, Grenoble, France</footnote>

<footnote n="2">²Université Côte d'Azur, Inria, INRAE, CNRS, Sorbonne Université, Macbes Team, Sophia Antipolis, France</footnote>

<footnote>†Corresponding author, rand.asswad@inria.fr</footnote>
to the total culture volume. We assume a strictly positive dilution rate $d > 0$. Consider in particular the dynamical system describing the growth of *Escherichia coli* (*E. coli*) bacteria along with synthesis of a fluorescent protein. We assume as in [4] that the latter takes part of the resources away from the synthesis of the self-replicative biomass. For $e(t)$, $f(t)$ and $s(t)$ the biomass, fluorescent protein, and substrate concentration, in the same order, and $x(t) = (s(t), e(t), f(t))^{\top} \in \mathbb{R}_+^3$, the model becomes [4], [15], [16]

$$ 
\begin{aligned}
\dot{s} &= -\frac{1}{\gamma}\mu(x,t)e + d(s_{\text{in}} - s) \\
\dot{e} &= (1 - \alpha)\mu(x,t)e - de \\
\dot{f} &= \alpha\mu(x,t)e - df
\end{aligned} 
\quad (2) $$

where $\mu(x,t) \in \mathbb{R}_+$, the growth rate per capita at time $t \ge t_0$, may depend on time due to quantities that are not part of the model state (pH, temperature) [1]. Constant $0 < \alpha < 1$ is the proportion of substrate import dedicated to fluorescent protein synthesis, while $s_{\text{in}}$ represents the inflow substrate concentration. Overall, this is a constant yield model where the substrate uptake rate per capita is given as $\mu(x,t)/\gamma$ with $\gamma$ the growth yield coefficient. The values of $\alpha$, $s_{\text{in}}$, and $\gamma$ are assumed known and strictly positive. The model is in the form of (1) with $N = (-1/\gamma, 1-\alpha, \alpha)^{\top} \in \mathbb{R}^3$, $r(x,t) = \mu(x,t)e(t) \in \mathbb{R}_+$, and $x_{\text{in}} = (s_{\text{in}}, 0, 0)^{\top} \in \mathbb{R}_+^3$.

For the problem of online state observation, in section IV we consider the specific growth rate known and strictly a function of the system's state. In section V, we will instead address the case where the functional form of $\mu(\cdot)$ is unknown, and thus treat it as an unknown function of time.

### A. Conservation relations

In a biochemical system, linear dependence of the rows of $N$ defines conservation relations. The rows of $N$ are linearly dependent if and only if the left null-space of $N$ has dimension $\ell = n - \text{rank}(N) > 0$, which corresponds to the number of conservation relations. A conservation matrix $\Pi \in \mathbb{R}^{\ell \times n}$ is defined as a matrix with rows forming a basis of the left-null space of $N$. Consequently, $\Pi N = 0$ and $\text{rank}(\Pi) = \ell$. The conservation matrix, which is not unique [17], [18], projects the system's state into an auxiliary state $\sigma = \Pi x$ of dimension $\ell$. We also refer to this definition as the conservation law. From (1),

$$ \dot{\sigma} = d(\Pi x_{\text{in}} - \sigma) \quad (3) $$

(in the “batch mode” $d \equiv 0$, $\dot{\sigma} = 0$, whence the name “conservation”). The conservation manifold (nullcline of the conservation law) is defined as $\Sigma = \{x \in \mathbb{R}^n, \Pi \dot{x} = 0\}$. From (3), $\sigma(t) = \Pi x_{\text{in}} + e^{-d(t-t_0)}(\sigma(t_0) - \Pi x_{\text{in}})$, which shows that $\Sigma$ is positively invariant. Since $d > 0$, $\Sigma$ is also exponentially attracting, thus it is asymptotic invariant [17].

System (2) has dimension $n = 3$ and $\text{rank}(N) = 1$. Hence there are $\ell = 2$ conservation relations. Notably, we choose

$$ \Pi = \begin{pmatrix} 1 & 0 & \frac{1}{\alpha\gamma} \\ 0 & 1 & -\frac{1-\alpha}{\alpha} \end{pmatrix} \implies \begin{cases} \sigma_1 = s + \frac{1}{\alpha\gamma}f \\ \sigma_2 = e - \frac{1-\alpha}{\alpha}f \end{cases} \quad (4) $$

### B. Steady states

Regardless of the form of $\mu(x,t)$, system (2) always admits an equilibrium $x_{\text{washout}} = (s_{\text{in}}, 0, 0)^{\top}$, which is globally asymptotically stable on $\Omega$ if the set of equilibria is a singleton and unstable otherwise. The existence and stability of other equilibria depends on $\mu(x,t)$ and on the system parameters, in this case $d$, $\alpha$, and $s_{\text{in}}$.

We focus on the Monod law [1]. In this case $\mu$ takes the form

$$ \mu(s) = \mu_{\text{max}} \frac{s}{k_s + s}, \quad s \ge 0, \quad (5) $$

with $k_s > 0$ known as half-saturation constant and $\mu_{\text{max}} > 0$ the maximal growth rate. Eq. (5) is strictly increasing, therefore it defines a bijection from $[0, \infty)$ to $[0, \mu_{\text{max}})$. If $d/(1-\alpha) < \mu(s_{\text{in}})$, (2) admits another equilibrium $x^* = (s^*, \gamma(1-\alpha)(s_{\text{in}}-s^*), \gamma\alpha(s_{\text{in}}-s^*))^{\top}$ with $s^* = \mu^{-1}(d/(1-\alpha))$. If this exists, it is globally asymptotically stable on any open subset of $\Omega \setminus \{e=0\}$. From now on we work under this assumption. These results extend similar known results for the subsystem $(s, e)$ without fluorescence synthesis [12]. They apply similarly to several other common growth models, such as Teissier law or Contois law [1].

## III. OBSERVABILITY AND STATE ESTIMATION

We review the system's observability and the definition of a state observer in the idealized case of noiseless measurements in continuous time. A dynamical system describes the time evolution of the system's state $x(t) \in \Omega \subset \mathbb{R}^n$. Typically, information about the state variables is obtained through measurements $y(t) = h(x(t)) \in \mathbb{R}^p$ (often referred to as the system's output), where $h$ is a function assumed $C^{\infty}(\Omega)$. The problem of interest in this paper is the construction of a state observer $\hat{x}(t) \in \mathbb{R}^n$ from the measurements $y(t)$ for all $t \ge t_0$. $\hat{x}(t)$ is an observer if (i) $\|x(t) - \hat{x}(t)\|$ converges to zero as $t$ tends to $\infty$, and (ii) $\hat{x}(t) = x(t) \implies \hat{x}(t+\tau) = x(t+\tau), \forall \tau \ge 0$. If the first condition holds for any $x(t_0)$ and $\hat{x}(t_0)$, the observer is global. If the convergence rate for (i) can be tuned, the observer is tunable. A tunable observer can be constructed if the system is observable, which is discussed in subsection III-A [7]. An open loop (non-tunable) observer (i.e. with no output feedback) is also referred to as a detector.

The available measurements for (2) are fluorescent protein concentrations $y(t) = f(t)$, thus $h$ is defined by $y(t) = Cx(t)$ with $C = (0, 0, 1)$.

### A. Observability

The local observability of a dynamical system $\dot{x} = f(x)$ is studied via its observation space $O(h)$, that is the vector field containing $h$ and its Lie derivatives along the vector field $f$. The Lie derivative along $f$ of a function $\varphi: \mathbb{R}^n \to \mathbb{R}^p$ is $L_f(\varphi) = \frac{\partial\varphi}{\partial x}f$. The system is locally observable if and only if there exists $n$ elements of $O(h)$ defining a diffeomorphism around $x$ for all $x \in \Omega$ [7]. By the chain rule, $\dot{y}(t) = \frac{\partial h}{\partial x}(x(t))\dot{x}(t) = L_f(h)(x(t))$. Ergo, the observation space is the space of time-derivatives of the outputs.
to the total culture volume. We assume a strictly positive dilution rate $d > 0$. Consider in particular the dynamical system describing the growth of *Escherichia coli* (*E. coli*) bacteria along with synthesis of a fluorescent protein. We assume as in [4] that the latter takes part of the resources away from the synthesis of the self-replicative biomass. For $e(t)$, $f(t)$ and $s(t)$ the biomass, fluorescent protein, and substrate concentration, in the same order, and $x(t) = (s(t), e(t), f(t))^{\top} \in \mathbb{R}_+^3$, the model becomes [4], [15], [16]

$$ 
\begin{aligned}
\dot{s} &= -\frac{1}{\gamma}\mu(x,t)e + d(s_{\text{in}} - s) \\
\dot{e} &= (1 - \alpha)\mu(x,t)e - de \\
\dot{f} &= \alpha\mu(x,t)e - df
\end{aligned} 
\quad (2) 
$$

where $\mu(x,t) \in \mathbb{R}_+$, the growth rate per capita at time $t \ge t_0$, may depend on time due to quantities that are not part of the model state (pH, temperature) [1]. Constant $0 < \alpha < 1$ is the proportion of substrate import dedicated to fluorescent protein synthesis, while $s_{\text{in}}$ represents the inflow substrate concentration. Overall, this is a constant yield model where the substrate uptake rate per capita is given as $\mu(x,t)/\gamma$ with $\gamma$ the growth yield coefficient. The values of $\alpha$, $s_{\text{in}}$, and $\gamma$ are assumed known and strictly positive. The model is in the form of (1) with $N = (-1/\gamma, 1-\alpha, \alpha)^{\top} \in \mathbb{R}^3$, $r(x,t) = \mu(x,t)e(t) \in \mathbb{R}_+$, and $x_{\text{in}} = (s_{\text{in}}, 0, 0)^{\top} \in \mathbb{R}_+^3$.

For the problem of online state observation, in section IV we consider the specific growth rate known and strictly a function of the system's state. In section V, we will instead address the case where the functional form of $\mu(\cdot)$ is unknown, and thus treat it as an unknown function of time.

### A. Conservation relations

In a biochemical system, linear dependence of the rows of $N$ defines conservation relations. The rows of $N$ are linearly dependent if and only if the left null-space of $N$ has dimension $l = n - \text{rank}(N) > 0$, which corresponds to the number of conservation relations. A conservation matrix $\Pi \in \mathbb{R}^{l \times n}$ is defined as a matrix with rows forming a basis of the left-null space of $N$. Consequently, $\Pi N = 0$ and $\text{rank}(\Pi) = l$. The conservation matrix, which is not unique [17], [18], projects the system's state into an auxiliary state $\sigma = \Pi x$ of dimension $l$. We also refer to this definition as the conservation law. From (1),

$$ \dot{\sigma} = d(\Pi x_{\text{in}} - \sigma) \quad (3) $$

(in the “batch mode” $d \equiv 0$, $\dot{\sigma} = 0$, whence the name “conservation”). The conservation manifold (nullcline of the conservation law) is defined as $\Sigma = \{x \in \mathbb{R}^n, \Pi \dot{x} = 0\}$. From (3), $\sigma(t) = \Pi x_{\text{in}} + e^{-d(t-t_0)}(\sigma(t_0) - \Pi x_{\text{in}})$, which shows that $\Sigma$ is positively invariant. Since $d > 0$, $\Sigma$ is also exponentially attracting, thus it is asymptotic invariant [17].

System (2) has dimension $n = 3$ and $\text{rank}(N) = 1$. Hence there are $l = 2$ conservation relations. Notably, we choose

$$ \Pi = \begin{pmatrix} 1 & 0 & \frac{1}{\alpha\gamma} \\ 0 & 1 & -\frac{1-\alpha}{\alpha} \end{pmatrix} \implies \left\{ \begin{array}{l} \sigma_1 = s + \frac{1}{\alpha\gamma}f \\ \sigma_2 = e - \frac{1-\alpha}{\alpha}f \end{array} \right. . \quad (4) $$

### B. Steady states

Regardless of the form of $\mu(x,t)$, system (2) always admits an equilibrium $x_{\text{washout}} = (s_{\text{in}}, 0, 0)^{\top}$, which is globally asymptotically stable on $\Omega$ if the set of equilibria is a singleton and unstable otherwise. The existence and stability of other equilibria depends on $\mu(x,t)$ and on the system parameters, in this case $d$, $\alpha$, and $s_{\text{in}}$.

We focus on the Monod law [1]. In this case $\mu$ takes the form

$$ \mu(s) = \mu_{\text{max}} \frac{s}{k_s + s}, \quad s \ge 0, \quad (5) $$

with $k_s > 0$ known as half-saturation constant and $\mu_{\text{max}} > 0$ the maximal growth rate. Eq. (5) is strictly increasing, therefore it defines a bijection from $[0, \infty)$ to $[0, \mu_{\text{max}})$. If $d/(1-\alpha) < \mu(s_{\text{in}})$, (2) admits another equilibrium $x^* = (s^*, \gamma(1-\alpha)(s_{\text{in}} - s^*), \gamma\alpha(s_{\text{in}} - s^*))^{\top}$ with $s^* = \mu^{-1}(d/(1-\alpha))$. If this exists, it is globally asymptotically stable on any open subset of $\Omega \setminus \{e=0\}$. From now on we work under this assumption. These results extend similar known results for the subsystem $(s, e)$ without fluorescence synthesis [12]. They apply similarly to several other common growth models, such as Teissier law or Contois law [1].

## III. OBSERVABILITY AND STATE ESTIMATION

We review the system's observability and the definition of a state observer in the idealized case of noiseless measurements in continuous time. A dynamical system describes the time evolution of the system's state $x(t) \in \Omega \subset \mathbb{R}^n$. Typically, information about the state variables is obtained through measurements $y(t) = h(x(t)) \in \mathbb{R}^p$ (often referred to as the system's output), where $h$ is a function assumed $C^{\infty}(\Omega)$. The problem of interest in this paper is the construction of a state observer $\hat{x}(t) \in \mathbb{R}^n$ from the measurements $y(t)$ for all $t \ge t_0$. $\hat{x}(t)$ is an observer if (i) $\|x(t) - \hat{x}(t)\|$ converges to zero as $t$ tends to $\infty$, and (ii) $\hat{x}(t) = x(t) \implies \hat{x}(t+\tau) = x(t+\tau), \forall \tau \ge 0$. If the first condition holds for any $x(t_0)$ and $\hat{x}(t_0)$, the observer is global. If the convergence rate for (i) can be tuned, the observer is tunable. A tunable observer can be constructed if the system is observable, which is discussed in subsection III-A [7]. An open loop (non-tunable) observer (i.e. with no output feedback) is also referred to as a detector.

The available measurements for (2) are fluorescent protein concentrations $y(t) = f(t)$, thus $h$ is defined by $y(t) = Cx(t)$ with $C = (0, 0, 1)$.

### A. Observability

The local observability of a dynamical system $\dot{x} = f(x)$ is studied via its observation space $O(h)$, that is the vector field containing $h$ and its Lie derivatives along the vector field $f$. The Lie derivative along $f$ of a function $\varphi: \mathbb{R}^n \to \mathbb{R}^p$ is $L_f(\varphi) = \frac{\partial\varphi}{\partial x}f$. The system is locally observable if and only if there exists $n$ elements of $O(h)$ defining a diffeomorphism around $x$ for all $x \in \Omega$ [7]. By the chain rule, $\dot{y}(t) = \frac{\partial h}{\partial x}(x(t))\dot{x}(t) = L_f(h)(x(t))$. Ergo, the observation space is the space of time-derivatives of the outputs.
For (2), $y = f$, therefore $\dot{y} = \alpha\mu(x, t)e - d \cdot y$. Setting $z = (d \cdot y + \dot{y})/\alpha = \mu(x_g, t)e$ defines a reduced dynamical system for the state $x_g = (s, e)^\top$ with output $z = h_g(x_g) = \mu(x_g, t)e$, provided $\mu(x, t)$ does not involve $f$. Consequently, $O(h_g) \subset O(h)$, and since a diffeomorphism exists from $f$ to $z$ by construction, the system is locally observable if the reduced system is observable.

Biologically, output $z = \mu(x_g, t)e$ corresponds to measurements of biogas flow rate in the anaerobic digestion process, which is studied and proven locally observable under suitable assumptions in [14]. Observability of (2) for the Monod growth law follows from the results in [14].

## B. Estimation problem and Kalman filtering

We stated the observer problem in continuous time, but in truth measurements are in discrete time. This is not an issue in practice if the time steps are sufficiently small (for example, for continuously operating sensors built on the bioreactor). In many experimental scenarios of interest, however, measurements are obtained by time-consuming procedures (see e.g. [6]), hence motivates us to explicitly account for the sampled nature of the data. Moreover, measurements are typically noisy due to biophysical uncertainty in the measurement process This gives rise to an estimation problem, as follows.

Let measurements at times $t_0 < t_1 < \dots < t_{M-1} < t_M$ be given by $y_k = Cx(t_k) + v_k$, where $v_k$ is an independent Gaussian noise with zero mean and known covariance matrix $R_k \in \mathbb{R}^{p \times p}$ (scalar variance if $p=1$). The problem is to estimate the continuous-time state $x(t)$ from the discrete-time measurements $y_k$. It can be addressed by a hybrid (continuous-discrete) Kalman filter [10]. This filter is optimal for linear state-space models, while our system dynamics are nonlinear. Though variants exist that cope with nonlinear dynamics, such as the Extended KF (EKF) or the Unscented KF (UKF) [10], [11], [19], convergence and performance of these filters are hardly proven. We instead propose to reformulate the estimation problem by treating the nonlinearities of the system dynamics as a time-varying input. This input is either pre-estimated (as explained in Sec. IV) or described by a probabilistic prior in the form of a conveniently chosen linear stochastic process (as explained in Sec. V). In both cases, this leads to a hybrid filtering problem on a linear (time-varying) stochastic system of the form

$$ \begin{aligned} \dot{x}(t) &= A(t)x(t) + b(t) + G(t)w(t), \\ y_k &= C_k x(t_k) + v_k, \end{aligned} \quad (6) $$

with $A(t) \in \mathbb{R}^{n \times n}$, $b(t) \in \mathbb{R}^n$, $C_k \in \mathbb{R}^{p \times n}$, $G(t) \in \mathbb{R}^{n \times m}$, $Q(t) \in \mathbb{R}^{m \times m}$ and $R_k \in \mathbb{R}^{p \times p}$ known. The process noise $\{w(t), t \ge t_0\}$ is formally the derivative of the Brownian motion [10], i.e. it is a zero-mean, stochastic process with $E\{w(t)w(\tau)^\top\} = Q(t)\delta(t - \tau)$, $\delta$ being the Dirac distribution.

The hybrid KF, which is optimal for this class of models, alternates prediction based on the dynamical model and update based on a new measurement. The prediction step

consists of integrating the system's dynamics for the state estimate $\hat{x}_{k|k-1}(t)$ and a Riccati equation for the state covariance matrix $P_{k|k-1}(t)$ over each time interval $[t_{k-1}, t_k]$, with $\hat{x}_{k|k-1}(t_{k-1}) = \hat{x}_{k-1|k-1}$ and $P_{k|k-1}(t_{k-1}) = P_{k-1|k-1}$. The update step applied to $\hat{x}_{k|k-1}(t_k)$ and $P_{k|k-1}(t_k)$ for a new measurement $y_k$ is identical to the discrete KF [10].

## C. Parameter tuning from data

In our estimator designs, Eq. (6) will comprise unknown parameters $\theta \in \Theta$. We propose (and later demonstrate in simulation) to estimate these parameters from a preliminary data set $Y_M = \{y_0, \dots, y_M\}$ through a Maximum Likelihood (ML) approach, which closely relates with the filter's prediction error methods for tuning filters [20].

Given $Y_M$, ML estimation is the problem of maximizing the value of the joint distribution $f(Y_M|\theta)$ over $\theta \in \Theta$. By recursive application of the Bayes law, one can write $f(Y_M|\theta)$ as the product of the conditional (Gaussian) densities $f(y_k|Y_{k-1}, \theta)$ for $k=0, \dots, M$. Mean and variance of these densities correspond to the Kalman output predictor $C_k \hat{x}_{k|k-1}$ and its covariance $S_k = C_k P_{k|k-1} C_k^\top + R_k$ [16], [20]. Writing ML as the equivalent negative log-likelihood minimization and working out the equations, the ML estimator of $\theta$ is

$$ \hat{\theta} = \underset{\theta \in \Theta}{\operatorname{argmin}} \left[ -\frac{1}{2} \sum_{k=0}^{M} p \log(2\pi) + \log |S_k| + \tilde{y}_k^\top S_k^{-1} \tilde{y}_k \right]. \quad (7) $$

For any candidate $\theta$, all of these quantities can be computed by the corresponding Kalman filtering iteration, which enables fast numerical solution of the optimization problem.

## IV. APPROACH FOR KNOWN $\mu(x)$

In this section, leveraging conservation laws, we propose a pipeline for state estimation for the case where $\mu$ is a known function of the state variables. We develop this for system (2) with Monod growth law, yet, the approach is viable for other growth laws and similar biochemical systems. The resulting estimator is proven to be a global observer.

As seen in section II, nonlinearities of CSTR dynamics are in the reaction rate vector $r(x, t)$ of (1). For system (2), it is however possible to compute a pre-estimate $\hat{\mu}(t)$ of the specific growth rate $\mu$. Replacing $\mu$ in (2) with $\hat{\mu}(t)$ yields linear time-varying dynamics, eligible for linear Kalman filtering. We next show how to compute $\hat{\mu}(t)$ in the case of the Monod growth rate $\mu(s)$.

### A. Proposed pipeline

The Monod law is a function of the unknown $s(t)$. By the conservation law (4), in the (exponentially attracting) conservation manifold, $s = \sigma_1 - \frac{1}{\alpha\gamma}f$. In view of this, we define a pre-estimate $\hat{s}_\sigma$ of $s$ by $\hat{s}_\sigma = \hat{\sigma}_1 - \frac{1}{\alpha\gamma}\hat{y}$, with $\hat{\sigma}_1$ calculated based on the dynamics of $\sigma_1$ and $\hat{y}$ a suitable smoothening of the measurements $y_k$ of $f$, as detailed next.

The conservation variable $\sigma_1$ is estimated through an open-loop detector $\hat{\sigma}_1(t)$ with dynamics as in (3), resulting in

$$ \hat{\sigma}_1(t) = s_{in} + e^{-d(t-t_0)} (\hat{\sigma}_1(t_0) - s_{in}). \quad (8) $$

This converges to $s_{in}$ exponentially at a rate $d$.
For (2), $y = f$, therefore $\dot{y} = \alpha\mu(x, t)e - d \cdot y$. Setting $z = (d \cdot y + \dot{y})/\alpha = \mu(x_g, t)e$ defines a reduced dynamical system for the state $x_g = (s, e)^\top$ with output $z = h_g(x_g) = \mu(x_g, t)e$, provided $\mu(x, t)$ does not involve $f$. Consequently, $O(h_g) \subset O(h)$, and since a diffeomorphism exists from $f$ to $z$ by construction, the system is locally observable if the reduced system is observable.

Biologically, output $z = \mu(x_g, t)e$ corresponds to measurements of biogas flow rate in the anaerobic digestion process, which is studied and proven locally observable under suitable assumptions in [14]. Observability of (2) for the Monod growth law follows from the results in [14].

## B. Estimation problem and Kalman filtering

We stated the observer problem in continuous time, but in truth measurements are in discrete time. This is not an issue in practice if the time steps are sufficiently small (for example, for continuously operating sensors built on the bioreactor). In many experimental scenarios of interest, however, measurements are obtained by time-consuming procedures (see e.g. [6]), hence motivates us to explicitly account for the sampled nature of the data. Moreover, measurements are typically noisy due to biophysical uncertainty in the measurement process This gives rise to an estimation problem, as follows.

Let measurements at times $t_0 < t_1 < \dots < t_{M-1} < t_M$ be given by $y_k = Cx(t_k) + v_k$, where $v_k$ is an independent Gaussian noise with zero mean and known covariance matrix $R_k \in \mathbb{R}^{p \times p}$ (scalar variance if $p = 1$). The problem is to estimate the continuous-time state $x(t)$ from the discrete-time measurements $y_k$. It can be addressed by a hybrid (continuous-discrete) Kalman filter [10]. This filter is optimal for linear state-space models, while our system dynamics are nonlinear. Though variants exist that cope with nonlinear dynamics, such as the Extended KF (EKF) or the Unscented KF (UKF) [10], [11], [19], convergence and performance of these filters are hardly proven. We instead propose to reformulate the estimation problem by treating the nonlinearities of the system dynamics as a time-varying input. This input is either pre-estimated (as explained in Sec. IV) or described by a probabilistic prior in the form of a conveniently chosen linear stochastic process (as explained in Sec. V). In both cases, this leads to a hybrid filtering problem on a linear (time-varying) stochastic system of the form

$$ \begin{aligned} \dot{x}(t) &= A(t)x(t) + b(t) + G(t)w(t), \\ y_k &= C_k x(t_k) + v_k, \end{aligned} \quad (6) $$

with $A(t) \in \mathbb{R}^{n \times n}$, $b(t) \in \mathbb{R}^n$, $C_k \in \mathbb{R}^{p \times n}$, $G(t) \in \mathbb{R}^{n \times m}$, $Q(t) \in \mathbb{R}^{m \times m}$ and $R_k \in \mathbb{R}^{p \times p}$ known. The process noise $\{w(t), t \ge t_0\}$ is formally the derivative of the Brownian motion [10], i.e. it is a zero-mean, stochastic process with $E\{w(t)w(\tau)^\top\} = Q(t)\delta(t - \tau)$, $\delta$ being the Dirac distribution.

The hybrid KF, which is optimal for this class of models, alternates prediction based on the dynamical model and update based on a new measurement. The prediction step

consists of integrating the system's dynamics for the state estimate $\hat{x}_{k|k-1}(t)$ and a Riccati equation for the state covariance matrix $P_{k|k-1}(t)$ over each time interval $[t_{k-1}, t_k]$, with $\hat{x}_{k|k-1}(t_{k-1}) = \hat{x}_{k-1|k-1}$ and $P_{k|k-1}(t_{k-1}) = P_{k-1|k-1}$. The update step applied to $\hat{x}_{k|k-1}(t_k)$ and $P_{k|k-1}(t_k)$ for a new measurement $y_k$ is identical to the discrete KF [10].

## C. Parameter tuning from data

In our estimator designs, Eq. (6) will comprise unknown parameters $\theta \in \Theta$. We propose (and later demonstrate in simulation) to estimate these parameters from a preliminary data set $Y_M = \{y_0, \dots, y_M\}$ through a Maximum Likelihood (ML) approach, which closely relates with the filter's prediction error methods for tuning filters [20].

Given $Y_M$, ML estimation is the problem of maximizing the value of the joint distribution $f(Y_M|\theta)$ over $\theta \in \Theta$. By recursive application of the Bayes law, one can write $f(Y_M|\theta)$ as the product of the conditional (Gaussian) densities $f(y_k|Y_{k-1}, \theta)$ for $k = 0, \dots, M$. Mean and variance of these densities correspond to the Kalman output predictor $C_k \hat{x}_{k|k-1}$ and its covariance $S_k = C_k P_{k|k-1} C_k^\top + R_k$ [16], [20]. Writing ML as the equivalent negative log-likelihood minimization and working out the equations, the ML estimator of $\theta$ is

$$ \hat{\theta} = \underset{\theta \in \Theta}{\operatorname{argmin}} \left[ -\frac{1}{2} \sum_{k=0}^{M} p \log(2\pi) + \log |S_k| + \tilde{y}_k^\top S_k^{-1} \tilde{y}_k \right]. \quad (7) $$

For any candidate $\theta$, all of these quantities can be computed by the corresponding Kalman filtering iteration, which enables fast numerical solution of the optimization problem.

## IV. APPROACH FOR KNOWN $\mu(x)$

In this section, leveraging conservation laws, we propose a pipeline for state estimation for the case where $\mu$ is a known function of the state variables. We develop this for system (2) with Monod growth law, yet, the approach is viable for other growth laws and similar biochemical systems. The resulting estimator is proven to be a global observer.

As seen in section II, nonlinearities of CSTR dynamics are in the reaction rate vector $r(x, t)$ of (1). For system (2), it is however possible to compute a pre-estimate $\hat{\mu}(t)$ of the specific growth rate $\mu$. Replacing $\mu$ in (2) with $\hat{\mu}(t)$ yields linear time-varying dynamics, eligible for linear Kalman filtering. We next show how to compute $\hat{\mu}(t)$ in the case of the Monod growth rate $\mu(s)$.

### A. Proposed pipeline

The Monod law is a function of the unknown $s(t)$. By the conservation law (4), in the (exponentially attracting) conservation manifold, $s = \sigma_1 - \frac{1}{\alpha\gamma}f$. In view of this, we define a pre-estimate $\hat{s}_\sigma$ of $s$ by $\hat{s}_\sigma = \hat{\sigma}_1 - \frac{1}{\alpha\gamma}\hat{y}$, with $\hat{\sigma}_1$ calculated based on the dynamics of $\sigma_1$ and $\hat{y}$ a suitable smoothening of the measurements $y_k$ of $f$, as detailed next.

The conservation variable $\sigma_1$ is estimated through an open-loop detector $\hat{\sigma}_1(t)$ with dynamics as in (3), resulting in

$$ \hat{\sigma}_1(t) = s_{\text{in}} + e^{-d(t-t_0)} (\hat{\sigma}_1(t_0) - s_{\text{in}}). \quad (8) $$

This converges to $s_{\text{in}}$ exponentially at a rate $d$.
To calculate $\hat{y}$, a Kalman-based regularization method is used. In the same spirit as [21], the method introduces a regularizing linear stochastic dynamical model for the time profile $f$ to estimate, and solves the problem by linear state-space estimation [22]. We consider a second-order linear stochastic differential equation $\ddot{y}(t) = \beta w(t)$ with $\{w(t)\}$ a scalar white noise process, and assume that the noisy measurements $y_k$ of $f$ are noisy measurements of process $y$. Such model of the data can obviously be expressed in the form of (6). An optimal parameter $\hat{\beta}$ is the solution to the ML problem discussed in subsection III-C. This provides a model of the data-generating process $f$ independent of the knowledge of $f$. Applying a Kalman filter based on this model to the measurements $y_k$ yields estimates of $y$, i.e. of $f$. The chosen form of process $y$ is known to penalize overly irregular solutions [16], [21], whereas the ML estimation of $\hat{\beta}$ ensures an appropriate tradeoff between noisy data interpolation and regularity of the solution. The filter produces piecewise-continuous predictions $\hat{y}_{k|k-1}(t)$ and filtered estimates $\hat{y}_{k|k}$. Both choices being equally viable as a definition of the estimate $\hat{y}(t)$ sought, we simply refer to $\hat{y}(t)$ from now on and add details only when necessary.

Armed with $\hat{\sigma}_1$ and $\hat{y}$, one can now compute $\hat{\mu}(t) = \mu(\hat{s}_\sigma(t))$ with $\hat{s}_\sigma(t) = \hat{\sigma}_1(t) - \frac{1}{\alpha\gamma}\hat{y}(t)$, and plug it into (2) to approach Kalman filtering estimation of the whole system state. Let $\hat{\mu}(t) = \mu(s(t)) - \hat{\mu}(t)$ be the pre-estimation error of $\mu$ at time $t$. Then

$$ \begin{aligned} \dot{x} &= N(\hat{\mu}(t) + \tilde{\mu}(t))e + d(x_{\text{in}} - x) \\ &= A_{\sigma}(t)x + b + Nw(t) \end{aligned} \quad (9) $$

with $w(t) = \tilde{\mu}(t)e(t)$ as the process error term, where

$$ A_{\sigma}(t) = \begin{pmatrix} -d & -\hat{\mu}(t)/\gamma & 0 \\ 0 & (1-\alpha)\hat{\mu}(t) - d & 0 \\ 0 & \alpha\hat{\mu}(t) & -d \end{pmatrix}, b = \begin{pmatrix} ds_{\text{in}} \\ 0 \\ 0 \end{pmatrix} $$

are now known. To run the linear Kalman filter of subsection III-B for this system (with $G = N$ and $C = (0, 0, 1)$ in Eq. (6)), we exploit the analytic expression of $\hat{\sigma}_1(t)$ in (8), for a given $\hat{\sigma}_1(t_0)$, while $\hat{y}$ and $\hat{\mu}$ can be obtained online sequentially for every interval $T_k = [t_k, t_{k+1})$ for $k = 0, \dots, M-1$. The process noise $w(t)$ remains unknown, and linearly dependent on the state variable $e$. In the next subsection, an approximate characterization of $w(t)$ is proposed, which is used to tune the KF. The specific choice of $Q(t)$ is especially important in determining performance in presence of measurement noise.

## B. Tuning of $Q(t)$

We now derive an approximate expression for the covariance $Q(t)$ of $w(t) = \tilde{\mu}(t)e(t)$. Consider the linear approximation of $\mu$ around $s(t)$. For sufficiently small $\tilde{s}_\sigma = s - \hat{s}_\sigma$, $\mu(\tilde{s}_\sigma) = \mu(s - \hat{s}_\sigma) \approx \mu(s) - \mu'(s)\tilde{s}_\sigma$. Then $\tilde{\mu}(t) = \mu'(s(t))\tilde{s}_\sigma(t)$ implies $w(t) = \mu'(s(t))e(t)\tilde{s}_\sigma(t)$. We choose to study the process noise around the system equilibrium. From (4), (8) and subsection II-B, it can be shown that for $t \in T_k$, $\tilde{s}_\sigma(t) = v_k/\alpha\gamma$ at the system's steady states. Consequently, $w(t) = \frac{\mu'(s^*)e^*}{\alpha\gamma}v_k$ for $t \in T_k$, near

the system's equilibria. It follows that $\{w(t)\}$ is a piecewise constant stochastic process. By integrating $E\{w(t)w(\tau)\} = Q(t)\delta(t - \tau)$ over $\mathbb{R}$ we obtain the piecewise-constant function $Q(t) = \left(\frac{\mu'(s^*)e^*}{\alpha\gamma}\right)^2 R_k(t_{k+1} - t_k)$ for $t \in T_k$.

## C. Asymptotic stability of the filter

We now prove that the proposed filtering pipeline defines a global observer in the sense of Section III. That is, we consider noiseless measurements and look at convergence at the discrete measurement times [7]. We first prove convergence of the observer error, then its positive invariance, i.e. that $x_k - \hat{x}_k = 0$ implies $x_{k+l} - \hat{x}_{k+l} = 0$ for $l > 0$.

As seen before, the KF is applied to the pre-estimated system (9) with $A_\sigma(t)x = N\hat{\mu}(t)e - dx$. In the absence of measurement noise, $\hat{s}_\sigma(t)$ converges to $s^*$ at the same rate as $\hat{\sigma}_1$ (globally, exponentially at rate $d$). Consequently, $\hat{\mu}(t) = \mu(\hat{s}_\sigma(t))$ converges to $\mu^* = \mu(s^*) = d/(1 - \alpha)$ and the system (9) is asymptotically time-invariant with $A_\sigma(t)$ converging globally to $A$.

$$ A = \begin{pmatrix} -d & -\mu^*/\gamma & 0 \\ 0 & 0 & 0 \\ 0 & \alpha\mu^* & -d \end{pmatrix}. $$

Integrating the linear time-invariant dynamics over $T_k$ gives the difference equation

$$ \begin{aligned} \dot{x}(t) &= Ax(t) + b + Nw(t) \\ x_{k+1} &= e^{A(t_{k+1}-t_k)}x_k + \int_{T_k} e^{A(t_{k+1}-t)}(b + Nw(t))dt \\ &= F_kx_k + b_k + Nw_k \end{aligned} $$

with $F_k = \exp(A(t_{k+1}-t_k))$ the state transition matrix [10], $b_k = (\int_{T_k} \exp(A(t_{k+1}-t_k))dt)_b$, and $w_k = \int_{T_k} w(t)dt$. Let us fix a time step $\Delta t = t_{k+1} - t_k, \forall k$. Then $F_k = \exp(A\Delta t) = F$, and one can verify that $e^{A\Delta t}N = N, \forall \tau$. For a fixed $Q_c$ for the continuous system, $E\{w_k w_l\} = Q\delta_{k,l}$ with $Q = Q_c\Delta t$.

The equivalent discrete Kalman predictor equation is

$$ \hat{x}_{k+1|k} = (F - K_k C)\hat{x}_{k|k-1} + K_k y_k + b_k \quad (10) $$

with $K_k \in \mathbb{R}^{n \times 1}$ the Kalman gain. If the pair $[F, C]$ is completely detectable and $[F, NG]$ is completely stabilizable for any $G$ such that $GG^\top = Q$, then (i) for any positive semidefinite initial matrix $P_0$, $P_{k+1|k}$ converges to the solution $P_\infty$ of the algebraic Riccati equation, which defines the associated gain $K_\infty = FP_\infty C^\top(CP_\infty C^\top + R)^{-1}$, (ii) all eigenvalues $\lambda$ of $F - K_\infty C$ verify $|\lambda| < 1$ [19]. The pair $[F, C]$ is completely detectable if $\exists K \in \mathbb{R}^3 : |\lambda| < 1, \forall \lambda \in \text{Sp}(F - KC)$ [19]. Let $K = (0, k_2, k_3)^\top$. For $k_2 = (1-\alpha)/4\alpha$, $k_3 = e^{-d\Delta t}$, $\text{Sp}(F - KC) = \{e^{-d\Delta t}, (1 \pm e^{-d\Delta t/2})/2\}$, which fulfills the condition for complete detectability. Similarly, the pair $[F, NG]$ is completely stabilizable if $\exists K \in \mathbb{R}^3 : |\lambda| < 1, \forall \lambda \in \text{Sp}(F - NGK^\top)$ [19]. Here $G = \sqrt{Q} > 0$. Taking $K = (0, \frac{1}{G(1-\alpha)\Delta t}, 0)^\top$, $\text{Sp}(F - NGK^\top) = \{e^{-d\Delta t}, 0\}$, which proves complete controllability. Hence $P_\infty$ and $K_\infty$ are well-defined.
To calculate $\hat{y}$, a Kalman-based regularization method is used. In the same spirit as [21], the method introduces a regularizing linear stochastic dynamical model for the time profile $f$ to estimate, and solves the problem by linear state-space estimation [22]. We consider a second-order linear stochastic differential equation $\ddot{y}(t) = \beta w(t)$ with $\{w(t)\}$ a scalar white noise process, and assume that the noisy measurements $y_k$ of $f$ are noisy measurements of process $y$. Such model of the data can obviously be expressed in the form of (6). An optimal parameter $\hat{\beta}$ is the solution to the ML problem discussed in subsection III-C. This provides a model of the data-generating process $f$ independent of the knowledge of $f$. Applying a Kalman filter based on this model to the measurements $y_k$ yields estimates of $y$, i.e. of $f$. The chosen form of process $y$ is known to penalize overly irregular solutions [16], [21], whereas the ML estimation of $\hat{\beta}$ ensures an appropriate tradeoff between noisy data interpolation and regularity of the solution. The filter produces piecewise-continuous predictions $\hat{y}_{k|k-1}(t)$ and filtered estimates $\hat{y}_{k|k}$. Both choices being equally viable as a definition of the estimate $\hat{y}(t)$ sought, we simply refer to $\hat{y}(t)$ from now on and add details only when necessary.

Armed with $\hat{\sigma}_1$ and $\hat{y}$, one can now compute $\hat{\mu}(t) = \mu(\hat{s}_\sigma(t))$ with $\hat{s}_\sigma(t) = \hat{\sigma}_1(t) - \frac{1}{\alpha\gamma}\hat{y}(t)$, and plug it into (2) to approach Kalman filtering estimation of the whole system state. Let $\hat{\mu}(t) = \mu(s(t)) - \hat{\mu}(t)$ be the pre-estimation error of $\mu$ at time $t$. Then

$$ 
\begin{aligned}
\dot{x} &= N(\hat{\mu}(t) + \tilde{\mu}(t))e + d(x_{\text{in}} - x) \\
&= A_{\sigma}(t)x + b + Nw(t)
\end{aligned}
\qquad (9) $$

with $w(t) = \tilde{\mu}(t)e(t)$ as the process error term, where

$$ A_{\sigma}(t) = \begin{pmatrix} -d & -\hat{\mu}(t)/\gamma & 0 \\ 0 & (1-\alpha)\hat{\mu}(t) - d & 0 \\ 0 & \alpha\hat{\mu}(t) & -d \end{pmatrix}, b = \begin{pmatrix} ds_{\text{in}} \\ 0 \\ 0 \end{pmatrix} $$

are now known. To run the linear Kalman filter of subsection III-B for this system (with $G = N$ and $C = (0, 0, 1)$ in Eq. (6)), we exploit the analytic expression of $\hat{\sigma}_1(t)$ in (8), for a given $\hat{\sigma}_1(t_0)$, while $\hat{y}$ and $\hat{\mu}$ can be obtained online sequentially for every interval $T_k = [t_k, t_{k+1})$ for $k = 0, \dots, M-1$. The process noise $w(t)$ remains unknown, and linearly dependent on the state variable $e$. In the next subsection, an approximate characterization of $w(t)$ is proposed, which is used to tune the KF. The specific choice of $Q(t)$ is especially important in determining performance in presence of measurement noise.

## B. Tuning of $Q(t)$

We now derive an approximate expression for the covariance $Q(t)$ of $w(t) = \tilde{\mu}(t)e(t)$. Consider the linear approximation of $\mu$ around $s(t)$. For sufficiently small $\tilde{s}_\sigma = s - \hat{s}_\sigma$, $\mu(\tilde{s}_\sigma) = \mu(s - \hat{s}_\sigma) \approx \mu(s) - \mu'(s)\tilde{s}_\sigma$. Then $\tilde{\mu}(t) = \mu'(s(t))\tilde{s}_\sigma(t)$ implies $w(t) = \mu'(s(t))e(t)\tilde{s}_\sigma(t)$. We choose to study the process noise around the system equilibrium. From (4), (8) and subsection II-B, it can be shown that for $t \in T_k$, $\tilde{s}_\sigma(t) = v_k/\alpha\gamma$ at the system's steady states. Consequently, $w(t) = \frac{\mu'(s^*)e^*}{\alpha\gamma}v_k$ for $t \in T_k$, near

the system's equilibria. It follows that $\{w(t)\}$ is a piecewise constant stochastic process. By integrating $E\{w(t)w(\tau)\} = Q(t)\delta(t - \tau)$ over $\mathbb{R}$ we obtain the piecewise-constant function $Q(t) = \left(\frac{\mu'(s^*)e^*}{\alpha\gamma}\right)^2 R_k(t_{k+1} - t_k)$ for $t \in T_k$.

## C. Asymptotic stability of the filter

We now prove that the proposed filtering pipeline defines a global observer in the sense of Section III. That is, we consider noiseless measurements and look at convergence at the discrete measurement times [7]. We first prove convergence of the observer error, then its positive invariance, i.e. that $x_k - \hat{x}_k = 0$ implies $x_{k+l} - \hat{x}_{k+l} = 0$ for $l > 0$.

As seen before, the KF is applied to the pre-estimated system (9) with $A_\sigma(t)x = N\hat{\mu}(t)e - dx$. In the absence of measurement noise, $\hat{s}_\sigma(t)$ converges to $s^*$ at the same rate as $\hat{\sigma}_1$ (globally, exponentially at rate $d$). Consequently, $\hat{\mu}(t) = \mu(\hat{s}_\sigma(t))$ converges to $\mu^* = \mu(s^*) = d/(1 - \alpha)$ and the system (9) is asymptotically time-invariant with $A_\sigma(t)$ converging globally to $A$.

$$ A = \begin{pmatrix} -d & -\mu^*/\gamma & 0 \\ 0 & 0 & 0 \\ 0 & \alpha\mu^* & -d \end{pmatrix}. $$

Integrating the linear time-invariant dynamics over $T_k$ gives the difference equation

$$ 
\begin{aligned}
\dot{x}(t) &= Ax(t) + b + Nw(t) \\
x_{k+1} &= e^{A(t_{k+1}-t_k)}x_k + \int_{T_k} e^{A(t_{k+1}-t)}(b + Nw(t))dt \\
&= F_kx_k + b_k + Nw_k
\end{aligned}
 $$

with $F_k = \exp(A(t_{k+1}-t_k))$ the state transition matrix [10], $b_k = (\int_{T_k} \exp(A(t_{k+1}-t_k))dt)_b$, and $w_k = \int_{T_k} w(t)dt$. Let us fix a time step $\Delta t = t_{k+1} - t_k, \forall k$. Then $F_k = \exp(A\Delta t) = F$, and one can verify that $e^{A\Delta t}N = N, \forall \tau$. For a fixed $Q_c$ for the continuous system, $E\{w_k w_l\} = Q\delta_{k,l}$ with $Q = Q_c\Delta t$.

The equivalent discrete Kalman predictor equation is

$$ \hat{x}_{k+1|k} = (F - K_k C)\hat{x}_{k|k-1} + K_k y_k + b_k \quad (10) $$

with $K_k \in \mathbb{R}^{n \times 1}$ the Kalman gain. If the pair $[F, C]$ is completely detectable and $[F, NG]$ is completely stabilizable for any $G$ such that $GG^\top = Q$, then (i) for any positive semidefinite initial matrix $P_0$, $P_{k+1|k}$ converges to the solution $P_\infty$ of the algebraic Riccati equation, which defines the associated gain $K_\infty = FP_\infty C^\top(CP_\infty C^\top + R)^{-1}$, (ii) all eigenvalues $\lambda$ of $F - K_\infty C$ verify $|\lambda| < 1$ [19]. The pair $[F, C]$ is completely detectable if $\exists K \in \mathbb{R}^3 : |\lambda| < 1, \forall \lambda \in \text{Sp}(F - KC)$ [19]. Let $K = (0, k_2, k_3)^\top$. For $k_2 = (1-\alpha)/4\alpha$, $k_3 = e^{-d\Delta t}$, $\text{Sp}(F - KC) = \{e^{-d\Delta t}, (1 \pm e^{-d\Delta t/2})/2\}$, which fulfills the condition for complete detectability. Similarly, the pair $[F, NG]$ is completely stabilizable if $\exists K \in \mathbb{R}^3 : |\lambda| < 1, \forall \lambda \in \text{Sp}(F - NGK^\top)$ [19]. Here $G = \sqrt{Q} > 0$. Taking $K = (0, \frac{1}{G(1-\alpha)\Delta t}, 0)^\top$, $\text{Sp}(F - NGK^\top) = \{e^{-d\Delta t}, 0\}$, which proves complete controllability. Hence $P_\infty$ and $K_\infty$ are well-defined.
From (10), we have $x_{k+1} = Fx_k + b_k$, therefore the estimation error for the asymptotic gain $K_\infty$ is

$$x_{k+1} - \hat{x}_{k+1|k} = (F - K_\infty C)(x_k - \hat{x}_{k|k-1}) - K_\infty v_k. \quad (11)$$

For $v_k = 0$, this converges to zero since $\|F - K_\infty C\| < 1$.

Finally, if $\exists k : x_k - \hat{x}_{k|k-1} = 0$, it follows from (11) that $x_{k+l} - \hat{x}_{k+l|k+l-1} = 0, \forall l > 0$. This concludes the proof.

## V. APPROACH FOR UNKNOWN $\mu(\cdot)$

The lack of knowledge of the specific growth rate function eliminates the possibility of classical approaches for state estimation. To remedy this issue, a Bayesian regularization approach is proposed for the online estimation of the reaction rate as an unknown function of time $r(t) = \mu(x(t), t)e(t)$. This is obtained by augmenting the system (2) with a probabilistic prior on $r(t)$ in the form of linear stochastic dynamics [16], [22], and directly solving the augmented linear estimation problem via the hybrid KF. Consider the following stochastic prior on $r$

$$\dot{r}(t) = -\theta r(t) + \kappa w(t) \quad (12)$$

where $\{w(t)\}$ is a standard Gaussian noise process. Together with the assumption that $r(0)$ is Gaussian with mean zero and variance $\kappa^2/(2\theta)$, this defines a stationary process. This convenience choice is discussed for a related application in [16]. In essence, $\theta$ defines the timescale of admissible fluctuations (the larger the $\theta$, the faster the fluctuations), and the standard deviation $\kappa/\sqrt{2\theta}$ defines the admissible magnitude of $r$. This already enables an educated guess for the choice of the design parameters $\theta$ and $\kappa$. Precise tuning of these parameters can be obtained from preliminary data via the ML procedure discussed in III-C.

For the ML-optimal values $\hat{\theta}$ and $\hat{\kappa}$, simultaneous online estimates of state $x$ and $r$ are obtained from measurements $y_k$ by the KF for the joint system (1)–(12),

$$\frac{d}{dt} \begin{bmatrix} x(t) \\ r(t) \end{bmatrix} = \begin{bmatrix} -dI & N \\ 0 & -\hat{\theta} \end{bmatrix} \begin{bmatrix} x(t) \\ r(t) \end{bmatrix} + \begin{bmatrix} dx_{in} \\ 0 \end{bmatrix} + \begin{bmatrix} 0 \\ \hat{\kappa} \end{bmatrix} w(t). \quad (13)$$

<figure><img src="image_1.png" /><figcaption>Fig. 1. Measurements and the reaction rate plotted against time. [Top] The fluorescent protein concentration $f(t)$ (black curve), noisy fluorescence measurements $y_k$ (green dots), and the pre-smoothed measurements $\hat{y}_k$ (red curve). [Bottom] The specific growth rate $\mu(s(t))$ and its pre-estimation via the conservation law $\mu_k$ (blue curves, solid and dotted, respectively), and the total growth rate $r(t) = \mu(s(t)e(t))$ and its BKF estimation $\hat{r}_k$ (orange curves, solid and dotted, respectively).</figcaption></figure>

Compared with the approach in Sec. IV, no further noise component appears in the dynamics of $x$. Performance of this approach will be demonstrated in simulation in Sec. VI.

## VI. SIMULATION RESULTS

In this section, numerical simulations of the proposed approaches are presented for system (2) with Monod growth rate function. The system is simulated over 36 hours, with measurements taken at a fixed interval of 5 minutes. Parameters are fixed arbitrarily but realistically based on the related system in [5]. We use $\gamma = 1$, $k_s = 0.2 gL^{-1}$, $\alpha = 0.3$, $s_{in} = 2 gL^{-1}$, the maximal growth rate $\mu_{max} = \ln 2 h^{-1}$ corresponds to a doubling time of 1 $h$, and we fix $d = 0.48 h^{-1}$. The initial conditions for the simulation are $x(0) = (0, 1, 0)^\top gL^{-1}$ whereas $\hat{x}(0) = (0.5, 1.5, 0)^\top gL^{-1}$.

For the approach presented in Sec. IV, $\hat{\sigma}_1(0) = 0.5 gL^{-1}$ is selected arbitrarily. The ML-optimal smoothing parameter is $\hat{\beta} \approx 0.013$. This method being based on the system's conservation law, we refer to it as the Conservation KF (CKF). As for the Bayesian regularization approach in Sec. V, we call it the Bayesian KF (BKF). Its optimal process parameters are $\hat{\theta} \approx 0.03$ and $\hat{\kappa} \approx 0.002$. The ML problems are solved through the optimizers provided in the Python library SciPy over numerically generated preliminary data, with the same initial conditions and the noise covariance $R$. Different data are later generated for running the filters.

<figure><img src="image_2.png" /><figcaption>Fig. 2. Comparison of KF for state estimation from fluorescence measurements. [Left] Reconstructed trajectories using CKF, BKF, and EKF using set of generated noisy data compared with the true state trajectory. [Right] Mean trajectory of estimated states for each filter over 500 data sets and a confidence interval given by twice the standard error around the mean trajectory.</figcaption></figure>
From (10), we have $x_{k+1} = Fx_k + b_k$, therefore the estimation error for the asymptotic gain $K_\infty$ is

$$x_{k+1} - \hat{x}_{k+1|k} = (F - K_\infty C)(x_k - \hat{x}_{k|k-1}) - K_\infty v_k. \quad (11)$$

For $v_k = 0$, this converges to zero since $\|F - K_\infty C\| < 1$.

Finally, if $\exists k: x_k - \hat{x}_{k|k-1} = 0$, it follows from (11) that $x_{k+l} - \hat{x}_{k+l|k+l-1} = 0, \forall l > 0$. This concludes the proof.

## V. APPROACH FOR UNKNOWN $\mu(\cdot)$

The lack of knowledge of the specific growth rate function eliminates the possibility of classical approaches for state estimation. To remedy this issue, a Bayesian regularization approach is proposed for the online estimation of the reaction rate as an unknown function of time $r(t) = \mu(x(t), t)e(t)$. This is obtained by augmenting the system (2) with a probabilistic prior on $r(t)$ in the form of linear stochastic dynamics [16], [22], and directly solving the augmented linear estimation problem via the hybrid KF. Consider the following stochastic prior on $r$

$$\dot{r}(t) = -\theta r(t) + \kappa w(t) \quad (12)$$

where $\{w(t)\}$ is a standard Gaussian noise process. Together with the assumption that $r(0)$ is Gaussian with mean zero and variance $\kappa^2/(2\theta)$, this defines a stationary process. This convenience choice is discussed for a related application in [16]. In essence, $\theta$ defines the timescale of admissible fluctuations (the larger the $\theta$, the faster the fluctuations), and the standard deviation $\kappa/\sqrt{2\theta}$ defines the admissible magnitude of $r$. This already enables an educated guess for the choice of the design parameters $\theta$ and $\kappa$. Precise tuning of these parameters can be obtained from preliminary data via the ML procedure discussed in III-C.

For the ML-optimal values $\hat{\theta}$ and $\hat{\kappa}$, simultaneous online estimates of state $x$ and $r$ are obtained from measurements $y_k$ by the KF for the joint system (1)–(12),

$$\frac{d}{dt} \begin{bmatrix} x(t) \\ r(t) \end{bmatrix} = \begin{bmatrix} -dI & N \\ 0 & -\hat{\theta} \end{bmatrix} \begin{bmatrix} x(t) \\ r(t) \end{bmatrix} + \begin{bmatrix} dx_{in} \\ 0 \end{bmatrix} + \begin{bmatrix} 0 \\ \hat{\kappa} \end{bmatrix} w(t). \quad (13)$$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Measurements and the reaction rate plotted against time. [Top] The fluorescent protein concentration $f(t)$ (black curve), noisy fluorescence measurements $y_k$ (green dots), and the pre-smoothed measurements $\hat{y}_k$ (red curve). [Bottom] The specific growth rate $\mu(s(t))$ and its pre-estimation via the conservation law $\mu_k$ (blue curves, solid and dotted, respectively), and the total growth rate $r(t) = \mu(s(t)e(t))$ and its BKF estimation $\hat{r}_k$ (orange curves, solid and dotted, respectively).</figcaption></figure>

Compared with the approach in Sec. IV, no further noise component appears in the dynamics of $x$. Performance of this approach will be demonstrated in simulation in Sec. VI.

## VI. SIMULATION RESULTS

In this section, numerical simulations of the proposed approaches are presented for system (2) with Monod growth rate function. The system is simulated over 36 hours, with measurements taken at a fixed interval of 5 minutes. Parameters are fixed arbitrarily but realistically based on the related system in [5]. We use $\gamma = 1, k_s = 0.2 \ gL^{-1}, \alpha = 0.3, s_{in} = 2 \ gL^{-1}$, the maximal growth rate $\mu_{max} = \ln 2 \ h^{-1}$ corresponds to a doubling time of 1 $h$, and we fix $d = 0.48 \ h^{-1}$. The initial conditions for the simulation are $x(0) = (0, 1, 0)^\top \ gL^{-1}$ whereas $\hat{x}(0) = (0.5, 1.5, 0)^\top \ gL^{-1}$.

For the approach presented in Sec. IV, $\hat{\sigma}_1(0) = 0.5 \ gL^{-1}$ is selected arbitrarily. The ML-optimal smoothing parameter is $\hat{\beta} \approx 0.013$. This method being based on the system's conservation law, we refer to it as the Conservation KF (CKF). As for the Bayesian regularization approach in Sec. V, we call it the Bayesian KF (BKF). Its optimal process parameters are $\hat{\theta} \approx 0.03$ and $\hat{\kappa} \approx 0.002$. The ML problems are solved through the optimizers provided in the Python library SciPy over numerically generated preliminary data, with the same initial conditions and the noise covariance $R$. Different data are later generated for running the filters.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Comparison of KF for state estimation from fluorescence measurements. [Left] Reconstructed trajectories using CKF, BKF, and EKF using set of generated noisy data compared with the true state trajectory. [Right] Mean trajectory of estimated states for each filter over 500 data sets and a confidence interval given by twice the standard error around the mean trajectory.</figcaption></figure>
The pre-estimates from the CKF pipeline are shown in Fig. 1. Estimate $\hat{\mu}_k$ converges to $\mu(s(t))$ at rate determined by $d$ after an onset resulting from the difference $\sigma_1(0) - \hat{\sigma}_1(0)$. The quality of the ML-estimates of $\beta$ (for CKF) and of $\theta$ and $\kappa$ (for BKF) is witnessed by the estimates $\hat{y}$, which duly provide a smoothened version of the data, and the performance of $\hat{r}$, which converges to $r$ after an onset. Reported in 1 in the interest of illustration, $\hat{r}$ is obtained in practice along with the estimation of the other state variables $(e, s, f)$, which is the focus of the next paragraph.

We then run each filter on a noisy dataset to get estimates $\hat{x}_k$ of $x(t)$. For comparison, we also run the well-known Extended KF (EKF). The $Q$ parameter for this filter is also fixed by solving the ML problem of subsection III-C on preliminary data, with 3 degrees of freedom, yielding $\hat{Q} \approx \text{diag}(0.04, 0.004, 0.00002)$. Estimation results on a single dataset are given in Fig. 2 (left). We also evaluate each filter statistically by Monte Carlo analysis, i.e. by estimating the same $x(t)$ from $N_{MC} = 500$ different noisy data sets with same noise intensity. A mean estimated trajectory is obtained for each filter as well as a 95% confidence band. Results from this statistical analysis are in Fig. 2 (right). Results show that all studied filters perform well generally. The estimated trajectories seem to converge to the true trajectory. However, both the CKF and BKF converge faster than EKF.

## VII. CONCLUSIONS

We proposed two online state estimators for a bacterial growth model in a bioreactor that includes synthesis of a fluorescent reporter protein. The estimators reconstruct the three-dimensional system's state from discrete-time noisy fluorescence measurements.

We analyzed the system observability showing the possibility to construct a tunable observer, which we achieved by our CKF. Using conservation laws, this solution pre-estimates the growth rate profile to turn state estimation into a linear time-varying problem. This allowed us to benefit from the KF convergence properties to also prove deterministic convergence in absence of noise. Our study focused on the Monod growth function, yet the CKF can be extended to other growth laws depending on both substrate and biomass concentration, such as the Contois law [1]. Indeed, for such a growth law, the a-priori unknown biomass concentration profile can be pre-estimated based on the second conservation relation in (4), obtaining $\hat{e}_\sigma$ from the fluorescence measurements sided by a detector for the conservation variable. In more generality, the approach is viable for any growth model depending on state variables that enter suitable conservation laws.

We also proposed what we called BKF as a solution that does not assume knowledge of the growth rate law, a relevant scenario for bioprocesses. The approach allows estimating the system's state efficiently when classical approaches are not applicable. By its nature, the BKF can be applied to a more general class of growth rate functions, provided a Bayesian prior can be fixed to describe the system's reaction rate behavior over time.

## REFERENCES

[1] G. Bastin and D. Dochain, Eds., *On-Line Estimation and Adaptive Control of Bioreactors*, ser. Process Measurement and Control. Amsterdam: Elsevier, Jan. 1990.

[2] F. Clementschitsch and K. Bayer, "Improvement of bioprocess monitoring: Development of novel concepts," *Microbial Cell Factories*, vol. 5, no. 1, p. 19, May 2006.

[3] S. Ben Said and D. Or, "Synthetic Microbial Ecology: Engineering Habitats for Modular Consortia," *Frontiers in microbiology*, vol. 8, Jun. 2017.

[4] A. dos Reis de Souza, D. Efimov, A. Polyakov, J.-L. Gouzé, and E. Cinquemani, "State observation in microbial consortia: A case study on a synthetic producer-cleaner consortium," *International Journal of Robust and Nonlinear Control*, vol. 33, no. 9, pp. 5011-5022, 2021.

[5] M. Mauri, J.-L. Gouzé, H. de Jong, and E. Cinquemani, "Enhanced production of heterologous proteins by a synthetic microbial community: Conditions and trade-offs," *PLOS Computational Biology*, vol. 16, no. 4, p. e1007795, Apr. 2020.

[6] F. Bertaux, S. Sosa-Carrillo, V. Gross, A. Fraisse, C. Aditya, M. Furstenheim, and G. Batt, "Enhancing bioreactor arrays for automated measurements and reactive control with ReacSight," *Nature Communications*, vol. 13, no. 1, p. 3363, Jun. 2022.

[7] G. Besançon, *Nonlinear Observers and Applications*, ser. Lecture Notes in Control and Information Sciences. Berlin, Heidelberg: Springer, 2007, vol. 363.

[8] J.-P. Gauthier and I. Kupka, "Deterministic Observation Theory and Applications." Cambridge University Press, Oct. 2001.

[9] R. V. Hogg, "An Introduction to Robust Estimation," in *Robustness in Statistics*, R. L. Launer and G. N. Wilkinson, Eds. Academic Press, Jan. 1979, pp. 1-17.

[10] A. Jazwinski, *Stochastic Processes and Filtering Theory*, 1st ed. Academic Press, Jan. 1970.

[11] S. Sarkka, "On Unscented Kalman Filtering for State Estimation of Continuous-Time Nonlinear Systems," *IEEE Transactions on Automatic Control*, vol. 52, no. 9, pp. 1631-1641, Sep. 2007.

[12] H. L. Smith and P. Waltman, *The Theory of the Chemostat: Dynamics of Microbial Competition*, ser. Cambridge Studies in Mathematical Biology. Cambridge: Cambridge University Press, 1995.

[13] R. Neeleman and A. Boxtel, "Estimation of specific growth rate from cell density measurements," *Bioprocess and Biosystems Engineering*, vol. 24, no. 3, pp. 179-185, Oct. 2001.

[14] I. Didi, A. Rapaport, H. Dib, and B. Cherki, "A positive observer for the chemostat model with biogas measurement," *Journal of Process Control*, vol. 122, pp. 49-58, Feb. 2023.

[15] H. de Jong, C. Ranquet, D. Ropers, C. Pinel, and J. Geiselmann, "Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria," *BMC Systems Biology*, vol. 4, no. 1, p. 55, Apr. 2010.

[16] A. Pavlou, E. Cinquemani, J. Geiselmann, and H. de Jong, "Maturation models of fluorescent proteins are necessary for unbiased estimates of promoter activity," *Biophysical Journal*, vol. 121, no. 21, pp. 4179-4188, Nov. 2022.

[17] M. Fjeld, O. A. Asbjørnsen, and K. J. Åström, "Reaction invariants and their importance in the analysis of eigenvectors, state observability and controllability of the continuous stirred tank reactor," *Chemical Engineering Science*, vol. 29, no. 9, pp. 1917-1926, Sep. 1974.

[18] R. Heinrich and S. Schuster, "Stoichiometric Analysis," in *The Regulation of Cellular Systems*, R. Heinrich and S. Schuster, Eds. Boston, MA: Springer US, 1996, pp. 75-111.

[19] B. Anderson and J. Moore, *Optimal Filtering*, ser. Dover Books on Electrical Engineering. Dover Publications, 1979.

[20] H. Lütkepohl, "State Space Models," in *Introduction to Multiple Time Series Analysis*, H. Lütkepohl, Ed. Berlin, Heidelberg: Springer, 1991, pp. 415-445.

[21] G. De Nicolao, G. Sparacino, and C. Cobelli, "Nonparametric input estimation in physiological systems: Problems, methods, and case studies," *Automatica*, vol. 33, no. 5, pp. 851-870, May 1997.

[22] G. De Nicolao and G. Ferrari-Trecate, "Regularization networks for inverse problems: A state-space approach," *Automatica*, vol. 39, no. 4, pp. 669-676, Apr. 2003.
The pre-estimates from the CKF pipeline are shown in Fig. 1. Estimate $\hat{\mu}_k$ converges to $\mu(s(t))$ at rate determined by $d$ after an onset resulting from the difference $\sigma_1(0) - \hat{\sigma}_1(0)$. The quality of the ML-estimates of $\beta$ (for CKF) and of $\theta$ and $\kappa$ (for BKF) is witnessed by the estimates $\hat{y}$, which duly provide a smoothened version of the data, and the performance of $\hat{r}$, which converges to $r$ after an onset. Reported in 1 in the interest of illustration, $\hat{r}$ is obtained in practice along with the estimation of the other state variables $(e, s, f)$, which is the focus of the next paragraph.

We then run each filter on a noisy dataset to get estimates $\hat{x}_k$ of $x(t)$. For comparison, we also run the well-known Extended KF (EKF). The $Q$ parameter for this filter is also fixed by solving the ML problem of subsection III-C on preliminary data, with 3 degrees of freedom, yielding $\hat{Q} \approx \text{diag}(0.04, 0.004, 0.00002)$. Estimation results on a single dataset are given in Fig. 2 (left). We also evaluate each filter statistically by Monte Carlo analysis, i.e. by estimating the same $x(t)$ from $N_{MC} = 500$ different noisy data sets with same noise intensity. A mean estimated trajectory is obtained for each filter as well as a 95% confidence band. Results from this statistical analysis are in Fig. 2 (right). Results show that all studied filters perform well generally. The estimated trajectories seem to converge to the true trajectory. However, both the CKF and BKF converge faster than EKF.

## VII. CONCLUSIONS

We proposed two online state estimators for a bacterial growth model in a bioreactor that includes synthesis of a fluorescent reporter protein. The estimators reconstruct the three-dimensional system's state from discrete-time noisy fluorescence measurements.

We analyzed the system observability showing the possibility to construct a tunable observer, which we achieved by our CKF. Using conservation laws, this solution pre-estimates the growth rate profile to turn state estimation into a linear time-varying problem. This allowed us to benefit from the KF convergence properties to also prove deterministic convergence in absence of noise. Our study focused on the Monod growth function, yet the CKF can be extended to other growth laws depending on both substrate and biomass concentration, such as the Contois law [1]. Indeed, for such a growth law, the a-priori unknown biomass concentration profile can be pre-estimated based on the second conservation relation in (4), obtaining $\hat{e}_\sigma$ from the fluorescence measurements sided by a detector for the conservation variable. In more generality, the approach is viable for any growth model depending on state variables that enter suitable conservation laws.

We also proposed what we called BKF as a solution that does not assume knowledge of the growth rate law, a relevant scenario for bioprocesses. The approach allows estimating the system's state efficiently when classical approaches are not applicable. By its nature, the BKF can be applied to a more general class of growth rate functions, provided a Bayesian prior can be fixed to describe the system's reaction rate behavior over time.

## REFERENCES

[1] G. Bastin and D. Dochain, Eds., *On-Line Estimation and Adaptive Control of Bioreactors*, ser. Process Measurement and Control. Amsterdam: Elsevier, Jan. 1990.

[2] F. Clementschitsch and K. Bayer, "Improvement of bioprocess monitoring: Development of novel concepts," *Microbial Cell Factories*, vol. 5, no. 1, p. 19, May 2006.

[3] S. Ben Said and D. Or, "Synthetic Microbial Ecology: Engineering Habitats for Modular Consortia," *Frontiers in microbiology*, vol. 8, Jun. 2017.

[4] A. dos Reis de Souza, D. Efimov, A. Polyakov, J.-L. Gouzé, and E. Cinquemani, "State observation in microbial consortia: A case study on a synthetic producer-cleaner consortium," *International Journal of Robust and Nonlinear Control*, vol. 33, no. 9, pp. 5011-5022, 2021.

[5] M. Mauri, J.-L. Gouzé, H. de Jong, and E. Cinquemani, "Enhanced production of heterologous proteins by a synthetic microbial community: Conditions and trade-offs," *PLOS Computational Biology*, vol. 16, no. 4, p. e1007795, Apr. 2020.

[6] F. Bertaux, S. Sosa-Carrillo, V. Gross, A. Fraisse, C. Aditya, M. Furstenheim, and G. Batt, "Enhancing bioreactor arrays for automated measurements and reactive control with ReacSight," *Nature Communications*, vol. 13, no. 1, p. 3363, Jun. 2022.

[7] G. Besançon, *Nonlinear Observers and Applications*, ser. Lecture Notes in Control and Information Sciences. Berlin, Heidelberg: Springer, 2007, vol. 363.

[8] J.-P. Gauthier and I. Kupka, "Deterministic Observation Theory and Applications." Cambridge University Press, Oct. 2001.

[9] R. V. Hogg, "An Introduction to Robust Estimation," in *Robustness in Statistics*, R. L. Launer and G. N. Wilkinson, Eds. Academic Press, Jan. 1979, pp. 1-17.

[10] A. Jazwinski, *Stochastic Processes and Filtering Theory*, 1st ed. Academic Press, Jan. 1970.

[11] S. Sarkka, "On Unscented Kalman Filtering for State Estimation of Continuous-Time Nonlinear Systems," *IEEE Transactions on Automatic Control*, vol. 52, no. 9, pp. 1631-1641, Sep. 2007.

[12] H. L. Smith and P. Waltman, *The Theory of the Chemostat: Dynamics of Microbial Competition*, ser. Cambridge Studies in Mathematical Biology. Cambridge: Cambridge University Press, 1995.

[13] R. Neeleman and A. Boxtel, "Estimation of specific growth rate from cell density measurements," *Bioprocess and Biosystems Engineering*, vol. 24, no. 3, pp. 179-185, Oct. 2001.

[14] I. Didi, A. Rapaport, H. Dib, and B. Cherki, "A positive observer for the chemostat model with biogas measurement," *Journal of Process Control*, vol. 122, pp. 49-58, Feb. 2023.

[15] H. de Jong, C. Ranquet, D. Ropers, C. Pinel, and J. Geiselmann, "Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria," *BMC Systems Biology*, vol. 4, no. 1, p. 55, Apr. 2010.

[16] A. Pavlou, E. Cinquemani, J. Geiselmann, and H. de Jong, "Maturation models of fluorescent proteins are necessary for unbiased estimates of promoter activity," *Biophysical Journal*, vol. 121, no. 21, pp. 4179-4188, Nov. 2022.

[17] M. Fjeld, O. A. Asbjørnsen, and K. J. Åström, "Reaction invariants and their importance in the analysis of eigenvectors, state observability and controllability of the continuous stirred tank reactor," *Chemical Engineering Science*, vol. 29, no. 9, pp. 1917-1926, Sep. 1974.

[18] R. Heinrich and S. Schuster, "Stoichiometric Analysis," in *The Regulation of Cellular Systems*, R. Heinrich and S. Schuster, Eds. Boston, MA: Springer US, 1996, pp. 75-111.

[19] B. Anderson and J. Moore, *Optimal Filtering*, ser. Dover Books on Electrical Engineering. Dover Publications, 1979.

[20] H. Lütkepohl, "State Space Models," in *Introduction to Multiple Time Series Analysis*, H. Lütkepohl, Ed. Berlin, Heidelberg: Springer, 1991, pp. 415-445.

[21] G. De Nicolao, G. Sparacino, and C. Cobelli, "Nonparametric input estimation in physiological systems: Problems, methods, and case studies," *Automatica*, vol. 33, no. 5, pp. 851-870, May 1997.

[22] G. De Nicolao and G. Ferrari-Trecate, "Regularization networks for inverse problems: A state-space approach," *Automatica*, vol. 39, no. 4, pp. 669-676, Apr. 2003.