Kalman-based approaches for online estimation of bioreactor dynamics from fluorescent reporter measurements
Rand Asswad Eugenio Cinquemani Jean-Luc Gouzé


An obstruction theory for strictly commutative algebras in positive characteristic
    Oisín Flynn-Connolly
    April 29, 2024
==================================================================================

empty
empty





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.





§ INTRODUCTION


Estimation and control of cellular growth in bioreactors have been dedicated a significant amount of work over several decades <cit.>.
Nowadays, these challenges receive renewed attention
in relation with new biotechnological
developments. Among the new research frontiers
is the feedback control of microbial consortia <cit.>. Motivated by the potential to outperform single species in biosynthetic
processes and other applications <cit.>, 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 <cit.>, we have explored 
(local) observability of the prototypical microbial consortium model from <cit.> 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 <cit.>)
have not been explored in detail. Given the nonlinear nature of these models, the problem is
nontrivial. General observer synthesis approaches exist 
<cit.>, 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 <cit.>, but they are generally not conceived for sampled data. Nonlinear generalizations of Kalman filtering exist that cope explicitly
with sampled noisy data <cit.>, 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 <cit.>)
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
<cit.>
or gaseous outflow <cit.>.
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. <ref>, we describe bioreactor growth models and illustrate the properties of our model. In Sec. <ref> we review concepts of observability
and Kalman filtering estimation that are used in the sequel.
In Sec. <ref> we introduce our estimation method 
for a perfectly known growth law,
and prove its convergence properties.
In Sec. <ref> we present our method to
address state estimation for unknown growth laws.
In Sec. <ref>, 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. <ref>.



§ 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

    = N r(x,t) + d( - x)

where
x(t)∈Ω=_+^n represents the reactants concentrations vector,
∈_+^n is the reactor feed concentrations vector,
N∈^n× q is the stoichiometric matrix,
r(x,t)∈_+^q is the reaction rates vector (),
and d∈_+^* is the dilution rate () <cit.>.
Concentrations (in ) are relative 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 <cit.> 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))∈_+^3, the model becomes <cit.>

    =    -1/γ   μ(x,t)e +d( - s)
       =    (1-α)   μ(x,t)e -de
       =    α   μ(x,t)e -df

where μ(x,t)∈_+, the growth rate per capita
at time t≥ t_0, 
may depend on time due to quantities that are not part of the model state (pH, temperature) <cit.>.
Constant 0<α<1 is the proportion of substrate import dedicated to
fluorescent protein synthesis, while  represents the inflow substrate concentration.
Overall, this is
a constant yield model where the substrate
uptake rate per capita is given as
μ(x,t)/γ with γ the 
growth yield coefficient.
The values of α, , and γ are assumed known and strictly
positive.
The model is in the form of (<ref>)
with N=(-1/γ,1-α,α)∈^3,
r(x,t)=μ(x,t)e(t)∈_+, and =(,0,0)∈_+^3.

For the problem of online state observation, in section <ref>
we consider the specific growth rate known and strictly a function of the
system's state.
In section <ref>, we will instead address the case where the functional
form of μ(·) is unknown, and thus treat it as an unknown function of
time.



 §.§ 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 ℓ = n - (N)>0,
which corresponds to the number of conservation relations.
A conservation matrix Π∈^ℓ× n is defined as a matrix 
with rows forming a basis of the left-null space of N.
Consequently, Π N=0 and (Π) = ℓ.
The conservation matrix, which is not unique
<cit.>,
projects the system's state into an auxiliary
state σ = Π x of dimension ℓ.
We also refer to this definition as the conservation law.
From (<ref>),

    = d(Π - σ)

(in the “batch mode” d≡0, =0, whence the name “conservation”).
The conservation manifold (nullcline of the conservation law) is defined as Σ=x∈^n,Π=0.
From (<ref>), 
σ(t) = Π + e^-d(t-t_0)(σ(t_0) - Π),
which shows that 
Σ is positively invariant.
Since d>0, Σ is also exponentially attracting, thus it  
is asymptotic invariant <cit.>.

