# Simultaneous design of fermentation and microbe

Anita L. Ziegler  Process Systems Engineering (AVT.SVT), RWTH Aachen
University, 52074 Aachen, Germany  Ashutosh Manchanda  Process Systems
Engineering (AVT.SVT), RWTH Aachen University, 52074 Aachen, Germany  Marc-
Daniel Stumm  Process Systems Engineering (AVT.SVT), RWTH Aachen University,
52074 Aachen, Germany  Lars M. Blank  Institute of Applied Microbiology
(iAMB), Aachen Biology and Biotechnology (ABBt), RWTH Aachen University, 52074
Aachen, Germany  Alexander Mitsos  JARA-ENERGY, 52056 Aachen, Germany  Process
Systems Engineering (AVT.SVT), RWTH Aachen University, 52074 Aachen, Germany
Institute of Energy and Climate Research, Energy Systems Engineering (IEK-10),
Forschungszentrum Jülich GmbH, 52425 Jülich, Germany  Corresponding author:
Alexander Mitsos<EMAIL_ADDRESS>
## Abstract

Constraint-based optimization of microbial strains and model-based bioprocess
design have been used extensively to enhance yields in biotechnological
processes. However, strain and process optimization are usually carried out in
sequential steps, causing underperformance of the biotechnological process
when scaling up to industrial fermentation conditions. Herein, we propose the
optimization formulation _SimulKnock_ that combines the optimization of a
fermentation process with metabolic network design in a bilevel optimization
program. The upper level maximizes space-time yield and includes mass balances
of a continuous fermentation, while the lower level is based on flux balance
analysis. SimulKnock predicts optimal gene deletions and finds the optimal
trade-off between growth rate and product yield. Results of a case study with
a genome-scale metabolic model of E. coli indicate higher space-time yields
than a sequential approach using OptKnock for almost all target products
considered. By leveraging SimulKnock, we reduce the gap between strain and
process optimization.

## Keywords

process optimization, computational strain design, constraint-based metabolic
modeling, metabolic engineering

##  1 Introduction

Industrial microbiology promises product synthesis from renewable feedstocks
representing a sustainable alternative to petro-chemical synthesis [1]. Often,
these biotechnological products are produced by high-performing microbial
strains that have been designed using metabolic engineering. Strain design is
supported by computational methods to reduce the experimental effort [2]. In
the first place, constraint-based metabolic modeling formulations predict the
influence of a given genetic modification on the microorganism [3]. The best-
known formulation is flux balance analysis (FBA) [4, 5]. Using linear
programming, FBA predicts the internal fluxes of a microorganism based on its
genome-scale metabolic model (GEM). Variations of FBA refine the cellular
objective of the organism [6] or better suit a genetically modified
microorganism [7, 8, 9]. Constraint-based strain optimization formulations go
one step further and predict targets for genetic modification [10, 3]. The
first formulation, OptKnock, was presented by Burgard et al. [11]. OptKnock
proposes optimal gene knockouts by solving a bilevel optimization problem. The
upper level represents the bioengineering perspective to maximize the product
yield on the substrate. The lower level, based on FBA, represents the
microorganism with the cellular objective of maximizing biomass production.
OptKnock has been extended, e.g., to account for worst-case predictions [12],
insertion of genes [13], and up and down regulation of genes [14]. It was also
modified to suit a genetically modified organism better [15]. For strain
design, the predicted modifications are experimentally tested and evaluated at
a laboratory scale, and the best-performing strain is chosen.

For industrial production, industrial-scale process conditions come into play.
Process design involves the analysis of mass and energy balances an analyzing
production costs, including up- and downstream processing [16]. In bioprocess
design, a suitable kinetic represents the microorganism, e.g., its growth [17,
18]. Computational process optimization supports process design by adjusting
the process parameters to achieve maximal space-time yield or minimal
production cost [18, 19, 20].

Typically, strain and process design are separate steps, performed in
sequence. Process conditions, however, influence the behavior and performance
of the microorganism. Hence, a high-performing strain in the laboratory may
underperform in an industrial process. For example, a strain is first designed
under small-scale batch fermentation conditions in the laboratory. Next, the
industrial process is designed for the designated strain, whereby industrial
processes are envisaged to run on a large scale in continuous mode, and
example processes already exist in the pharmaceutical and food sectors. This
switch in conditions entails a significant adjustment for the organism.
Reasons for the under-performance can be steep gradients in the spatial
distribution of substrates (e.g., oxygen) in the fermenter, high downstream
processing costs, or lacking genetic and phenotypic stability by the microbe
[21, 22]. Strain and process design coupling is needed to overcome these
scale-up difficulties [22, 23]. There exist different approaches to do so. In
their attempt to capture spatial changes in a stirred bioreactor, Lapin et al.
[24] coupled computational fluid dynamics and metabolic modeling. To model
diauxic growth, Mahadevan et al. [25] connected the mass balance equations of
batch fermentation with metabolic modeling in their dynamic flux balance
analysis (dFBA). Structurally, dFBA are differential‐algebraic equations with
embedded optimization criteria. For fed-batch fermentation, de Oliveira et al.
[26] replaced the metabolic network in dFBA with a surrogate model. In an
iterative process, Zhuang et al. [27] connected dFBA with computational strain
optimization techniques to suit the organism to the future process conditions
of batch fermentation. Ploch et al. [28] developed a differential-algebraic
equation system with embedded dFBA optimization to model a biorefinery process
under changing conditions. Jabarivelisdeh and Waldherr [29] introduced a
bilevel dynamic optimization framework, where a previously found genetic
modification can be switched on during batch fermentation. In a step-wise
approach, Tafur Rangel et al. [30] selected a microbial strain regarding the
required downstream process units. Konstantinos and Antonis [31] considered
downstream process optimization simultaneous to strain optimization, namely in
a superstructure optimization program. However, the implications of the switch
from laboratory (batch or fed-batch) fermentation conditions to industrial
(continuous) fermentation conditions have yet to be considered.

We propose a new formulation combining strain optimization with process
optimization. Our simultaneous approach _SimulKnock_ , suggests optimal gene
knockouts that suit a strain to its future industrial-scale conditions, and at
the same time, the process conditions are customized to the strain. SimulKnock
is a bilevel optimization program where the upper level maximizes the space-
time yield subject to the mass balances of fermentation, and the lower level
maximizes biomass production based on the FBA formulation. In this paper, we
consider continuous fermentation, an envisaged fermentation mode in industry.
SimulKnock, however, can readily be extended to account for other fermentation
modes, too. Certainly, extensions may result in different optimization
classes, e.g., simultaneous optimization of batch operation and strain design
would be a dynamic bilevel optimization problem, which is challenging [28,
32]. Michaelis-Menten or Monod kinetics are employed to connect the two
levels. Compared to existing strain optimization formulations, SimulKnock can
be interpreted as an extension of OptKnock by process equations or as a
variation of dFBA with a switch to continuous fermentation and implementation
of a process objective and knockout predictions. Similarly, SimulKnock can be
seen as an extension of process optimization by the strain optimization level
(lower level). For our case studies, we use strong dualization to reformulate
the bilevel optimization program into a single-level program and solve it
globally. We apply our formulation to a GEM of E. coli for the production of
ethanol, succinate, acetate, formate, fumarate, and lactate and compare the
results with sequential optimization results and experimental results for
succinate production from the literature. The proposed formulation SimulKnock
will help close the gap in profitability between petrochemical and biochemical
processes.

##  2 Method: SimulKnock

###  2.1 Mass balances for a continuous process optimization

For optimal process design, the fermentation process was modeled with
mathematical equations. These equations may include mass and energy balances
of different reactor types and process units. In SimulKnock, we chose the mass
balance equations of a continuous stirred-tank reactor. In principle, however,
SimulKnock can easily be extended to account for other reactor types or even
full processes. We assumed that the reactor is ideally mixed and contains only
one homogeneous phase, that only one substrate limits growth and production,
and that only one population of organisms is in the tank. Furthermore, we
assumed that the volumetric inflow equals the volumetric outflow. The biomass,
the substrate, and the product were chosen as the representative compounds of
the fermentation. Similarly to standard literature, mass conservation leads to

$\displaystyle\frac{\mathrm{d}c_{bio}}{\mathrm{d}t}$
$\displaystyle=(\mu-D)c_{bio}$ (1)
$\displaystyle\frac{\mathrm{d}c_{S}}{\mathrm{d}t}$
$\displaystyle=-\frac{1}{Y_{bio/S}}\mu
c_{bio}-\frac{1}{Y_{P/S}}q_{p}c_{bio}-m_{S}c_{bio}+D(c_{S,Feed}-c_{S})$
$\displaystyle\frac{\mathrm{d}c_{P}}{\mathrm{d}t}$
$\displaystyle=q_{P}c_{bio}-c_{P}\cdot D,$

where $D$ is the dilution rate in ${\mathrm{h}}^{-1}$; $c_{bio}$ is the
biomass concentration in gram Cell Dry Weight (CDW) per liter
($\mathrm{g_{CDW}}\text{\,}{\mathrm{L}}^{-1}$); $c_{S,Feed}$, $c_{S}$ and
$c_{P}$ are the substrate feed concentration, substrate concentration, and
product concentration, respectively (all in
$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$). $Y_{bio/S}$ and $Y_{P/S}$, both with
the unit $\mathrm{g}\text{\,}{\mathrm{g}}^{-1}$, are the yields of biomass and
product per substrate, respectively. The maintenance factor is denoted as
$m_{S}$ in $\mathrm{g}\text{\,}{\mathrm{g}}^{-1}\text{\,}{\mathrm{h}}^{-1}$
and describes the required substrate uptake rate for cellular maintenance. Two
kinetics were used in the equations: the product kinetics $q_{P}$ in
$\mathrm{g}\text{\,}{\mathrm{g}}^{-1}\text{\,}{\mathrm{h}}^{-1}$, describing
the rate of formation of the product, and the growth rate kinetics $\mu$ in
${\mathrm{h}}^{-1}$, describing the rate of formation of the biomass. The
dilution rate $D$ is the quotient of the volumetric inflow/outflow and the
culture volume within the bioreactor. The dilution rate is also the inverse of
the residence time. As is typical, we assumed that a flow equilibrium is
achieved after a sufficiently long period of operation at a constant dilution
rate. It follows that the continuous fermentation is at a steady state, such
that $\frac{\mathrm{d}c_{bio}}{\mathrm{d}t}=0$,
$\frac{\mathrm{d}c_{S}}{\mathrm{d}t}=0$, and
$\frac{\mathrm{d}c_{P}}{\mathrm{d}t}=0$ apply. Thereby, it directly follows
that $D=\mu$, which means that biomass is neither washed out nor accumulated
in the process.

###  2.2 Flux balance analysis

