# Hybrid physics-informed metabolic cybergenetics: process rates augmented
with machine-learning surrogates informed by flux balance analysis

Sebastián Espinel-Ríos<EMAIL_ADDRESS>Department of Chemical and
Biological Engineering, Princeton University, Princeton NJ, USA  José L.
Avalos<EMAIL_ADDRESS>Department of Chemical and Biological
Engineering, Princeton University, Princeton NJ, USA

###### Abstract

Metabolic cybergenetics is a promising concept that interfaces gene expression
and cellular metabolism with computers for real-time dynamic metabolic
control. The focus is on control at the transcriptional level, serving as a
means to modulate intracellular metabolic fluxes. Recent strategies in this
field have employed constraint-based dynamic models for process optimization,
control, and estimation. However, this results in bilevel dynamic optimization
problems, which pose considerable numerical and conceptual challenges. In this
study, we present an alternative hybrid physics-informed dynamic modeling
framework for metabolic cybergenetics, aimed at simplifying optimization,
control, and estimation tasks. By utilizing machine-learning surrogates, our
approach effectively embeds the physics of metabolic networks into the process
rates of structurally simpler macro-kinetic models coupled with gene
expression. These surrogates, informed by flux balance analysis, link the
domains of manipulatable intracellular enzymes to metabolic exchange fluxes.
This ensures that critical knowledge captured by the system’s metabolic
network is preserved. The resulting models can be integrated into metabolic
cybergenetic schemes involving single-level optimizations. Additionally, the
hybrid modeling approach maintains the number of system states at a necessary
minimum, easing the burden of process monitoring and estimation. Our hybrid
physics-informed metabolic cybergenetic framework is demonstrated using a
computational case study on the optogenetically-assisted production of
itaconate by Escherichia coli.

###### keywords:

cybergenetics, metabolic engineering, machine learning, predictive control,
estimation, optogenetics

Omenn-Darling Bioengineering Institute, Princeton University, Princeton NJ,
USA The Andlinger Center for Energy and the Environment, Princeton NJ, USA
High Meadows Environmental Institute, Princeton University, Princeton NJ, USA

##  1 Introduction

Microbial biotechnology seeks to exploit the metabolic machinery of cells
through genetic and metabolic engineering to produce valuable products from
renewable resources 1, 2. In this context, dynamic metabolic engineering
strategies can be implemented to enhance the production efficiency of
(heterologous) metabolites and proteins 3, 4. Compared to static metabolic
engineering, dynamic metabolic engineering aims to adjust the cell’s metabolic
flux distribution over time, offering more optimization degrees of freedom and
flexibility for effectively addressing intrinsic metabolic trade-offs in
bioprocesses.

A typical metabolic trade-off occurs when an increase in product yield diverts
resources from biomass synthesis, consequently reducing biomass yields and
volumetric productivity rates 5. To address this issue, dynamic transitions
between metabolic states, such as those from growth to production, can lead to
improved volumetric productivity in bioprocesses 6, 3, 4, 5. Dynamic metabolic
engineering can also be considered for modulating the flux through production
pathways toward achieving a desired product composition 7.

There are different strategies to enable dynamic metabolic engineering 6, 3,
4, 5. At the transcriptional level, the most widely implemented strategy, the
gene transcription of key metabolic enzymes is modulated over time. With this,
the intracellular concentration of these enzymes is adjusted dynamically,
influencing the achievable metabolic fluxes in the cell. In this study, such
key enzymes are referred to as manipulatable enzymes, and the metabolic fluxes
they catalyze as manipulatable fluxes. At the post-translational level,
metabolic modulation targets enzymes after protein synthesis, often
constitutively expressed. For instance, modifying the physical clustering of
constitutively-expressed enzymes can impact metabolic fluxes 8. However,
constitutive expression of enzymes can impose a metabolic burden in cells 9.

Transcriptional gene modulation will be the focus of this paper. External
input signals used to modulate gene transcription range from chemicals such as
IPTG 10 and methanol 11, to process conditions, including temperature 12, pH
13, and oxygen supply 14. Many of these inducers, nevertheless, are
irreversible, i.e., it is challenging to remove them or reduce their levels
once triggered. Their effective delivery can also be hindered by mass- or
heat-transfer limitation in bioreactors 15. In addition, chemical inducers
might also induce cellular stress and toxicity 16, 17. Moreover, there may be
unintended cellular effects due to the non-orthogonality of the inducers. For
example, pH, temperature, and oxygen availability may have overreaching
effects on the overall physiological and metabolic state of cells 18, 19, 20.
These issues have motivated the rise of optogenetics in dynamic metabolic
engineering, employing light-responsive transcription factors to regulate gene
expression 21, 22, 23, 24, 25, 26, 27. Light is particularly appealing as a
control input signal due to its high tunability, reversibility, orthogonality,
and cost-effectiveness.

Determining the optimal inputs in dynamic metabolic engineering is not
trivial; for example, the optimal light input profiles in optogenetically
controlled systems. Model-based dynamic optimization has been shown to be a
powerful tool in predicting dynamic inputs for optimally steering the
expression of manipulatable enzymes in bioprocesses 28, 29. Yet, the
application of open-loop model-based optimization often fails to consider
system uncertainties, such as disturbances or model-plant mismatch. The
concept of metabolic cybergenetics (see Fig. 1), which is the focus of this
paper, has emerged as a way to address the limitations of open-loop control by
integrating state feedback 22, 29. These systems adopt closed-loop control,
such as model predictive control, by iteratively updating and solving the
open-loop optimization problem with the current system state 29. Other control
configurations are also possible, e.g., by merging external control via open-
loop optimization and in-cell feedback encoded by genetic circuits 30.

Figure 1: Scheme of a metabolic cybergenetic control system. Optimal system
inputs are obtained through computational methods such as dynamic
optimization. These external input signals modulate the expression of
metabolic enzymes and thereby metabolic fluxes. Closed-loop control is
implemented by incorporating state feedback and re-optimizing the system.
Optionally, if not all the state information can be measured, soft sensors in
the form of state estimators can be employed to reconstruct the system state.
The dynamic model can contain prior knowledge such as physics (e.g., metabolic
networks) and phenomenological relations (e.g., kinetic functions), which can
either be modeled explicitly or indirectly via machine-learning surrogates.
Process data can also be used to capture certain parts of the model using
machine learning.

Implementing model-based optimization for metabolic cybergenetics requires a
suitable mathematical model that captures the complex dynamics of relevant
metabolic processes. The model should integrate, for example, the dynamics of
the expression of manipulatable enzymes (micro-scale phenomena) and the
corresponding effect on process exchange rates (macro-scale phenomena).
Dynamic enzyme-cost flux balance analysis (deFBA) and its variations offer a
constraint-based modeling platform for coupling the dynamics of metabolism
with gene expression 31, 32, 33, 34, 35, 28, 29. These models are typically
underdetermined and thus formulated as dynamic optimization problems subject
to constraints on mass balances, resource allocation, and regulatory
mechanisms. They assume that cells are driven to optimize a dynamic biological
function, requiring for its formulation a significant understanding of the
short-term and long-term cellular evolutionary driving forces.

Structured models such as deFBA integrate the dynamics of key extracellular
species and intracellular biomass components, such as catalytic enzymes,
ribosomes, and storage compounds, among others. These models give valuable
insight into cellular metabolism and its underlying processes, and have been
considered in previous metabolic cybergenetic frameworks 29. Nevertheless,
using these models in process optimization leads to bilevel optimization
problems, as the optimal control problems are constrained by the constraint-
based metabolic models 32, 34, 35, 29. Solving bilevel optimizations is
inherently challenging both numerically and conceptually, often requiring
assumptions about the objectives’ collaborative interaction (e.g., optimistic
or pessimistic), which is not trivial. In an optimistic scenario, the inner
optimization, the constraint-based metabolic model, can be replaced by its
Karush-Kuhn-Tucker (KKT) conditions (see Refs. 36, 37, 38 for more details).
This substitution transforms the bilevel optimization into a single-level
optimization. Unfortunately, it also introduces non-convexity due to the
Lagrangian and the complementary slackness constraints of the KKT conditions.
Additionally, the Lagrange multipliers from the KKT conditions must be
included as additional optimization variables. Moreover, deFBA and similar
models usually involve a large number of states. In feedback control schemes
such as model predictive control, where one needs to update the system state
at different sampling times 39, using deFBA-type models imposes significant
hassles for implementation due to the large numbers of states to be monitored.

In the context of dynamic metabolic engineering, we recently introduced an
alternative simpler hybrid dynamic metabolic modeling technique augmented with
machine-learning surrogates 7 (submitted to the 22nd European Control
Conference). There, we represent metabolic exchange fluxes as functions of
manipulatable intracellular metabolic fluxes through neural-network surrogates
informed by flux balance analysis (FBA). Classical FBA operates under steady-
state assumptions while considering a static cell’s objective function 40, 41,
42. These surrogates for metabolic exchange fluxes, which link intra- and
extra-cellular metabolic domains, can be plugged into the process rates of
macro-kinetic dynamic models. Using a toy metabolic network as an example in
our previous work 7, we demonstrated that these hybrid models facilitate the
determination of optimal flux profiles through single-level dynamic
optimization problems, with the manipulatable intracellular fluxes being the
dynamic degrees of freedom.

A limitation of the previous hybrid metabolic modeling approach 7 is the
assumption that manipulatable intracellular fluxes can be adjusted directly
and without any delay. This simplification, however, is not applicable to
transcriptional-level metabolic flux modulation, the focus of this paper.
Here, we build upon our initial modeling framework by integrating the dynamics
of manipulatable enzyme expression. In this study, we differentiate between
gene expression in prokaryotic and eukaryotic cells to offer a more
comprehensive framework. We apply enzyme-capacity relationships that map the
concentration of manipulatable enzymes to their corresponding manipulatable
metabolic fluxes. After systematically exploring the solution space of FBA,
constrained by manipulatable enzyme concentrations, we develop neural-network
surrogates with the manipulatable enzyme concentrations as input features and
the FBA-predicted metabolic exchange fluxes as output labels. These surrogates
informed by the considered FBA model can then be used within the process rates
of macro-kinetic dynamic models augmented with manipulatable enzyme
expression. Our hybrid modeling can be straightforwardly integrated into
simplified single-level optimization problems for control and estimation in
metabolic cybergenetics. Moreover, due to the macro-kinetic nature of the
proposed models, the number of states involved is comparatively smaller than
it would be with for example deFBA-type models, simplifying process
monitoring.

Previous hybrid bioprocess models aim to learn or correct process rates using
process data 43, 44, 45, 46, 47. Our proposed modeling strategy differs by
embedding the physics of FBA into the process rates, making it physics-
informed. However, our approach does not exclude the use of machine learning
to refine uncertain parts of the model using process data, potentially in a
synergistic manner. Furthermore, this work extends beyond previous polynomial-
based functions for mapping sets of exchange metabolic fluxes predicted by FBA
48. Here we explicitly consider and exploit the intracellular domain for
cybergenetic control applications in dynamic metabolic engineering.

The remainder of this paper is structured as follows. In Section 2, we present
our hybrid physics-informed metabolic modeling framework with manipulatable
enzyme expression. Utilizing this hybrid machine-learning-supported model, we
formulate in Section 3 a feedback controller based on model predictive control
aimed at maximizing the production efficiency of metabolic cybergenetic
systems. We also outline an optimization-based estimator to reconstruct
unmeasured intracellular states. Remark that model-based optimization,
control, and estimation tasks are formulated as single-level optimization
problems. Finally, the proposed modeling and metabolic cybergenetic control
framework is demonstrated with a computational case dealing with the
optogenetically-assisted production of itaconate in batch fermentation by
Escherichia coli. Open-loop optimizations are used to aid in the in silico
design of the optogenetic gene expression system. More details on the case
study are given in Section 4.