System (<ref>) has dimension n=3  and (N)=1. Hence
there are ℓ=2 conservation relations. Notably, we choose

    Π = 
    1     0    1αγ
    
    0     1     -1-αα{σ_1    = s + 1/αγf,
    σ_2    = e - 1-α/αf.
    .




 §.§ Steady states


Regardless of the form of μ(x,t), system
(<ref>) always admits an equilibrium
x_washout = (,0,0), which is globally asymptotically
stable on Ω if the set of equilibria is a singleton and unstable
otherwise.
The existence and stability of other equilibria depends on μ(x,t)
and on the system parameters, in this case d, α, and .

We focus on the  Monod law <cit.>. In this case μ
takes the form

    μ(s) = s/k_s + s,    s≥ 0,

with k_s>0 known as half-saturation constant and
>0 the maximal growth rate.
Eq. (<ref>) is strictly increasing,
therefore it defines a bijection from [0,∞) to [0,).
If d/(1-α) < μ(), (<ref>) admits another equilibrium
x^*=(s^*, γ(1-α)(-s^*), γα(-s^*))
with s^*=μ^-1d/(1-α).
If this exists, it is globally asymptotically stable
on any open subset of Ω∖e=0.
From now on we work under this assumption.
These results extend similar known results for the subsystem
(s,e) without fluorescence synthesis <cit.>.
They apply similarly to several other common growth models,
such as Teissier law or Contois law <cit.>.



§ 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)∈Ω⊂^n. Typically, information about the state variables
is obtained through measurements y(t)=h(x(t))∈^p
(often referred to as the system's output),
where h is a function assumed ^∞(Ω).
The problem of interest in this paper is the construction of a state
observer (t)∈^n from the measurements y(t) for all t≥ t_0.
(t) is an observer if
(i) x(t)-(t) converges to zero as t tends to ∞,
and (ii) (t)=x(t)(t+τ)=x(t+τ),∀τ≥0.
If the first condition holds for any x(t_0) and (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 <ref> <cit.>.
An open loop (non-tunable) observer (i.e. with no output feedback)
is also referred to as a detector.

The available measurements for (<ref>) are
fluorescent protein concentrations y(t)=f(t), thus
h is defined by 
y(t)=C x(t) with C=(0,0,1). 



 §.§ Observability


The local observability of a dynamical system =f(x) is studied
via its observation space (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 ϕ:^n→^p
is L_f(ϕ) = ϕxf.
The system is locally observable if and only if there exists n elements
of (h) defining a diffeomorphism around x for all x∈Ω
<cit.>.
By the chain rule, (t)=hx(x(t))(t)=L_f(h)(x(t)).
Ergo, the observation space is the space of time-derivatives of the outputs.

For (<ref>), y=f, therefore = αμ(x,t)e - d· y.
Setting z = (d· y + )/α = μ(x_g,t)e defines a reduced
dynamical system for the state x_g=(s,e) with 
output
z=h_g(x_g)=μ(x_g,t)e, provided μ(x,t) does not involve f.
Consequently, (h_g)⊂(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=μ(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 <cit.>.
Observability of (<ref>) for the Monod growth law follows from
the results in <cit.>.



 §.§ 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. <cit.>)), 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<⋯<t_M-1<t_M 
be given by
y_k = C x(t_k) + v_k,
where v_k is an independent Gaussian noise with zero mean and known
covariance matrix R_k∈^p× 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 <cit.>.
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)
<cit.>,
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. <ref>) or described
by a probabilistic prior in the form of a conveniently chosen linear stochastic
process (as explained in Sec. <ref>).
In both cases, this leads to a hybrid filtering problem on a linear
(time-varying) stochastic system of the form

    (t)    = A(t)x(t) + b(t) + G(t) w(t),
    
        y_k    = C_k x(t_k) + v_k,

with 
A(t)∈^n× n,b(t)∈^n,C_k∈^p× n,G(t)∈^n× m,
Q(t)∈^m× m and R_k∈^p× p known. The process noise 
{w(t),t≥ t_0} is formally the derivative of the Brownian motion
<cit.>, i.e. it is a zero-mean, stochastic process with
w(t)w(τ)=Q(t)δ(t-τ), δ 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 
the system's dynamics for the state estimate _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 _k|k-1(t_k-1)=_k-1|k-1 and
P_k|k-1(t_k-1)=P_k-1|k-1.
The update step applied to _k|k-1(t_k) and
P_k|k-1(t_k) for a new measurement y_k  
is identical to the discrete KF <cit.>. 



 §.§ Parameter tuning from data


In our estimator designs, Eq. (<ref>) will comprise unknown parameters
θ∈Θ. We propose (and later demonstrate in simulation) to estimate
these parameters from a preliminary data set Y_M=y_0,…,y_M
through a Maximum Likelihood (ML) approach,
which closely relates with the filter's prediction error methods
for tuning filters <cit.>.

Given Y_M, ML estimation is the problem of maximizing
the value of the joint distribution
f(Y_M|θ) over θ∈Θ.
From the Bayes law,
f(y_M, Y_M-1|θ) = f(y_M|Y_M-1,θ)f(Y_M-1|θ).
By recursive application one can write
f(Y_M|θ) as the product of the conditional
(Gaussian) densities f(y_k|Y_k-1,θ) for k=0,…,M.
Mean and covariance of these densities correspond to the Kalman prediction
C_k_k|k-1 and its covariance S_k=C_k P_k|k-1C_k + R_k
<cit.>.
Writing ML as the equivalent negative log-likelihood minimization and working
out the equations, the ML estimator of θ is

    = _θ∈Θ f(Y_M|θ)
       = _θ∈Θ-∑_k=0^Mlog f(y_k|Y_k-1,θ)
       = _θ∈Θ-1/2∑_k=0^M plog(2π) + logS_k
                + _kS_k^-1_k.

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



§ APPROACH FOR KNOWN Μ(X)
 

In this section, leveraging conservation laws, we propose a pipeline for state
estimation for the case where μ is a known function of the state variables.
We develop this for system (<ref>) 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 <ref>, nonlinearities of CSTR dynamics are
in the reaction rate vector r(x,t) of (<ref>).
For system (<ref>), it is however possible to compute a pre-estimate
(t) of the specific growth rate μ(x,t).
Replacing μ(x,t) in (<ref>) with (t) yields linear
time-varying dynamics,
eligible for linear Kalman filtering.
We next show how to compute (t) in the case of the
Monod growth rate μ(s).



 §.§ Proposed pipeline


The Monod law is a function of the unknown s(t).
By the conservation law (<ref>),
in the (exponentially attracting) conservation manifold, s=σ_1 - 1/αγf.
In view of this, we define a pre-estimate _σ of s by _σ=σ̂_1 - 1/αγŷ, with σ̂_1 calculated based on the dynamics of σ_1 and ŷ a suitable smoothening of the measurements y_k of f, as detailed next.

The conservation variable σ_1 is estimated through an open-loop detector
_1(t) with dynamics as in (<ref>), resulting in 

    _1(t) =  + e^-d(t-t_0)(_1(t_0) - ).

This converges to  exponentially at a rate d.

To calculate , a Kalman-based regularization method is used. 
In the same spirit as <cit.>,  
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 <cit.>.
We consider a second-order linear stochastic differential equation
ÿ(t) = β 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
(<ref>) as

    y(t)
    (t) =
        0     1
    0     0y(t)
    (t) + 0
    βw(t)

with C=(1,0).
An optimal parameter β is the solution to the ML problem discussed
in subsection <ref>. 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 measurments y_k yields estimates of y, i.e. of f. The chosen form of process y 
is known to penalize overly irregular solutions <cit.>, whereas the ML estimation of β ensures an appropriate tradeoff between noisy data interpolation and regularity of the solution. The filter produces 
piecewise-continuous predictions
_k|k-1(t) and filtered estimates _k|k. 
Both choices being equally viable as a definition of the estimate (t) sought,
we simply refer to (t) from now on and add details only when necessary.

Armed with _1 and , one can now compute 
(t)=μ(_σ(t)) 
with _σ(t) = _1(t) - 1/αγ(t),
and plug it into (<ref>) to approach Kalman filtering estimation of the whole system state.
Let (t) = μ(s(t)) - (t) be the pre-estimation error
of μ at time t. Then

    = N((t) + (t))e + d( - x)
       = A_σ(t)x + b + N w(t)

with w(t)=(t)e(t) as the process error term, where 

    A_σ(t) = 
        -d     -(t)/γ    0
    
        0      (1-α)(t) -d     0
    
        0     α(t)            -d
    , b=d·
    0
    0

are now known.
To run the linear Kalman filter of subsection <ref> for this system
(with G=N and C=(0,0,1) in Eq. (<ref>)), we exploit the analytic
expression of _1(t) in (<ref>), for a given _1(t_0),
while  and  can be obtained online sequentially
for every interval =[t_k,t_k+1) for k=0,…,M-1,
as shown in Figure <ref>.


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.





 §.§ Tuning of Q(t)


We now derive an approximate expression for Q(t) that characterizes
w(t)=(t)e(t).
Consider the linear approximation of μ around s(t).
For sufficiently small _σ=s - _σ, 

    μ(_σ)=μ(s - _σ)≈μ(s) - μ'(s)_σ.

Then (t) = μ'(s(t))_σ(t) implies 
w(t)=μ'(s(t))e(t)_σ(t).
We choose to study the process noise around the system equilibrium.
From (<ref>), (<ref>) and subsection <ref>,
it can be shown that for t∈, _σ(t) = v_k/αγ
at the system's steady states. Consequently,
w(t) = μ'(s^*)e^*/αγ v_k for t∈,
near the system's equilibria.
It follows that w(t) is a piecewise constant stochastic process.
By integrating w(t)w(τ) = Q(t)δ(t-τ) over 
we obtain the piecewise-constant function
Q(t)=μ'(s^*)e^*/αγ^2 R_k (t_k+1 - t_k) for t∈.



 §.§ Asymptotic stability of the filter


We now prove that the proposed filtering pipeline defines a global observer in the sense of Section <ref>. That is, we consider 
noiseless measurements and look at convergence at the discrete measurement times <cit.>.
We first prove convergence of the observer error, then
its positive invariance, i.e.
that x_k - _k = 0 implies x_k+l - _k+l = 0 for l>0.

As seen before, the KF is applied to the pre-estimated system
(<ref>) 
with A_σ(t)x=N(t)e-dx.
In the absence of measurement noise, _σ(t) converges to s^*
at the same rate as _1 (globally, exponentially at rate d).
Consequently, (t)=μ(_σ(t)) converges to
μ^*=μ(s^*)=d/(1-α)
and the system (<ref>) is asymptotically time-invariant
with A_σ(t) converging globally to

    A = 
        -d     -μ^*/γ    0
    
        0     0     0
    
        0    αμ^*     -d
    .

By integrating the linear time-invariant dynamics over ,

    (t)    = A x(t) + b + N w(t)
    
        x_k+1   = e^A(t_k+1-t_k) x_k + ∫_ e^A(t_k+1-t) (b + N w(t)).

One can verify that here, e^AτN=N,∀τ.
We obtain the difference equation, for all k, 

    x_k+1 = F_k x_k + b_k + N w_k

with F_k=exp(A(t_k+1-t_k)) the state transition matrix <cit.>,
b_k=∫_exp(A(t_k+1-t_k))b,
and w_k=∫_ w(t).
Let us fix a time step Δ t=t_k+1-t_k,∀ k. Let F=F_k.
For a fixed Q_c for the continuous system, w_k w_l=Q δ_k,l
with Q = Q_cΔ t. 

The equivalent discrete Kalman equation is

    _k+1|k = (F - K_k C) _k|k-1 + K_k y_k + b_k

with K_k∈^n× 1 the Kalman gain.
If the pair [F,C] is completely detectable and [F,NG]
is completely stabilizable for any G such that GG=Q,
then (i) for any positive semidefinite initial matrix P_0,
P_k+1|k converges to the solution P_∞ of the algebraic
Riccati equation, which defines the associated gain
K_∞ = F P_∞C(C P_∞C + R)^-1,
(ii) all eigenvalues λ of F-K_∞ C verify λ<1
<cit.>. The pair [F,C] is completely detectable if
∃ K∈^3: λ<1, ∀λ∈(F-KC) <cit.>.
Let K=(0, k_2, k_3). For
k_2=(1-α)/4α,k_3=e^-dΔ t, (F-KC)=e^-dΔ t, 1± e^-dΔ t/2/2,
which fulfills the condition for complete detectability.
Similarly, the pair [F,NG] is completely stabilizable if 
∃ K∈^3: λ<1, 
∀λ∈(F-NGK) <cit.>.
Here G=√(Q)>0. Taking K=(0,1/G(1-α)Δ t,0),
(F-NGK)=e^-dΔ t,0, which proves complete
controllability. Hence P_∞ and K_∞ are well-defined.

From (<ref>), we have x_k+1=F x_k + b_k, therefore 
the estimation error for the asymptotic gain K_∞ is

    x_k+1 - _k+1|k = (F - K_∞ C)(x_k - _k|k-1) - K_∞ v_k.

For v_k=0, this converges to zero since F-K_∞ C<1.

Finally, if ∃ k: x_k-_k|k-1=0, it follows from
(<ref>) that x_k+l-_k+l|k+l-1=0,∀ l>0.
This concludes the proof.



§ APPROACH FOR UNKNOWN Μ(·)


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)=μ(x(t),t)e(t).
This is obtained by augmenting the system (<ref>) with a
probabilistic prior on r(t) in the form of linear stochastic dynamics
<cit.>,
and directly solving the augmented linear estimation problem via the hybrid KF.
Consider the following stochastic prior on r

    (t) = -θ r(t) + κ w(t)