FBA analyzes the internal fluxes within a cell using linear programming [5,
4]. The analysis is based on the metabolic network of an organism that
includes the metabolites and the stoichiometry of the internal reactions. In
FBA, based on the steady state assumption, the metabolites’ mass balances are
constrained to zero, and the reaction fluxes are constrained to upper and
lower bounds. One cellular objective is defined, e.g., the maximization of the
biomass flux. The mathematical formulation of FBA reads

$\displaystyle\max_{\bm{v}\in\mathbb{R}^{n}}$ $\displaystyle v_{bio}$ s.t.
$\displaystyle\bm{Sv}=\bm{0}$
$\displaystyle\bm{v^{lower}}\leq\bm{v}\leq\bm{v^{upper}}$

where $\bm{v}$ denotes the vector of reaction fluxes in
$\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$ and
where $n$ is the number of reactions in a metabolic network, where all
reversible reactions were split into one forward and one backward reaction,
further called irreversible network. The biomass flux is denoted by $v_{bio}$
in ${\mathrm{h}}^{-1}$ and is an element of $\bm{v}$; $\bm{S}$ is the
stoichiometric matrix of the irreversible network, and $lower$ and $upper$ are
lower and upper bounds of the flux values, respectively. These bounds may
include thresholds on the biomass flux or the ATP maintenance reaction.

###  2.3 Proposed combined optimization: SimulKnock

We designed SimulKnock by embedding the FBA into the optimization of a
continuous fermentation accounting for cellular degrees of freedom, resulting
in a bilevel optimization formulation (Figure 1), following the paradigm of
OptKnock [11]. Using the steady-state mass balances of metabolites stemming
from FBA, we assume that the cellular metabolism adapts infinitesimally
quickly to shifting environmental conditions [33].

Figure 1: The bilevel optimization formulation of SimulKnock. Note that either
Monod or Michaelis-Menten kinetics are applied at the upper level or the lower
level.

SimulKnock maximizes the space-time yield in the upper-level program like
other bioprocess optimization examples [34, 35]. The connection between upper-
and lower-level program is achieved via the expression of process parameters
with metabolic variables, based on Ploch et al. [28] and Mahadevan et al.
[25]. The growth rate directly transforms to $\mu=v_{bio}$, and due to steady-
state continuous fermentation conditions, it holds $D=v_{bio}$. Moreover,
instead of considering that everything taken up goes into product formation,
biomass formation, or maintenance, the substrate uptake is now directly
expressed by $-\frac{1}{Y_{X/S}}\mu-\frac{1}{Y_{P/S}}q_{p}-m_{S}=-v_{S}\cdot
M_{S}$, where $v_{S}$ is the substrate uptake flux in
$\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{L}}^{-1}$ and
$M_{S}$ denotes the molar mass of the substrate in
$\mathrm{g}\text{\,}{\mathrm{mmol}}^{-1}$. Thus, the maintenance is only
directly considered at the lower level in SimulKnock by setting a threshold on
the ATP maintenance reaction; there is no doubling at the upper level. The
rate of product formation can also be directly considered now by setting
$q_{P}=v_{P}\cdot M_{P}$, where $v_{P}$ is the product flux in
$\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{L}}^{-1}$ and
is an element of the flux vector $\bm{v}$ and $M_{P}$ is the molar mass of the
product in $\mathrm{g}\text{\,}{\mathrm{mmol}}^{-1}$. Hence, (1) transforms to

$\displaystyle D$ $\displaystyle=v_{bio}$ $\displaystyle 0$
$\displaystyle=-v_{S}\cdot M_{S}\cdot c_{bio}+D\cdot(c_{S,f}-c_{S})$
$\displaystyle 0$ $\displaystyle=v_{P}\cdot M_{P}\cdot c_{bio}-c_{P}\cdot D.$

Note that instead of fixing the dilution rate in the upper level and, thereby,
fixing the growth rate, the dilution rate $D$ will be set according to the
optimal value of $v_{bio}$ after the lower-level optimization. Hence, in the
following, $D$ will be replaced with $v_{bio}$.

Furthermore, we considered two alternative kinetics to connect the upper and
lower levels: Monod and Michaelis-Menten. In Section A.1 in the Supporting
Information, we discuss why including both kinetics simultaneously is not
advisable.

Monod is an empirical, widely known, and easy-to-measure model of microbial
growth. It links the substrate concentration with the growth rate and reads

$\displaystyle v_{bio}$
$\displaystyle=v_{bio}^{max}\frac{c_{S}}{c_{S}+K_{S}}\xRightarrow{}c_{S}=K_{S}\frac{v_{bio}}{v_{bio}^{max}-v_{bio}},$
$\displaystyle(Monod)$ (2)

where $v_{bio}^{max}$ denotes the maximum growth rate in ${\mathrm{h}}^{-1}$,
and $K_{S}$ denotes the Monod affinity constant in
$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$ for the specified substrate. We assumed
that gene deletions would not affect the kinetic parameters, i.e., we assumed
that the kinetic parameters were constants for a given organism and substrate.
Thus, both $v_{bio}^{max}$ and $K_{S}$ are parameters taken from literature
or, in the case of $v_{bio}^{max}$, can be approximated with FBA. The growth
rate results from the lower-level program, and thus, the Monod kinetics need
to be imposed on the upper-level program (reactor level). As depicted on the
right-hand side in (2), we, therefore, reformulated the kinetics such that the
substrate concentration is a function of the growth rate: $c_{S}=f(v_{bio})$.
Then, we implemented the kinetics in the upper level by replacing $c_{S}$ with
the kinetic term. If we would keep $c_{S}$, the growth rate would be fixed by
the kinetic term (equality constraint) instead of being optimized at the lower
level. SimulKnock with included Monod kinetics reads

$\displaystyle\underset{\bm{y},c_{S,Feed},c_{P},c_{bio},v_{S}^{\prime},v_{bio}^{\prime},v_{P}^{\prime}}{\max}$
$\displaystyle\quad c_{P}\cdot v_{bio}^{\prime}$ (3) s.t.
$\displaystyle\sum_{i=1}^{r}(1-y_{i})$ $\displaystyle\leq K$ $\displaystyle 0$
$\displaystyle=-v_{S}^{\prime}\cdot M_{S}\cdot
c_{bio}+v_{bio}^{\prime}\cdot(c_{S,Feed}-K_{S}\frac{v_{bio}^{\prime}}{v_{bio}^{max}-v_{bio}^{\prime}})$
$\displaystyle 0$ $\displaystyle=v_{P}^{\prime}\cdot M_{P}\cdot
c_{bio}-c_{P}\cdot v_{bio}^{\prime}$ $\displaystyle\bm{v}^{\prime}$
$\displaystyle\in\text{arg
}\underset{\bm{v}\in\mathbb{R}^{n}}{\text{max}}\quad v_{bio}$
$\displaystyle\begin{aligned} \text{s.t.}\quad\bm{S}\bm{v}&=\bm{0}\\\
v_{bio}&\geq f\cdot v_{bio,WT}\\\
\bm{v}&\geq\bm{v^{lower}}\circ(\bm{B}\bm{y})\\\
\bm{v}&\leq\bm{v^{upper}}\circ(\bm{B}\bm{y}),\end{aligned}$

where $\circ$ denotes the element-wise product, $\bm{y}\in\\{0,1\\}^{r}$
denotes the binary knockout vector, and $r$ denotes the number of reactions in
a metabolic network, which contains irreversible and reversible reactions,
further referred to as a reversible network. The parameter $K$ is the number
of maximal allowed knockouts, and the parameter $f$ denotes a value between
zero and one. The matrix $\bm{B}$ maps the reactions in the irreversible
metabolic network with those in the reversible network. The apostrophe ′
denotes the cellular variables in the upper level. Note that, in the
algorithm, they are fixed by the lower-level program.

As described above, we chose to eliminate $D$ and $c_{S}$ from the
formulation. In theory, the process variables $c_{bio}$ and $c_{P}$ could also
be eliminated, as the equality constraints at the upper level are fixing them.
Nevertheless, we decided on this formulation because preliminary tests and our
intuition suggest that it results in simpler convergence. Note that our recent
research indicates that an intermediate between full-space (keeping all
variables) and reduced-space formulation (eliminating all variables except the
degrees of freedom) is most beneficial in some cases [36, 37]. With $c_{bio}$
and $c_{P}$ being fixed by equality constraints, the substrate concentration
in the feed $c_{S,Feed}$ is the only degree of freedom in the process
conditions. At the cellular level, $v_{S}$ will reach its upper bound based on
an experimental value, as it is standard in other strain optimization examples
(e.g., [11]). Even if the Monod kinetics constrain $v_{bio}$ and $c_{S,Feed}$
is fixed during optimization; spare carbon atoms will go into product
formation.

In our study, the Monod constant was set to
$0.044\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$, and $v_{bio}^{max}$ was
set to $0.73\text{\,}{\mathrm{h}}^{-1}$, according to Wick et al. [38]. The
wild-type biomass flux $v_{bio,WT}$ was determined by a previously performed
FBA, and the parameter $f$ was set to 0.1. Furthermore, the ATP maintenance
reaction flux threshold was set to
$6.86\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$,
according to the default value in iML1515 [39]. The upper bound of the glucose
exchange reaction was set to
$10\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$.

Alternatively, we used Michaelis-Menten, the simplest form of enzyme kinetics.
Enzymes catalyze reactions within a cell, with metabolites as reactants. An
enzyme-catalyzed reaction can be described using Michaelis-Menten kinetics. In
alignment with Mahadevan et al. [25] and Ploch et al. [28], we chose to apply
Michaelis-Menten kinetics to the substrate uptake reaction. Thus, the kinetics
link the substrate concentration in the bioreactor with the substrate uptake
rate. More precisely, the substrate uptake rate is a function of the substrate
concentration: $v_{S}=g(c_{S})$. The kinetics read

$\displaystyle v_{S}$ $\displaystyle=v_{S}^{max}\frac{c_{S}}{c_{S}+K_{S,MM}}.$
$\displaystyle(Michaelis-Menten)$

The maximal substrate uptake is denoted with $v_{S}^{max}$ in
$\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$, and
the Michaelis constant is $K_{S,MM}$ in
$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$. Again, we assumed that the kinetic
parameters were constants and, thus, can be retrieved from the literature. As
$v_{S}$ is a cellular variable, the kinetics were implemented in the lower
level. SimulKnock with embedded Michaelis-Menten kinetics reads