##  2 Hybrid physics-informed dynamic modeling for metabolic cybergenetics

We introduce our hybrid physic-informed dynamic modeling framework for
metabolic cybergenetics distinguishing between prokaryotic and eukaryotic
cells.

###  2.1 Model structure for prokaryotic cells

Let the external states, encompassing substrates, products, and biomass dry
weight, be denoted as $\bm{z}\in\mathbb{R}^{n_{z}}$. The biomass dry weight is
indicated with $b\in\mathbb{R}$. The manipulatable intracellular enzymes are
represented by $\bm{e}\in\mathbb{R}^{n_{e}}$, while the external inputs for
inducing gene expression are denoted by $\bm{u}\in\mathbb{R}^{n_{u}}$.

We lump up transcription and translation processes in prokaryotic cells. In
prokaryotes such as E. coli, we deem this a valid assumption since
transcription and translation occur simultaneously in the cytoplasm, i.e.,
they are coupled in space and time 49. Therefore, for metabolic cybergenetics
using prokaryotic cells, we model a batch process as follows (cf. Fig. 2-a):

$\displaystyle\odv{\bm{z}(t)}{t}=b(t)\cdot\bm{q_{z}}(\bm{z}(t),\bm{e}(t),\bm{\theta},\bm{\Theta}),$
(1a)
$\displaystyle\odv{\bm{e}(t)}{t}=\bm{q_{e_{p}}}(\bm{u}(t),\bm{\theta})-\left(\mu(\bm{z}(t),\bm{e}(t),\bm{\theta},\bm{\Theta})\cdot\bm{1_{n}}+\bm{d_{e}}(\bm{e}(t),\bm{\theta})\right)\circ\bm{e}(t),$
(1b) $\displaystyle\bm{z}(t_{0})=\bm{z_{0}},\,\bm{e}(t_{0})=\bm{e_{0}}.$ (1c)

$\bm{\theta}\in\mathbb{R}^{n_{\theta}}$ encompasses constant parameters of the
explicitly-modeled prior knowledge/physics.
$\bm{\Theta}\in\mathbb{R}^{n_{\Theta}}$ represents constant parameters of the
neural-network surrogates (weights and biases) informed by FBA; these
correspond to the physics-informed machine-learning component of the model.
The accumulation of external states depends on biomass-specific process
exchange rates
$\bm{q_{z}}:\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{e}}\times\mathbb{R}^{n_{\theta}}\times\mathbb{R}^{n_{\Theta}}\rightarrow\mathbb{R}^{n_{z}}$.
The dynamics of intracellular manipulatable enzymes in prokaryotic cells
depend on the input-dependent production rate (lumped transcription and
translation)
$\bm{q_{e_{p}}}:\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}^{n_{e}}$,
cell dilution given by the growth rate
$\mu:\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{e}}\times\mathbb{R}^{n_{\theta}}\times\mathbb{R}^{n_{\Theta}}\rightarrow\mathbb{R}$,
and the degradation rate
$\bm{d_{e}}:\mathbb{R}^{n_{e}}\times\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}^{n_{e}}$.
Note that $\mu\in\bm{q_{z}}$, $\bm{1_{n}}\in\mathbb{R}^{n_{e}}$ is a vector of
ones, and $\circ$ indicates the Hadamard product. The initial process time is
denoted as $t_{0}$. We assume that each external input uniquely modulates the
expression of one corresponding manipulatable enzyme, thus $n_{u}=n_{e}$.

(a) Prokaryotic cell

(b) Eukaryotic cell

Figure 2: Flow diagram showing the interdependence of the external inputs and
the dynamic states of the proposed hybrid physics-informed cybergenetic model
for (a) prokaryotic and (b) eukaryotic cells. Refer to the text for details on
the notation.

It is worth mentioning that the external inputs affect the rates
$\bm{q_{e_{p}}}$, influencing enzyme expression. In contrast, the process
exchange rates are affected by the manipulatable enzyme concentrations, but
not directly by the external inputs. This helps to capture an input-output
delay in terms of the effect of the external input on the process exchange
rates. In an optogenetically controlled system, the external inputs could be,
for example, the light intensity or the frequency of light pulses.

###  2.2 Model structure for eukaryotic cells

Transcription and translation in eukaryotic cells are uncoupled in space and
time. Transcription takes place in the nucleus, followed by translation in the
cytoplasm 50. Therefore, to properly account for time delays in metabolic
cybergenetics using eukaryotic cells, we need to consider the dynamics of
transcription separate from translation. With this in mind, we model the batch
process as follows (cf. Fig. 2-b):

$\displaystyle\odv{\bm{z}(t)}{t}=b(t)\cdot\bm{q_{z}}(\bm{z}(t),\bm{e}(t),\bm{\theta},\bm{\Theta}),$
(2a)
$\displaystyle\odv{\bm{e}(t)}{t}=\bm{q_{e_{e}}}(\bm{p}(t),\bm{\theta})-\left(\mu(\bm{z}(t),\bm{e}(t),\bm{\theta},\bm{\Theta})\cdot\bm{1_{n}}+\bm{d_{e}}(\bm{e}(t),\bm{\theta})\right)\circ\bm{e}(t),$
(2b)
$\displaystyle\odv{\bm{p}(t)}{t}=\bm{q_{p}}(\bm{u}(t),\bm{\theta})-\left(\mu(\bm{z}(t),\bm{e}(t),\bm{\theta},\bm{\Theta})\cdot\bm{1_{n}}+\bm{d_{p}}(\bm{p}(t),\bm{\theta})\right)\circ\bm{p}(t),$
(2c)
$\displaystyle\bm{z}(t_{0})=\bm{z_{0}},\,\bm{e}(t_{0})=\bm{e_{0}},\,\bm{p}(t_{0})=\bm{p_{0}}.$
(2d)

In this case, $\bm{p}\in\mathbb{R}^{n_{e}}$ represents the product of
transcription, i.e., the mRNA transcripts. The dynamics of the latter state is
affected by the rates of transcription
$\bm{q_{p}}:\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}^{n_{e}}$,
cell dilution (growth rate), and mRNA degradation
$\bm{d_{p}}:\mathbb{R}^{n_{p}}\times\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}^{n_{e}}$.
In contrast to the model for prokaryotic cells, the dynamics of enzyme
production in eukaryotic cells are affected by the translation rate
$\bm{q_{e_{e}}}:\mathbb{R}^{n_{e}}\times\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}^{n_{e}}$,
which does not directly depend on the external input, but rather on the
availability of mRNA transcripts.

Throughout this work, $\bm{z}$ and $b$ are expressed in $\mathrm{mmol/L}$ and
$\mathrm{g_{b}/L}$, respectively. $\mathrm{g_{b}}$ refers to grams of cell dry
weight. The intracellular states $\bm{e}$ and $\bm{p}$ are both expressed in
$\mathrm{mmol/g_{b}}$. While the models above are written for batch processes
(the focus of this work), they can be expanded to other modes of operation,
e.g., fed-batch or continuous processes, by incorporating proper dilution
rates. From now on, we will omit the time dependency of the variables unless
unclear from the context. In addition, we consistently use bold-notation to
represent vectors and matrices and non-bold notation to represent scalars and
unidimensional variables and functions.

###  2.3 Macro-scale process exchange rates

Defining the rates $\bm{q_{z}}$ requires elucidating the relationship between
manipulatable enzyme concentrations and the change in the cell’s metabolic
flux distribution. One could try to use phenomenological relationships, but
these often tend to overlook important intracellular aspects such as intrinsic
metabolic trade-offs and energy/redox balancing 51. Therefore, we consider a
hybrid physics-informed approach where:

$\displaystyle q_{z_{i}}=v_{\mathrm{ext},{i}}(\bm{e},\bm{\Theta})\cdot
h(\bm{z},\bm{\theta}),\,\forall i\in\\{1,2,...,n_{z}\\}.$ (3)

In the above equation,
$v_{\mathrm{ext},{i}}:\mathbb{R}^{n_{e}}\times\mathbb{R}^{n_{\Theta}}\rightarrow\mathbb{R}$
represents the metabolic exchange flux for an extracellular state $i$
predicted by a neural-network surrogate of FBA. $v_{\mathrm{ext},{i}}$ can be
regarded as a maximum theoretical rate value as classical FBA typically
neglects limitations due to external concentrations or conditions. The
function
$h:\mathbb{R}^{n_{z}}\times\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}$
explicitly captures possible prior knowledge on aspects such as substrate-
limitation or product-inhibition that reduce the flux $v_{\mathrm{ext},{i}}$
–not captured by the FBA model, thus also not by the corresponding surrogate–.
The function $h$ can in principle also integrate, partly or fully, a machine-
learning component, e.g., a correction or missing term learned from process
data. For clarity of notation, remark that we refer to $q_{z_{i}}$ as process
(exchange) rates and $v_{\mathrm{ext},{i}}$ as metabolic (exchange) fluxes
throughout this work.

We collect all the metabolic exchange fluxes in
$\bm{v_{\mathrm{ext}}}\in\mathbb{R}^{n_{\mathrm{ext}}}$. Then, the neural-
network surrogates can be formulated as:

$\bm{v_{\mathrm{ext}}}=f_{\mathrm{NN}}(\bm{e},\bm{\Theta}),$ (4)

where
$\bm{v_{\mathrm{ext}}}=\\{v_{\mathrm{ext},i}\,|\,i\in\\{1,2,...,n_{z}\\}\\}$.
The function
$f_{\mathrm{NN}}:\mathbb{R}^{n_{e}}\times\mathbb{R}^{n_{\Theta}}\rightarrow\mathbb{R}^{n_{z}}$
represents the neural-network surrogate that maps the manipulatable
intracellular enzyme concentrations to the metabolic exchange fluxes.

In this work, we consider feedforward neural networks, which consist of
interconnected nodes organized into $n_{l}$ layers. The output vector
$\bm{a_{l}}\in\mathbb{R}^{n_{a_{l}}}$ of a layer $l$ in a trained feedforward
neural network is represented as follows 52:

$\bm{a_{l}}=\bm{\sigma_{l}}(\mathbf{W}_{l}\cdot\bm{a_{l-1}}+\bm{b_{l}}),\,\forall
l\in\\{1,2,...,n_{l}\\}.$ (5)

The features are the inputs to the first layer, i.e., $\bm{a_{0}}=\bm{e}$, and
the labels are the outputs of the last layer, i.e.,
$\bm{a_{n_{l}}}=\bm{v_{\mathrm{ext}}}$.
$\bm{\sigma_{l}}:\mathbb{R}^{n_{a_{l-1}}}\rightarrow\mathbb{R}^{n_{a_{l}}}$ is
the layer’s activation function. $\mathbf{W}_{l}\in\mathbb{R}^{n_{a_{l}}\times
n_{a_{l-1}}}$ and $\bm{b_{l}}\in\mathbb{R}^{n_{a_{l}}}$ are weight matrices
and bias vectors, respectively. Note that different machine-learning methods,
besides neural networks, can in principle also be used. However, the benefit
of feedforward neural networks is their feature of being universal function
approximations 53.

###  2.4 Micro-scale rates of the model

