Modeling 2: composition and balances

Writing every balance of a bioprocess by hand is where errors hide: a yield that does not conserve electrons, a dilution term forgotten in one equation. In this tutorial you describe an aerobic batch culture by its species and its reaction instead, Labinetix writes the balances, derives the yields the chemistry determines, and refuses a reaction that would create atoms or electrons.

It continues Your first model; the commands are the same.

The idea

Every species follows the same balance (Bastin & Dochain’s general dynamical model), per liquid volume:

dc_i/dt = Σ_feeds D_k (c_in,ik − c_i) + Σ_reactions K_ij r_j + T_i

with D_k the dilution rate of feed k, K_ij the stoichiometric coefficient of species i in reaction j, r_j the reaction’s rate and T_i the gas–liquid transfer. Everything but the rates is generic. What is specific to your process is which reactions exist and how fast they run; the rest follows from the species’ formulas.

1. The model

labinetix define new model_scaffold aerobic_batch

Replace models/aerobic_batch/authoring/aerobic_batch.toml with:

[model]
name = "aerobic_batch"
contract_version = "0.1.0"
kind = "ode"
time_unit = "h"
doc = "Aerobic batch growth on glucose, composed."

[parameters.mu_max]
unit = "1/h"
default = 0.4
doc = "Maximum specific growth rate."

[parameters.k_s]
unit = "g/L"
default = 0.1
doc = "Monod constant of glucose."

[parameters.y_xs]
unit = "g/g"
default = 0.5
doc = "Biomass yield on glucose."

[parameters.y_os]
unit = "g/g"
default = 0.4
doc = "Oxygen taken per glucose."

[parameters.v_r]
unit = "L"
default = 1.0
doc = "Liquid volume."

[parameters.kla]
unit = "1/h"
default = 300.0
doc = "Oxygen transfer coefficient."

[parameters.c_o2_sat]
unit = "g/L"
default = 0.007
doc = "Oxygen saturation concentration."

[states.c_x]
unit = "g/L"
min = 0.0
doc = "Biomass."

[states.c_s]
unit = "g/L"
min = 0.0
doc = "Glucose."

[states.c_o2]
unit = "g/L"
min = 0.0
doc = "Dissolved oxygen."

[liquid]
volume = "v_r"

[species.biomass]
formula = "CH1.8O0.5N0.2"
state = "c_x"

[species.glucose]
formula = "C6H12O6"
state = "c_s"

[species.oxygen]
formula = "O2"
state = "c_o2"

[species.co2]
formula = "CO2"

[species.water]
formula = "H2O"

[species.ammonia]
formula = "NH3"

[[reactions]]
name = "growth"
catalyst = "biomass"
rate = "mu_max / y_xs * monod(c_s_pos, k_s)"
stoichiometry = { glucose = -1, biomass = "y_xs", oxygen = "-y_os" }
balance = ["co2", "water", "ammonia"]

[[transfers]]
species = "oxygen"
kla = "kla"
saturation = "c_o2_sat"

[outputs.c_x_out]
unit = "g/L"
expr = "c_x_pos"
doc = "Biomass."

What is new:

The rate here is per gram of biomass (catalyst = "biomass") in grams of glucose: glucose has coefficient −1, so the reaction consumes mu_max / y_xs · s/(k_s + s) g glucose per g biomass per hour and makes y_xs g biomass from each. A given coefficient is a mass ratio (g_X / g_S), and the unit check holds the generated balances to the states’ units like any equation.

You write no [equations]: the composition writes the balance of every species state. The generated locals c_x_pos, c_s_pos, … (a state read as zero when negative), r_growth (the reaction’s rate per litre) and tr_oxygen are available to your own locals and outputs.

2. A reaction that does not balance

Sync:

labinetix define sync model_scaffold/aerobic_batch

labinetix: define sync: invalid entity state: …/aerobic_batch.toml: [reactions] `growth`: does not balance for every parameter value: per g of `biomass`, element residuals [C +2.2740e-2, H +5.6850e-3, N -1.7055e-2, O -1.1370e-2] mol/g; degree of reduction +1.7055e-1 mol e⁻/g. The elements fix 4 of its 6 coefficients and `balance` derives 3: derive 1 more, one of `biomass`, `oxygen` (a coefficient the elements determine is never a parameter)

Six species take part in the reaction and four elements constrain them, so two coefficients are free, and the other four follow. You gave three (glucose, the biomass yield and the oxygen yield) and derived three. With all three fixed, the elements C, H, O and N cannot all close for every value of the yields, and the electrons (degree of reduction) do not either. A fit could move y_os away from the one value that conserves them, and the model would quietly create electrons: the error that a soft sensor reading the off-gas turns into a wrong biomass.

The fix is to let the chemistry decide: keep the biomass yield, the one quantity the culture chooses, and derive oxygen. Delete the y_os parameter and change the reaction:

stoichiometry = { glucose = -1, biomass = "y_xs" }
balance = ["oxygen", "co2", "water", "ammonia"]