$\displaystyle\underset{\bm{y},c_{S,Feed},c_{S},c_{P},c_{bio},v_{S}^{\prime},v_{bio}^{\prime},v_{P}^{\prime}}{\max}$
$\displaystyle\quad c_{P}\cdot v_{bio}^{\prime}$ (4) s.t.
$\displaystyle\sum_{i=1}^{r}(1-y_{i})$ $\displaystyle\leq K$ $\displaystyle 0$
$\displaystyle=-v_{S}^{\prime}\cdot M_{S}\cdot
c_{bio}+v_{bio}^{\prime}\cdot(c_{S,Feed}-c_{S})$ $\displaystyle 0$
$\displaystyle=v_{P}^{\prime}\cdot M_{P}\cdot c_{bio}-c_{P}\cdot
v_{bio}^{\prime}$ $\displaystyle\bm{v}^{\prime}$ $\displaystyle\in\text{arg
}\underset{\bm{v}\in\mathbb{R}^{n}}{\text{max}}\quad v_{bio}$
$\displaystyle\begin{aligned} \text{s.t.}\quad\bm{S}\bm{v}&=\bm{0}\\\
v_{bio}&\geq f\cdot v_{bio,WT}\\\
v_{S}&=v_{S}^{max}\frac{c_{S}}{c_{S}+K_{S,MM}}\\\
\bm{v}&\geq\bm{v^{lower}}\circ(\bm{B}\bm{y})\\\
\bm{v}&\leq\bm{v^{upper}}\circ(\bm{B}\bm{y}).\\\ \end{aligned}$

Again, $c_{bio}$ and $c_{P}$ are being fixed by the equality constraints in
the upper level. However, unlike the formulation including Monod, $c_{S}$ is a
degree of freedom at the process level now. It is included in the upper-level
equality constraint in the summand $v_{bio}\cdot(c_{S,Feed}-c_{S})$. With
$v_{bio}$ being part of the upper-level objective, this term also has an
optimum, and $c_{S}$ is determined during the optimization. Another difference
is that $v_{S}$ is now constrained by the kinetics instead of reaching the
upper bound set by the user.

In our study, the Michaelis constant was calculated with
$0.53\text{\,}\mathrm{mmol}\text{\,}{\mathrm{L}}^{-1}$ multiplied by the molar
weight of glucose, and $v_{S}^{max}$ was set to
$10\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$,
both values according to Meadows et al. [40].

To simplify the comparison of Monod and Michaelis-Menten, Table 1 summarizes
the different aspects of both kinetics applied in SimulKnock.

Table 1: Comparison of Monod and Michaelis-Menten in SimulKnock.
Abbreviations: $c_{S}$: substrate concentration; $v_{bio}$: growth rate;
$v_{S}$: substrate uptake rate; $f$, $g$: notation for functions.

| Monod | Michaelis-Menten
---|---|---
functionality | $c_{S}=f(v_{bio})$ | $v_{S}=g(c_{S})$
level | whole-cell kinetics $\rightarrow$ process level | enzyme kinetics $\rightarrow$ cellular level
$v_{S}$ | unrestricted, reaches its upper bound | restricted by kinetics and $c_{S}$
$c_{S}$ | calculated by kinetics | determined by equality constraint and objective in process level
assumption | gene deletion does not change growth behavior | gene deletion does not change substrate uptake behavior
parameter availability | higher | lower
possible extensions | inclusion of inhibitions | consideration of multiple substrates

Apart from the aspects already discussed above, we found that, in general, the
parameter availability for Monod is higher than for Michaelis-Menten.
Furthermore, the inclusion of inhibitions is supposedly a straightforward
extension with Monod, whereas including multiple substrates would potentially
be easier with Michaelis-Menten (cf. [28]).

###  2.4 Reformulation of the bilevel program to a single-level quadratically-
constrained quadratic program

We used the software packages libALE [41] and libDIPS [42] to implement
SimulKnock in both versions—Monod (3) and Michaelis-Menten (4)—and
successfully solved the bilevel programs for the E. coli core network [43]
using the algorithm [44] without KKT-based tightening and with Gurobi v10.0
[45] as the subsolver. However, this algorithm was too expensive for the
genome-scale metabolic model iML1515. [39]. Consequently, we reformulated the
bilevel programs to single-level mixed-integer quadratically-constrained
quadratic programs (MIQCQP). We conducted two reformulation steps for
numerical reasons: first, we reformulated SimulKnock as a single-level
program. Second, we eradicated nonlinear terms to achieve an MIQCQP
formulation of SimulKnock.

To the first reformulation: In both versions of SimulKnock (Monod or
Michaelis-Menten embedded), the cellular level is linear (in the lower-level
variables, as the process level fixes the process variables). Thus, we can
apply strong duality and reformulate to a single-level program. The
reformulations for SimulKnock with Monod kinetics embedded are identical to
those presented in OptKnock [11], as the Monod kinetics are in the upper-level
program. In the following, we will show the reformulation of SimulKnock with
Michaelis-Menten kinetics embedded. For ease of notation, we introduced

$\bm{\tilde{C}}=\begin{bmatrix}\bm{I}^{|n|}\\\ \bm{I}^{|n|}\\\
\bm{c}^{T}_{kin}\end{bmatrix},\;\bm{\tilde{b}}=\begin{bmatrix}\bm{v^{lower}}\circ(\bm{B}\bm{y})\\\
-\bm{v^{upper}}\circ(\bm{B}\bm{y})\\\
d\end{bmatrix},\;c_{i}=\begin{cases}0,\text{if}\;i\neq bio\\\
1,\text{if}\;i=bio,\end{cases}\forall i\in{1,...,n},$

where $\bm{c}$ describes the position of the biomass reaction and $\bm{I}$ is
the unity matrix. The vector $\bm{c}_{kin}\in\\{0,1\\}^{|n|}$ has exactly one
nonzero entry at the index of the glucose uptake reaction. The scalar $d$
describes the right-hand side of the kinetics, i.e.,
$d=v_{S}^{max}\frac{c_{S}}{c_{S}+K_{S,MM}}$. By applying the definitions
above, the lower-level program of (4) can be rewritten as

$\displaystyle\underset{\bm{v}}{\text{max}}$ $\displaystyle\bm{c}^{T}\bm{v}$
(5) s.t. $\displaystyle\bm{S}\cdot\bm{v}=\bm{0}$
$\displaystyle\bm{\tilde{C}}\cdot\bm{v}\leq\bm{\tilde{b}}$
$\displaystyle\bm{v}\geq 0.$

Converting program 5 to its dual and applying the strong duality theorem
yields the system of equations

$\displaystyle\bm{c}^{T}\bm{v}$
$\displaystyle=\bm{\tilde{b}}^{T}\cdot\bm{\tilde{\mu}}$ (6)
$\displaystyle\bm{S}\cdot\bm{v}$ $\displaystyle=\bm{0}$
$\displaystyle\bm{\tilde{C}}\cdot\bm{v}$ $\displaystyle\leq\bm{\tilde{b}}$
$\displaystyle\bm{c}$
$\displaystyle=\bm{S}^{T}\cdot\bm{\lambda}+\bm{\tilde{C}}^{T}\cdot\bm{\tilde{\mu}}$
$\displaystyle\bm{\tilde{\mu}}\geq\bm{0},$

with $\bm{\lambda}\in\mathbb{R}^{|m|}$ being the dual variables corresponding
to the mass balances for the metabolites and
$\bm{\tilde{\mu}}\in\mathbb{R}_{+}^{|2n+1|}$ being the dual variables
corresponding to the bounds set on the reaction fluxes through the inequality
constraints and the kinetics. The set of equations 6 can readily replace the
cellular level in (4). The resulting program is a single-level mixed-integer
nonlinear program. The nonlinearity is introduced by the fact that process and
cellular variables are at the same level now. The nonlinear term is the
kinetic term in $d$. These considerations about nonlinearity are also
applicable to the reformulated SimulKnock with Monod embedded.

In the second reformulation step, we reformulated the nonlinear kinetic term.
We introduced an additional variable for the fraction on the right-hand side
of the kinetics, thereby introducing an additional equality constraint in the
problem. For example, in the Michaelis-Menten kinetics, the fraction is
replaced with the optimization variable $\sigma_{MM}$, and the constraints

$\displaystyle v_{S}=v_{S}^{max}\sigma_{MM}$ (7)
$\displaystyle\sigma_{MM}(c_{S}+K_{S,MM})=c_{S}$

are introduced, which results in the single-level optimization problem
becoming an MIQCQP. This reformulation is beneficial because MIQCQPs are
solved directly by commercial solvers, such as Gurobi [45].

Replacing the lower level in (4) with the set of constraints of (6) and
applying the reformulation of the kinetic term suggested in (7) yields the
single-level mixed-integer quadratic form of SimulKnock with Michaelis-Menten
embedded:

$\displaystyle\underset{\bm{y},c_{S},c_{S,Feed},\sigma}{\max}$
$\displaystyle\quad c_{P}\cdot v_{bio}$ (8) s.t.:
$\displaystyle\sum(1-y_{i})\leq 1$ $\displaystyle 0=-v_{S}\cdot M_{S}\cdot
c_{bio}+v_{bio}\cdot(c_{S,Feed}-c_{S})$ $\displaystyle 0=v_{P}\cdot M_{P}\cdot
c_{bio}-c_{P}\cdot v_{bio}$ $\displaystyle v_{S}=v_{S}^{max}\sigma_{MM}$
$\displaystyle c_{S}=\sigma_{MM}(c_{S}+K_{S,MM})$ $\displaystyle
v_{bio}=\left[-\bm{v^{lower}}\circ\left(\bm{B}\bm{y}\right)\quad\bm{v^{upper}}\circ\left(\bm{B}\bm{y}\right)\quad
v_{S}^{max}\sigma\right]\cdot\bm{\mu}$ $\displaystyle\bm{S}\cdot\bm{v}=\bm{0}$
$\displaystyle\bm{v}\geq\bm{v^{lower}}\circ\left(\bm{B}\bm{y}\right)$
$\displaystyle\bm{v}\leq\bm{v^{upper}}\circ\left(\bm{B}\bm{y}\right)$
$\displaystyle\bm{S}^{T}\cdot\bm{\lambda}+\left[\bm{I}^{|n|}\quad\bm{I}^{|n|}\quad\bm{c}_{S}\right]\cdot\bm{\mu}=\bm{c}$
$\displaystyle v_{bio}\geq f\cdot v_{bio}^{max}$
$\displaystyle\bm{v},\bm{\mu}\geq\bm{0}$

Similar reformulations also reduce SimulKnock with Monod kinetics to an
MIQCQP. The reformulated programs were implemented using Pyomo [46, 47] and
solved with Gurobi v10.0 [45].

##  3 Case studies

We performed three case studies with SimulKnock. The first case study
illustrates SimulKnock’s mode of action compared to OptKnock and a sequential
optimization approach. The second case study employs the genome-scale
metabolic model iML1515 [39] of E. coli for one to three knockout predictions
and elaborates on the differences between Monod and Michaelis-Menten kinetics.
The third case study compares the results of SimulKnock with the published
results of an experimental study on E. coli continuous fermentation [48]. The
first case study was performed on up to 8 Intel Xeon E5-2640 CPU threads and
was solved within less than a second. Case studies two and three were solved
within five hours on the RWTH high-performance computing cluster using 48
threads and up to 4 GB of memory per thread.