Regarding the rates involved in the dynamics of the expression of
manipulatable enzymes, with or without lumped transcription and translation,
one can use, e.g., the widely-known Hill function for either input-dependent
activation or repression 54, 55, 56. Degradation rates can be modeled, e.g.,
following commonly employed first-order degradation kinetics 57, 58, 59. An
example of these rate functions will be shown in the itaconate computational
case study in Section 4.

###  2.5 Training the neural-network surrogates for enzyme-dependent metabolic
exchange fluxes

We now outline how to build the enzyme-dependent surrogates for the metabolic
exchange fluxes based on the physics of FBA. To do so, we first link
manipulatable intracellular enzymes to manipulatable intracellular fluxes
$\bm{v_{\mathrm{man}}}\in\mathbb{R}^{n_{\mathrm{man}}}$ via enzyme-capacity
relationships. For a manipulatable enzyme $i$:

$e_{i}=\left|\frac{v_{\mathrm{man},{i}}}{k_{\mathrm{cat},i}}\right|,\,\forall
i\in\\{1,2,...,n_{e}\\},$ (6)

where
$\bm{k_{\mathrm{cat}}}=\\{k_{\mathrm{cat},i}\,|\,i\in\\{1,2,...,n_{e}\\}\\}$
contains the catalytic constants of the manipulatable enzymes. Eq. (6) assumes
that the manipulatable enzymes operate under substrate saturation conditions,
close to their maximum reaction rates. In fact, in microorganisms such as E.
coli, most enzymes typically work under substrate saturation 60. This
simplification avoids incorporating intracellular metabolite concentrations in
Eq. (6), thereby keeping the model simple. Note that $|\cdot|$ indicates the
absolute value operation, which is necessary to handle possible negative
metabolic fluxes in case of reversible reactions.

The next step is to populate a dataset
$\mathcal{D}\in\mathbb{R}^{n_{g}\times(n_{\mathrm{ext}}+n_{e})}$ with input
features (manipulatable enzyme concentrations) and output labels (metabolic
exchange fluxes) to train our surrogate model of FBA. Using a systematic grid
search approach, we run simulations of a given FBA model across feasible
combinations of intracellular manipulatable fluxes. The exploration space is
determined by the Cartesian product $\mathbb{G}$ of the sets
$\mathbb{V}_{\mathrm{man},1},\mathbb{V}_{\mathrm{man},2},\ldots,\mathbb{V}_{\mathrm{man},n}$:

$\mathbb{G}=\mathbb{V}_{\mathrm{man},1}\times\mathbb{V}_{\mathrm{man},2}\times\ldots\times\mathbb{V}_{\mathrm{man},n},$
(7)

where the set $\mathbb{V}_{\mathrm{man},n}$ contains user-defined values in
the grid for the $n$-th manipulatable flux. Note that the total number of
combinations of manipulatable fluxes is indicated by
$n_{g}=\text{size}(\mathbb{G})$ and each element of $\mathbb{G}$ is an
$n_{\mathrm{man}}$-tuple.

Let us refer to the $j$-th element of $\mathbb{G}$ as $\bm{g_{j}}$. To fill the dataset $\mathcal{D}$, the FBA model constrained by the manipulatable enzyme concentrations is systematically solved as: v,eF(v,e), ˙m=Sv = 0, v_min ≤v ≤v_max, v_man = g_j, v_man ⊆v, e_i = | vman,ikcat,i |, ∀i ∈{1,2,…,n_e}, ∀j ∈{1,2,…,n_g}.

$F:\mathbb{R}^{n_{v}}\times\mathbb{R}^{n_{e}}\rightarrow\mathbb{R}$ captures
the assumed static objective function that the cell maximizes. The decision
variables of the optimization are the (intracellular and extracellular)
metabolic fluxes of the metabolic network $\bm{v}\in\mathbb{R}^{n_{v}}$ and
the intracellular manipulatable enzymes. Furthermore,
$\bm{m}\in\mathbb{R}^{n_{m}}$ represents intracellular metabolites and
$\bm{S}\in\mathbb{R}^{n_{m}\times n_{v}}$ is the stoichiometric matrix of the
intracellular metabolites. Eq. (2.5) is equal to a zero vector 0 of
appropriate dimension, consistent with the steady-state assumption of FBA.
Feasible lower ($\bm{v_{\mathrm{min}}}\in\mathbb{R}^{n_{v}}$) and upper bounds
($\bm{v_{\mathrm{max}}}\in\mathbb{R}^{n_{v}}$) constrain the metabolic fluxes
in Eq. (2.5). The different combinations of intracellular manipulatable fluxes
are integrated into the FBA optimization via the constraint in Eq. (2.5). The
enzyme-capacity relationships between the manipulatable enzymes and
manipulatable fluxes are represented by Eq. (2.5). Note that, if necessary,
the FBA model’s predictive accuracy can be in principle enhanced by
considering additional constraints dealing with resource allocation,
thermodynamics, and regulation, among other aspects. For example, one can
consider enzyme-capacity constraints as in Eq. (2.5) on other (non-
manipulatable) fluxes as well. For simplicity, we only consider in this work
enzyme-capacity constraints on the manipulatable fluxes.

The solution of the FBA is iteratively appended to the dataset as follows:

$\mathcal{D}\leftarrow\mathcal{D}\cup\\{(\bm{v_{\mathrm{ext}}}^{(j)},\bm{e}^{(j)})\\},\,\forall
j\in\\{1,2,\ldots,n_{g}\\},$ (8)

where $\bm{v_{\mathrm{ext}}}^{(j)}$ and $\bm{e}^{(j)}$ represent the solution
vectors of $\bm{v_{\mathrm{ext}}}$ and $\bm{e}$. The populated dataset,
containing all the solutions from the FBA model under the different
combinations of intracellular manipulatable fluxes in $\mathbb{G}$, is then
used to train the neural-network surrogate in Eq. (4).

The training of the neural network is performed by iterative adaptation of the
weight matrices and bias vectors (cf. Eq. (5)) so that the loss function is
minimized. Here we use the mean squared error (MSE) as the loss function:

$\mathrm{MSE}=\frac{1}{n_{g}\times n_{z}}\sum_{i=1}^{n_{g}\times
n_{z}}(\mathfrak{y}_{i}-\hat{\mathfrak{y}}_{i})^{2},$ (9)

where $\hat{\mathfrak{y}}_{i}\in\mathbb{R}$ is the predicted value of the
output label by the neural network for the $i$-th data entry, and
$\mathfrak{y}_{i}\in\mathbb{R}$ is the actual output label value for the
$i$-th data entry as predicted by FBA. In total, we have $n_{g}\times n_{z}$
data entries, which corresponds to having $n_{g}$ FBA solutions, each with
$n_{z}$ associated output labels.

The proposed surrogate approach assumes that an appropriate metabolic network
is a priori available. Note that there are several genome-scale models at hand
for different microorganisms derived from gene-enzyme-reaction associations
61. Furthermore, the surrogate strategy presented in this section can be in
principle expanded to learn other key metabolic fluxes, even if they are not
involved in the differential equations of the hybrid models. This can help to
achieve continued insight into metabolism, e.g., for process monitoring,
without the need of solving an FBA problem at every sampling point. In this
work, however, we focus only on obtaining machine-learning surrogates of FBA
for $\bm{v_{\mathrm{ext}}}$.

##  3 Optimization-based feedback control and estimation using the hybrid
physics-informed model

With the hybrid physics-informed dynamic model outlined in the previous
section, we now aim to find the external input trajectories that maximize
bioprocess production efficiency. To do so, we formulate a model-based optimal
controller following a metabolic cybergenetic approach. The overall
methodology is presented in Fig. 3 and will be outlined in detail in this
section.

Figure 3: Flow diagram of our integrated hybrid physics-informed dynamic
modeling and control strategy for metabolic cybergenetics. Details on steps 1
and 2 are provided in Section 2, while details on steps 3 and 4 are provided
in Section 3.

###  3.1 Model predictive control

Let us assume for simplicity equidistant sampling points throughout the
process. The sampling time is denoted as $t_{k}:=kh_{s}$,
$k\in\mathbb{N}_{0}$, with a sampling interval $h_{s}\in\mathbb{R}_{>0}$.
Also, $\bm{x}\in\mathbb{R}^{n_{x}}$ is the state vector, i.e.,
$\bm{x}:=[\bm{z}^{\mkern-1.5mu\mathsf{T}},\bm{e}^{\mkern-1.5mu\mathsf{T}}]^{\mkern-1.5mu\mathsf{T}}$
for prokaryotes and
$\bm{x}:=[\bm{z}^{\mkern-1.5mu\mathsf{T}},\bm{e}^{\mkern-1.5mu\mathsf{T}},\bm{p}^{\mkern-1.5mu\mathsf{T}}]^{\mkern-1.5mu\mathsf{T}}$
for eukaryotes. Following a model predictive control approach 39, we
repetitively solve an optimal control problem with the updated system state
$\bm{x}(t_{k})$, thereby incorporating state feedback. The state information
can be obtained from hardware or soft sensors 62, 63. The model predictive
controller at time $t_{k}$ is formulated as the following optimization problem
constrained by the hybrid physics-informed dynamic model: uJ(⋅), Eqs.
(1a)-(1b), if prokaryote, Eqs. (2a)-(2c), if eukaryote, Eqs. (3)-(4), 0
≤g_p(x,u,θ,Θ), u = {{u(⋅)|t0th, x0}, if  tk=t0,u(⋅)|tkth, if  tk>t0, x(t_k) =
x_k.

The efficiency of the process is captured by the function $J(\cdot)$. In a
production context, the optimal control problem can be formulated to maximize
a metric such as volumetric productivity, economic profit, environmental
impact, or a weighted combination of those. One could also think of minimizing
the error in set-point- or trajectory-following-type problems, in case an
optimal set-point or trajectory is a priori known. An important consideration
is made regarding the decision variables of the controller $\bm{\mathfrak{u}}$
(cf. Eq. (3.1)). The dynamic decision variable of the optimization is the
input function over $t\in[t_{k},t_{h}]$, where $t_{h}$ is the time at the end
of the prediction horizon. In batch setups, we often set $t_{h}=t_{f}$, where
$t_{f}$ is the final batch time, hence a shrinking prediction horizon. In
other circumstances, such as in continuous processes, one may consider that
$t_{h}$ moves alongside the prediction horizon.

Within the optimization, at $t_{k}=t_{0}$, one could also determine the
optimal initial state of the system $\bm{x_{0}}$ (cf. Eq. (3.1)). For example,
one could start the process at specific concentrations of biomass, substrate,
and manipulatable enzymes (via previous induction). However, remark that once
the process has started, i.e., at $t_{k}>t_{0}$, we are only allowed to adjust
the dynamic input. Possible (nonlinear) state or input constraints are also
considered in the optimization problem. These are represented by the function
$\bm{g_{p}}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{\theta}}\times\mathbb{R}^{n_{\Theta}}\rightarrow\mathbb{R}^{n_{g_{p}}}$
and can include economic, technical, environmental, and safety constraints. It
is also worth mentioning that if the optimization in (3.1) is only solved at
$t_{k}=t_{0}$ and does not incorporate state feedback, then it is regarded as
an open-loop optimization.

###  3.2 Estimation of mRNA transcripts and manipulatable enzymes