labinetix define sync model_scaffold/aerobic_batch

Labinetix solves the element balances for the four listed species. Every derived coefficient is a generated local nu_<reaction>_<species>, linear in the yields it depends on, which your rates may read too:

#![allow(unused)]
fn main() {
let l_nu_growth_oxygen = -1.0656764137747288 + 1.3643097014167802 * p[2];
}

That is −1.066 g O2 per g glucose for full combustion, less 1.364 g for every g of biomass built: the electrons that end in biomass need no oxygen.

3. Read the structure

labinetix define balance model_scaffold/aerobic_batch --measured oxygen

model_scaffold/aerobic_batch — composition at the parameter defaults

species (atoms per formula unit, γ the degree of reduction):
              M g/mol      C      H      N      O       γ  state
  biomass     24.6263  1.000  1.800  0.200  0.500   4.200  c_x
  glucose    180.1560  6.000 12.000  0.000  6.000  24.000  c_s
  oxygen      31.9980  0.000  0.000  0.000  2.000  -4.000  c_o2
  co2         44.0090  1.000  0.000  0.000  2.000   0.000  -
  water       18.0150  0.000  2.000  0.000  1.000   0.000  -
  ammonia     17.0310  0.000  3.000  1.000  0.000   0.000  -

stoichiometry (g per unit extent; * derived from the elements), electron residual (mol e⁻ per unit extent), degrees of freedom and the free yields (given coefficients that read parameters):
               biomass     glucose      oxygen         co2       water     ammonia          e⁻  dof  free yields
  growth       0.5000     -1.0000     -0.3835*     0.5722*     0.3805*    -0.0692*    -1.2e-16    2  biomass (y_xs)

reaction invariants of the states (2 of 3 states; concentrations, g/L):
  +2.0000·c_x +1.0000·c_s
  +0.7670·c_x +1.0000·c_o2

rates from the measured conversion rates of oxygen:
  rank 1 of 1 reactions, redundancy 0: every reaction rate is determined
  determined directions: 1; balance tests on the measured rates: 0
  calculable conversion rates: biomass, glucose, co2, water, ammonia
  not calculable: none

What it tells you:

• The stoichiometry, with the derived coefficients marked: growth takes 0.38 g O2 and makes 0.57 g CO2 per g glucose. The electron residual is zero to rounding.
• Degrees of freedom of each reaction: what you may choose (here two: the reference glucose and the biomass yield). The free yields are the given coefficients that read parameters: what an estimate fits; every derived one follows from them.
• Reaction invariants: combinations of states the reaction cannot change. 2 c_x + c_s stays constant in a batch: every gram of glucose consumed makes half a gram of biomass. These are what a kinetics-free observer can track without knowing any rate.
• Which measurements determine the rates. With the oxygen uptake measured, the one reaction’s rate is determined; a model with several reactions tells you which further species you must measure (try models/overflow_yeast, where oxygen and CO2 are not enough). Measuring more than the rates need leaves balance tests (the redundancy), and every species whose rate the measured ones fix is calculable.
• An estimator’s instruments. --estimator <name> reads the measured outputs of models/<id>/estimators/<name>/ (From a model to a released package) instead: each output that reads one species, an off-gas fraction, a concentration, the dissolved oxygen, counts as that species’ measured rate.

4. Gas phase, feeds and more reactions

The larger models in models/ show the rest:

• models/aerobic_fed_batch: a growing volume fed with glucose ([liquid] volume = "v", [[feeds]]), growth and maintenance as two reactions, and a quasi-steady [gas] whose outputs are the off-gas O2 and CO2 fractions. Every species a reaction makes or takes without a state is listed either as a gas species (inlet) or as untracked (water, ammonia).

• models/overflow_yeast: the Sonnleitner–Käppeli bottleneck, three reactions on one catalyst; the fermentation’s oxygen is given as 0, so its ethanol and CO2 yields are derived. Its rates read the derived oxygen coefficients to turn a respiratory capacity into rates of glucose and ethanol.

• models/e4b_family_r: a model in modules. Its [model] names the organism organisms/e4b.toml (E. coli with acetate overflow and re-uptake, four reactions) and the vessel vessels/pilot_fed_batch.toml (feed, titrant, oxygen transfer, [gas], and the [headspace] and [storage] declarations). The sections merge into one model with one hash.

Simulate them as in the first tutorial; their scenarios are in models/<id>/scenarios/. The first two have an independent cross-check in crates/studio/analyze-model/tests/: the model written a second time by hand, with closed-form yields, integrated by an adaptive solver. e4b_family_r’s right-hand side is held to v1’s soft-sensor benchmark at sampled states (crates/studio/analyze-fixtures/goldens/softsensor/family_r_rhs.json).

Next

Estimation and identifiability asks which parameters a measurement can determine, draws synthetic data, and estimates the parameters with their uncertainty. From a model to a released package takes a composed model from the store to a package taktwerk runs, and back from its record.