###  3.1 SimulKnock predicts different knockouts than OptKnock with embedded
illustrative network

To show SimulKnock’s mode of action compared to sequential optimization using
OptKnock, we constructed an illustrative network (Figure 2).

Figure 2: Illustrative network. The grey circles A to G denote metabolites,
$v_{O}$ and $v_{S}$ denote the flux of oxygen and the substrate uptake
reaction, respectively; $v_{P}$ and $v_{bio}$ are the target chemical and
biomass flux, respectively. All stoichiometric coefficients are one, except
where indicated with numbers. For one maximum allowable knockout, the wrenches
indicate the knockout prediction of SimulKnock and OptKnock.

We performed four optimizations with the illustrative network embedded: $(i)$
an FBA with maximization of biomass as the objective (referred to as wild-
type), $(ii)$ an OptKnock prediction, $(iii)$ a sequential approach using
OptKnock (see below), and $(iv)$ SimulKnock. Table 2 shows the results of the
case study.

Table 2: Results of the case study with the illustrative network (Figure 2) embedded. The substrate is glucose, and the maximum substrate flux allowed is set to $v_{S}^{upper}=$10\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$$; the upper bounds on A-B and A-F are set to $7\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$ and $3\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$, respectively; the maximum oxygen flux allowed is $v_{O}^{upper}=$3\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$$. Michaelis-Menten kinetics were applied in SimulKnock and the sequential approach, with the Michaelis constant $K_{S,MM}=$0.53\text{\,}\mathrm{mmol}\text{\,}{\mathrm{L}}^{-1}$$ [40]. In the sequential approach, 1. OptKnock, and 2. A process optimization is performed with embedded knockouts from OptKnock. Abbreviations: $v_{S}$: substrate uptake flux, $v_{O}$: oxygen flux, $v_{bio}$: growth rate, $c_{bio}$: biomass concentration, $c_{P}$: product concentration, STY: space-time yield | Knockout | $v_{S}$ | $v_{O}$ | $v_{bio}$ | $v_{P}$ | $c_{bio}$ | $c_{P}$ | STY
---|---|---|---|---|---|---|---|---
| [-] | [$\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$] | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$] | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}\text{\,}{\mathrm{h}}^{-1}$]
Wild-type | - | 10 | 3 | 13 | 0 | - | - | -
OptKnock | B-C | 10 | 3 | 9.5 | 3.5 | - | - | -
Sequential | B-C | 7.8 | 3 | 8.4 | 2.4 | 8.7 | 2.5 | 21.0
SimulKnock | F-C | 3 | 3 | 3 | 3 | 9.8 | 9.8 | 29.3

The sequential approach denotes a two-step procedure. First, OptKnock is
performed. Second, the predicted knockouts are applied to the network, and a
continuous process optimization is performed. The continuous process
optimization looks like the formulation of SimulKnock but with fixed
knockouts, which were determined from OptKnock. Thus, OptKnock is included in
the sequential approach. Still, we added the results of OptKnock for
completeness to allow for a direct comparison of functionalities with
SimulKnock. Michaelis-Menten kinetics were applied in the sequential approach
and SimulKnock. The FBA and OptKnock predictions do not imply process
variables, e.g., the substrate or product concentration, which is also why the
kinetics were not applied in OptKnock or FBA.

The results indicate that SimulKnock and OptKnock predict different optimal
knockouts for one maximum allowable knockout. The reason for this difference
is that OptKnock optimizes for the target chemical production, whereas
SimulKnock optimizes for space-time yield. Indeed, the product flux $v_{P}$ is
higher for OptKnock than for SimulKnock. The applied kinetics in SimulKnock
and the sequential approach force the substrate uptake lower than the maximum
allowed value. Namely, in the sequential approach, the substrate uptake is set
to
$8\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$
by the process optimization, whereas OptKnock predicted it to be
$10\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$.
In SimulKnock, only the reactions in the upper part of the network, i.e., the
reactions going via the metabolites A, F, G, E, and C, are active and used for
biomass and target chemical production. In contrast, in OptKnock and the
sequential approach, the reactions going via B, D, and E are used for target
chemical production, while the route via F and C produces additional biomass.
The space-time yield, which is the product of the growth rate $v_{bio}$ and
the product concentration $c_{P}$, is higher for SimulKnock than for the
sequential approach. Note that both formulations predict the same two
knockouts, B-C and F-C, and the same space-time yield for two maximum
allowable knockouts. However, this does not generally indicate that with more
allowed knockouts, the results of OptKnock and SimulKnock would be similar.
Instead, the illustrative network is constructed so small that no other
options exist.

###  3.2 SimulKnock can achieve higher space-time yields than sequential
optimization with embedded iML1515

To further investigate the computational tractability of the SimulKnock
approach, we chose a genome-scale metabolic model, iML1515 [39], which
describes the metabolism of E. coli. The metabolic network includes 1516
genes, resulting in 1877 metabolites and 2712 metabolic reactions. The
objective was set to maximize the space-time yield of six target chemicals for
a glucose-limited medium.
Furthermore, to highlight the superiority of the simultaneous strain and
process optimization, SimulKnock was compared against the sequential approach
on the space-time yield of the target chemical.

####  3.2.1 One knockout

For one knockout prediction, SimulKnock and the sequential approach furnished
identical results. The knockout predictions and space-time yields are in the
Supporting Information (Figure A.1 and Table A.1).

####  3.2.2 Two knockouts: Monod vs. Michaelis Menten

In this case study, we allowed for two gene knockouts. The comparison results
of Michaelis-Menten and Monod kinetics are plotted in Figure 3.

(a) $Michaelis-Menten$

(b) $Monod$

Figure 3: Comparison of space-time yield using Michaelis-Menten and Monod
kinetics for two reaction eliminations on the genome-scale metabolic model
iML1515 [39]. Ethanol, succinate, and lactate were produced via anaerobic
pathways; acetate, formate, and fumarate were produced via aerobic pathways.

Note that for E. coli, the growth rate is higher with oxygen than without. The
lower-level optimization problem always activates the oxygen uptake to
maximize growth. Hence, in SimulKnock, the oxygen supply must be by the
modeler. We performed runs with and without an oxygen supply. Figure 3 shows
the results where SimulKnock achieved the higher absolute value.

SimulKnock predicts significantly higher space-time yields for four of the six
considered target chemicals with embedded Michaelis-Menten kinetics. The
difference in the space-time yield from the two approaches is due to different
reaction knockouts predicted (see the knockouts table A.2 in the Supporting
Information for reference). For example, SimulKnock targets the
6-phosphogluconolactonase and the acetate reversible transport for higher
ethanol space-time yield instead of the ATP synthase and the triose-phosphate
isomerase. The most significant difference becomes visible for formate.
Instead of the phosphoglycerate mutase, SimulKnock predicts the
phosphofructokinase will be knocked out. Even if SimulKnock and the sequential
approach do not predict the same knockouts for fumarate, they reach the same
result with Michaelis-Menten.

The kinetics also influence the knockout predictions. SimulKnock with embedded
Michaelis-Menten kinetics suggests blocking Acetaldehyde dehydrogenase and
D-alanine-D-alanine dipeptidase for succinate production. With Monod embedded
dytosine deaminase is knocked out instead of the dipeptidase,. Interestingly,
SimulKnock achieves a lower space-time yield for fumarate than the sequential
approach with embedded Monod kinetics.

The sequential formulation maximizes the target chemical concentration at the
process level, with the biomass growth rate being used only in the metabolism
level problem. Therefore, the targeted reaction knockouts decrease the growth
rate to achieve the maximum target chemical concentration. In contrast, the
SimulKnock formulation recognizes the trade-off between a higher growth rate
and a higher target chemical flux. It is, therefore, not surprising that the
sequential approach predicts a higher target chemical concentration and lower
biomass concentrations than the SimulKnock approach. The different modes of
action become even more visible when comparing the molar product yield on the
substrate. We computed the yield by dividing the target chemical flux $v_{P}$
by the substrate flux $v_{S}$ (cf. Table A.2 in the Supporting Information for
the data). In all cases, the molar product yield was lower with SimulKnock
than with the sequential approach. At maximum, in the case of acetate
production, the yield was $45\text{\,}\%$ lower, whereas the space-time yield
was 2.6 times higher.

####  3.2.3 Three knockouts

Upon increasing the possible reaction eliminations to three, the SimulKnock
formulation with Michaelis-Menten kinetics could not be solved within a
feasible time limit. Thus, only the Monod kinetics results are presented in
Figure 4.

Figure 4: Space-time yield using Monod kinetics for three reaction
eliminations on the genome-scale metabolic model iML1515 Monk et al. [39]. All
six target chemicals were produced via aerobic pathways.

Note that, again, we performed runs with and without an oxygen supply. Figure
4 shows the results where SimulKnock achieved the higher absolute value, all
achieved using aerobic pathways.

Acetate and formate space-time yield is increased by about 100% when using
SimulKnock, whereas ethanol, succinate, and fumarate increase marginally, and
lactate remains the same. Once again, the reason for this difference is the
different reaction knockouts predicted by the two approaches (cf. Table A.3 in
the Supporting Information). The knockouts for lactate and fumarate are
similar, often targeting neighboring reactions, thereby predicting similar
yields. When ethanol is considered, two of the three reaction eliminations
target the same reactions with SimulKnock when compared to the sequential
approach. When the target chemical is acetate, all three knockouts are
different; with formate, two of the three knockouts are different. This
observation hints towards the possibility of missing possible knockout
strategies with the sequential approach.

It is worth noting that the assumption made in this study about the biomass
reaction being the rate-limiting step has yet to be experimentally validated.
Therefore, reaction eliminations could, theoretically, target the rate-
limiting step on which the Monod kinetic parameters were fitted. In that case,
a new parameter fitting using the modified organism would be required. In
contrast, the Michaelis-Menten kinetics do not require parameter refitting due
to reaction eliminations, as the substrate uptake reaction is always active.

When we line up space-time yield and growth rate with an increasing number of
knockouts, we observe different behaviors: For ethanol, succinate, lactate,
and formate, space-time yield and growth rate with SimulKnock stay more or
less constant over one to three knockouts. Exemplary, Figure 5(a) shows the
course for succinate.

(a) $Succinate$

(b) $Acetate$

Figure 5: Change of space-time yield and growth rate over the number of
knockouts for iML1515 [39] with embedded Monod kinetics. Succinate was
produced anaerobically for 1 and 2 knockouts and aerobically for 3 knockouts.
Acetate was produced via an aerobic pathway. STY: space-time yield, Seq.
Appr.: Sequential approach.