One of the challenges of (metabolic) cybergenetics is how to measure changes
in intracellular states in real-time for updating the model predictive
controller. For example, accurately measuring concentrations of mRNA
transcripts and intracellular enzymes often relies on offline technologies
such as transcriptomics 64 and proteomics 65. The latter technologies are
time-consuming and thus impractical to be employed for real-time monitoring.
One possibility for real-time monitoring of intracellular enzymes is to use
fluorescent protein tags as reporters or biosensors 66, 67. However, the read-
out can be challenging if the spectra of multiple reporters overlap.
Additionally, the biosensor machinery might impose a resource burden on the
cells. Also, biosensors are not without delays in their activation and
inactivation. Consequently, in scenarios where measuring the intracellular
states is challenging or impractical, we propose to estimate them via soft
sensors.

Specifically, we propose an estimator based on moving horizon estimation and
full information estimation that considers the dynamic model of the system and
a history of past measurements, to infer the current state of the system 68,
69, 70, 71, 72, 73, 29, 45. Our soft sensor at time $t_{k}$ is formulated in a
discrete manner as follows 39, 52: x_k-N, w∥x_k-N \- ¯x_k-N ∥_P^2 + ∑_i=k-N^k
∥^y_i \- y_i ∥_R^2 + ∑_i=k-N^k ∥w_i ∥_Q^2, x_i + h_s= x_i \+
∫_t_i^t_i+h_sf_d(x,^u,θ,Θ) dt + w_i, y_i = η(x_i,^u_i,θ,Θ), 0
≤g_p(x_i,^u_i,θ,Θ), ∀i ∈{k-N,…,k}.

Here
$\bm{f_{d}}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{\theta}}\times\mathbb{R}^{n_{\Theta}}\rightarrow\mathbb{R}^{n_{x}}$
collects all the dynamic equations of the hybrid physics-informed dynamic
model. $\bm{y}\in\mathbb{R}^{n_{y}}$ is the measurement vector and
$\bm{\eta}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\times\mathbb{R}^{n_{\theta}}\times\mathbb{R}^{n_{\Theta}}\rightarrow\mathbb{R}^{n_{y}}$
is the measurement function. $N\in\mathbb{N}$ is the length of the estimation
window. $\bm{\hat{u}_{i}}$ and $\bm{\hat{y}_{i}}$ are the input and output
measurements at time $t_{i}$, respectively.
$\bm{{w}_{i}}\in\mathbb{R}^{n_{x}}$ is the state noise vector at time $t_{i}$,
thus $\bm{w}:=\\{\bm{w_{i}},\,\forall i\in\\{k-N,...,k\\}\\}$.

Furthermore, $\bm{\bar{x}_{k-N}}$ is the best guess of the state at the
beginning of the estimation window. $\mathbf{P}$, $\mathbf{R}$ and
$\mathbf{Q}$ are weighting matrices of appropriate dimensions.
$\|\bm{a}\|_{\mathbf{A}}^{2}:=\bm{a}^{\mkern-1.5mu\mathsf{T}}\mathbf{A}\bm{a}$
is the weighted squared norm of a vector $\bm{a}$ with respect to a weight
matrix $\mathbf{A}$. The first term of the objective function in (3.2)
represents the arrival cost at the beginning of the estimation window. This
serves as a kind of "memory" of past points left behind by the estimation
window. The second term of the objective function represents the measurement
error, while the third term represents the process noise or state
disturbances. The decision variable of the estimator is the initial state at
the beginning of the estimation window $\bm{x_{k-N}}$ and the state noise
$\bm{w}$, which can be used to regenerate the filtered/smoothed system
dynamics via Eq. (3.2). Note that even if the states can be measured through
biosensors or other means, the estimator can still be beneficial to filter
possible noisy measurements.

Note that uncertain parameters can in principle also be estimated with the
described soft sensor. However, for simplicity, here we assume constant model
parameters in the estimation. It is also worth noting that the estimation
window grows in time, i.e., the number of data points within the estimation
window increases, until the length $N$ is reached. Afterward, the estimation
window only moves, and the estimator is regarded as a moving horizon
estimator. If we choose $N$ such that the estimation window always starts from
$t_{0}$, the estimation window always grows, and we refer to this estimator as
a full information estimator. We employ a full information estimator in the
case study presented in the next section. Although here we consider
equidistant sampling points, multi-rate estimators are also possible 69.

##  4 Computational case study: itaconate biosynthesis by E. coli with
optogenetic control

As a starting point for our metabolic cybergenetic case study, we consider E.
coli ita36A, a strain engineered for itaconate biosynthesis from glucose under
aerobic conditions 12. E. coli ita36A features deletions of the genes aceA
(isocitrate lyase), pta (phosphate acetyltransferase), pykF (pyruvate kinase
I), and pykA (pyruvate kinase II). Additionally, the original strain harbors a
plasmid that encodes the cadA gene (cis-aconitate decarboxylase) from
Aspergillus terreus. Cis-aconitate decarboxylase enables the synthesis of
itaconate from cis-aconitate, an intermediate metabolite of the tricarboxylic
acid (TCA) cycle, through a decarboxylation reaction. The plasmid also encodes
gene gltA (citrate synthase) from Corynebacterium glutamicum, enhancing the
flux through the upper branch of the TCA cycle. It is important to note that
both cadA and gltA are constitutively expressed in E. coli ita36A.

Furthermore, the native icd (isocitrate dehydrogenase) in E. coli ita36A has
been replaced by a synthetic genetic switch that regulates icd transcription
via temperature control. The enzyme isocitrate dehydrogenase modulates the
flux of the TCA cycle downstream isocitrate. Therefore, as demonstrated by
Harder et al., the temperature-regulation of the expression of icd can enhance
the volumetric productivity of itaconate in two-stage fermentation processes
12. This can be done by performing a growth phase (with icd expression off at
37 °C) followed by a production phase (with icd expression on at 28 °C). The
flux through the TCA cycle is intrinsically linked to growth, supplying
precursor metabolites, energy cofactors, and reducing equivalents essential
for biomass synthesis. Thus, reducing the flux through the TCA cycle by
suppressing icd expression indirectly increases the flux through itaconate
biosynthesis, although at the expense of growth due to the less active TCA
cycle. The temperature-inducible control of icd showed promising results for
enhancing the process volumetric productivity; however, fine-tuning the
bioreactor’s temperature in large-scale setups poses challenges, particularly
due to heat-transfer limitations 15. Additionally, temperature is a non-
orthogonal control input, potentially affecting the transcription of
untargeted genes in E. coli 19.

In our computational case study, we introduce several modifications to the
original E. coli ita36A strain 12. First, we consider the dynamic regulation
of cadA as an alternative and more direct mechanism to balance flux between
the TCA cycle (growth mode) and itaconate biosynthesis (production mode). An
increase in the flux of cis-aconitate decarboxylase is expected to divert
carbon flux from the TCA cycle toward itaconate biosynthesis, as predicted by
FBA simulations (cf. Section 4.2.2). In this modified strain, icd is
constitutively expressed. For experimental implementations, however, it is
crucial to evaluate whether the cell will indeed reduce the flux through the
TCA cycle in favor of itaconate synthesis with increased expression of cadA.
There is a well-grounded possibility that the evolutionary driving force for
cellular growth will make the TCA cycle preferable, despite increased
expression of cadA. If such a phenotype is observed, one could further
consider simultaneously inducing cadA externally and repressing icd expression
to achieve tighter control. Having said that, for simplicity of demonstration
of the hybrid physics-informed modeling and control cybergenetic framework, we
assume that regulating cadA alone is sufficient to balance the flux through
the TCA cycle and itaconate biosynthesis.

Secondly, this computational study incorporates an optogenetic gene expression
system to regulate the expression of cadA. Utilizing light as a gene
expression control input offers enhanced tunability and orthogonality compared
to other inducers such as temperature. A simplified diagram illustrating the
metabolism of this E. coli strain for itaconate biosynthesis is presented in
Fig. 4. Throughout this study, our adapted E. coli strain with optogenetic
regulation is referred to as E. coli vOpt, in which cadA expression is
activated by light in the fermentation.

Figure 4: Scheme of E. coli vOpt for itaconate biosynthesis. We represent gene
deletions with red crosses. We indicate in italics genes encoding relevant
enzymes. Orange dotted arrows denote metabolic exchange reactions. The green
arrow with a resistor shape represents the reaction associated with the light-
inducible cis-aconitate decarboxylase. $v_{i}\in\mathbb{R}$ denotes the
metabolic exchange flux of the external species $i$.

###  4.1 Assumptions and remarks on the hybrid physics-informed model

Since E. coli vOpt is a prokaryotic organism, the hybrid physics-informed
metabolic model, as described by Eqs. (1a)-(1c), is applicable to our case
study. We consider four external states: the concentrations of glucose,
itaconate, acetate, and biomass; hence, the state vector is denoted as
$\bm{z}:=[z_{\mathrm{glc}},z_{\mathrm{ita}},z_{\mathrm{ace}},b]^{\top}$. For
simplicity, gas exchange is not included in the model and is assumed not to be
a rate-limiting factor. Enzyme cadA is the only enzyme under external control
of gene expression, thus $\bm{e}:=e_{\mathrm{cadA}}$,
$\bm{q_{e_{p}}}:=q_{e_{\mathrm{cadA}}}$, and
$\bm{d_{e}}:=d_{e_{\mathrm{cadA}}}$. The input-dependent expression rate of
cadA, $q_{e_{\mathrm{cadA}}}$, is modeled using the Hill activation function
for lumped transcription and translation:

$\displaystyle
q_{e_{\mathrm{cadA}}}=\theta_{1}+\theta_{2}\frac{u_{\mathrm{cadA}}^{\theta_{3}}}{\theta_{4}^{\theta_{3}}+u_{\mathrm{cadA}}^{\theta_{3}}},$
(10)

with parameters $\theta_{3}=2.0780\,\mathrm{[dimensionless]}$ and
$\theta_{4}=0.3799\,\mathrm{W/m^{2}}$ based on the values reported for the
CcaS-CcaR 2.0 two-component optogenetic system 74. Here the input
$u_{\mathrm{cadA}}\in\mathbb{R}_{\geq 0}$ is the intensity of green light
($\lambda_{\mathrm{max}}=520\,\mathrm{nm}$) in $\mathrm{W/m^{2}}$, triggering
cadA transcription, hence $\bm{u}:=u_{\mathrm{cadA}}$. See Ref. 74 for
additional details on the mechanism of the CcaS-CcaR 2.0 optogenetic system.
For simplicity, we assume no constitutive expression of the target enzyme,
setting $\theta_{1}=0$. The parameter $\theta_{2}\in\mathbb{R}_{\geq 0}$ is
the maximal input-dependent rate of expression in $\mathrm{mmol/g_{b}/h}$. In
later sections, we will compare different values of $\theta_{2}$, which
represent various design parameter values of the optogenetic gene expression
system. In this work, we assume homogeneous light penetration inside the
bioreactor, however, non-homogeneous light penetration can be also considered
as in Ref. 29.

The degradation rate of enzyme cadA is modeled using first-order kinetics.
Therefore:

$\displaystyle d_{e_{\mathrm{cadA}}}=\theta_{5},$ (11)

where $\theta_{5}=0.6931\,\mathrm{h^{-1}}$ assuming an average protein half-
life time of 1 $\mathrm{h}$ in E. coli 75.

The function $h$ in Eq. (3) considers limitation by the substrate glucose,
modeled using a hyperbolic (Monod) function, and inhibition by the by-product
acetate, modeled using a hyperbolic decay function. Therefore,

$\displaystyle
h=\left(\frac{z_{\mathrm{glc}}}{\theta_{6}+z_{\mathrm{glc}}}\right)\left(\frac{\theta_{7}}{\theta_{7}+z_{\mathrm{ace}}}\right),$
(12)