where w(t) is a standard Gaussian noise process.
Together with the assumption that r(0) is Gaussian with mean zero and
variance κ^2/(2θ), this defines a stationary process.
This convenience choice is discussed for a related application in
<cit.>.
In essence, θ defines the timescale of admissible fluctuations
(the larger the θ, the faster the fluctuations),
and the standard deviation κ/√(2θ) defines the admissible
magnitude of r. This already enables an educated guess for the choice of the
design parameters θ and κ.
Precise tuning of these parameters can be obtained from preliminary data via
the ML procedure discussed in <ref>.

For the ML-optimal values  and , simultaneous online estimates
of state x and r are obtained from measurements y_k by the KF
for the joint system (<ref>)–(<ref>),

    x(t)
    r(t) =
        -dI     N
    0     -x(t)
    r(t) + d
    0
        + 0
    w(t).




Compared with the approach in Sec. <ref>, no further noise
component appears in the dynamics of x. Performance of this approach
will be demonstrated in simulation in Sec. <ref>.



§ SIMULATION RESULTS


In this section, numerical simulations of the proposed approaches are presented
for system (<ref>) 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 <cit.>.
We use γ=1, k_s=0.2 , α=0.3, =2, the maximal growth
rate =ln2 corresponds to a doubling time of 1 h, and 
we fix d=0.48.
The initial conditions for the simulation are
x(0)=(0,1,0) whereas (0)=(0.5,1.5,0).

For the approach presented in Sec. <ref>, _1(0)=0.5
is selected arbitrarily. The ML-optimal smoothing parameter is
β≈ 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. <ref>,
we call it the Bayesian KF (BKF). Its optimal process parameters are
≈0.03 and ≈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.

The pre-estimates from the CKF pipeline are shown in Fig. <ref>.
Estimate _k converges to μ(s(t)) at rate determined by d after an
onset resulting from the difference σ_1(0)-_1(0).
The quality of the ML-estimates of β (for CKF) and of θ and
κ (for BKF) is witnessed by the estimates , which duly provide
a smoothened version of the data, and the performance of r,
which converges to r after an onset. 
Reported in <ref> in the interest of illustration, 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 _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 <ref> on preliminary
data, with 3 degrees of freedom, yielding
Q≈diag(0.04, 0.004, 0.00002).
Estimation results on a single dataset are given in Fig. <ref> (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. <ref> (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.



§ 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 <cit.>.
Indeed, for such a growth law, the a-priori unknown biomass concentration
profile can be pre-estimated based on the second conservation relation in
(<ref>), obtaining _σ 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.




IEEEtran