In the sequential approach, space-time yield decreases with an increasing
number of knockouts, as depicted at two gene knockouts in Figure 5(a). This
decrease in space-time yield comes due to the mathematical structure of
OptKnock. In OptKnock, high product yields are achieved at the expense of
reduced growth (cf. Table A.3 in the Supporting Information for data). In
turn, space-time yield (equals product concentration times growth rate) is
linearly dependent on the growth rate. The figure depicts the aerobic pathway
for three gene knockouts, following our decision to show the results for the
highest SimulKnock space-time yield. Hence, the fermentation conditions change
compared to the anaerobic conditions that apply for one and two gene
knockouts. While the space-time yield of the sequential approach would decline
for anaerobic conditions (not shown in the figure), it reaches similar results
as SimulKnock for aerobic conditions.

Acetate and fumarate exhibit a different behavior, as depicted in an example
in Figure 5(b). With SimulKnock, the space-time yield rises significantly at
three knockouts, in alignment with the increase in growth. With the sequential
approach, the space-time yield is lowest for two knockouts and then increases
again, decreasing the growth rate. This interesting behavior stems from the
two-step procedure. Only in the second step are the process conditions
optimized, but they can only influence one of the two factors, namely the
product concentration $c_{P}$. Thus, there is no monotonous behavior with the
sequential approach. This finding displays the problems often occurring during
scale-up when sequentially optimizing the microbial strain and the process
conditions. On the other hand, acetate case underpins the potential of
SimulKnock.

###  3.3 SimulKnock can furnish meaningful knockout predictions but also
exhibits model-experiment mismatch

In the last case study, we took the experimental study of van Heerden and
Nicol [48] as a test case to elaborate on how far SimulKnock reproduces their
results regarding space-time yield, dilution rate, and knockouts. We aimed to
see whether SimulKnock furnishes experimentally meaningful results. In their
experimental study, van Heerden and Nicol [48] produced succinic acid from E.
coli KJ134 in a continuous fermentation with two different glucose feed
concentrations. We applied SimulKnock with embedded Monod kinetics to the E.
coli GEMs iML1515 [39] as well as iEC1349_Crooks [49] with three maximum
allowable knockouts. While van Heerden and Nicol [48] tested different
dilution rates and measured the space-time yield for each, SimulKnock
predicted the optimal space-time yield and the corresponding dilution rate.
Figure 6 shows the case study results, depicted as the space-time yield over
the dilution rate. Note that, in continuous fermentation, the dilution rate
equals the growth rate $\mu$ and the biomass flux $v_{bio}$.

Figure 6: Comparison of experimental laboratory data for succinic acid
production using E. coli KJ134 [48] with results from SimulKnock, based on the
metabolic networks iML1515 [39] and iEC1349_Crooks [49]. Two glucose feed
concentrations were studied: $20\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$
and $50\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$.

In the experimental study, 13 genetic modifications were applied to the
organism to reach a maximum space-time yield of around
$1\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}\text{\,}{\mathrm{h}}^{-1}$ on
$20\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$ glucose feed. SimulKnock
predicted a maximum space-time yield of about 4 and
$6\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}\text{\,}{\mathrm{h}}^{-1}$
with iML1515 and iEC1349_Crooks, respectively, on
$50\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$. The space-time yields were
reached with three knockouts. The reactions to be eliminated, predicted with
iML1515 embedded, were acetate kinase, ATP synthase, and fumarase. With
iEC1349_Crooks embedded, the predicted reaction eliminations were
phosphotransacetylase, ATP synthase, and succinyl-CoA synthetase. Amongst
them, acetate kinase and phosphotransacetylase were also eliminated in the
experimental strain. Hence, SimulKnock proves to predict experimentally
meaningful results. However, the figure reveals a substantial model-experiment
mismatch. Especially the results of iEC1349_Crooks are higher than any
experimental value. In iEC1349_Crooks, the high growth rate might result from
the biomass function difference. For illustration, for an FBA with a glucose
limitation of $10\text{\,}\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}$,
iML1515 has a growth rate of $0.87\text{\,}{\mathrm{h}}^{-1}$, whereas
iEC1349_Crooks has a growth rate of $1.02\text{\,}{\mathrm{h}}^{-1}$. The
mismatch in space-time yield that comes visible with both metabolic networks
might result from SimulKnock not considering inhibition and repression effects
and being overly optimistic by not considering byproduct formation unless the
byproduct formation hinders biomass formation. The mismatch could be reduced
by including experimental flux data of the organism in question, for example,
from Fischer et al. [50] for E. coli, or by using additional information from
a protein allocation model [51]. In both cases, the bounds of specific fluxes
would be fixed in a preprocessing step, which can readily be done in the
SimulKnock code. These refinements of the metabolic network and adaptations at
the cellular level are interesting in case of application but were not in the
scope of this work.

##  4 Conclusion

We presented SimulKnock, an optimization formulation combining continuous
fermentation optimization with strain optimization via gene knockouts. Such an
optimization formulation allows us to consider industrial fermentation
conditions already in the strain design, thereby overcoming current scale-up
difficulties. SimulKnock is a bilevel optimization formulation, which we
transformed into a single-level mixed-integer nonlinear optimization program
and solved globally for an illustrative example network as well as for the E.
coli GEM iML1515 [39]. We applied Monod or Michaelis-Menten kinetics to
connect the process level to the cellular level in SimulKnock. SimulKnock with
applied Monod kinetics showed similar results as Michaelis-Menten kinetics,
with Monod being less computationally intensive. SimulKnock predicted
different knockouts than OptKnock [11]. Also, the space-time yield was
significantly higher with SimulKnock compared to a sequential approach of
OptKnock plus process optimization for aerobic and anaerobic conditions.
Compared to experimental data [48], SimulKnock indicated that higher space-
time yields could be achieved with fewer knockouts.

The computations for SimulKnock with embedded GEM were performed on a high-
performance computing cluster. A decrease in the computation demand and
reduction of computation time could be achieved by applying further network
reduction methods [52] or by considering essential genes [53]. Future work
should consider that kinetics would have to be adjusted when knockouts are
applied to the organism because Michaelis-Menten, describing the substrate
uptake, and Monod, describing the growth behavior, are affected by genetic
modifications. Especially for altered growth behavior, the update of the
kinetics should be considered. Our case study with experimental data indicated
that SimulKnock furnishes overly optimistic results. The conversion to a
robust optimization formulation, as suggested in RobustKnock [12], could
tackle this observation. Moreover, experimental flux data (for E. coli, see
Pharkya and Maranas [14] and Fischer et al. [50]) could reduce the model-
experiment mismatch by fixing the bounds of specific fluxes in a preprocessing
step. To make a more accurate prediction, SimulKnock could be extended both at
the process level and at the cellular level, e.g., by more elaborate microbial
formulations [8, 7, 2], as depicted by OptKnock [11] and its extensions [54,
55]. These formulations include a reference flux distribution, which could be
furnished using a protein allocation model [51] or experimental in vivo flux
data [56]. The extension could also comprise other genetic modifications,
i.e., gene insertion [13] and regulation [14], as well as other process
operation modes, e.g., continuous mode with cell retention, batch, and fed-
batch.

## Author contributions

ALZ designed the SimulKnock formulation, the reformulation, the
implementation, and the case studies. MS performed the formulation setup, the
reformulation, the implementation, and the case studies under the supervision
of ALZ and AsMa. ALZ and MS visualized the data. ALZ, MS, and AsMa analyzed
and discussed the data. ALZ and AsMa wrote the manuscript draft. AlMi and LMB
had the initial idea of SimulKnock, discussed the data, and reviewed the
draft. AlMi secured funding. All authors read and approved the final
manuscript.

CRediT:
ALZ: Conceptualization, Methodology, Software, Validation, Investigation,
Writing - Original Draft (Lead), Writing - Review and Editing, Visualization
AsMa: Methodology, Software, Investigation, Writing - Review and Editing
MS: Methodology, Software, Validation, Investigation, Visualization, Writing -
Review and Editing
LMB: Conceptualization, Writing - Review and Editing
AlMi: Conceptualization, Supervision, Methodology, Writing- Review and
Editing, Funding acquisition

## Acknowledgments

This project was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany´s Excellence Strategy – Cluster of
Excellence 2186 “The Fuel Science Center” – ID: 390919832. This project
benefited from the work of Clemens Kortmann, who improved the modularity and
the user-friendliness of our git repository. Computations were performed with
computing resources granted by RWTH Aachen University under project thes1376.

## Data Availability and Reproducibility Statement

The implementation of the SimulKnock approach, the data pre-processing, and
the interface with the solver Gurobi are openly available in our GitLab
repository “SimulKnock: Simultaneous design of fermentation and microbe” at
https://git.rwth-aachen.de/avt-svt/public/simulknock. Data from the
manuscript’s figures are tabulated in the Supporting Information.

## Declaration of conflicts

No conflicts to declare.

## Graphical Abstract

Figure 7: Microbial network and process optimizations are combined in our
optimization formulation SimulKnock.