where the parameter $\theta_{6}=2.964\times 10^{-04}\,\mathrm{mmol/L}$ is the
assumed substrate affinity constant in the fermentation of glucose by E. coli
from previous determinations 76. The parameter
$\theta_{7}=134.63\,\mathrm{mmol/L}$ was fitted using SciPy’s curve_fit module
77 with experimental data of acetate inhibition on E. coli’s growth at pH 7.4
78.

To relate the intracellular concentration of the enzyme cadA to the attainable
metabolic flux through the itaconate biosynthesis reaction (cf. Eq. (6)), we
consider a catalytic constant $k_{\mathrm{cadA}}=66,240\,\mathrm{h^{-1}}$
reported for the Aspergillus version of the enzyme 79.

###  4.2 Neural-network surrogates informed by FBA

Here, we explain how we built the physics-informed component of our hybrid
dynamic model (cf. Eq. (4)) for E. coli vOpt. First, we need a suitable
metabolic network such that we can run FBA simulations under different levels
of manipulatable flux values as described in the optimization problem (2.5).
Thus, we employed a version of the metabolic network EColiCore2 containing the
heterologous reaction catalyzed by cadA provided in Ref. 80. The metabolic
network comprises 115 reactions and 93 metabolites. Note that the original
EColiCore2 81, derived from the genome-scale model iJO1366 of E. coli 82, is a
reduced metabolic network representing the cell’s central metabolism.

####  4.2.1 Preliminary accuracy check of the constraint-based model

To validate the accuracy of the metabolic network to be used in our case study
with E. coli vOpt, we first ran FBA simulations using COBRApy 83 constrained
by experimentally-reported metabolic exchange fluxes for E. coli ita36A with
temperature-dependent dynamic control of the TCA cycle 12. The objective
function of the FBA was set to maximize the biomass metabolic flux
$v_{\mathrm{bio}}$. As shown in Table 1, the biomass metabolic fluxes
predicted by FBA closely matched the experimentally-reported values during
both growth and production phases by E. coli ita36A. Therefore, we deemed this
version of EColiCore2 suitable for testing our hybrid metabolic cybergenetic
framework with E. coli vOpt. Note that only constraints represented by Eqs.
(2.5) and (2.5) were considered in the preliminary accuracy check of the FBA
model outlined in Table 1, i.e., no enzyme-related constraints were taken into
account.

Table 1: Preliminary accuracy check of the proposed constraint-based model to
be used in our case study with E. coli vOpt. FBA simulations were performed
considering E. coli ita36A, with temperature-dependent dynamic control of the
TCA cycle, constrained by experimentally-measured metabolic exchange fluxes.
Exp.: experimental; pred.: predicted.

Exp. fluxesa [$\mathrm{mmol/g_{b}/h}$] | Exp. $v_{\mathrm{bio}}$a [$\mathrm{h^{-1}}$] | Pred. $v_{\mathrm{bio}}$b [$\mathrm{h^{-1}}$] | Reported phasea
---|---|---|---
$v_{\mathrm{glc}}=3.48$, $v_{\mathrm{ita}}=0.17$ | 0.25 $\pm$ 0.01 | 0.27 | Growth
$v_{\mathrm{glc}}=0.85$, $v_{\mathrm{ita}}=0.43$ | 0.03 $\pm$ 0.00 | 0.01 | Production

  * a

Based on the experimental data for E. coli ita36A provided in Ref. 12.

  * b

For these FBA simulations, we assumed an ATP maintenance demand of 8.4
$\mathrm{mmol/g_{b}/h}$ and a formate excretion rate equal to zero based on
the phenotype of previous iterations of E. coli ita36A 80. We also set the
exchange flux of glutamate, a possible way to replenish $\alpha$-ketoglutarate
in the TCA, equal to zero. Gene knock-outs were simulated by setting to zero
the metabolic fluxes associated with the genes aceA, pta, pykF, and pykA.

####  4.2.2 Systematic exploration of cis-aconitate decarboxylase flux
manipulation

Prior to systematically exploring manipulatable intracellular flux values in
E. coli vOpt, we constrained the FBA model to match the expected phenotype. We
constrained the FBA model for E. coli vOpt with an ATP demand of 8.4
$\mathrm{mmol/g_{b}/h}$ and a formate exchange flux equal to zero, following a
similar rationale as in the FBA simulations with E. coli ita36A in Table 1.
The upper bound of the glucose exchange flux was set to 3.48
$\mathrm{mmol/g_{b}/h}$, matching the experimental growth scenario in Table 1,
while the fluxes associated with genes knock-outs (aceA, pta, pykF, and pykA)
were set equal to zero.

We consider only one manipulatable flux $v_{\mathrm{cadA}}$, namely the flux
catalyzed by enzyme cadA. The set of explored manipulatable flux values for
this enzyme is denoted as $\mathbb{V}_{\mathrm{man,cadA}}$, hence
$\mathbb{G}=\mathbb{V}_{\mathrm{man,cadA}}$ (cf. Eq. (7)). The set
$\mathbb{V}_{\mathrm{man,cadA}}$ was defined as
$\mathbb{V}_{\mathrm{man,cadA}}:=\\{0.000,0.007,0.014,\ldots,3.476\\}$,
covering the feasible space of growth-production yield trade-offs. With
$\mathbb{G}$ defined, we systematically solved the FBA problem, as detailed in
the optimization problem (2.5), to fill our dataset represented by Eq. (8).

Upon systematic exploration, Figure 5-a illustrates the predicted impact of
intracellular concentration of enzyme cadA on the metabolic exchange fluxes
predicted by the FBA model. It is observed that the exchange flux of itaconate
increases linearly with the concentration of enzyme cadA, aligning with the
assumed enzyme-capacity relationship described by Eq. (6). As expected, this
comes at the expense of a reduced biomass flux. Similarly, a decrease in the
acetate exchange flux is observed with increasing concentrations of enzyme
cadA, despite the fact that the acetate exchange flux is already minimal due
to the pta knock-out (the residual acetate in the FBA model originates from
the synthesis reactions of the amino acids arginine and cysteine). In all
scenarios, the FBA model predicts a glucose exchange flux of 3.48
$\mathrm{mmol/g_{b}/h}$, consistent with the flux’s upper limit. The explored
metabolic trade-off between growth and production is presented in Fig. 5-b. It
becomes clear that the defined minimum and maximum values in
$\mathbb{V}_{\mathrm{man,cadA}}$ were chosen to cover the range from the
maximum biomass flux ($v_{\mathrm{bio}}=0.277\,\mathrm{h^{-1}}$), with no
itaconate exchange flux ($v_{\mathrm{ita}}=0.000\,\mathrm{mmol/g_{b}/h}$ at
$v_{\mathrm{cadA}}=0.000\,\mathrm{mmol/g_{b}/h}$), to the maximum itaconate
exchange flux ($v_{\mathrm{ita}}=3.476\,\mathrm{mmol/g_{b}/h}$), with no
growth ($v_{\mathrm{bio}}=0.000\,\mathrm{h^{-1}}$ at
$v_{\mathrm{cadA}}=3.476\,\mathrm{mmol/g_{b}/h}$).

(a) Relevant metabolic exchange fluxes (b) Yield trade-off between growth and
production

Figure 5: (a) Changes of metabolic exchange fluxes with varying intracellular
cis-aconitate decarboxylase concentration predicted by the FBA model. (b)
Solution space of the yield trade-off between itaconate and biomass synthesis
from the exploration with FBA. glc: glucose, ita: itaconate, bio: biomass,
ace: acetate, cadA: cis-aconitate decarboxylase.

With the data generated during the systematic exploration, we trained a fully-
connected feedforward neural network using the MSE loss function represented
by Eq. (9). To do so, we used the Python toolbox HILO-MPC 52 interfaced with
PyTorch 84. The goal was to generate Eq. (4), the physics-informed machine-
learning part of the model. As the input feature, we considered the
intracellular concentration of enzyme cadA, and as output labels, the
metabolic exchange fluxes of itaconate, acetate, and biomass. The metabolic
exchange flux of glucose was not included as a label because it remained
constant throughout the FBA exploration (see Fig. 5-a), and was therefore
assumed constant in Eq. (3). We shuffled the dataset, populated with the FBA
solutions from the exploration phase, and reserved 15 % of the data for
testing the model’s quality post-training. Of the remaining 85 % of the data,
80 % was used for training and 20 % for validation, employing early stopping
with a patience setting of 30. To find the most suitable neural network
hyperparameters, we performed a grid search considering different activation
functions, number of layers, number of neurons per layer, and learning rates.
Using the independent testing dataset to discriminate between trained neural
networks, we selected an architecture of 2 hidden layers, 3 neurons per layer,
and the ReLU (rectified linear unit) activation function, trained with a
learning rate of 0.01. For this configuration, the parity plots comparing the
independent dataset with the neural network’s predictions yielded coefficients
of determination $R^{2}$ of 1.00 (Fig. 6). This indicates that the neural-
network surrogate was able to effectively learn the physics of the FBA model
for the metabolic exchange fluxes as a function of the manipulatable enzyme
concentration without overfitting.

(a) Biomass (b) Acetate

(c) Itaconate

Figure 6: Parity plots of relevant metabolic exchange flux values as predicted
by FBA against the values as predicted by the neural network (NN).

###  4.3 Evaluating design parameters of the optogenetic gene expression
system

With the neural-network surrogate now trained, our hybrid physics-informed
dynamic metabolic model is almost complete. The model still has the parameter
$\theta_{2}$ as a free design parameter of the optogenetic gene expression
system (cf. Eq. (10)). Therefore, to aid in the in silico design of the gene
expression system, we assessed the performance of various $\theta_{2}$ values
in open-loop optimization problems. That is, the optimization problem in (3.1)
is solved just once at $t_{0}$, without state feedback, and assuming no
mismatch between the plant’s and controller’s model. Specifically, we
considered $\theta_{2}$ values of $3.674\times 10^{-05}$, $3.306\times
10^{-05}$, $2.939\times 10^{-05}$, $2.572\times 10^{-05}$, $2.204\times
10^{-05}$, and $1.837\times 10^{-05}$ $\mathrm{mmol/g_{b}/h}$. The latter
values are within the range of orders of magnitude used in previous simulation
studies of gene expression 29. The Hill activation function
$q_{e_{\mathrm{cadA}}}$ for these different scenarios is plotted in Fig. 7.
The primary distinction among these scenarios is the maximum achievable
expression rate, effectively the strength of the gene expression system. The
objective function was set to maximize itaconate concentration at the end of
the batch, hence $J(\cdot)=z_{\mathrm{ita}}(t_{f})$. We also set the upper
bound of the input to $5\,\mathrm{W/m^{2}}$ as the function
$q_{e_{\mathrm{cadA}}}$ plateaus at this level for all values of $\theta_{2}$.
Note that we only optimized the dynamic light input, while the initial
conditions of the plant were a priori given. We considered 24 equidistant
inputs with an interval size of 1 h. Such discretization of the input function
aids in approximating the otherwise infinite-dimensional function 85. The
dynamic optimization using the hybrid physics-informed model was implemented
in HILO-MPC 52 interfaced with CasADi 86.

Figure 7: Hill activation function for the expression of the enzyme cis-
aconitate decarboxylase across different design parameter values of
$\theta_{2}$ in the considered optogenetic system (cf. Eq. (10)). $\theta_{2}$
values; OLO_1: $3.674\times 10^{-05}$, OLO_2: $3.306\times 10^{-05}$, OLO_3:
$2.939\times 10^{-05}$, OLO_4: $2.572\times 10^{-05}$, OLO_5: $2.204\times
10^{-05}$, and OLO_6: $1.837\times 10^{-05}$ $\mathrm{mmol/g_{b}/h}$.