## References

  * Clomburg et al. [2017] Clomburg JM, Crumbley AM, Gonzalez R.  Industrial biomanufacturing: The future of chemical production.  Science (New York, NY) 2017;355(6320).
  * Landon et al. [2019] Landon S, Rees-Garbutt J, Marucci L, Grierson C.  Genome-driven cell engineering review: in vivo and in silico metabolic and genome engineering.  Essays in biochemistry 2019;63(2):267–284.
  * Maia et al. [2016] Maia P, Rocha M, Rocha I.  In Silico Constraint-Based Strain Optimization Methods: the Quest for Optimal Cell Factories.  Microbiology and molecular biology reviews : MMBR 2016;80(1):45–67.
  * Varma and Palsson [1994] Varma A, Palsson BO.  Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110.  Applied and environmental microbiology 1994;60(10):3724–3731.
  * Orth et al. [2010] Orth JD, Thiele I, Palsson BØ.  What is flux balance analysis?  Nature biotechnology 2010;28(3):245–248.
  * Schuetz et al. [2012] Schuetz R, Zamboni N, Zampieri M, Heinemann M, Sauer U.  Multidimensional optimality of microbial metabolism.  Science (New York, NY) 2012;336(6081):601–604.
  * Segrè et al. [2002] Segrè D, Vitkup D, Church GM.  Analysis of optimality in natural and perturbed metabolic networks.  Proceedings of the National Academy of Sciences of the United States of America 2002;99(23):15112–15117.
  * Brochado et al. [2012] Brochado AR, Andrejev S, Maranas CD, Patil KR.  Impact of stoichiometry representation on simulation of genotype-phenotype relationships in metabolic networks.  PLoS computational biology 2012;8(11):e1002758.
  * Shlomi et al. [2005] Shlomi T, Berkman O, Ruppin E.  Regulatory on/off minimization of metabolic flux changes after genetic perturbations.  Proceedings of the National Academy of Sciences of the United States of America 2005;102(21):7695–7700.
  * Valderrama-Gomez et al. [2017] Valderrama-Gomez MA, Kreitmayer D, Wolf S, Marin-Sanguino A, Kremling A.  Application of theoretical methods to increase succinate production in engineered strains.  Bioprocess and biosystems engineering 2017;40(4):479–497.
  * Burgard et al. [2003] Burgard AP, Pharkya P, Maranas CD.  Optknock: a bilevel programming framework for identifying gene knockout strategies for microbial strain optimization.  Biotechnology and bioengineering 2003;84(6):647–657.
  * Tepper and Shlomi [2010] Tepper N, Shlomi T.  Predicting metabolic engineering knockout strategies for chemical production: accounting for competing pathways.  Bioinformatics (Oxford, England) 2010;26(4):536–543.
  * Pharkya et al. [2004] Pharkya P, Burgard AP, Maranas CD.  OptStrain: a computational framework for redesign of microbial production systems.  Genome research 2004;14(11):2367–2376.
  * Pharkya and Maranas [2006] Pharkya P, Maranas CD.  An optimization framework for identifying reaction activation/inhibition or elimination candidates for overproduction in microbial systems.  Metabolic engineering 2006;8(1):1–13.
  * Kim et al. [2011] Kim J, Reed JL, Maravelias CT.  Large-scale bi-level strain design approaches and mixed-integer programming solution techniques.  PloS one 2011;6(9):e24162.
  * Biegler [1997] Biegler LT.  Systematic methods of chemical process design.  Prentice Hall international series in the physical and chemical engineering sciences, Upper Saddle River, N.J.: Prentice Hall PTR; 1997.
  * Chmiel et al. [2018] Chmiel H, Takors R, Weuster-Botz D, editors.  Bioprozesstechnik.  4\. auflage ed. Berlin and Heidelberg: Springer Spektrum; 2018.
  * Villadsen et al. [2011] Villadsen J, Nielsen J, Lidén G.  Bioreaction Engineering Principles.  3\. ed. ed. Boston, MA: Springer Science+Business Media LLC; 2011.
  * Edgar et al. [2001] Edgar TF, Himmelblau DM, Lasdon LS.  Optimization of chemical processes.  2nd ed. ed. McGraw-Hill chemical engineering series, New York, NY: McGraw-Hill; 2001.
  * Gordeeva and Gordeev [2015] Gordeeva YL, Gordeev LS.  Optimization of continuous microbiological synthesis processes with nonlinear microbial growth kinetics.  Theoretical Foundations of Chemical Engineering 2015;49(6):829–835.
  * Wehrs et al. [2019] Wehrs M, Tanjore D, Eng T, Lievense J, Pray TR, Mukhopadhyay A.  Engineering Robust Production Microbes for Large-Scale Cultivation.  Trends in microbiology 2019;27(6):524–537.
  * Olsson et al. [2022] Olsson L, Rugbjerg P, Torello Pianale L, Trivellin C.  Robustness: linking strain design to viable bioprocesses.  Trends in biotechnology 2022;40(8):918–931.
  * Richelle et al. [2020] Richelle A, David B, Demaegd D, Dewerchin M, Kinet R, Morreale A, et al.  Towards a widespread adoption of metabolic modeling tools in biopharmaceutical industry: a process systems biology engineering perspective.  NPJ systems biology and applications 2020;6(1):6.
  * Lapin et al. [2004] Lapin A, Müller D, Reuss M.  Dynamic Behavior of Microbial Populations in Stirred Bioreactors Simulated with Euler$-$Lagrange Methods: Traveling along the Lifelines of Single Cells.  Industrial & Engineering Chemistry Research 2004;43(16):4647–4656.
  * Mahadevan et al. [2002] Mahadevan R, Edwards JS, Doyle FJ.  Dynamic Flux Balance Analysis of Diauxic Growth in Escherichia coli.  Biophysical Journal 2002;83(3):1331–1340.
  * de Oliveira et al. [2021] de Oliveira RD, Guedes MN, Matias J, Le Roux GAC.  Nonlinear Predictive Control of a Bioreactor by Surrogate Model Approximation of Flux Balance Analysis.  Industrial & Engineering Chemistry Research 2021;60(40):14464–14475.
  * Zhuang et al. [2013] Zhuang K, Yang L, Cluett WR, Mahadevan R.  Dynamic strain scanning optimization: an efficient strain design strategy for balanced yield, titer, and productivity. DySScO strategy for strain design.  BMC biotechnology 2013;13:8.
  * Ploch et al. [2019] Ploch T, Zhao X, Hüser J, von Lieres E, Hannemann-Tamás R, Naumann U, et al.  Multiscale dynamic modeling and simulation of a biorefinery.  Biotechnology and bioengineering 2019;116(10):2561–2574.
  * Jabarivelisdeh and Waldherr [2018] Jabarivelisdeh B, Waldherr S.  Optimization of bioprocess productivity based on metabolic-genetic network models with bilevel dynamic programming.  Biotechnology and bioengineering 2018;115(7):1829–1841.
  * Tafur Rangel et al. [2022] Tafur Rangel AE, Oviedo AG, Mojica FC, Gómez JM, Gónzalez Barrios AF.  Development of an integrating systems metabolic engineering and bioprocess modeling approach for rational strain improvement.  Biochemical Engineering Journal 2022;178:108268.
  * Konstantinos and Antonis [2022] Konstantinos D, Antonis K.  Simultaneous Synthesis of Metabolic and Process Engineering for the Production of Muconic Acid.  Proceedings of the 32nd European Symposium on Computer Aided Process Engineering (ESCAPE32) 2022;51:889–894.
  * Mitsos et al. [2009] Mitsos A, Chachuat B, Barton PI.  Towards global bilevel dynamic optimization.  Journal of Global Optimization 2009;45(1):63–93.
  * Stephanopoulos et al. [1998] Stephanopoulos G, Aristidou AA, Nielsen J.  Metabolic engineering: Principles and methodologies.  [nachdr.] ed. San Diego, Calif.: Acad. Press; 1998.
  * Gordeeva et al. [2021] Gordeeva EL, Ravichev LV, Borodkin AG, Gordeeva YL.  Mathematical Modeling of a Biotechnological Continuous Fermentation Process for Lactic Acid Production: A Review.  Theoretical Foundations of Chemical Engineering 2021;55(6):1192–1203.
  * Sinner et al. [2019] Sinner P, Kager J, Daume S, Herwig C.  Model-based Analysis and Optimisation of a Continuous Corynebacterium glutamicum Bioprocess Utilizing Lignocellulosic Waste.  IFAC-PapersOnLine 2019;52(26):181–186.
  * Bongartz and Mitsos [2019] Bongartz D, Mitsos A.  Deterministic global flowsheet optimization: Between equation–oriented and sequential–modular methods.  AIChE Journal 2019;65(3):1022–1034.
  * Najman et al. [2021] Najman J, Bongartz D, Mitsos A.  Linearization of McCormick relaxations and hybridization with the auxiliary variable method.  Journal of Global Optimization 2021;80(4):731–756.
  * Wick et al. [2001] Wick LM, Quadroni M, Egli T.  Short- and long-term changes in proteome composition and kinetic properties in a culture of Escherichia coli during transition from glucose-excess to glucose-limited growth conditions in continuous culture and vice versa.  Environmental microbiology 2001;3(9):588–599.
  * Monk et al. [2017] Monk JM, Lloyd CJ, Brunk E, Mih N, Sastry A, King Z, et al.  iML1515, a knowledgebase that computes Escherichia coli traits.  Nature biotechnology 2017;35(10):904–908.
  * Meadows et al. [2010] Meadows AL, Karnik R, Lam H, Forestell S, Snedecor B.  Application of dynamic flux balance analysis to an industrial Escherichia coli fermentation.  Metabolic engineering 2010;12(2):150–160.
  * Djelassi and Mitsos [2019] Djelassi H, Mitsos A, libALE – a library for algebraic-logical expression trees; 2019.  https://git.rwth-aachen.de/avt.svt/public/libale.git.
  * Jungen et al. [2023] Jungen D, Zingler A, Djelassi H, Mitsos A, libDIPS - a library for Discretization-Based Semi-Infinite Programming Solvers; 2023.  https://git.rwth-aachen.de/avt-svt/public/libdips/.
  * Orth et al. [2010] Orth JD, Fleming RMT, Palsson BØ.  Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide.  EcoSal Plus 2010;4(1).
  * Mitsos et al. [2008] Mitsos A, Lemonidis P, Barton PI.  Global solution of bilevel programs with a nonconvex inner program.  Journal of Global Optimization 2008;42(4):475–513.
  * Gurobi Optimization, LLC [2023] Gurobi Optimization, LLC, Gurobi Optimizer Reference Manual; 2023.  https://www.gurobi.com.
  * Bynum et al. [2021] Bynum ML, Hackebeil GA, Hart WE, Laird CD, Nicholson BL, Siirola JD, et al.  Pyomo–optimization modeling in python, vol. 67.  Third ed. Springer Science & Business Media; 2021.
  * Hart et al. [2011] Hart WE, Watson JP, Woodruff DL.  Pyomo: modeling and solving mathematical programs in Python.  Mathematical Programming Computation 2011;3(3):219–260.
  * van Heerden and Nicol [2013] van Heerden CD, Nicol W.  Continuous and batch cultures of Escherichia coli KJ134 for succinic acid fermentation: metabolic flux distributions and production characteristics.  Microbial cell factories 2013;12:80.
  * Monk et al. [2016] Monk JM, Koza A, Campodonico MA, Machado D, Seoane JM, Palsson BO, et al.  Multi-omics Quantification of Species Variation of Escherichia coli Links Molecular Features with Strain Phenotypes.  Cell systems 2016;3(3):238–251.e12.
  * Fischer et al. [2004] Fischer E, Zamboni N, Sauer U.  High-throughput metabolic flux analysis based on gas chromatography-mass spectrometry derived 13C constraints.  Analytical biochemistry 2004;325(2):308–316.
  * Alter et al. [2021] Alter TB, Blank LM, Ebert BE.  Proteome Regulation Patterns Determine Escherichia coli Wild-Type and Mutant Phenotypes.  mSystems 2021;6(2).
  * Erdrich et al. [2015] Erdrich P, Steuer R, Klamt S.  An algorithm for the reduction of genome-scale metabolic network models to meaningful core models.  BMC systems biology 2015;9:48.
  * Goodall et al. [2018] Goodall ECA, Robinson A, Johnston IG, Jabbari S, Turner KA, Cunningham AF, et al.  The Essential Genome of Escherichia coli K-12.  mBio 2018;9(1).
  * Apaydin et al. [2017] Apaydin M, Xu L, Zeng B, Qian X.  Robust mutant strain design by pessimistic optimization.  BMC genomics 2017;18(Suppl 6):677.
  * Kim and Reed [2010] Kim J, Reed JL.  OptORF: Optimal metabolic and regulatory perturbations for metabolic engineering of microbial strains.  BMC systems biology 2010;4:53.
  * Kuepfer et al. [2005] Kuepfer L, Sauer U, Blank LM.  Metabolic functions of duplicate genes in Saccharomyces cerevisiae.  Genome research 2005;15(10):1421–1430.