The results of the open-loop optimizations are shown in Fig. 8. In all
scenarios, the light input follows a two-stage trajectory. Initially, the
light input is off, completely repressing the expression of cadA.
Subsequently, the light input increases to its maximum allowable value. This
aligns with an initial growth phase marked by maximal flux through the TCA
cycle, thus no itaconate production. It is followed by a gradual increase in
the intracellular concentration of enzyme cadA, which enhances itaconate
biosynthesis, at the expense of a reduced flux through the TCA cycle and
biomass synthesis. The difference between input profiles is the triggering
time for the second phase. The higher the value of $\theta_{2}$, the more
delayed the triggering of the production phase. That is, with the enzyme being
expressed at higher intracellular levels, an extended growth phase is
affordable with increasing $\theta_{2}$, allowing for significant biomass
accumulation for enhanced volumetric productivity. Conversely, with lower
values of $\theta_{2}$, less manipulatable enzyme accumulates intracellularly,
necessitating a longer production phase to somewhat compensate for the reduced
itaconate yields.

It is worth noting that the manipulatable enzyme concentration eventually
plateaus, approaching a pseudo-steady-state intracellular level, influenced by
the value of $q_{e_{\mathrm{cadA}}}$ (Fig. 7). To prove this point, let us
consider for simplicity a scenario where the manipulatable enzyme
concentration is such that it leads to the maximum biologically-feasible flux
of $v_{\mathrm{cadA}}$, hence $v_{\mathrm{bio}}=0$ and $\mu=0$. In this
specific scenario, the pseudo-steady-state concentration of the enzyme
$e_{\mathrm{cadA}_{\mathrm{ss}}}$ is:

$\displaystyle
e_{\mathrm{cadA}_{\mathrm{ss}}}:=\frac{q_{e_{\mathrm{cadA}}}}{d_{e_{\mathrm{cadA}}}}.$
(13)

Therefore, one can say that the higher the value of $q_{e_{\mathrm{cadA}}}$,
the higher the intracellular pseudo-steady-state concentration of the enzyme.
Note that the case with the maximum value of $\theta_{2}=3.674\times
10^{-05}\,\mathrm{mmol/g_{b}/h}$ closely aligns with the scenario described in
the derivation of Eq. (13) at prolonged induction levels. To illustrate this,
we have represented $e_{\mathrm{cadA}_{\mathrm{ss}}}$ with a dotted blue line
in Fig. 8-b for the scenario with the highest value of $\theta_{2}$. There,
the intracellular concentration of enzyme cadA approaches the calculated value
of $e_{\mathrm{cadA}_{\mathrm{ss}}}$ when the enzyme is fully induced.

(a) Light input (b) Manipulatable enzyme (c) Biomass (d) Glucose (e) Acetate
(f) Itaconate

Figure 8: Open-loop optimizations using the hybrid physics-informed dynamic
model with varying design parameter values of $\theta_{2}$ in the optogenetic
system (cf. Eq. (10)). It is assumed there is no mismatch between the
controller’s model and the plant’s model. $\theta_{2}$ values; OLO_1:
$3.674\times 10^{-05}$, OLO_2: $3.306\times 10^{-05}$, OLO_3: $2.939\times
10^{-05}$, OLO_4: $2.572\times 10^{-05}$, OLO_5: $2.204\times 10^{-05}$, and
OLO_6: $1.837\times 10^{-05}$ $\mathrm{mmol/g_{b}/h}$. Initial conditions:
$z_{\mathrm{glc}}(t_{0})=120\,\mathrm{mmol/L}$,
$z_{\mathrm{ita}}(t_{0})=0.00\,\mathrm{mmol/L}$,
$z_{\mathrm{ace}}(t_{0})=0.00\,\mathrm{mmol/L}$,
$b(t_{0})=0.05\,\mathrm{g_{b}/L}$, $e_{\mathrm{cadA}}(t_{0})=0.00\times
10^{-05}\,\mathrm{mmol/g_{b}}$. glc: glucose, ita: itaconate, bio: biomass,
ace: acetate, cadA: cis-aconitate decarboxylase.

In terms of biomass growth, only the scenario with the highest $\theta_{2}$
value exhibits a stagnation in biomass concentration towards the mid-to-end of
the batch. This is explained by the fact that the enzyme cadA accumulates to
levels such that $v_{\mathrm{bio}}\approx 0$ (cf. Fig. 5-a). In the other
scenarios, featuring decreasing values of $\theta_{2}$, the growth rate
becomes correspondingly more exponential, a trend explained by a lower
intracellular concentration of enzyme cadA and a more active TCA cycle. Since
the acetate metabolic flux is growth-coupled (see Fig. 5-a), the dynamic
behavior of acetate mirrors that of biomass. In terms of itaconate, the
product of interest, production only starts after the light input is switched
on. As expected, the rate of itaconate production correlates with the
intracellular concentration of enzyme cadA. Therefore, higher itaconate titers
are observed for the scenarios with higher enzyme cadA accumulation.
Naturally, enhanced titers are also enabled by the two-stage fermentation
process, where the initial growth-dominant stage facilitates sufficient
biomass accumulation to support the subsequent production-dominant phase.

Based on the predicted performance of the open-loop optimizations, we selected
the highest value of $\theta_{2}=3.674\times 10^{-05}\,\mathrm{mmol/g_{b}/h}$
as the most suitable design parameter value for the optogenetic gene
expression system. This value was chosen because it enables the coverage of
the entire space of growth-production trade-offs, as presented in Fig. 5, a
feature desirable for efficient dynamic metabolic control. We refer to the
described hybrid physics-informed model, with the specified parameters, as the
nominal model. It should be noted that the results in this section do not
account for system uncertainty, and therefore, the controller’s model is
assumed to precisely describe the plant’s dynamics. That is, the nominal model
applies to both the controller and the plant. In reality, bioprocesses are
frequently subject to uncertainties, including model-plant mismatch and
disturbances that need to be addressed with closed-loop control, which will be
the focus of the next section.

###  4.4 Putting hybrid physics-informed metabolic cybergenetics into action

We implemented a metabolic cybergenetic scheme to evaluate the performance of
the bioprocess under system uncertainty. This scheme follows the optimization
formulation presented in (3.1), which spans in the time domain from $t_{k}$ to
$t_{f}$. Samples are taken hourly to re-optimize the system using the current
state information at time $t_{k}$. As a source model-plant mismatch, we scaled
the function $h$ in Eq. (12) by a factor of 1.04. This model is henceforth
denoted as the modified model.

In the first feedback control case, denoted as MPC_1, we provided the
controller with the modified model and used the nominal model to simulate the
plant behavior. We also assumed full and accurate state information; that is,
we considered all states measurable, including the intracellular enzyme cadA,
without noise at every sampling instance. Control simulation results for MPC_1
are presented in Fig. 9. Note that we also plotted the open-loop control
performance, denoted as OLO_mis, for comparison; remark that this represents
the case without state feedback. Although MPC_1 is a very optimistic scenario,
it provides a benchmark for the performance of the model predictive
controller.

(a) Light input (b) Manipulatable enzyme (c) Biomass (d) Glucose (e) Acetate
(f) Itaconate

Figure 9: Model predictive control scenarios using the hybrid physics-informed
dynamic model. It is assumed that there is a mismatch between the controller’s
model and the plant’s model. To ensure a fair comparison between scenarios,
exact plant values are plotted. Initial conditions:
$z_{\mathrm{glc}}(t_{0})=120\,\mathrm{mmol/L}$,
$z_{\mathrm{ita}}(t_{0})=0.00\,\mathrm{mmol/L}$,
$z_{\mathrm{ace}}(t_{0})=0.00\,\mathrm{mmol/L}$,
$b(t_{0})=0.05\,\mathrm{g_{b}/L}$, $e_{\mathrm{cadA}}(t_{0})=0.00\times
10^{-05}\,\mathrm{mmol/g_{b}}$.

When analyzing MPC_1 against OLO_mis, the feedback controller effectively
adapts the dynamic light input to compensate for the model-plant mismatch. The
adjusted input facilitates greater biomass accumulation explained by a
slightly lower intracellular concentration of enzyme cadA after the onset of
the induction. As a result, there is a higher volumetric productivity,
resulting in a final itaconate titer of $86.64\,\mathrm{mmol/L}$ in MPC_1,
which is $14.32\,\%$ higher than the achieved with OLO_mis. Furthermore, MPC_1
depletes up to $98.2\,\%$ of the initial glucose, compared to only $83.2\,\%$
achieved with OLO_mis.

In a second feedback control case, denoted as MPC_2, we also provided the
controller with the modified model, while we used the nominal model to
simulate the plant behavior. However, unlike the external states, we
considered that the intracellular concentration of enzyme cadA could not be
measured. Therefore, we used a full information estimator, as described in the
optimization problem (3.2), to reconstruct the system state. We incorporated
measurement noise by sampling the measurements from a Gaussian distribution
with a 1.5 % standard deviation with respect to the true plant value. We set
$\mathbf{P}$ as a 5 $\times$ 5 diagonal matrix with diagonal elements equal to
10. Similarly, we set $\mathbf{R}$ as a 4 $\times$ 4 diagonal matrix with
diagonal elements equal to 1,000. These weights were tuned to provide
convergence and (to our judgment) good estimation quality. We omitted the term
associated with state disturbances from the objective function of the state
estimator as it did not make a difference in terms of the quality of the
outcome. The measured external states and the reconstructed manipulatable
enzyme were provided to the MPC controller at each sampling instance. Note
that, like the controller, the estimator was provided with the modified model.
The control simulation results for MPC_2 are also plotted in Fig. 9.

When comparing MPC_2 against OLO_mis, the feedback controller is also able to
compensate for the model-plant mismatch by adjusting the dynamic input. Like
MPC_1, MPC_2 enables more biomass accumulation due to the slightly less
intracellular concentration of enzyme cadA after the onset of the induction.
In fact, the performance of MPC_2 is comparable to MPC_1, yielding a final
itaconate concentration of $87.07\,\mathrm{mmol/L}$, which is $14.88\,\%$
higher than that achieved with OLO_mis. Unlike OLO_mis, MPC_2 is able to
deplete up to $98.9\,\%$ of the initial glucose. Regarding the state
estimator’s performance, as shown in Fig. 9-b, the estimator accurately
reconstructs the manipulatable enzyme concentration, closely matching the
exact plant values. Only the first two estimates are slightly off, which can
be explained by the few data points considered in those estimation instances.

##  5 Conclusion

This study outlines a hybrid physics-informed dynamic modeling framework for
predictive control and estimation in metabolic cybergenetics. We embed the
physics of metabolic networks into the process rates of macro-kinetic models
augmented with gene expression, which helps to account for intrinsic delays
related to transcription and translation of manipulatable enzymes. Such models
are built by employing machine-learning surrogates informed by FBA, wherein
manipulatable intracellular enzymes serve as input features and relevant
metabolic exchange fluxes as output labels. This approach ensures that the
information from the system’s metabolic network is not lost; rather, we
extract valuable physics from it to enhance structurally simpler dynamic
models. In comparison to previous metabolic cybergenetic frameworks,
particularly those utilizing dynamic constraint-based metabolic modeling 29,
our hybrid machine-learning-supported modeling approach circumvents complex
bilevel optimization problems. That is, process control and estimation tasks
can be formulated as single-level optimizations, thereby simplifying the
implementation of metabolic cybergenetics. Additionally, our framework reduces
the number of states that need to be measured or estimated in real-time during
the operation of the bioprocess due to the smaller number of states involved
in the models.