##  Appendix A Supporting Information

###  A.1 Analysis of the formulation with both kinetics embedded

Including both kinetics, Michaelis-Menten and Monod, may be tempting to tackle
the drawback of Monod, which is that the substrate uptake flux $v_{S}$ is
unrestricted. Regarding the parameter availability in the literature, the
inclusion of two kinetics is impractical, as it typically is difficult to
attain the parameters for both kinetics. Moreover, including both kinetics may
result in an inconsistency, which we analyze in the following.

Looking at the mathematical functionalities, the inclusion of both kinetics
leads to

$\displaystyle\text{Monod:\,}c_{S}=f(v_{bio})\text{\,and Michaelis-
Menten:\,}v_{S}=g(c_{S})$ $\displaystyle\xRightarrow{}v_{S}=g(f(v_{bio}))$
$\displaystyle\text{metabolic\,network:\,}v_{bio}\in\mathcal{H}(v_{S}),$

where the set $\mathcal{H}$ indicates that $v_{bio}$ depends on $v_{S}$
through their connections in the metabolic network and that there is no unique
solution. These relations could lead to inconsistency. Alternatively, the
kinetics could constrain $v_{S}$ and $v_{bio}$ to such an extent that the
metabolic network would only simulate the product and side-product fluxes
instead of optimizing $v_{bio}$. This case is not desirable either because we
assumed that the metabolic network is more accurate than the kinetics, which
is why we developed SimulKnock in the first place.

Based on these considerations, we tested a formulation where we inserted Monod
into Michaelis-Menten kinetics such that, as described above,
$v_{S}=g(f(v_{bio}))$. We conducted exemplary runs using iML1515. The
formulation was infeasible for anaerobic production, exemplarily shown for
succinate (1 and 2 knockouts) and acetate (1 knockout). The formulation is
feasible only when ignoring the maintenance threshold, with the resulting
space-time yield being zero. For aerobic production, the formulation did not
converge within the time limit of $20.000\text{\,}\mathrm{s}$ for formate
production. For succinate (2 knockouts) and acetate production (1 knockout),
the achieved space-time yields were far lower
($0.16\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}\text{\,}{\mathrm{h}}^{-1}$
and
$0.013\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}\text{\,}{\mathrm{h}}^{-1}$,
respectively) than the results with only one kinetics embedded. At the same
time, we observed elevated growth rates. For succinate production (2
knockouts), for example, the growth rate was $0.68\text{\,}{\mathrm{h}}^{-1}$,
which, set in relation to the experimental results of van Heerden and Nicol
[48], is not realistic.

From our mathematical analysis and the exemplary study, which aligned with the
mathematical analysis, we conclude that it is not beneficial to include both
kinetics.

###  A.2 Data for one reaction elimination

(a) $Michaelis-Menten$

(b) $Monod$

Figure A.1: Comparison of space-time yield using Michaelis-Menten and Monod
kinetics for one reaction elimination on the genome-scale metabolic model
iML1515 [39]. Ethanol, succinate, and lactate were produced via anaerobic
pathways; acetate, formate, and fumarate were produced via aerobic pathways.

Table A.1: One reaction elimination: Comparison of space-time yield, product
concentration, substrate concentration, biomass concentration, growth rate,
substrate uptake flux $v_{S}$, product flux $v_{P}$, and molar yield using
Michaelis-Menten and Monod kinetics for the genome-scale metabolic model
iML1515 [39]. Note that, for all results, the substrate feed concentration
reached its upper bound of $10\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$.
STY: Space-time yield, MM: Michaelis-Menten

Chemical | Kinetic | Optimization | Knockouts | STY | Product | Substrate | Biomass | Growth | vs | vp | Yield
---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}\text{\,}{\mathrm{h}}^{-1}$] | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$] | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$] | [$\mathrm{g_{CDW}}\text{\,}{\mathrm{L}}^{-1}$] | [${\mathrm{h}}^{-1}$] | [$\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$] | [$\mathrm{mol}\text{\,}{\mathrm{mol}}^{-1}$]
Ethanol | MM | Sequential | Acetate reversible transport via proton symport (periplasm) | 0.43 | 4.11 | 1.09 | 0.56 | 0.1 | 9.19 | 16.57 | 1.8
|  | SimulKnock | Acetate reversible transport via proton symport (periplasm) | 0.43 | 4.12 | 1.07 | 0.56 | 0.1 | 9.18 | 16.57 | 1.8
| Monod | Sequential | Acetate reversible transport via proton symport (periplasm) | 0.55 | 4.58 | 0.01 | 0.66 | 0.12 | 10.0 | 17.93 | 1.79
|  | SimulKnock | Acetate reversible transport via proton symport (periplasm) | 0.55 | 4.58 | 0.01 | 0.66 | 0.12 | 10.0 | 17.93 | 1.79
Succinate | MM | Sequential | Ethanol transport via diffusion (extracellular to periplasm) | 0.54 | 5.23 | 1.1 | 0.56 | 0.1 | 9.2 | 8.39 | 0.91
|  | SimulKnock | Ethanol transport via diffusion (extracellular to periplasm) | 0.54 | 5.25 | 1.1 | 0.56 | 0.1 | 9.2 | 8.42 | 0.91
| Monod | Sequential | Ethanol transport via diffusion (extracellular to periplasm) | 0.69 | 5.84 | 0.01 | 0.66 | 0.12 | 10.0 | 9.07 | 0.91
|  | SimulKnock | Ethanol transport via diffusion (extracellular to periplasm) | 0.69 | 5.84 | 0.01 | 0.66 | 0.12 | 10.0 | 9.07 | 0.91
Lactate | MM | Sequential | Ethanol transport via diffusion (extracellular to periplasm) | 0.83 | 8.01 | 1.1 | 0.56 | 0.1 | 9.2 | 16.74 | 1.82
|  | SimulKnock | Ethanol transport via diffusion (extracellular to periplasm) | 0.83 | 8.04 | 1.09 | 0.56 | 0.1 | 9.19 | 16.78 | 1.83
| Monod | Sequential | Ethanol transport via diffusion (extracellular to periplasm) | 1.06 | 8.95 | 0.01 | 0.66 | 0.12 | 10.0 | 18.12 | 1.81
|  | SimulKnock | Ethanol transport via diffusion (extracellular to periplasm) | 1.06 | 8.95 | 0.01 | 0.66 | 0.12 | 10.0 | 18.12 | 1.81
Acetate | MM | Sequential | ATP synthase (four protons for one ATP) (periplasm) | 1.19 | 4.68 | 0.96 | 1.4 | 0.25 | 9.09 | 14.37 | 1.58
|  | SimulKnock | ATP synthase (four protons for one ATP) (periplasm) | 1.19 | 4.69 | 0.96 | 1.4 | 0.25 | 9.09 | 14.38 | 1.58
| Monod | Sequential | ATP synthase (four protons for one ATP) (periplasm) | 1.46 | 5.14 | 0.02 | 1.58 | 0.28 | 10.0 | 15.71 | 1.57
|  | SimulKnock | ATP synthase (four protons for one ATP) (periplasm) | 1.46 | 5.13 | 0.03 | 1.58 | 0.28 | 10.0 | 15.71 | 1.57
Formate | MM | Sequential | ATP synthase (four protons for one ATP) (periplasm) | 3.9 | 15.34 | 0.96 | 1.4 | 0.25 | 9.1 | 61.8 | 6.79
|  | SimulKnock | ATP synthase (four protons for one ATP) (periplasm) | 3.9 | 15.35 | 0.96 | 1.4 | 0.25 | 9.09 | 61.79 | 6.79
| Monod | Sequential | ATP synthase (four protons for one ATP) (periplasm) | 4.8 | 16.85 | 0.02 | 1.58 | 0.28 | 10.0 | 67.56 | 6.76
|  | SimulKnock | ATP synthase (four protonsfor one ATP) (periplasm) | 4.8 | 16.83 | 0.03 | 1.58 | 0.28 | 10.0 | 67.56 | 6.76
Fumarate | MM | Sequential | Fumarase | 0.35 | 0.48 | 0.91 | 4.04 | 0.72 | 9.05 | 0.76 | 0.08
|  | SimulKnock | Fumarase | 0.35 | 0.47 | 1.15 | 3.93 | 0.74 | 9.24 | 0.78 | 0.08
| Monod | Sequential | Tryptophanase (L-tryptophan) | 0.0 | 0.0 | 0.01 | 0.0 | 0.16 | 10.0 | 0.0 | 0.0
|  | SimulKnock | 6-phosphogluconolactonase | 0.0 | 0.0 | 0.01 | 0.0 | 0.16 | 10.0 | 0.0 | 0.0

###  A.3 Data for two reaction eliminations

Table A.2: Two reaction eliminations: Comparison of space-time yield, product
concentration, substrate concentration, biomass concentration, growth rate,
substrate uptake flux $v_{S}$, product flux $v_{P}$, and molar yield using
Michaelis-Menten and Monod kinetics for the genome-scale metabolic model
iML1515 [39]. Note that, for all results, the substrate feed concentration
reached its upper bound of $10\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$.
STY: Space-time yield, MM: Michaelis-Menten