This work was supported by the U.S. Department of Energy, Office of Science,
Office of Biological and Environmental Research (Award Number DE-SC0022155).

## References

  * Liu and Nielsen 2019 Liu, Y.; Nielsen, J. Recent trends in metabolic engineering of microbial chemical factories. _Curr Opin Biotechnol_ 2019, _60_ , 188–197
  * Cho et al. 2022 Cho, J. S.; Kim, G. B.; Eun, H.; Moon, C. W.; Lee, S. Y. Designing microbial cell factories for the production of chemicals. _JACS Au_ 2022, _2_ , 1781–1799
  * Hartline et al. 2021 Hartline, C. J.; Schmitz, A. C.; Han, Y.; Zhang, F. Dynamic control in metabolic engineering: theories, tools, and applications. _Metab Eng_ 2021, _63_ , 126–140
  * Ni et al. 2021 Ni, C.; Dinh, C. V.; Prather, K. L. Dynamic control of metabolism. _Annu Rev Chem Biomol Eng_ 2021, _12_ , 519–541
  * Banerjee and Mukhopadhyay 2023 Banerjee, D.; Mukhopadhyay, A. Perspectives in growth production trade-off in microbial bioproduction. _RSC Sustain_ 2023, _1_ , 224–233
  * Lalwani et al. 2018 Lalwani, M. A.; Zhao, E. M.; Avalos, J. L. Current and future modalities of dynamic control in metabolic engineering. _Curr Opin Biotechnol_ 2018, _52_ , 56–65
  * Espinel-Ríos and Avalos 2023 Espinel-Ríos, S.; Avalos, J. L. Linking intra- and extra-cellular metabolic domains via neural-network surrogates for dynamic metabolic control. 2023, arXiv:2310.17179
  * Walls et al. 2023 Walls, M. T.; Xu, K.; Brangwynne, C. P.; Avalos, J. L. A modular design for synthetic membraneless organelles enables compositional and functional control. 2023, bioRxiv:10.1101/2023.10.03.560789
  * Wu et al. 2016 Wu, G.; Yan, Q.; Jones, J. A.; Tang, Y. J.; Fong, S. S.; Koffas, M. A. Metabolic burden: cornerstones in synthetic biology and metabolic engineering applications. _Trends Biotechnol_ 2016, _34_ , 652–664
  * Xu et al. 2023 Xu, J.; Wu, Z.; Zhao, K.; Xi, Z.; Wang, L.; Cheng, F.; Xue, Y.; Zheng, Y. IPTG-induced high protein expression for whole-cell biosynthesis of L-phosphinothricin. _Biotechnol J_ 2023, _18_ , 2300027
  * Vogl and Glieder 2013 Vogl, T.; Glieder, A. Regulation of Pichia pastoris promoters and its consequences for protein production. _New Biotechnol_ 2013, _30_ , 385–404
  * Harder et al. 2018 Harder, B.; Bettenbrock, K.; Klamt, S. Temperature-dependent dynamic control of the TCA cycle increases volumetric productivity of itaconic acid production by Escherichia coli. _Biotechnol Bioeng_ 2018, _115_ , 156–164
  * Hoynes-O’Connor et al. 2017 Hoynes-O’Connor, A.; Shopera, T.; Hinman, K.; Creamer, J. P.; Moon, T. S. Enabling complex genetic circuits to respond to extrinsic environmental signals. _Biotechnol Bioeng_ 2017, _114_ , 1626–1631
  * Wichmann et al. 2023 Wichmann, J.; Behrendt, G.; Boecker, S.; Klamt, S. Characterizing and utilizing oxygen-dependent promoters for efficient dynamic metabolic engineering. _Metab Eng_ 2023, _77_ , 199–207
  * Chisti 2006 Chisti, Y. In _Basic Biotechnology_ , 3rd ed.; Kristiansen, B., Ratledge, C., Eds.; Cambridge University Press: Cambridge, 2006; pp 181–200
  * Dvorak et al. 2015 Dvorak, P.; Chrast, L.; Nikel, P. I.; Fedr, R.; Soucek, K.; Sedlackova, M.; Chaloupkova, R.; de Lorenzo, V.; Prokop, Z.; Damborsky, J. Exacerbation of substrate toxicity by IPTG in Escherichia coli BL21(DE3) carrying a synthetic metabolic pathway. _Microb Cell Fact_ 2015, _14_ , 201
  * Lin et al. 2021 Lin, N.-X.; He, R.-Z.; Xu, Y.; Yu, X.-W. Oxidative stress tolerance contributes to heterologous protein production in Pichia pastoris. _Biotechnol Biofuels_ 2021, _14_ , 160
  * Maurer et al. 2005 Maurer, L. M.; Yohannes, E.; Bondurant, S. S.; Radmacher, M.; Slonczewski, J. L. pH regulates genes for flagellar motility, catabolism, and oxidative stress in Escherichia coli K-12. _J Bacteriol_ 2005, _187_ , 304–319
  * Gadgil et al. 2008 Gadgil, M.; Kapur, V.; Hu, W.-S. Transcriptional response of Escherichia coli to temperature shift. _Biotechnol Prog_ 2008, _21_ , 689–699
  * Bettenbrock et al. 2014 Bettenbrock, K.; Bai, H.; Ederer, M.; Green, J.; Hellingwerf, K. J.; Holcombe, M.; Kunz, S.; Rolfe, M. D.; Sanguinetti, G.; Sawodny, O.; Sharma, P.; Steinsiek, S.; Poole, R. K. _Advances in Microbial Physiology_ ; Elsevier, 2014; Vol. 64; pp 65–114
  * Tandar et al. 2019 Tandar, S. T.; Senoo, S.; Toya, Y.; Shimizu, H. Optogenetic switch for controlling the central metabolic flux of Escherichia coli. _Metab Eng_ 2019, _55_ , 68–75
  * Carrasco-López et al. 2020 Carrasco-López, C.; García-Echauri, S. A.; Kichuk, T.; Avalos, J. L. Optogenetics and biosensors set the stage for metabolic cybergenetics. _Curr Opin Biotechnol_ 2020, _65_ , 296–309
  * Pouzet et al. 2020 Pouzet, S.; Banderas, A.; Le Bec, M.; Lautier, T.; Truan, G.; Hersen, P. The promise of optogenetics for bioproduction: dynamic control strategies and scale-up instruments. _Bioengineering_ 2020, _7_ , 151
  * Zhao et al. 2021 Zhao, E. M.; Lalwani, M. A.; Chen, J.-M.; Orillac, P.; Toettcher, J. E.; Avalos, J. L. Optogenetic amplification circuits for light-induced metabolic control. _ACS Synth Biol_ 2021, _10_ , 1143–1154
  * Lalwani et al. 2021 Lalwani, M. A.; Ip, S. S.; Carrasco-López, C.; Day, C.; Zhao, E. M.; Kawabe, H.; Avalos, J. L. Optogenetic control of the lac operon for bacterial chemical and protein production. _Nat Chem Biol_ 2021, _17_ , 71–79
  * Hoffman et al. 2022 Hoffman, S. M.; Tang, A. Y.; Avalos, J. L. Optogenetics illuminates applications in microbial engineering. _Annu Rev Chem Biomol_ 2022, _13_ , 373–403
  * Wegner et al. 2022 Wegner, S. A.; Barocio-Galindo, R. M.; Avalos, J. L. The bright frontiers of microbial metabolic optogenetics. _Curr Opin Chem Biol_ 2022, _71_ , 102207
  * Espinel-Ríos et al. 2022 Espinel-Ríos, S.; Morabito, B.; Pohlodek, J.; Bettenbrock, K.; Klamt, S.; Findeisen, R. Optimal control and dynamic modulation of the ATPase gene expression for enforced ATP wasting in batch fermentations. _IFAC-PapersOnLine_ 2022, _55_ , 174–180
  * Espinel-Ríos et al. 2024 Espinel-Ríos, S.; Morabito, B.; Pohlodek, J.; Bettenbrock, K.; Klamt, S.; Findeisen, R. Toward a modeling, optimization, and predictive control framework for fed-batch metabolic cybergenetics. _Biotechnol Bioeng_ 2024, _121_ , 366–379
  * Ohkubo et al. 2023 Ohkubo, T.; Soma, Y.; Sakumura, Y.; Hanai, T.; Kunida, K. A hybrid in silico/in-cell controller for microbial bioprocesses with process-model mismatch. _Sci Rep_ 2023, _13_ , 13608
  * Waldherr et al. 2015 Waldherr, S.; Oyarzún, D. A.; Bockmayr, A. Dynamic optimization of metabolic networks coupled with gene expression. _J Theor Biol_ 2015, _365_ , 469–485
  * Jabarivelisdeh and Waldherr 2018 Jabarivelisdeh, B.; Waldherr, S. Optimization of bioprocess productivity based on metabolic-genetic network models with bilevel dynamic programming. _Biotechnol Bioeng_ 2018, _115_ , 1829–1841
  * Liu and Bockmayr 2020 Liu, L.; Bockmayr, A. Regulatory dynamic enzyme-cost flux balance analysis: a unifying framework for constraint-based modeling. _J Theor Biol_ 2020, _501_ , 110317
  * Jabarivelisdeh et al. 2020 Jabarivelisdeh, B.; Carius, L.; Findeisen, R.; Waldherr, S. Adaptive predictive control of bioprocesses with constraint-based modeling and estimation. _Comput Chem Eng_ 2020, _135_ , 106744
  * Espinel-Ríos et al. 2022 Espinel-Ríos, S.; Bettenbrock, K.; Klamt, S.; Findeisen, R. Maximizing batch fermentation efficiency by constrained model-based optimization and predictive control of adenosine triphosphate turnover. _AIChE J_ 2022, _68_ , e17555
  * Dempe et al. 2015 Dempe, S.; Kalashnikov, V.; Pérez-Valdés, G. A.; Kalashnykova, N. _Bilevel Programming Problems_ ; Springer Berlin Heidelberg: Berlin, Heidelberg, 2015; pp 41–115, Series Title: Energy Systems
  * Dempe and Franke 2019 Dempe, S.; Franke, S. Solution of bilevel optimization problems using the KKT approach. _Optimization_ 2019, _68_ , 1471–1489
  * Sinha et al. 2019 Sinha, A.; Soun, T.; Deb, K. Using Karush-Kuhn-Tucker proximity measure for solving bilevel optimization problems. _Swarm Evol Comput_ 2019, _44_ , 496–510
  * Rawlings et al. 2020 Rawlings, J.; Mayne, D.; Diehl, M. _Model predictive control: theory, computation and design_ , 2nd ed.; Nob Hill Publishing, LLC: Santa Barbara, 2020
  * Orth et al. 2010 Orth, J. D.; Thiele, I.; Palsson, B. Ø. What is flux balance analysis? _Nat Biotechnol_ 2010, _28_ , 245–248
  * Klamt et al. 2014 Klamt, S.; Hädicke, O.; Von Kamp, A. In _Large-Scale Networks in Engineering and Life Sciences_ ; Benner, P., Findeisen, R., Flockerzi, D., Reichl, U., Sundmacher, K., Eds.; Springer International Publishing: Cham, 2014; pp 263–316
  * Klamt et al. 2018 Klamt, S.; Müller, S.; Regensburger, G.; Zanghellini, J. A mathematical framework for yield (vs. rate) optimization in constraint-based modeling and applications in metabolic engineering. _Metab Eng_ 2018, _47_ , 153–169
  * Teixeira et al. 2005 Teixeira, A.; Cunha, A.; Clemente, J.; Moreira, J.; Cruz, H.; Alves, P.; Carrondo, M.; Oliveira, R. Modelling and optimization of a recombinant BHK-21 cultivation process using hybrid grey-box systems. _J Biotechnol_ 2005, _118_ , 290–303
  * Marques et al. 2017 Marques, R.; Von Stosch, M.; Portela, R. M.; Torres, C. A.; Antunes, S.; Freitas, F.; Reis, M. A.; Oliveira, R. Hybrid modeling of microbial exopolysaccharide (EPS) production: the case of Enterobacter A47. _J Biotechnol_ 2017, _246_ , 61–70
  * Espinel-Ríos et al. 2023 Espinel-Ríos, S.; Bettenbrock, K.; Klamt, S.; Avalos, J. L.; Findeisen, R. _Comput Aided Chem Eng_ ; Elsevier, 2023; Vol. 52; pp 2601–2606
  * Rogers et al. 2023 Rogers, A. W.; Song, Z.; Ramon, F. V.; Jing, K.; Zhang, D. Investigating ‘greyness’ of hybrid model for bioprocess predictive modelling. _Biochem Eng J_ 2023, _190_ , 108761
  * Pinto et al. 2023 Pinto, J.; Ramos, J. R. C.; Costa, R. S.; Rossell, S.; Dumas, P.; Oliveira, R. Hybrid deep modeling of a CHO-K1 fed-batch process: combining first-principles with deep neural networks. _Front Bioeng Biotechnol_ 2023, _11_ , 1237963
  * De Oliveira et al. 2021 De Oliveira, R. D.; Guedes, M. N.; Matias, J.; Le Roux, G. A. C. Nonlinear predictive control of a bioreactor by surrogate model approximation of flux balance analysis. _Ind Eng Chem Res_ 2021, _60_ , 14464–14475
  * Proshkin et al. 2010 Proshkin, S.; Rahmouni, A. R.; Mironov, A.; Nudler, E. Cooperation between translating ribosomes and RNA polymerase in transcription elongation. _Science_ 2010, _328_ , 504–508
  * Rodríguez-Molina et al. 2023 Rodríguez-Molina, J. B.; West, S.; Passmore, L. A. Knowing when to stop: transcription termination on protein-coding genes by eukaryotic RNAPII. _Mol Cell_ 2023, _83_ , 404–415
  * Zhang et al. 2019 Zhang, D.; Del Rio-Chanona, E. A.; Petsagkourakis, P.; Wagner, J. Hybrid physics-based and data-driven modeling for bioprocess online simulation and optimization. _Biotechnol Bioeng_ 2019, _116_ , 2919–2930
  * Pohlodek et al. 2022 Pohlodek, J.; Morabito, B.; Schlauch, C.; Zometa, P.; Findeisen, R. Flexible development and evaluation of machine-learning-supported optimal control and estimation methods via HILO-MPC. 2022, arXiv:2203.13671
  * Hornik et al. 1989 Hornik, K.; Stinchcombe, M.; White, H. Multilayer feedforward networks are universal approximators. _Neural Netw_ 1989, _2_ , 359–366
  * Widder et al. 2007 Widder, S.; Schicho, J.; Schuster, P. Dynamic patterns of gene regulation I: Simple two-gene systems. _J Theor Biol_ 2007, _246_ , 395–419
  * Santillán 2008 Santillán, M. On the Use of the Hill Functions in Mathematical Models of Gene Regulatory Networks. _Math Model Nat Phenom_ 2008, _3_ , 85–97
  * Ang et al. 2013 Ang, J.; Harris, E.; Hussey, B. J.; Kil, R.; McMillen, D. R. Tuning response curves for synthetic biology. _ACS Synth Biol_ 2013, _2_ , 547–567
  * Mosteller and Goldstein 1984 Mosteller, R. D.; Goldstein, B. E. A mathematical model that applies to protein degradation and post-translational processing of proteins and to analogous processes for other molecules in non-growing and exponentially growing cells. _J Theor Biol_ 1984, _108_ , 597–621
  * Subramanian and Srienc 1996 Subramanian, S.; Srienc, F. Quantitative analysis of transient gene expression in mammalian cells using the green fluorescent protein. _J Biotechnol_ 1996, _49_ , 137–151
  * Yin and Redovich 2018 Yin, J.; Redovich, J. Kinetic modeling of virus growth in cells. _Microbiol Mol Biol_ 2018, _82_ , e00066–17
  * Bennett et al. 2009 Bennett, B. D.; Kimball, E. H.; Gao, M.; Osterhout, R.; Van Dien, S. J.; Rabinowitz, J. D. Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli. _Nat Chem Biol_ 2009, _5_ , 593–599
  * Passi et al. 2021 Passi, A.; Tibocha-Bonilla, J. D.; Kumar, M.; Tec-Campos, D.; Zengler, K.; Zuniga, C. Genome-scale metabolic modeling enables in-depth understanding of big data. _Metabolites_ 2021, _12_ , 14
  * Biechele et al. 2015 Biechele, P.; Busse, C.; Solle, D.; Scheper, T.; Reardon, K. Sensor systems for bioprocess monitoring. _Eng Life Sci_ 2015, _15_ , 469–488
  * Noll and Henkel 2020 Noll, P.; Henkel, M. History and evolution of modeling in biotechnology: modeling & simulation, application and hardware performance. _Comput Struct Biotechnol J_ 2020, _18_ , 3309–3323
  * Lowe et al. 2017 Lowe, R.; Shirley, N.; Bleackley, M.; Dolan, S.; Shafee, T. Transcriptomics technologies. _PLoS Comput Biol_ 2017, _13_ , e1005457
  * Shuken 2023 Shuken, S. R. An introduction to mass spectrometry-based proteomics. _J Proteome Res_ 2023, _22_ , 2151–2171
  * Su 2005 Su, W. W. Fluorescent proteins as tools to aid protein production. _Microb Cell Fact_ 2005, _4_ , 12
  * Lee et al. 2013 Lee, S.; Lim, W. A.; Thorn, K. S. Improved blue, green, and red fluorescent protein tagging vectors for S. cerevisiae. _PLoS ONE_ 2013, _8_ , e67902
  * Hernández Rodríguez et al. 2021 Hernández Rodríguez, T.; Posch, C.; Pörtner, R.; Frahm, B. Dynamic parameter estimation and prediction over consecutive scales, based on moving horizon estimation: applied to an industrial cell culture seed train. _Bioprocess Biosyst Eng_ 2021, _44_ , 793–808
  * Elsheikh et al. 2021 Elsheikh, M.; Hille, R.; Tatulea-Codrean, A.; Krämer, S. A comparative review of multi-rate moving horizon estimation schemes for bioprocess applications. _Comput Chem Eng_ 2021, _146_ , 107219
  * Espinel-Ríos et al. 2022 Espinel-Ríos, S.; Huber, N.; Alcalá-Orozco, E. A.; Morabito, B.; Rexer, T. F.; Reichl, U.; Klamt, S.; Findeisen, R. Cell-free biosynthesis meets dynamic optimization and control: a fed-batch framework. _IFAC-PapersOnLine_ 2022, _55_ , 92–97
  * Espinel-Ríos et al. 2022 Espinel-Ríos, S.; Morabito, B.; Bettenbrock, K.; Klamt, S.; Findeisen, R. Soft sensor for monitoring dynamic changes in cell composition. _IFAC-PapersOnLine_ 2022, _55_ , 98–103
  * Tuveri et al. 2023 Tuveri, A.; Nakama, C. S.; Matias, J.; Holck, H. E.; Jäschke, J.; Imsland, L.; Bar, N. A regularized moving horizon estimator for combined state and parameter estimation in a bioprocess experimental application. _Comput Chem Eng_ 2023, _172_ , 108183
  * Espinel-Ríos et al. 2023 Espinel-Ríos, S.; Slaviero, G.; Bettenbrock, K.; Klamt, S.; Findeisen, R. Monitoring intracellular metabolite concentrations by moving horizon estimation based on kinetic modeling. _IFAC-PapersOnLine_ 2023, _56_ , 4608–4613
  * Schmidl et al. 2014 Schmidl, S. R.; Sheth, R. U.; Wu, A.; Tabor, J. J. Refactoring and optimization of light-switchable Escherichia coli two-component systems. _ACS Synth Biol_ 2014, _3_ , 820–831
  * Milo and Phillips 2016 Milo, R.; Phillips, R. _Cell biology by the numbers_ ; Garland science: New York, 2016
  * Senn et al. 1994 Senn, H.; Lendenmann, U.; Snozzi, M.; Hamer, G.; Egli, T. The growth of Escherichia coli in glucose-limited chemostat cultures: a re-examination of the kinetics. _Biochim Biophys Acta Gen Subj_ 1994, _1201_ , 424–436
  * Virtanen et al. 2020 Virtanen, P. et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. _Nat Methods_ 2020, _17_ , 261–272
  * Pinhal et al. 2019 Pinhal, S.; Ropers, D.; Geiselmann, J.; De Jong, H. Acetate metabolism and the inhibition of bacterial growth by acetate. _J Bacteriol_ 2019, _201_
  * Chen et al. 2019 Chen, F.; Lukat, P.; Iqbal, A. A.; Saile, K.; Kaever, V.; Van Den Heuvel, J.; Blankenfeldt, W.; Büssow, K.; Pessler, F. Crystal structure of cis -aconitate decarboxylase reveals the impact of naturally occurring human mutations on itaconate synthesis. _Proc Natl Acad Sci USA_ 2019, _116_ , 20644–20654
  * Harder et al. 2016 Harder, B.-J.; Bettenbrock, K.; Klamt, S. Model-based metabolic engineering enables high yield itaconic acid production by Escherichia coli. _Metab Eng_ 2016, _38_ , 29–37
  * Hädicke and Klamt 2017 Hädicke, O.; Klamt, S. EColiCore2: a reference network model of the central metabolism of Escherichia coli and relationships to its genome-scale parent model. _Sci Rep._ 2017, _7_ , 39647
  * Orth et al. 2011 Orth, J. D.; Conrad, T. M.; Na, J.; Lerman, J. A.; Nam, H.; Feist, A. M.; Palsson, B. O. A comprehensive genome-scale reconstruction of Escherichia coli metabolism—2011. _Mol Syst Biol_ 2011, _7_ , 535
  * Ebrahim et al. 2013 Ebrahim, A.; Lerman, J. A.; Palsson, B. O.; Hyduke, D. R. COBRApy: COnstraints-Based Reconstruction and Analysis for Python. _BMC Syst Biol_ 2013, _7_ , 74
  * Paszke et al. 2019 Paszke, A. et al. _Proceedings of the 33rd International Conference on Neural Information Processing Systems_ ; Curran Associates Inc.: Red Hook, NY, USA, 2019; pp 8026–8037
  * Findeisen and Allgöwer 2002 Findeisen, R.; Allgöwer, F. An introduction to nonlinear model predictive control. 21st Benelux Meeting on Systems and Control. 2002; pp 119–141
  * Andersson et al. 2019 Andersson, J. A. E.; Gillis, J.; Horn, G.; Rawlings, J. B.; Diehl, M. CasADi: a software framework for nonlinear optimization and optimal control. _Math Program Comput_ 2019, _11_ , 1–36