Chemical | Kinetic | Optimization | Knockouts | STY | Product | Substrate | Biomass | Growth | vs | vp | Yield
---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}\text{\,}{\mathrm{h}}^{-1}$] | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$] | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$] | [$\mathrm{g_{CDW}}\text{\,}{\mathrm{L}}^{-1}$] | [${\mathrm{h}}^{-1}$] | [$\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$] | [$\mathrm{mol}\text{\,}{\mathrm{mol}}^{-1}$]
Ethanol | MM | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 0.09 | 4.2 | 1.64 | 0.11 | 0.02 | 9.45 | 18.58 | 1.97
(anaerob) |  |  | 2\. Triose-phosphate isomerase |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. 6-phosphogluconolactonase | 0.43 | 4.16 | 1.1 | 0.56 | 0.1 | 9.2 | 16.81 | 1.83
|  |  | 2\. Acetate reversible transport via proton symport (periplasm) |  |  |  |  |  |  |  |
| Monod | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 0.13 | 5.02 | 0.0 | 0.14 | 0.03 | 10.0 | 19.62 | 1.96
|  |  | 2\. Triose-phosphate isomerase |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. 6-phosphogluconolactonase | 0.55 | 4.65 | 0.01 | 0.65 | 0.12 | 10.0 | 18.19 | 1.82
|  |  | 2\. Acetate reversible transport via proton symport (periplasm) |  |  |  |  |  |  |  |
Succinate | MM | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 0.42 | 7.0 | 1.23 | 0.32 | 0.06 | 9.28 | 11.51 | 1.24
(anaerob) |  |  | 2\. D-glucose transport via PEP:Pyr PTS (periplasm) |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. Acetaldehyde dehydrogenase (acetylating) | 0.54 | 5.25 | 1.09 | 0.56 | 0.1 | 9.2 | 8.41 | 0.91
|  |  | 2\. D-alanine-D-alanine dipeptidase |  |  |  |  |  |  |  |
| Monod | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 0.55 | 7.98 | 0.0 | 0.39 | 0.07 | 10.0 | 12.39 | 1.24
|  |  | 2\. D-glucose transport via PEP:Pyr PTS (periplasm) |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. Acetaldehyde dehydrogenase (acetylating) | 0.69 | 5.84 | 0.01 | 0.66 | 0.12 | 10.0 | 9.07 | 0.91
|  |  | 2\. Cytosine deaminase |  |  |  |  |  |  |  |
Lactate | MM | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 0.17 | 8.12 | 1.65 | 0.11 | 0.02 | 9.45 | 18.59 | 1.97
(anaerob) |  |  | 2\. Triose-phosphate isomerase |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. 6-phosphogluconolactonase | 0.83 | 8.04 | 1.1 | 0.56 | 0.1 | 9.2 | 16.8 | 1.83
|  |  | 2\. Acetate reversible transport via proton symport (periplasm) |  |  |  |  |  |  |  |
| Monod | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 0.25 | 9.7 | 0.0 | 0.14 | 0.03 | 10.0 | 19.62 | 1.96
|  |  | 2\. Triose-phosphate isomerase |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. 6-phosphogluconolactonase | 1.06 | 8.98 | 0.01 | 0.65 | 0.12 | 10.0 | 18.19 | 1.82
|  |  | 2\. Acetate reversible transport via proton symport (periplasm) |  |  |  |  |  |  |  |
Acetate | MM | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 0.55 | 5.34 | 1.1 | 0.55 | 0.1 | 9.2 | 16.85 | 1.83
(aerob) |  |  | 2\. Phosphoglycerate mutase |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. Phosphoserine transaminase | 1.2 | 4.85 | 0.96 | 1.36 | 0.25 | 9.1 | 14.89 | 1.64
|  |  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
| Monod | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 0.7 | 5.97 | 0.01 | 0.65 | 0.12 | 10.0 | 18.23 | 1.82
|  |  | 2\. Phosphoglycerate mutase |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. ATP synthase (four protons for one ATP)(periplasm) | 1.48 | 5.32 | 0.03 | 1.53 | 0.28 | 10.0 | 16.29 | 1.63
|  |  | 2\. Phosphoglycerate dehydrogenase |  |  |  |  |  |  |  |
Formate | MM | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 1.85 | 17.96 | 1.1 | 0.55 | 0.1 | 9.2 | 74.23 | 8.07
(aerob) |  |  | 2\. Phosphoglycerate mutase |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. Phosphofructokinase (s7p) | 3.9 | 15.35 | 0.96 | 1.4 | 0.25 | 9.09 | 61.79 | 6.79
|  |  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
| Monod | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 2.35 | 20.06 | 0.01 | 0.65 | 0.12 | 10.0 | 80.35 | 8.04
|  |  | 2\. Phosphoglycerate mutase |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. Phosphogluconate dehydrogenase | 4.8 | 16.83 | 0.03 | 1.58 | 0.28 | 10.0 | 67.56 | 6.76
|  |  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
Fumarate | MM | Sequential | 1\. Fumarase | 0.35 | 0.48 | 0.9 | 4.04 | 0.72 | 9.04 | 0.75 | 0.08
(aerob) |  |  | 2\. Cytochrome oxidase bd (menaquinol-8: 2 protons) (periplasm) |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. Fumarase | 0.35 | 0.49 | 0.9 | 4.04 | 0.72 | 9.04 | 0.77 | 0.08
|  |  | 2\. 1-tetradecanoyl-sn-glycerol 3-phosphate O-acyltransferase (n-C14:0) |  |  |  |  |  |  |  |
| Monod | Sequential | 1\. Fumarase | 0.41 | 0.51 | 0.41 | 4.28 | 0.8 | 10.0 | 0.84 | 0.08
|  |  | 2\. Cytochrome oxidase bd (menaquinol-8: 2 protons) (periplasm) |  |  |  |  |  |  |  |
|  | SimulKnock | 1\. Fumarase | 0.29 | 0.42 | 0.78 | 3.54 | 0.69 | 10.0 | 0.72 | 0.07
|  |  | 2\. Cytochrome oxidase bo3 (ubiquinol-8: 4 protons) (periplasm) |  |  |  |  |  |  |  |

###  A.4 Data for three reaction eliminations

Table A.3: Three reaction eliminations: Comparison of space-time yield,
product concentration, substrate concentration, biomass concentration, growth
rate, substrate uptake flux $v_{S}$, product flux $v_{P}$, and molar yield
using Monod kinetics for the genome-scale metabolic model iML1515 [39] for
aerobic production. Note that, for all results, the substrate feed
concentration reached its upper bound of
$10\text{\,}\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$. STY: Space-time yield

Chemical | Optimization | Knockouts | STY | Product | Substrate | Biomass | Growth | vs | vp | Yield
---|---|---|---|---|---|---|---|---|---|---
|  |  | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}\text{\,}{\mathrm{h}}^{-1}$] | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$] | [$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$] | [$\mathrm{g_{CDW}}\text{\,}{\mathrm{L}}^{-1}$] | [${\mathrm{h}}^{-1}$] | [$\mathrm{mmol}\text{\,}{\mathrm{g_{CDW}}}^{-1}\text{\,}{\mathrm{h}}^{-1}$] | [$\mathrm{mol}\text{\,}{\mathrm{mol}}^{-1}$]
Ethanol | Sequential | 1\. ATP synthase (four protons for one ATP) (periplasm) | 0.6 | 4.6 | 0.01 | 0.72 | 0.13 | 10.0 | 17.99 | 1.8
|  | 2\. Acetate reversible transport via proton symport (periplasm) |  |  |  |  |  |  |  |
|  | 3\. Succinate dehydrogenase (irreversible) |  |  |  |  |  |  |  |
| SimulKnock | 1\. Phosphotransacetylase | 0.6 | 4.45 | 0.01 | 0.75 | 0.14 | 10.0 | 17.41 | 1.74
|  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
|  | 3\. Succinyl-CoA synthetase (ADP-forming) |  |  |  |  |  |  |  |
Succinate | Sequential | 1\. Fumarase | 0.75 | 5.86 | 0.01 | 0.71 | 0.13 | 10.0 | 9.1 | 0.91
|  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
|  | 3\. Acetate reversible transport via proton symport (periplasm) |  |  |  |  |  |  |  |
| SimulKnock | 1\. Phosphotransacetylase | 0.76 | 5.6 | 0.01 | 0.75 | 0.14 | 10.0 | 8.7 | 0.87
|  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
|  | 3\. Succinyl-CoA synthetase (ADP-forming) |  |  |  |  |  |  |  |
Lactate | Sequential | 1\. 2-Oxogluterate dehydrogenase | 1.16 | 8.69 | 0.01 | 0.74 | 0.13 | 10.0 | 17.58 | 1.76
|  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
|  | 3\. Acetate reversible transport via proton symport (periplasm) |  |  |  |  |  |  |  |
| SimulKnock | 1\. Phosphotransacetylase | 1.16 | 8.6 | 0.01 | 0.75 | 0.14 | 10.0 | 17.41 | 1.74
|  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
|  | 3\. Succinyl-CoA synthetase (ADP-forming) |  |  |  |  |  |  |  |
Acetate | Sequential | 1\. 6-phosphogluconate dehydratase | 0.87 | 7.6 | 0.01 | 0.64 | 0.11 | 10.0 | 23.21 | 2.32
|  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
|  | 3\. Triose-phosphate isomerase |  |  |  |  |  |  |  |
| SimulKnock | 1\. Phosphoserine transaminase | 2.24 | 4.1 | 0.13 | 3.0 | 0.55 | 10.0 | 12.68 | 1.27
|  | 2\. Phosphogluconate dehydrogenase |  |  |  |  |  |  |  |
|  | 3\. Succinate dehydrogenase (irreversible) |  |  |  |  |  |  |  |
Formate | Sequential | 1\. 2-dehydro-3-deoxy-phosphogluconate aldolase | 2.57 | 22.32 | 0.01 | 0.64 | 0.11 | 10.0 | 89.39 | 8.94
|  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
|  | 3\. Triose-phosphate isomerase |  |  |  |  |  |  |  |
| SimulKnock | 1\. ATP synthase (four protons for one ATP) (periplasm) | 4.8 | 16.83 | 0.03 | 1.58 | 0.28 | 10.0 | 67.56 | 6.76
|  | 2\. Glucosyltransferase II (LPS core synthesis) |  |  |  |  |  |  |  |
|  | 3\. Glycerol-3-phosphate acyltransferase (C16:0) |  |  |  |  |  |  |  |
Fumarate | Sequential | 1\. Fumarase | 0.74 | 5.75 | 0.01 | 0.71 | 0.13 | 10.0 | 9.1 | 0.91
|  | 2\. ATP synthase (four protons for one ATP) (periplasm) |  |  |  |  |  |  |  |
|  | 3\. Acetate reversible transport via proton symport (periplasm) |  |  |  |  |  |  |  |
| SimulKnock | 1\. Fumarase | 2.27 | 4.12 | 0.14 | 3.02 | 0.55 | 10.0 | 6.6 | 0.66
|  | 2\. Phosphoserine transaminase |  |  |  |  |  |  |  |
|  | 3\. Phosphogluconate dehydrogenase |  |  |  |  |  |  |  |

###  A.5 Data for Figure 6: experimental data compared with SimulKnock results

Table A.4: Comparison of SimulKnock results with lab data from van Heerden and Nicol [48]. Glucose feed concentration | Dilution rate | space-time yield | source
---|---|---|---
$\mathrm{g}\text{\,}{\mathrm{L}}^{-1}$ | ${\mathrm{h}}^{-1}$ | $\mathrm{g}\text{\,}{\mathrm{L}}^{-1}\text{\,}{\mathrm{h}}^{-1}$ | -
20 | 0.06 | 0.72 |
| 0.15 | 0.73 |
| 0.15 | 0.68 |
| 0.09 | 0.85 | van Heerden and Nicol [48]
| 0.09 | 0.87 |
| 0.19 | 0.46 |
| 0.04 | 0.52 |
| 0.13 | 1.51 | SimulKnock, iML1515
| 0.22 | 2.42 | SimulKnock, iEC1349_Crooks
50 | 0.06 | 0.71 |
| 0.06 | 0.69 |
| 0.09 | 0.65 | van Heerden and Nicol [48]
| 0.09 | 0.59 |
| 0.04 | 0.58 |
| 0.04 | 0.57 |
| 0.02 | 0.38 |
| 0.13 | 3.79 | SimulKnock, iML1515
| 0.22 | 6.05 | SimulKnock, iEC1349_Crooks