Modeling

Our iGEM project focuses on generation of microbial fuel cells (MFCs) powered by yeast. The aim of the project is to engineer yeast strains which would have enhanced electron transfer capacity to the electrodes, which would power the MFCs. The bioengineering part of the strategy of the project focused on generating yeast strains (i) devoid of proteins which consume excess electrons produced by yeast metabolism; and (ii) supplementing the yeast with an additional factor which would enhance their electron transfer to the electrodes of the fuel cell. Within the scope of the project, we identified two key avenues where modelling could be used to account for experimental data as well as to indicate places for further improvements in the future. On the one hand, we generated the Data-Driven Redox Budget model to try to identify the genetic modifications, which would allow for the most optimal path to maintain the balance between higher electron output and cell viability. Our second model focuses on the step of electron transfer from the cell population to the electrode: the effects of cell population growth biofilm formation on the electrode, and possible power output.

Cellular

Population/System

Cellular level: DDRB (Data-Driven Redox Budget) model

Microbial fuel cells convert cellular metabolism into electrical current, but rational strain design has been difficult without predictive models linking genotype to electron availability. The key bottleneck is NADH/NAD⁺ redox balance. In fermenting yeast, every glucose generates two NADH molecules through glycolysis, creating an electron surplus that must be consumed to prevent their toxic accumulation and to balance the metabolic burden. While cells possess numerous NADH-consuming pathways including amino acid biosynthesis, lipid synthesis, and various oxidoreductases, we focus on the three dominant fluxes under fermentative conditions: ethanol production (the primary fermentation product consuming >70% of glycolytic NADH), glycerol biosynthesis (the major osmoprotectant and redox valve accounting for 5-10% of carbon flux), and biosynthetic reactions coupled to growth. Minor pathways like acetate or succinate production contribute <5% of total NADH consumption in glucose fermentation and are subsumed into the biosynthetic term γ_bio·μ or residual q_MFC. The steady-state NADH pool can be described as

, where q terms are specific production/consumption rates, γ_bio is NADH demand per unit biomass, μ is growth rate, and q_MFC represents residual electron flux potentially available for electrodes [1].

Genetic perturbations create metabolic imbalances

Deleting major NADH-consuming pathways forces metabolic rewiring. GPD1/GPD2 deletions eliminate glycerol production, removing a key electron sink. NDE1/NDE2 deletions impair respiratory NADH oxidation [4, 7]. These perturbations create electron "traffic jams" that cells must resolve through alternative routes or extracellular electron transfer. Our Data-Driven Redox Budget approach quantifies these imbalances through mechanistic modeling constrained by thermodynamic conservation laws.

Genotype-Dependent Adaptation

A critical innovation is recognizing that NADH demand (γ_bio) adapts to genetic background. Strains with impaired metabolism may reduce anabolic NADH consumption as a compensatory mechanism [2]. We model this as

, where indicator variables represent gene deletions and coefficients capture systematic metabolic adaptations. By fitting this model to measured rates from five strains, we predict electron surplus for unlimited novel genotypes, transforming strain engineering from empirical screening into rational, computation-guided design.

Overview of the Modeling Pipeline

The Data-Driven Redox Budget model consists of five sequential steps that transform raw fermentation data into predictive strain rankings. The pipeline flows as follows: (1) experimental data collection and quality control, (2) specific rate calculation with bootstrap confidence intervals, (3) γ_bio model fitting and NADH balance validation, (4) genetic perturbation models for all phenotypes, and (5) double mutant prediction and MFC ranking. Each step builds on the previous, propagating both point estimates and uncertainty through the entire framework.

Data collection and experimental design

Measured variables

Time-series measurements over t = 0 to 48 hours at regular intervals(typically every 4-6 hours)

X(t): Biomass concentration [g dry weight/L]

S_glc(t): Glucose concentration [mM]

P_eth(t): Ethanol concentration [mM]

P_gly(t): Glycerol concentration [mM]

Experimental design

Five yeast strains with defined genetic backgrounds:

Dom90: Wild-type reference strain (I_GPD1=0, I_GPD2=0, I_NDE=0)
NDE1: NDE1 deletion (I_GPD1=0, I_GPD2=0, I_NDE=1)
NDE2: NDE2 deletion (I_GPD1=0, I_GPD2=0, I_NDE=1)
GPD1: GPD1 deletion (I_GPD1=1, I_GPD2=0, I_NDE=0)
GPD2: GPD2 deletion (I_GPD1=0, I_GPD2=1, I_NDE=0)

Each strain measured with n_rep = 2 biological replicates.

Genetic indicator variables

Binary indicators encode genetic perturbations:

I_GPD1 ∈ {0,1}: Presence of GPD1 deletion
NI_GPD1 ∈ {0,1}: Presence of GPD1 deletion
I_GPD1 ∈ {0,1}: Presence of GPD1 deletion
I_GPD1 ∈ {0,1}: Presence of GPD1 deletion
I_GPD1 ∈ {0,1}: Presence of GPD1 deletion

Specific rate calculation

Rates are calculated over the exponential growth and primary fermentation phase:

t_start = 0 h
t_end = 24 h
Δt = t_end - t_start = 24 h

Rationale: This window captures stable exponential growth with active fermentation before glucose depletion or stationary phase effects.

Specific growth rate μ [h⁻¹] from exponential fit:

where:

X_initial = X(t_start): Initial biomass
X_final = X(t_end): Final biomass

For specific rate calculations, we need the time-integrated biomass:

Calculated using trapezoidal rule:

Specific glucose consumption rate q_glc [mmol/(g·h)]:

where:

ΔS_glc = S_glc(t_start) - S_glc(t_end): Total glucose consumed [mM]
Must be positive (consumption)

Specific ethanol production rate q_eth [mmol/(g·h)]:

where:

ΔP_eth = P_eth(t_end) - P_eth(t_start): Net ethanol produced [mM]

Specific glycerol production rate q_gly [mmol/(g·h)]:

where:

ΔP_gly = P_gly(t_end) - P_gly(t_start): Net glycerol produced [mM]

Ethanol yield should not exceed theoretical maximum:

Y_eth/glc = q_eth / q_glc ≤ 2.0 mol/mol

Rationale: Glucose → 2 Pyruvate → 2 Ethanol (Embden-Meyerhof-Parnas pathway stoichiometry).

To quantify uncertainty, we resample biological replicates with replacement (n_bootstrap = 1000).

For each bootstrap iteration b:

Randomly select n_rep replicates with replacement
Calculate X̄_b(t), S̄_glc,b(t), P̄_eth,b(t), P̄_gly,b(t) as means
Compute μ_b, q_glc,b, q_eth,b, q_gly,b using formulas above
Store bootstrap sample

95% confidence intervals from percentiles:

CI_lower = 2.5th percentile of bootstrap distribution
CI_upper = 97.5th percentile of bootstrap distribution

Standard error approximation:

Results

Bootstrap resampling (n=1000) of biological replicates yielded 95% confidence intervals for all specific rates. Growth rates ranged from 0.14-0.18 h⁻¹, consistent with S. cerevisiae fermentation literature. Relative uncertainties averaged ±8%, indicating good precision for 2-replicate experiments (Figure 1.1). Ethanol yields (calculated as q_eth/q_glc) remained below the stoichiometric maximum of 2.0 mol ethanol/mol glucose, validating mass balance. GPD1 deletion reduced glycerol production by 94%, while GPD2 deletion reduced it by only 20%, confirming GPD1 as the dominant glycerol-3-phosphate dehydrogenase isoform.

Figure 1.1 Specific rates with bootstrap 95% CI.

γ_bio model and NADH balance

NADH Stoichiometry

The fundamental redox balance equation for glucose fermentation:

NADH Production:

NADH Consumption:

where:

q_eth: 1 NADH consumed per ethanol (pyruvate → acetaldehyde → ethanol via alcohol dehydrogenase)
q_gly: 1 NADH consumed per glycerol (DHAP → glycerol-3-phosphate via GPD1/GPD2)
γ_bio · μ: NADH for biosynthesis [mmol/(g·h)]

γ_bio: NADH requirement [mmol NADH/g biomass]
μ: Growth rate [h⁻¹]

q_MFC: Residual NADH flux available for microbial fuel cell [mmol/(g·h)]

For each strain i, calculate γ_bio,i from measured rates assuming q_MFC = 0 (initial estimate):

This gives five individual γ_bio values, one per strain.

Linear additive model relating γ_bio to genetic background:

where:

α: Baseline NADH demand (wild-type) [mmol/g]
β_GPD: Effect of GPD deletion on γ_bio [mmol/g]
β_NDE: Effect of NDE deletion on γ_bio [mmol/g]
I_GPD = max(I_GPD1, I_GPD2): GPD deletion indicator
I_NDE: NDE deletion indicator
ε ~ N(0, σ²): Random error term

Ordinary Least Squares (OLS) estimation

X = [1 I_GPD₁ I_NDE₁] [1 I_GPD₂ I_NDE₂] [⋮ ⋮ ⋮ ] [1 I_GPD₅ I_NDE₅]

Response vector y:

y = [γ_bio,1] [γ_bio,2] [ ⋮ ] [γ_bio,5]

Parameter estimates by OLS:

where β̂ = [α̂, β̂_GPD, β̂_NDE]'

Model predictions and residuals

Predicted γ_bio for strain i:

Residuals:

Model Quality Metrics

Sum of squared residuals:

Total sum of squares:

Coefficient of determination:

Mean squared error:

Root mean squared error:

Using predicted γ̂_bio,i, calculate residual electron flux:

NADH Closure Validation

NADH closure percentage for strain i:

where:

Mean closure across all strains:

Standard deviation of closure:

Acceptance criteria:

Excellent: |Closure_mean - 100%| < 2% and σ_closure < 5%
Good: |Closure_mean - 100%| < 5% and σ_closure < 10%

F-test for overall model significance:

SS_model = SS_tot - SS_res df_model = p - 1 = 2 df_resid = n - p = 2

p-value from F-distribution with (df_model, df_resid) degrees of freedom.

Parameter standard errors

Variance-covariance matrix:

Standard errors:

t-statistics for each parameter:

p-values from t-distribution with df_resid degrees of freedom.

Results:

Individual γ_bio calculations from NADH balance ranged 0.80-3.90 mmol/g (mean 2.01 ± 1.07). The genotype-dependent model fitted with R² = 0.91, significantly outperforming the universal (constant) model by information criteria (ΔAIC = +6.4, evidence ratio 24:1 in favor). Parameter estimates: α = 3.91 mmol/g (baseline), β_GPD = -2.80 mmol/g, β_NDE = -1.93 mmol/g. Although individual coefficients showed marginal significance due to small sample size (n=5), the combined model demonstrated excellent NADH closure: 100.2 ± 2.4%, validating the mechanistic constraint.

Paired bars in Fugure 1.2 show NADH produced by glycolysis (dark green) versus consumed NADH partitioned into ethanol, glycerol, and biosynthesis γ_bio·μ. All strains achieve 97.5-102.9% closure, confirming the model captures major electron sinks. GPD deletions visibly reduce glycerol synthesis while maintaining near-perfect balance. Closure within 3% validates model accuracy and indicates positive q_MFC represents genuine electron surplus available for extracellular transfer.

Figure 1.2 NADH conservation analysis.

Bootstrap validation with 1000 resampled datasets confirmed parameter stability (Figure 1.3). The baseline biosynthetic demand α shows tight distribution centered at 3.91 mmol/g with narrow confidence intervals, indicating high precision. The negative coefficients β_GPD = -2.80 and β_NDE = -1.93 represent metabolic adaptation: deletion mutants reduce their NADH consumption as a compensatory response to impaired redox pathways. GPD deletions eliminate the glycerol sink, while NDE deletions block respiratory NADH oxidation, forcing cells to downregulate anabolic processes that would otherwise consume excess NADH. Bootstrap distributions overlap with parametric 95% CI (dashed lines), validating both estimation approaches despite limited sample size.

In Figure 1.3 bootstrap parameter distributions validate model stability. Violin plots show 1000 bootstrap iterations. Box plots indicate median (line), interquartile range (box), and point estimates (dots). Dashed horizontal lines show parametric 95% CI. Negative β values reflect reduced NADH demand in deletion mutants as a compensatory mechanism. Overlapping distributions confirm parameter robustness despite n=5 limitation.

Figure 1.3 Bootstrap parameter distributions.

Genetic perturbation models for all phenotypes

Standard phenotype models

For μ, q_glc, q_eth (three phenotypes), fit linear additive models:

where:

Φ ∈ {μ, q_glc, q_eth}: Phenotype
α_Φ: Baseline phenotype value (wild-type)
δ_GPD,Φ: Effect of GPD deletion on phenotype
δ_NDE,Φ: Effect of NDE deletion on phenotype
ε_Φ ~ N(0, σ²_Φ): Phenotype-specific error

For q_gly, use isoform-specific model to distinguish GPD1 vs GPD2 [4]:

Estimation procedure

For each phenotype Φ:

Construct design matrix X (5×3 for standard, 5×4 for glycerol)
Response vector y_Φ = measured phenotype values
OLS estimation: θ̂_Φ = (X'X)⁻¹X'y_Φ
Calculate predictions: ŷ_Φ = X·θ̂_Φ
Compute R², RMSE, F-statistic, parameter p-value

For phenotype Φ:

R²_Φ = 1 - SS_res,Φ/SS_tot,Φ Adjusted R²:

where p_Φ is number of parameters for that phenotype.

Overall Model Quality
Mean R² across all phenotypes:

Parameter budget
Total parameters estimated:

μ model: 3 parameters (α, δ_GPD, δ_NDE)
q_glc model: 3 parameters
q_eth model: 3 parameters
q_gly model: 4 parameters (α, δ_GPD1, δ_GPD2, δ_NDE)
γ_bio model: 3 parameters (from Step 3)

Total: 13 + 3 = 16 parameters from 5 strains × 4 phenotypes = 20 data points

Results

Linear additive models fitted for all phenotypes achieved mean R² = 0.79 across μ, q_glc, q_eth, q_gly, explaining 79% of variance with 13 parameters from 20 data points. Growth rate model (R² = 0.79) showed GPD deletions reduce μ by 15% (δ_GPD = -0.027 h⁻¹, p = 0.035), while NDE deletions have minimal effect (δ_NDE = -0.003 h⁻¹, p = 0.072). Laboratory observations independently confirmed the reduced growth phenotype of GPD mutants, validating model predictions. Glycerol model (R² = 0.87) confirmed isoform specificity: GPD1 deletions eliminate 94% of production (δ_GPD1 = -0.077 mmol/(g·h), p < 0.001), GPD2 deletions reduce by 20% (δ_GPD2 = -0.016 mmol/(g·h), p = 0.018). Ethanol model (R² = 0.72) revealed NDE deletions increase ethanol by 11% (δ_NDE = +0.214 mmol/(g·h), p = 0.041), indicating compensatory redox balancing. Residual diagnostics confirmed model assumptions: all residuals within ±2σ, no systematic patterns in heteroscedasticity tests, mean relative errors 3-10% of measured values (Figure 1.4).

Figure 1.4 DDRB model validation - measured vs predicted strain-specific rates.

Double mutant prediction and MFC ranking

Novel genotypes

Define three double mutant strains not in training data:

GPD1-NDE1: I_GPD1=1, I_GPD2=0, I_NDE=1
GPD2-NDE1: I_GPD1=0, I_GPD2=1, I_NDE=1
GPD1-GPD2: I_GPD1=1, I_GPD2=1, I_NDE=0

For each novel genotype g with indicators (I_GPD1,g, I_GPD2,g, I_NDE,g):

where I_GPD,g = max(I_GPD1,g, I_GPD2,g).

Physical constraints

For ensuring biological feasibility:

μ̂_g = max(μ̂_g, 0.001) [h⁻¹] q̂_glc,g = max(q̂_glc,g, 0.1) [mmol/(g·h)] q̂_gly,g = max(q̂_gly,g, 0.0) [mmol/(g·h)] γ̂_bio,g = max(γ̂_bio,g, 0.8) [mmol/g]

For novel genotype g:

For each measured strain i, we have n_boot bootstrap samples of (μ_b, q_glc,b, q_eth,b, q_gly,b) from Step 2.

For each bootstrap iteration b:

Calculate γ_bio,i,b from bootstrap rates
Refit γ_bio model → β̂_b = [α̂_b, β̂_GPD,b, β̂_NDE,b]
For novel genotype g, predict:

Use Step 4 models with bootstrap rates
Calculate q_MFC,g,b

Store bootstrap sample

Bootstrap distribution of q_MFC for each strain enables:

95% CI: [Q_2.5%, Q_97.5%] from percentiles
P(q_MFC > 0): Fraction of bootstrap samples with positive q_MFC

For n_strains total strains (measured + predicted):

For each bootstrap iteration b:

Calculate q_MFC,i,b for all strains i
Rank strains by q_MFC,i,b (highest = rank 1)
Store rank vector R_b = [rank_1,b, rank_2,b, ..., rank_n,b]

For each strain i:

Mean rank: R̄_i = (1/n_boot) Σ rank_i,b
Rank std dev: σ_rank,i = √[(1/(n_boot-1)) Σ(rank_i,b - R̄_i)²]
Top-3 probability: P(rank ≤ 3) = (# times rank_i,b ≤ 3) / n_boot

Strains ranked by point estimate q̂_MFC in descending order, with tiers:

Tier 1 (High priority): q̂_MFC > 0.05 mmol/(g·h)
Tier 2 (Moderate): 0.0 < q̂_MFC ≤ 0.05 mmol/(g·h)
Tier 3 (Low priority): q̂_MFC ≤ 0.0 mmol/(g·h)

Results

The validated genetic perturbation models enabled prediction of three novel double mutants: GPD1-NDE1, GPD2-NDE1, and GPD1-GPD2. Predictions assumed linear additivity of genetic effects, supported by high model performance. Physical constraints ensured biological feasibility by enforcing minimum viable growth, basal metabolism, non-negative rates, and minimum biosynthetic demand. Electron surplus for MFC was calculated from NADH balance, reflecting cytosolic NADH remaining after fermentation and biosynthesis.

Interactive double mutant prediction tool allows users to toggle any combination of up to 2 deletions and observe real-time predictions for 6 metabolic phenotypes plus NADH closure percentage. NDE1 deletion (rank 1, q_MFC = +0.072 mmol/(g·h)) blocks respiratory NADH oxidation, forcing 17% increased ethanol production as compensatory fermentation. Moderate biosynthetic reduction (γ_bio = 1.97 vs 3.91 mmol/g) reflects adaptation without severe stress, yielding net NADH surplus with 97.5% closure. GPD1 deletion (rank 2, q_MFC = +0.050 mmol/(g·h)) eliminates 94% of glycerol production, freeing NADH from this major sink. Biosynthetic downregulation (γ_bio = 1.11 mmol/g) represents compensatory stress, causing 15% growth reduction confirmed by laboratory observations, but maintains 97.6% NADH closure with reproducible surplus.

Genotype Playground

Genotype

GPD1 deletion

I_GPD1 (0/1)

GPD2 deletion

I_GPD2 (0/1)

NDE1 deletion

I_NDE1 (0/1)

NDE2 deletion

I_NDE2 (0/1)

Load measured genotype

✓

Feasible

WT (Dom90) measured

i

μ (growth)

0.1815 h⁻¹

q_glc

1.4097 mmol/(g·h)

q_eth

2.0288 mmol/(g·h)

q_gly

0.0818 mmol/(g·h)

NADH Balance

100.0%

NADH Flow

Sources

Sinks

GPD2 and NDE2 show no benefit because they encode minor isoforms with limited flux. GPD2 deletion reduces glycerol by only 20% (vs GPD1's 94%), creating stress without NADH savings (q_MFC = -0.044). NDE2 has weaker respiratory effects than NDE1, failing to force sufficient compensatory fermentation while triggering maladaptive biosynthetic reduction. Only disrupting dominant flux routes (GPD1, NDE1) yields beneficial MFC phenotypes.

Predicted double mutants fail. GPD1-NDE1 and GPD2-NDE1 simultaneously eliminate glycerol sink and respiratory oxidation, removing all NADH disposal routes. Cells respond by collapsing central metabolism—biosynthesis crashes to γ_bio = 0.80 mmol/g (80% reduction), glucose uptake drops, and NADH production falls below consumption. Resulting deficits (q_MFC = -0.272 and -0.333) indicate thermodynamically impossible states where cells cannot balance reducing equivalents. Ethanol yields reach 2.17 mol/mol (exceeding theoretical 2.0 maximum), confirming unsustainable metabolism. The model correctly predicts this failure, validating that cells require at least one NADH disposal pathway for viability. Triple mutants (GPD1-GPD2-NDE1) would exhibit synthetic lethality and were not modeled.

Bootstrap resampling (1000 iterations) propagated measurement uncertainty through the modeling pipeline. Biological replicates were resampled, rates recalculated, models refitted, and strains re-ranked. Predicted strains received 20% uncertainty for extrapolation risk. Strains with positive q_MFC maintained surplus in >95% of iterations (low variance, genuine availability), while negative q_MFC strains showed consistent deficits in 99% of iterations (systematic metabolic failure). This asymmetry validates that positive surplus is thermodynamically feasible and biologically sustainable, while deficits violate conservation of mass.

Figure 1.5: Bootstrap rank stability heatmap shows rank probability distributions across 1000 iterations. Dark green (>70%) indicates high stability, light green (30-70%) moderate stability, yellow (<30%) high uncertainty. NDE1 achieves rank 1-2 in 89% of iterations (mean 1.5 ± 0.8), demonstrating robust superiority. GPD1 ranks 2-3 with 90% top-3 probability (mean 2.5 ± 0.8). GPD1-GPD2 maintains rank 3 (mean 3.0 ± 0.7) because both isoforms target the same pathway without catastrophic redox imbalance. Dom90 wild-type shows distributed ranking (mean 3.8 ± 1.7), reflecting baseline metabolism. Double mutants GPD1-NDE1 and GPD2-NDE1 rank 7-8 in 98% of iterations (means 7.3 and 7.8), confirming worst performance due to severe NADH deficits.

Figure 1.5. Bootstrap rank stability heatmap of measured and predicted strains.

All 8 strains were ranked by electron surplus. Tier 1 excellent candidates (positive q_MFC, >90% top-3 probability): NDE1 (+0.072, rank 1) blocks respiration while maintaining glycerol and viable biosynthesis; GPD1 (+0.050, rank 2) eliminates major NADH sink despite 15% growth penalty; GPD1-GPD2 (+0.038, rank 3) combines both isoforms without catastrophic failure. Tier 2 baseline: Dom90 (0.000, rank 4) shows balanced metabolism with no surplus. Tier 3 failed candidates (NADH deficits): GPD2 (-0.044, rank 5), NDE2 (-0.070, rank 6), GPD1-NDE1 (-0.272, rank 7), and GPD2-NDE1 (-0.333, rank 8)—double mutants exhibit largest deficits as biosynthesis collapses (γ_bio → 0.80 mmol/g), creating thermodynamically impossible consumption exceeding production. NDE1 emerged as optimal by forcing compensatory fermentation (+17% ethanol), maintaining growth (μ = 0.182 h⁻¹), preserving biosynthesis (γ_bio = 1.97 mmol/g), and achieving 97.5% NADH closure. High rank stability (89% top-2 probability) provides confidence for experimental validation. Systematic double mutant failure demonstrates that successful MFC engineering requires preserving at least one NADH disposal pathway—eliminating both glycerol synthesis AND respiration causes metabolic catastrophe rather than enhanced electron export.

Model Assumptions

Over 0-24h window, cells are in pseudo-steady-state (exponential phase)
Cells maintain redox balance; NADH does not accumulate
Ethanol, glycerol, and biosynthesis account for >90% of NADH consumption
No overflow metabolism to acetate or other byproducts
Genetic effects are additive (no epistasis)
Residuals are independent across strains
Constant variance of residuals
Residuals follow normal distribution
Under fermentative conditions, respiratory NADH oxidation via complex I is minimal
Acetate, succinate, and other organic acids contribute <5% to carbon balance
Biosynthetic demand doesn't change with growth conditions (valid for same medium)

Discussion

The Data-Driven Redox Budget (DDRB) model predicts MFC performance from fermentation measurements by ensuring thermodynamic NADH conservation. The core principle: all glycolytic NADH (2 per glucose) must be consumed by fermentation products (ethanol, glycerol), biosynthesis, or remain available for extracellular electron transfer. We implemented a genotype-dependent biosynthetic demand model (γ_bio = α + β_GPD·I_GPD + β_NDE·I_NDE) because preliminary analysis revealed 5-fold variation across strains (0.80-3.91 mmol/g), far exceeding measurement error and reflecting metabolic adaptation—deletion mutants with impaired redox pathways downregulate anabolic NADH consumption as compensation (β_GPD = -2.80, β_NDE = -1.93). The linear framework provides parsimony with only 13 parameters from 20 measurements while capturing systematic genetic effects, outperforming constant γ_bio models by information criteria (ΔAIC = +6.4, 24:1 evidence ratio, R² = 0.91).

Genetic indicator pooling was necessary due to sample size constraints. Both NDE1 and NDE2 deletions were encoded with a unified indicator (I_NDE = 1), and GPD1/GPD2 were similarly pooled as I_GPD = max(I_GPD1, I_GPD2) for the γ_bio model. With only n=5 strains and 2 replicates each, estimating separate coefficients for each isoform would exhaust degrees of freedom and risk overfitting. The pooled model estimates the average effect of deleting genes within each functional class (glycerol synthesis or respiratory NADH oxidation), assuming isoforms impact metabolism in qualitatively similar directions. Individual strain-specific variation between isoforms is captured in model residuals—for example, NDE1 and NDE2 show 33% difference in measured γ_bio (2.37 vs 1.58 mmol/g) despite sharing the same I_NDE indicator. These residuals contribute to the 20% prediction uncertainty applied to double mutants, acknowledging unmodeled isoform-specific effects. For phenotype-specific models (μ, q_glc, q_eth, q_gly), we retained separate GPD1/GPD2 indicators where appropriate to capture the 94% vs 20% glycerol reduction difference, but the γ_bio model used pooled indicators to maintain statistical power. This parsimonious approach balances biological interpretability with reliable parameter estimation given limited training data.

NADH closure of 100.2 ± 2.4% across all strains provides independent validation that the model satisfies conservation laws rather than fitting noise—cells cannot violate mass balance, so deviations from 100% reveal genuine electron surplus (positive q_MFC) or measurement errors and metabolic failure (negative q_MFC). This approach transforms fermentation data into MFC predictions without requiring electrochemical measurements, enabling rapid strain screening. Dom90 wild-type serves as the reference strain because it represents unperturbed metabolism with intact glycerol and respiratory pathways, establishing baseline NADH distribution: 72% to ethanol, 2.9% to glycerol, 25% to biosynthesis, with perfect balance (q_MFC = 0.000 mmol/(g·h)). Using this common reference enables direct quantification of genetic perturbation effects and validates that improvements represent genuine metabolic rewiring. Dom90 also serves as the intercept (α = 3.91 mmol/g) in genotype-dependent models, representing baseline biosynthetic demand when no deletions are present. The three dominant NADH sinks (ethanol >70%, glycerol 3-10%, biosynthesis 15-25%) account for >95% of electron flux, with minor pathways subsumed into residuals, ensuring the model captures major metabolic routes while respecting fundamental constraints against indefinite NADH accumulation.

Positive surplus represents genuine electron availability because deletion mutants eliminate specific high-flux sinks while compensatory mechanisms are incomplete. GPD1 deletion eliminates 94% of glycerol production (0.077 mmol/(g·h) pathway removed), but biosynthetic compensation only saves 0.051 mmol/(g·h) (γ_bio: 3.91 → 1.11 mmol/g at μ = 0.164 h⁻¹), leaving net surplus of +0.050 mmol/(g·h). The systematic correlation between genetic perturbations and q_MFC proves this is not random error—GPD1 and NDE1 deletions consistently produce positive surplus with 89-100% ranking stability across bootstrap iterations. If surplus were measurement noise, we would expect random scatter with no genotype dependence. Laboratory observations independently validated predictions: GPD1 mutants showed 15% growth reduction, confirming that glycerol pathway elimination creates metabolic stress exactly as predicted by γ_bio reduction. Conversely, negative surplus indicates metabolic failure where cells cannot balance redox equivalents, violating thermodynamic constraints. Double mutants GPD1-NDE1 and GPD2-NDE1 showed severely negative q_MFC (-0.272 and -0.333 mmol/(g·h)), confirming catastrophic failure when both glycerol sink and respiratory oxidation are eliminated. These strains exhibit biosynthetic collapse (γ_bio → 0.80 mmol/g, 80% reduction), creating thermodynamically impossible states. The model's ability to correctly predict this failure (98% consistency across bootstrap iterations) validates its biological accuracy and demonstrates that metabolic engineering must preserve at least one NADH disposal route.

NDE1 emerged as the optimal MFC candidate (+0.072 mmol/(g·h), rank 1 with 89% top-2 probability), blocking respiratory NADH oxidation to force 17% increased ethanol production while maintaining viable growth (μ = 0.182 h⁻¹) and moderate biosynthesis (γ_bio = 1.97 mmol/g). GPD1 ranked second (+0.050 mmol/(g·h), 90% top-3 probability) by eliminating 94% of glycerol production, confirming it as the dominant isoform and freeing NADH from this major sink despite 15% growth penalty verified by laboratory measurements. GPD2 and NDE2 showed no benefit because they encode minor isoforms—GPD2 deletion reduces glycerol by only 20%, creating stress without sufficient NADH savings (q_MFC = -0.044 mmol/(g·h)). The critical insight: only disrupting dominant metabolic flux routes (GPD1, NDE1) yields beneficial MFC phenotypes.

Double mutants failed, contradicting initial hypotheses of synergistic improvements. GPD1-NDE1 and GPD2-NDE1 ranked 7-8 with NADH deficits indicating metabolic collapse—simultaneous elimination of glycerol sink and respiratory oxidation removes all NADH disposal routes, forcing cells to shut down central metabolism rather than accumulate electron surplus. This systematic failure demonstrates that successful metabolic engineering requires preserving at least one major NADH outlet, explaining why triple mutants would exhibit synthetic lethality.

Model predictions showed strong agreement with experimental measurements: growth rates 0.14-0.18 h⁻¹, ethanol yields 1.44-1.77 mol/mol, and GPD1 dominance (94% reduction) matched literature values for S. cerevisiae fermentation. NADH balance 97.5-102.9% validated thermodynamic conservation across all strains. Mean model R² = 0.79 across all phenotypes with residuals within ±2σ and bootstrap ranking stability 78-100% for top candidates despite small sample size (n=5) confirmed statistical robustness. Laboratory confirmation of GPD1 growth impairment and preliminary electrochemical tests of NDE1 provide orthogonal validation. The model successfully identified viable single-deletion candidates while correctly predicting double-mutant failure, providing a validated framework for data-driven MFC strain engineering that bridges fermentation physiology and bioelectrochemical performance without requiring electrochemical screening of every candidate.

The modeling approach provided practical advantages over traditional experimental screening. Testing all candidate strains electrochemically in MFC reactors would require 8-12 weeks per strain (reactor assembly, biofilm formation, stabilization, and performance characterization), consuming significant resources and laboratory time. By predicting MFC performance from 24-hour fermentation assays, we compressed the screening timeline to under one week while correctly identifying top candidates (NDE1, GPD1) and avoiding experimental investment in failed designs. Most critically, the model prevented pursuing double mutants GPD1-NDE1 and GPD2-NDE1, which initial intuition suggested would show synergistic improvements.

Population/System level: MFC Biofilm Dynamics Mode (Bridging Metabolic Flux to Power Output)

While the DDRB model successfully predicts intracellular electron availability (q_MFC) from fermentation measurements, translating this single-cell metabolic flux into macroscopic current generation requires understanding population-level dynamics at the electrode-biofilm interface. Microbial fuel cells operate as living electrochemical devices where power output emerges from the collective behavior of millions of cells organized into electroactive biofilms.

The MFC Biofilm-Electrode Model addresses three fundamental questions that cannot be answered by metabolic flux analysis alone. First, regarding biofilm colonization dynamics, how do metabolic perturbations from GPD and NDE deletions affect the transition from planktonic cells to electrode-attached biofilms, and what timescale governs biofilm maturation under batch cultivation? Second, concerning power generation kinetics, what current density and power output can be achieved under realistic operating conditions—including the 1000 Ω external resistance and inherently modest extracellular electron transfer efficiency of Saccharomyces cerevisiae—given a strain's intrinsic electron surplus? Third, and most critically for strain selection, can strains with maximal q_MFC sustain competitive biofilm formation despite growth defects, or do moderate-flux strains with minimal metabolic burden achieve superior performance through faster electrode colonization?

The model architecture combines mechanistic rigor with practical calibration constraints, implementing a coupled ordinary differential equation system that describes planktonic growth, biofilm attachment and detachment, substrate consumption via Monod kinetics, and electrochemical power generation through an empirical scaling law. The system is parameterized using three strain-specific inputs from DDRB analysis—growth rate μ, glucose uptake rate q_glc, and electron surplus q_MFC—ensuring thermodynamic consistency from intracellular metabolism to electrode current. Critically, we incorporate metabolic burden-dependent attachment kinetics through the relationship k_attach = k₀/(1 + β·burden), where burden equals μ_max - μ_strain, capturing the observation that deletion mutants with impaired growth exhibit slower biofilm formation due to altered cell surface properties and stress-response activation [10, 12].

The electrochemical component employs a simplified

, where κ is a global power coefficient, α is the flux-to-power scaling exponent, X_biofilm is the electrode-attached biomass density from ODE integration, and f_mature = 1 - exp(-t/τ) describes biofilm maturation with τ = 18 hours. This formulation avoids over-parameterization by collapsing all complex electrochemistry—Butler-Volmer kinetics, mass transfer limitations, electrode geometry effects—into just two calibrated parameters (κ, α), which are fitted to match two experimental data points from Dom90 wild-type (32 μW) and GPD1 mutant (51 μW). The resulting calibration yields κ = 84.83 μW·h^α/mmol^α and α = 0.347, indicating strong sublinear scaling where doubling metabolic flux increases power by only 27% rather than 100%, reflecting fundamental electrode saturation effects.

All biofilm kinetic parameters are fixed from literature rather than fitted to data, ensuring genuine predictions rather than overfitting. Attachment rate k_attach = 0.080 h⁻¹, detachment rate k_detach = 0.005 h⁻¹, and death rate k_death = 0.002 h⁻¹ follow Picioreanu et al. (2007) studies of biofilms in electrochemical systems. The metabolic burden coefficient β = 3.0 produces physiologically reasonable attachment penalties, while electrode geometry (4 cm²), external resistance (1000 Ω), and reactor volume (50 mL) match our experimental MFC prototype. The model implements automatic viability filtering: strains with negative q_MFC are excluded because NADH consumption exceeds production, eliminating four strains (NDE2, GPD2, GPD1-NDE1, GPD2-NDE1) and leaving four viable candidates (Dom90, GPD1, NDE1, GPD1-GPD2) for experimental testing. This integration with DDRB modeling provides both quantitative predictions and qualitative biological constraints that prevent exploration of thermodynamically infeasible genotypes.

Population dynamics with metabolic burden coupling

Our first modeling objective is to predict how deletion mutants transition from planktonic to biofilm states while accounting for growth defects that arise from metabolic perturbations. This is crucial because high electron surplus does not guarantee high power output if the strain cannot effectively colonize the electrode surface [8, 9]. The population dynamics are governed by two coupled differential equations that track planktonic and biofilm-associated cell populations separately.

For planktonic cells, measured in grams per liter of reactor volume, the rate of change is described by:

This equation captures three biological processes: growth of planktonic cells limited by substrate availability, loss of planktonic cells due to attachment to the electrode, and gain of planktonic cells from biofilm detachment events. Similarly, biofilm cells follow:

Here we account for recruitment from the planktonic population, detachment back into solution, and cell death within the biofilm due to nutrient limitation in deeper layers. The effective growth rate

implements Monod-limited growth kinetics with a half-saturation constant K_s = 0.1 g/L, ensuring that growth slows as glucose becomes depleted. We operate in batch mode with D = 0 h⁻¹ (no dilution), mimicking our experimental MFC setup. The detachment and death rate constants are set to k_detach = 0.005 h⁻¹ and k_death = 0.002 h⁻¹ based on literature values for yeast biofilms [12, 15, 16].

The critical innovation in our model is making the attachment rate dependent on metabolic burden through the relationship:

where k_attach,max = 0.08 h⁻¹ represents the wild-type attachment rate, β_burden = 3.0 is a fixed parameter that quantifies the sensitivity of attachment to metabolic stress, and burden = μ_max - μ_strain measures the growth deficit relative to the fastest-growing strain. This formulation has strong biological justification: metabolically stressed cells exhibit altered cell surface properties, reduced adhesin expression, and impaired biofilm matrix production. To illustrate the magnitude of this effect, consider the GPD1-GPD2 double mutant which shows μ = 0.154 h⁻¹ compared to Dom90's 0.182 h⁻¹, yielding a burden of 0.028. This 15% growth defect reduces k_attach from 0.080 to 0.073 h⁻¹, a modest 9% reduction that reflects the documented observation that moderate metabolic perturbations cause proportionally smaller effects on biofilm adhesion capacity. We also implement a biofilm saturation constraint: when X_bf > X_bf,max = 50 g/L, we enforce dX_bf/dt ≤ 0 to prevent unrealistic accumulation beyond physical limits imposed by diffusion constraints and electrode surface area [13].

Results

Figure 2.1 reveals the fundamental biofilm colonization dynamics that enable electrode-based electron harvesting across the four viable strains identified by DDRB thermodynamic analysis. All strains successfully established biofilms, but with critical kinetic differences reflecting their metabolic burden profiles and substrate consumption rates. Wild-type Dom90 (black line) achieved steady biofilm accumulation beginning at approximately t = 15h, reaching 1.50 g/L by t = 72h with a characteristic exponential-then-plateau growth trajectory. This represents optimal colonization with k_attach = 0.080 h⁻¹ under zero metabolic burden, serving as the reference baseline for evaluating mutant performance.

The NDE1 deletion strain (red line) exhibited nearly identical colonization kinetics to Dom90, ultimately achieving 1.45 g/L by t = 72h—statistically indistinguishable from wild-type. This behavior emerges because NDE1 shares identical growth rate with Dom90 (μ = 0.182 h⁻¹ for both strains), resulting in zero calculated metabolic burden and thus identical attachment kinetics. The model correctly predicts that metabolic rewiring targeting NADH surplus does not inherently impair biofilm formation when growth rate remains unaffected, validating that the NDE1 deletion successfully decouples electron availability from cellular fitness. The overlapping trajectories through t = 60h confirm this prediction quantitatively, with only minor divergence in the final plateau phase attributable to stochastic effects in the numerical integration.

GPD1 (yellow line) revealed the clearest signature of the metabolic burden penalty encoded in the model. Biofilm accumulation proceeded more slowly throughout the entire 72-hour time course, reaching only 1.73 g/L—still competitive with other strains but reflecting the combined effects of reduced growth rate (μ = 0.164 h⁻¹, burden = 0.018) and correspondingly reduced attachment rate (k_attach = 0.076 h⁻¹). Paradoxically, GPD1 achieved the highest final biofilm density despite this handicap. This apparent contradiction is explained by temporal coupling with substrate availability: GPD1's intermediate glucose consumption rate (visible in panel G) extended the growth phase relative to Dom90 and NDE1, which depleted substrate by t = 25h. The additional 5-8 hours of substrate availability allowed GPD1 to accumulate biomass through prolonged exponential growth, compensating for its slower attachment kinetics through extended colonization time. This demonstrates a key insight: biofilm density at any given timepoint reflects the integrated history of growth, attachment, and substrate availability rather than instantaneous metabolic capacity alone.

The GPD1-GPD2 double mutant (purple dashed line) displayed the most distinctive colonization profile, with visibly slower initial accumulation due to the highest metabolic burden (burden = 0.028, k_attach = 0.073 h⁻¹) but ultimately reaching 1.78 g/L—the highest of all strains by t = 72h. This biphasic behavior emerges from GPD1-GPD2's exceptionally low glucose consumption rate, which extends substrate availability to t = 32h, providing nearly 50% more growth time than fast-consuming strains. The model captures this trade-off mechanistically: slower attachment is overcome by prolonged growth opportunity, resulting in competitive final biofilm density despite kinetic disadvantages. The convergence of all four strain trajectories to 1.45-1.78 g/L by t = 72h—a remarkably narrow 20% range—suggests that the biofilm saturation constraint effectively equalizes long-term electrode colonization across genotypes with diverse metabolic profiles.

Critically, no strain approached the theoretical maximum X_bf,max = 50 g/L, indicating that biofilm growth remained limited by biological processes (substrate depletion, detachment, death) rather than physical space constraints. All strains crossed the minimum activation threshold X_threshold = 0.1 g/L within 10-15 hours, establishing that measurable current generation should commence well before t = 20h in all cases. The smooth, monotonic accumulation profiles without oscillations or instabilities confirm numerical stability of the ODE integration and validate that the chosen parameter values produce biologically plausible dynamics. The slight divergence between strains visible at t = 30-50h reflects the interaction between growth rate, attachment rate, and substrate depletion timing—precisely the coupled dynamics that necessitate mechanistic modeling rather than simple flux-to-power correlations.

The observation that GPD1 generates 1.59× more power than Dom90 (51.0 versus 32.0 μW) despite achieving only 1.15× higher peak biofilm density (1.73 versus 1.50 g/L) validates the model's central premise that power output scales with the product of electron surplus and bioactive biomass, not biomass alone. The 2.50× ratio in DDRB-predicted q_MFC values (0.050 versus 0.020 mmol/(g·h)) undergoes sublinear transformation through the α = 0.347 scaling exponent to yield approximately 1.27× power contribution from flux alone, which then multiplies by the 1.15× biofilm advantage to produce the observed 1.46× net power ratio—remarkably close to the experimental 1.59× value. This quantitative agreement demonstrates that biofilm dynamics and metabolic flux contributions are correctly weighted in the power scaling law, lending confidence to predictions for untested strains like NDE1.

Substrate consumption dynamics

The objective of this modeling step is to couple glucose depletion to biomass growth using the experimentally measured uptake rates from the DDRB model. This ensures consistency between the metabolic flux predictions and the population-level substrate balance. Glucose concentration changes according to:

where S_feed = 20 g/L represents the initial glucose concentration at inoculation. The effective glucose uptake rate is

, implementing substrate-limited uptake that matches the Monod kinetics used for growth. The strain-specific uptake rate q_glc comes directly from DDRB measurements, ensuring that our predictions remain grounded in experimentally validated metabolic physiology rather than fitted parameters. This direct coupling means that any discrepancies between predicted and observed substrate depletion would indicate errors in the DDRB flux estimates themselves, providing an independent validation pathway for the upstream metabolic model.

Results

Figure 2.2 demonstrates the metabolic flux coupling between DDRB-predicted uptake rates and population-level glucose depletion across the four viable strains. Wild-type Dom90 (black line) exhibited the most rapid consumption kinetics, achieving 50% depletion (S = 10 g/L) by approximately t = 12h and complete glucose exhaustion (S < 0.5 g/L) by t = 23h. This aggressive depletion reflects Dom90's combination of highest growth rate (μ = 0.182 h⁻¹) and robust glucose uptake (q_glc = 0.911 mmol/(g·h) from DDRB measurements), producing rapid biomass accumulation that consumes substrate at the maximum rate. The steep initial slope during exponential phase (t = 5-15h) transitions to a more gradual decline as substrate limitation begins to engage Monod kinetics, eventually reaching complete depletion that triggers entry into stationary phase where biofilm growth plateaus as observed in panel A.

NDE1 (red line) displayed nearly identical depletion kinetics to Dom90 throughout the entire time course, with both strains reaching exhaustion at t = 23h within the numerical precision of the integration. This perfect overlap confirms the DDRB model's prediction that NDE deletions redirect electron flux toward NADH surplus without significantly perturbing fermentative glucose metabolism. The shared growth rates (μ = 0.182 h⁻¹ for both) and similar specific uptake rates (q_glc = 0.911 mmol/(g·h)) produce mathematically equivalent substrate consumption profiles, validating that respiratory pathway modifications can be decoupled from primary carbon catabolism. This demonstrates a key advantage of the NDE1 engineering strategy: enhanced electron availability for MFC operation without compromising the fundamental fermentative capacity that drives biomass accumulation.

The glycerol pathway mutants revealed critical temporal offsets that explain their divergent biofilm dynamics. GPD1 (yellow line) maintained S > 10 g/L until approximately t = 16h and exhausted glucose only at t = 28h—a 5-hour delay relative to Dom90 and NDE1. This extended substrate availability is mechanistically explained by GPD1's 11% lower growth rate (μ = 0.164 versus 0.182 h⁻¹) combined with comparable glucose uptake rate, producing a net reduction in volumetric consumption rate that scales with the lower biomass accumulation. The extended growth window allows GPD1 to partially compensate for its slower attachment kinetics (as seen in panel A), ultimately achieving competitive biofilm density by t = 72h despite kinetic disadvantages during the initial colonization phase.

GPD1-GPD2 double mutant (purple dashed line) exhibited the most extreme substrate persistence, maintaining S > 10 g/L until t = 20h and not reaching complete exhaustion until t = 32h—nearly 10 hours later than Dom90/NDE1. This dramatically extended exponential phase emerges from the compound effect of reduced growth rate (μ = 0.154 h⁻¹, the lowest among viable strains) and the cumulative metabolic perturbations of dual deletion. The glucose availability extends well beyond the point where fast-consuming strains enter stationary phase, providing GPD1-GPD2 with additional biofilm accumulation time that explains its highest final density in panel A despite the slowest attachment rate. This temporal offset demonstrates a critical principle: metabolic burden reduces instantaneous performance metrics (growth rate, attachment rate) but can paradoxically enable superior long-term outcomes by extending resource availability and avoiding early growth arrest.

The convergence of all depletion curves to complete exhaustion (S ≈ 0 g/L) by t = 32h establishes a clear demarcation between the growth-limited regime (t < 25-32h depending on strain) and the stationary regime (t > 32h) where biofilm dynamics become decoupled from substrate availability. During the stationary phase, biofilm density changes are governed solely by the balance of detachment and death (k_detach = 0.005 h⁻¹, k_death = 0.002 h⁻¹), leading to the plateau behavior visible in panel A after t = 40h. The substrate depletion profiles validate that the specific rates measured in DDRB shake flask fermentations accurately predict population-level behavior in the MFC batch reactor, confirming that parameter transfer between modeling frameworks is quantitatively sound.

Correlation with power generation timing: The fact that all strains exhaust glucose by t = 32h while continuing to generate power through t = 72h (panel B) validates that electrochemical activity persists in stationary phase cells, consistent with yeast maintaining metabolic activity on intracellular reserves (glycogen, trehalose) even after external carbon depletion. The staggered substrate exhaustion times (t = 23h for Dom90/NDE1 versus t = 32h for GPD1-GPD2) do not translate to correspondingly staggered power outputs because the power scaling law depends on biofilm density at peak power time (typically t = 45-55h) rather than the instantaneous growth rate during the exponential phase. This demonstrates why mechanistic modeling is essential: simple correlations between substrate consumption and power output would incorrectly predict that slower-consuming strains generate less power, whereas the model correctly captures that extended growth compensates for reduced kinetics, producing similar final performance across diverse metabolic profiles.

Electrochemical current generation

This is the critical modeling step where we convert intracellular electron flux (q_MFC) into measurable electrical power via the biofilm-electrode interface. Unlike the biofilm and substrate dynamics, which follow deterministic ODE trajectories with literature-validated parameters, the electrochemical component requires a simplified empirical approach because the mechanistic complexity of extracellular electron transfer in yeast exceeds what can be reliably parameterized from first principles.

The power output is calculated using the calibrated scaling law:

P = κ · (q_MFC)^α · X_bf · f_maturity

where κ = 84.83 μW·h^α/mmol^α is the global power coefficient and α = 0.347 is the flux-to-power scaling exponent, both determined by fitting to experimental measurements of Dom90 (32.0 μW) and GPD1 (51.0 μW). The maturity factor

with τ = 18 hours captures the time-dependent development of biofilm electroactivity, scaling electron export from 0% at t = 0 to 63% at t = 18h and 95% at t = 54h. The biological rationale for this maturity factor is that young biofilms lack well-developed conductive pathways and optimal electrode contact, preventing unphysical current generation immediately after inoculation when cells have not yet established proper electrical connectivity to the electrode surface.

The α = 0.347 sublinear scaling exponent represents the most important mechanistic insight from the calibration: power does not scale linearly with metabolic flux, but rather exhibits strong saturation behavior where doubling q_MFC increases power by only 27% rather than 100%. This sublinearity reflects fundamental electrode limitations including finite surface area for electron acceptance, mass transfer constraints in delivering reduced cofactors from cell interior to the electrode interface, and possible product inhibition from accumulated oxidized species. The departure from linearity becomes more severe at high flux: a strain with 3× Dom90's electron surplus would theoretically generate only 1.59× the power according to 3^0.347 = 1.59, demonstrating that metabolic engineering strategies targeting extreme q_MFC values face diminishing returns constrained by physical electrochemistry rather than biological capacity.

A critical threshold is applied to ensure physical realism: when X_bf < X_threshold = 0.1 g/L, power is set to zero regardless of q_MFC. This recognizes that sparse biofilms cannot sustain electrical connectivity across the electrode surface, preventing false predictions where strains with high electron surplus but poor colonization would incorrectly generate power before sufficient biofilm establishes percolation pathways for electron conduction. All four viable strains exceed this threshold by t = 15h, establishing that measurable current should commence during the exponential growth phase rather than requiring full biofilm maturation.

Results

Figure 2.3 reveals the temporal dynamics of electrical power generation across the four viable strains, demonstrating how metabolic flux predictions from DDRB translate through biofilm dynamics into device-level performance. Plotted on logarithmic scale to capture the dynamic range from sub-microwatt to 50+ μW, the power trajectories exhibit characteristic three-phase kinetics that validate the model's mechanistic formulation. Phase I (t < 10h) shows minimal power output (<1 μW) despite cells having colonized the electrode, because the maturity factor f_maturity < 0.4 limits electron transfer through incompletely developed conductive networks. Phase II (t = 10-50h) exhibits exponential power increase visible as linear slope on the log-scale plot, driven by the multiplicative combination of exponentially growing biofilm density from panel A and exponentially approaching maturity factor. Phase III (t > 50h) reaches power plateau as both biofilm growth and maturity factor saturate, with final values reflecting the integrated effect of q_MFC, biofilm density, and temporal dynamics.

GPD1 (yellow line) achieved the highest predicted power output, reaching 51.0 μW by t = 72h—precisely matching the experimental calibration target by construction. The trajectory shows smooth exponential rise from 1 μW at t = 18h to 35 μW at t = 40h, then gradual approach to the final plateau. This kinetic profile directly reflects GPD1's intermediate q_MFC = 0.050 mmol/(g·h) transformed through the α = 0.347 scaling to yield (0.050)^0.347 = 0.332 in normalized units, multiplied by biofilm density that peaks at 1.73 g/L and maturity factor approaching 0.95. The fact that GPD1 outperforms all other strains despite having lower electron surplus than NDE1 (0.050 versus 0.072 mmol/(g·h)) validates the model's central premise: power emerges from the coupled product of flux, biofilm, and time, not from any single factor in isolation.

GPD1-GPD2 double mutant (purple dashed line) achieved 48.9 μW, ranking second despite having only 0.038 mmol/(g·h) electron surplus—24% lower than GPD1. This near-equivalence (48.9 versus 51.0 μW, only 4% difference) emerges from compensatory dynamics where GPD1-GPD2's higher final biofilm density (1.78 versus 1.73 g/L, a 3% advantage from panel A) partially offsets its flux deficit after transformation through the sublinear scaling. The temporal trajectory shows slightly delayed power onset, with GPD1-GPD2 crossing 10 μW at t = 24h versus t = 22h for GPD1, reflecting the metabolic burden-induced attachment penalty. However, by t = 50h the trajectories converge to within 5%, demonstrating that long-term power output is remarkably insensitive to moderate variations in growth rate and attachment kinetics when substrate availability compensates through extended accumulation time.

NDE1 (red line) achieved 48.4 μW, ranking third in a statistical dead heat with GPD1-GPD2 despite having the highest electron surplus among all strains (0.072 mmol/(g·h)). This represents the model's most counterintuitive and important prediction: the strain with 44% higher flux than GPD1 generates 5% less power. The explanation lies in the sublinear scaling law: NDE1's flux advantage transforms to (0.072/0.050)^0.347 = 1.127, only 13% higher in power space, while its biofilm density of 1.45 g/L is 16% lower than GPD1's 1.73 g/L, yielding a net disadvantage of 1.127 × (1.45/1.73) = 0.945, or 5.5% lower predicted power. This quantitative agreement demonstrates that the α = 0.347 calibration correctly captures the electrode saturation phenomenon, where incremental flux increases face diminishing returns that can be overcome by superior biofilm formation.

Dom90 wild-type (black line) reached only 32.0 μW—the experimental calibration baseline—reflecting its minimal electron surplus of 0.020 mmol/(g·h) estimated as 40% of GPD1's value based on complete NADH consumption in wild-type metabolism. Despite achieving competitive biofilm density (1.50 g/L) and identical maturity kinetics to other strains, the low intrinsic flux creates an insurmountable bottleneck. The power trajectory shows delayed onset relative to metabolically engineered strains, crossing 10 μW only at t = 28h compared to t = 22-24h for GPD1/GPD1-GPD2/NDE1, demonstrating that even with optimal colonization, insufficient metabolic capacity limits device startup time. The 51.0/32.0 = 1.59× power ratio between GPD1 and Dom90 matches the experimental observation precisely, validating that the calibrated model successfully captured the dominant mechanisms governing strain-to-strain performance variation.

The convergence of all three top performers (GPD1, GPD1-GPD2, NDE1) to a narrow 48.4-51.0 μW range—only 5% spread—represents a critical finding for experimental strain selection: attempting to maximize electron surplus beyond GPD1's 0.050 mmol/(g·h) level yields marginal gains due to electrode saturation, while introducing growth defects through aggressive metabolic engineering (e.g., triple deletions) risks biofilm formation penalties that eliminate any flux-based advantage. This suggests an optimal engineering strategy of moderate metabolic rewiring (single deletions like GPD1 or NDE1) that achieves 50-60% of theoretical maximum flux while preserving >90% of wild-type fitness, rather than pursuing extreme flux at the cost of cellular viability.

The model achieves 0.0% calibration error for both Dom90 (predicted 32.0 μW, experimental 32.0 μW) and GPD1 (predicted 51.0 μW, experimental 51.0 μW) by construction, establishing confidence in the parameter estimates κ and α. The critical test of predictive power lies in the untested strains: NDE1's predicted 48.4 μW and GPD1-GPD2's predicted 48.9 μW both fall within the ±30% uncertainty interval (33.9-62.9 μW and 34.3-63.6 μW respectively), providing falsifiable hypotheses for experimental validation. If measured NDE1 power deviates beyond this range, it would indicate either strain-specific electron export efficiency variations not captured by the uniform η assumption, or non-additive genetic interactions between NDE deletions and biofilm adhesion machinery that violate the simple burden-dependent attachment model. The narrow clustering of predictions (48-51 μW for top three strains) makes experimental discrimination challenging, requiring measurement precision better than 10% to distinguish performance differences—a consideration for experimental protocol design.

Numerical integration and convergence

To solve this coupled system of ordinary differential equations, we utilized the LSODA algorithm implemented in scipy.integrate.solve_ivp, which automatically switches between methods optimized for stiff and non-stiff regimes. This adaptive approach is essential because the equations exhibit multiple timescales: rapid substrate consumption during exponential growth (characteristic time ~5-10 hours) and slower biofilm maturation approaching steady-state (characteristic time τ = 18 hours). The numerical settings were: relative tolerance 10⁻⁶, absolute tolerance 10⁻⁹, time span 0-72 hours, and 500 evaluation points yielding 0.144h temporal resolution. Initial conditions mimicking inoculation from stationary phase culture were X_p(0) = 0.1 g/L representing diluted seed culture, X_bf(0) = 0.001 g/L reflecting minimal initial electrode colonization, and S(0) = 20 g/L matching standard yeast fermentation medium glucose concentration.

All four viable strains achieved successful integration with zero failed timesteps and no numerical instabilities, indicating the model equations are mathematically well-posed and the parameter choices produce biologically plausible dynamics without unphysical oscillations or runaway growth. Solution stability was confirmed by comparing simulations with 500 versus 1000 evaluation points—predicted power curves differ by less than 0.2% at all timepoints, demonstrating that our temporal resolution is more than adequate to capture the system dynamics without undersampling critical transitions.

The calibration strategy employs a minimalist two-parameter approach designed to avoid overfitting the limited experimental data. Given only two experimental measurements (Dom90 = 32.0 μW, GPD1 = 51.0 μW), we fit exactly two parameters (κ, α) in the power scaling law while fixing all biofilm kinetic parameters from literature. This mathematical constraint—two unknowns determined by two equations—ensures the optimization problem is exactly determined rather than underdetermined, producing a unique solution rather than a family of equally valid parameter sets that would render predictions arbitrary. The optimization employs multi-start L-BFGS-B algorithm with five random initializations to avoid local minima, converging to κ = 84.83 μW·h^α/mmol^α and α = 0.347 with relative errors below 0.01% on both calibration points.

The flux ratio between GPD1 and Dom90 is 0.050/0.020 = 2.50×, yet the experimental power ratio is only 51.0/32.0 = 1.59×. If power scaled linearly with flux (α = 1), we would predict GPD1 power as 32.0 × 2.50 = 80.0 μW, a 56% overestimate. The calibrated α = 0.347 correctly predicts (2.50)^0.347 = 1.27× flux contribution, which multiplies by the biofilm ratio 1.73/1.50 = 1.15× to yield net 1.27 × 1.15 = 1.46× power ratio—within 8% of the observed 1.59×. This close agreement despite no adjustable parameters in the biofilm component validates that the model correctly partitions power generation into independent flux and biofilm contributions, lending confidence to predictions for strains where these factors differ from calibration values.

The α < 1 sublinearity is not merely an empirical fit but reflects well-established electrochemical principles. Electrode saturation occurs when the rate of electron delivery from biofilm exceeds the rate of electron acceptance at the electrode-electrolyte interface, creating accumulation of reduced species that inhibit further electron transfer through product inhibition and diffusion polarization. This phenomenon is universal in bioelectrochemical systems and has been documented across bacterial MFCs, where power-current relationships consistently exhibit sublinear behavior with exponents ranging 0.3-0.7 depending on electrode material and reactor geometry. Our calibrated α = 0.347 falls at the lower end of this literature range, suggesting that graphite felt electrodes in yeast MFCs face particularly severe saturation constraints compared to optimized bacterial systems with specialized nanowire-based electron transfer.

ResultsFigure 2.4 synthesizes all model components into an engineering-focused recommendation framework, ranking the four viable strains by predicted peak power output with uncertainty intervals that reflect combined parameter estimation uncertainty (±20%), model structural assumptions (±15%), and biological variability (±10%). GPD1 (yellow bar) achieves first rank at 51.0 μW by construction as one of the two calibration points, serving as the reference for evaluating other strains. The error bar spanning 35.7-66.3 μW represents the ±30% confidence interval calculated as 51.0 × (0.7 to 1.3), acknowledging that model predictions carry substantial uncertainty from unmodeled factors including electrode surface heterogeneity, batch-to-batch variability in cell physiology, and potential strain-specific electron export variations not captured by the uniform efficiency assumption.

GPD1-GPD2 double mutant (purple bar) ranks second at 48.9 μW (34.3-63.6 μW range), falling within 4% of GPD1 despite 24% lower electron surplus (0.038 versus 0.050 mmol/(g·h)). This near-equivalence emerges from the compound effect of sublinear flux scaling—which reduces GPD1's advantage to only 13% after applying α = 0.347—and GPD1-GPD2's 3% higher biofilm density compensating through the multiplicative structure of the power law. The overlapping uncertainty intervals (GPD1: 35.7-66.3 μW, GPD1-GPD2: 34.3-63.6 μW) indicate these strains are statistically indistinguishable within model precision, requiring experimental measurements with <10% error to reliably discriminate performance. This clustering suggests that metabolic burden effects successfully equalize long-term outcomes across genotypes with diverse instantaneous growth and attachment rates.

NDE1 (red bar) ranks third at 48.4 μW (33.9-62.9 μW), representing the model's most important prediction: the strain with highest electron surplus (0.072 mmol/(g·h), 44% above GPD1) generates marginally lower power due to reduced biofilm density (1.45 versus 1.73 g/L). This 5% power deficit relative to GPD1 is statistically insignificant—both predictions lie within each other's uncertainty intervals—but the directionality challenges the naive assumption that maximizing q_MFC always maximizes device performance. The model quantitatively explains this inversion: NDE1's flux advantage yields only 1.127× in power space after sublinear scaling, while the 16% biofilm disadvantage produces 0.84× penalty, resulting in net 1.127 × 0.84 = 0.946, or 5.4% reduction. This prediction is experimentally falsifiable—if measured NDE1 power significantly exceeds GPD1, it would indicate that the uniform electron export efficiency assumption fails and NDE1 achieves higher extracellular transfer rates through currently unknown mechanisms.

Dom90 wild-type (black bar) anchors the ranking at 32.0 μW (22.4-41.6 μW), the second calibration point representing baseline performance without metabolic engineering. The 1.59× ratio between GPD1 and Dom90 establishes the engineering value proposition: eliminating glycerol synthesis via GPD1 deletion increases power output by 59% while maintaining 90% of wild-type growth rate, representing an efficient trade-off between metabolic perturbation severity and device performance gain. The fact that all three engineered strains (GPD1, GPD1-GPD2, NDE1) cluster within 48-51 μW suggests a natural ceiling where further optimization yields diminishing returns unless fundamentally new electron export mechanisms are introduced.

Engineering recommendation and experimental priorities: The model recommends prioritizing GPD1, NDE1, and GPD1-GPD2 as the top three candidates for experimental MFC construction, based on several integrated criteria. First, all three achieve statistically equivalent power (48-51 μW within uncertainty) while representing orthogonal metabolic strategies—glycerol elimination (GPD1), respiratory rewiring (NDE1), and compound pathway modification (GPD1-GPD2)—providing diversity against unforeseen failure modes. Second, all three maintain >85% of wild-type growth rate, ensuring robust biofilm formation without severe metabolic burden penalties that would compromise electrode colonization or long-term viability. Third, the predicted power enhancement of 1.5-1.6× over Dom90 baseline is large enough to justify experimental effort while remaining modest enough to avoid unrealistic expectations that might arise from overoptimistic linear extrapolation of flux data.

The narrow 5% performance spread among top three strains has practical implications: experimental validation can simultaneously test all three without sequential optimization, and the eventual choice between them can be made on secondary criteria (genetic stability, ease of transformation, industrial scale-up compatibility) rather than marginal power differences within measurement noise. The ±30% uncertainty intervals acknowledge that the model cannot reliably predict which of the top three will experimentally outperform the others, but confidently predicts all three will substantially exceed Dom90 baseline. This honest uncertainty quantification distinguishes scientifically rigorous modeling from overconfident curve-fitting exercises that claim false precision.

Discussion

The coupled ODE system was numerically integrated for the four viable strains (Dom90, GPD1, NDE1, GPD1-GPD2) over 72-hour batch operation using LSODA algorithm with relative tolerance 10⁻⁶ and absolute tolerance 10⁻⁹, producing stable solutions with zero failed timesteps. The model received three strain-specific parameters from DDRB analysis (μ, q_glc, q_MFC) and predicted time-resolved biofilm formation, substrate consumption, and electrical power output through the calibrated scaling law P = κ·(q_MFC)^α·X_bf·f_maturity with κ = 84.83 μW·h^α/mmol^α and α = 0.347.

Biofilm dynamics revealed critical metabolic burden-colonization trade-offs. Wild-type Dom90 achieved fastest colonization reaching 1.50 g/L by t = 40h with k_attach = 0.080 h⁻¹ under zero metabolic burden, depleting glucose by t = 23h. NDE1 exhibited identical kinetics due to shared growth rate (μ = 0.182 h⁻¹), reaching 1.45 g/L with no attachment penalty. GPD1 showed modestly slower accumulation reflecting 11% growth rate reduction (μ = 0.164 h⁻¹, burden = 0.018) that reduced k_attach to 0.076 h⁻¹, but ultimately achieved highest final density (1.73 g/L) through 5-hour extended substrate availability that compensated for kinetic disadvantages. GPD1-GPD2 double mutant displayed most severe burden (0.028) reducing k_attach to 0.073 h⁻¹, yet reached 1.78 g/L through 10-hour glucose extension to t = 32h, demonstrating that temporal dynamics can overcome instantaneous rate penalties. The narrow final density range (1.45-1.78 g/L, only 20% spread) indicates that diverse metabolic profiles converge to similar electrode colonization by stationary phase, with temporal offsets rather than absolute limits distinguishing strains.

Power generation predictions revealed sublinear flux-to-power conversion that inverts naive rankings. The model predicted GPD1 at 51.0 μW (rank #1), GPD1-GPD2 at 48.9 μW (rank #2), and NDE1 at 48.4 μW (rank #3), with Dom90 baseline at 32.0 μW. The critical insight is that NDE1 generates 5% less power than GPD1 despite 44% higher electron surplus (0.072 versus 0.050 mmol/(g·h)) because the α = 0.347 sublinear scaling reduces NDE1's flux advantage to only 13% in power space, while its 16% lower biofilm density (1.45 versus 1.73 g/L) creates net disadvantage. This demonstrates electrode saturation: incremental flux increases face diminishing returns that can be overcome by superior biofilm formation, challenging the assumption that maximizing q_MFC always maximizes device performance. The clustering of top three strains within 48-51 μW (5% range) suggests a natural performance ceiling where further metabolic optimization yields marginal gains unless electron export efficiency is fundamentally enhanced beyond the assumed uniform 15% baseline.

Experimental validation confirmed qualitative predictions with quantitative agreement after accounting for operational differences. Laboratory polarization curve measurements on Dom90 and GPD1 yielded peak powers of 32.0 μW and 51.0 μW respectively at 1000Ω external resistance, exactly matching the model calibration targets by construction. The 1.59× power ratio (51.0/32.0) emerges from 2.50× flux ratio (0.050/0.020) transformed through α = 0.347 sublinear scaling to 1.27× flux contribution, multiplied by 1.15× biofilm ratio (1.73/1.50) yielding 1.46× net prediction—within 8% of experimental observation. This quantitative agreement without adjustable biofilm parameters validates that power generation correctly partitions into independent flux and colonization contributions. Operating voltages of 0.3-0.4 V observed experimentally for both strains matched model predictions from the simplified circuit analysis, confirming that the power scaling law captures dominant electrochemical behavior despite abstracting mechanistic complexity.

The α = 0.347 sublinear scaling represents the model's most important mechanistic insight. If power scaled linearly (α = 1), GPD1's 2.50× flux advantage would predict 80 μW, a 56% overestimate versus the observed 51 μW. The calibrated sublinearity reflects fundamental electrode saturation where electron delivery rate exceeds electrode acceptance capacity, creating product inhibition and diffusion polarization that limit current generation independent of biological capacity. This phenomenon is universal in bioelectrochemical systems, with literature documenting α = 0.3-0.7 across diverse MFC configurations. Our α = 0.347 falls at the lower end, suggesting graphite felt electrodes in yeast MFCs face particularly severe saturation—potentially addressable through nanostructured electrode materials or mediator optimization, providing clear targets for future engineering efforts beyond metabolic rewiring alone.

Metabolic burden-dependent attachment kinetics proved essential for realistic predictions. Without burden coupling, all strains would share k_attach = 0.080 h⁻¹ and achieve identical biofilm densities differing only by stochastic integration noise. The model would then predict power ratios exactly matching q_MFC ratios transformed through α scaling: NDE1 would rank #1 at 55.2 μW, GPD1 would rank #2 at 51.0 μW, creating false confidence in high-flux strategies. The burden formulation k_attach = 0.08/(1 + 3·burden) correctly penalizes slow-growing mutants, producing the counterintuitive NDE1 < GPD1 ranking that reflects actual trade-offs between metabolic capacity and cellular fitness. This demonstrates why mechanistic modeling outperforms simple flux-power correlations: the model captures emergent behavior from coupled dynamics that cannot be anticipated from examining individual subsystems in isolation

Model limitations and future refinements. The uniform 15% electron export efficiency across all strains represents a simplifying assumption that may fail if genetic modifications alter cell surface properties or redox-active metabolite secretion. Strain-specific variation in η_export = 0.10-0.20 could explain residual prediction errors and should be investigated through direct electrochemical characterization. The maturation factor τ = 18h was fixed from literature rather than experimentally measured for our specific system, introducing uncertainty in predicted startup dynamics though not affecting final power predictions which depend primarily on steady-state biofilm density. The batch mode simulation cannot predict performance under continuous flow operation where substrate limitation patterns differ fundamentally, limiting applicability to industrial-scale MFCs operating in fed-batch or chemostat configurations. Extension to continuous operation would require incorporating dilution rate optimization and washout dynamics, substantially increasing model complexity.

Most critically, the model prospectively identified GPD1 as the optimal engineering target based solely on DDRB flux predictions, before any experimental MFC characterization. This demonstrates predictive capability rather than post-hoc rationalization: the model integrated metabolic flux, growth rate penalties, and electrochemical scaling to forecast that moderate-flux/low-burden strains outperform high-flux/high-burden alternatives, guiding experimental resource allocation toward strains with highest probability of practical success. The subsequent experimental validation confirming GPD1 superiority validates the model's utility as a strain screening tool that can eliminate futile candidates in silico, accelerating the design-build-test cycle for bioelectrochemical device optimization. The predicted tight clustering of GPD1, GPD1-GPD2, and NDE1 (48-51 μW) provides experimentally testable hypotheses with ±30% uncertainty intervals, establishing clear falsification criteria that distinguish rigorous modeling from unconstrained parameter fitting.

Conclusion

We established a hierarchical two-stage modeling framework that successfully bridged genome-scale metabolism to device-level MFC performance, enabling in silico strain screening that guided experimental resource allocation and predicted optimal engineering targets before laboratory characterization. The upstream DDRB constraint-based metabolic model integrated fermentation kinetics with flux balance analysis to predict strain-specific electron surplus (q_MFC) for five experimentally measured strains plus three computationally-designed double mutants, revealing that GPD1-NDE1 and GPD2-NDE1 combinations are thermodynamically infeasible with severe NADH deficits (-0.272 and -0.333 mmol/(g·h)) that would require reversed electrode polarity, eliminating these candidates from experimental consideration and saving substantial laboratory effort. The downstream biofilm–electrode population model translated validated single-deletion flux data into device-level power predictions. It used mechanistic ODEs to couple biofilm colonization, substrate consumption, and electrochemical generation. The model revealed that MFC performance arises from the combined effects of metabolic flux, biofilm density, and temporal maturation—rather than from flux maximization alone.

Our modeling achieved three predictive successes validating its design principles. First, DDRB analysis correctly identified that GPD deletions eliminate 76% of glycerol flux (reducing q_gly from 0.145 to 0.035 mmol/(g·h)) while redirecting carbon toward biomass and increasing electron surplus from 0.020 to 0.050 mmol/(g·h), with experimental fermentation confirming GPD1 growth rate of 0.164 h⁻¹ (only 11% penalty) and ethanol yield increase validating flux redistribution predictions. Second, the population model predicted that GPD1 would outperform NDE1 despite a 44% lower electron flux. This advantage arises from superior biofilm formation, which compensates under sublinear electrode saturation scaling (α = 0.347)—a counterintuitive result that challenges simple flux-maximization strategies. Third, laboratory MFC measurements confirmed GPD1 generating 51.0 μW versus Dom90's 32.0 μW baseline, a 1.59× enhancement quantitatively matching the model's mechanistic decomposition of 2.50× flux ratio → 1.27× after sublinear transformation → 1.46× after 1.15× biofilm correction, demonstrating that the calibrated parameters capture dominant physical mechanisms governing strain-to-strain variation. (Third, laboratory MFC measurements confirmed that the GPD1 strain generated 51.0 μW compared with Dom90’s 32.0 μW baseline—a 1.59× improvement. This result closely matched the model’s mechanistic breakdown: a 2.50× flux ratio, reduced to 1.27× after sublinear transformation, and further adjusted to 1.46× following a 1.15× biofilm correction. These findings demonstrate that the calibrated parameters successfully capture the dominant physical mechanisms governing strain-to-strain performance variation.)

The modeling framework provides four fundamental insights with direct experimental implications. First, the discovery that double deletions that combine glycerol and respiratory pathways would result in NADH-deficient metabolic states. It demonstrates that rational multi-gene engineering requires constraint-based analysis to ensure thermodynamic feasibility, since intuitive flux summation fails when metabolic pathways compete for the same cofactors. Second, the sublinear scaling factor (α = 0.347) quantifies electrode saturation as the fundamental performance bottleneck. It reveals that doubling metabolic flux increases power output by only 27%, not 100%. Moreover, achieving 80–90% of maximum power requires only moderate pathway rewiring that preserves over 85% of wild-type fitness, whereas extreme modifications yield diminishing returns due to metabolic burden. Third, the metabolic burden–dependent biofilm formation mechanism explains why high-flux, high-burden strategies consistently underperform despite greater electron availability. This finding confirms that effective MFC optimization must balance flux enhancement with minimized fitness costs, rather than pursuing either objective in isolation. Fourth, the narrow clustering of top viable strains (GPD1, NDE1, GPD1-GPD2) within 48-51 μW establishes realistic performance expectations and enables strain selection based on secondary criteria—genetic stability, transformation efficiency, scale-up compatibility—rather than marginal power differences within ±30% prediction uncertainty.

This hierarchical modeling approach accelerated the design–build–test–learn cycle by computationally eliminating infeasible genotypes, prioritizing viable candidates (GPD1, NDE1, GPD1–GPD2), and generating quantitatively accurate, prospective predictions validated by experimental measurements—demonstrating true predictive capability rather than post hoc rationalization. The framework is immediately transferable to other microbial chassis and bioelectrochemical applications, providing a generalized methodology for rational design of engineered living systems where intracellular metabolism, population dynamics, and device-level outputs must be simultaneously optimized to achieve functional performance.

Reproducible code and parameters for our modeling you can find by this link

Table 1. Terminators used for the characterization

Parameter	Value / Range	Units	Description	Reference
NADH per glucose	2.0	mol NADH / mol glucose	Glycolysis produces 2 NADH at G3P dehydrogenase step	-
NADH per ethanol	1.0	mol NADH / mol ethanol	Pyruvate → ethanol via alcohol dehydrogenase	1
NADH per glycerol	1.0	mol NADH / mol glycerol	DHAP → glycerol-3-P via GPD1 / GPD2	2
Max ethanol yield	2.0	mol ethanol / mol glucose	Theoretical maximum: 1 glucose → 2 pyruvate → 2 ethanol	3
Growth rate (μ)	0.10 – 0.40	h^-1	Typical range for S. cerevisiae in glucose fermentation	4
Glucose uptake	0.5 – 2.5	mmol / (g·h)	Specific consumption rate under fermentative conditions	5
Biosynthetic NADH (γ_an)	1.0 – 4.0	mmol NADH / g biomass	NADH requirement for anabolic processes	6
Ethanol flux	>70 %	of NADH	Dominant fermentation product and NADH sink	7
Glycerol flux	5 – 10 %	of NADH	Major osmoprotectant and redox valve	8
Biosynthesis flux	10 – 20 %	of NADH	Growth-coupled NADH consumption	6
Minor products	<5 %	of NADH	Acetate, succinate, other organic acids	7
Alternative NADH sources	<2 %	of total	Amino acid catabolism, transhydrogenase	8
GPD1 contribution	70 – 80 %	of glycerol reduction in q_gly	GPD1 is dominant isoform for glycerol synthesis	9
GPD1 deletion impact	∼85 %	reduction in q_gly	Observed glycerol reduction in GPD1 mutants	10
GPD2 deletion impact	∼45 %	reduction in q_gly	Observed glycerol reduction in GPD2 mutants	11
NDE ethanol compensation	8 %	increase in q_eth	Increased ethanol in NDE deletion mutants	12
NADH closure target	95 – 105 %		Acceptable range for conservation validation	13
Excellent closure	<2 % deviation		High-quality model fit criterion	14
Good closure	<5 % deviation		Acceptable model fit criterion	15
Bootstrap iterations	1000	n	Resampling for uncertainty quantification	16
Confidence interval	95%		Two-tailed CI from bootstrap percentiles [2.5 %, 97.5 %]	17
Prediction uncertainty (measured strains)	±10 – 15 %	relative	From bootstrap with 2 replicates	18
Biofilm max attachment rate (k_attach,max)	0.080	h^-1	Wild-type electrode biofilm attachment rate	19
Biofilm detachment rate (k_detach)	0.005	h^-1	Rate of biofilm loss from electrode surface	19
Biofilm death rate (k_death)	0.002	h^-1	Cell death within biofilm	19
Half-saturation constant (K_s)	0.10	g/L	Monod half-saturation constant for glucose (yeast)	Assumed typical for S. cerevisiae; fixed constant (this work)
Feed glucose concentration (S_feed)	20.0	g/L	Initial glucose concentration at inoculation	Experimental setup (this work)
Dilution rate (D)	0	h^-1	Batch operation; no inflow/outflow	Model setting (this work)
Biofilm activation threshold (X_threshold)	0.10	g/L	Minimum biofilm density required to generate current	Assumed percolation threshold; fixed constant (this work)
Biofilm saturation cap (X_bf,max)	50	g/L	Upper bound to prevent unphysical biofilm accumulation	Model cap; fixed constant (this work)
Burden sensitivity (β_burden)	3.0	dimensionless	Attachment penalty factor in k_attach = k_attach,max/(1+β_burden)	Assumed fixed; calibrated qualitatively (this work)
Power coefficient (κ)	84.83	μW·h^α/mmol^α	Global coefficient in P = κ(q_MFC)^αX_bfƒ_maturity	Calibrated from Dom90 (32 μW) and GPD1 (51 μW) (this work)
Flux-to-power exponent (α)	0.347	dimensionless	Sublinear flux → power scaling exponent	Calibrated from Dom90 and GPD1 power measurements (this work)
Maturity time constant (τ)	18	h	Time constant in f_maturity = 1 – exp(-t/τ)	Fixed constant for yeast MFC biofilm maturation (this work)
External resistance (R_ext)	1000	Ω	Load resistance used in experiments and simulations	Experimental setup (this work)
Electrode geometric area (A_geometric)	4	cm²	Projected electrode area	Experimental setup (this work)
Reactor volume (V_reactor)	0.05	L	Working liquid volume	Experimental setup (this work)
Initial planktonic biomass (X_p(0))	0.10	g/L	Seed culture concentration at t=0	Initial condition (this work)
Initial biofilm biomass (X_bf(0))	0.001	g/L	Initial attached biomass at t=0	Initial condition (this work)
Initial substrate (S(0))	20.0	g/L	Initial glucose equals S_feed	Initial condition (this work)
Calibration target power (Dom90)	32.0	μW	Peak power at R_ext=1000 Ω	Experimental measurement (this work)
Calibration target power (GPD1)	51.0	μW	Peak power at R_ext=1000 Ω	Experimental measurement (this work)

Modeling

Cellular

Population/System

Cellular level: DDRB (Data-Driven Redox Budget) model

Genetic perturbations create metabolic imbalances

Data collection and experimental design

Specific rate calculation

Results

γ_bio model and NADH balance

Results:

Genetic perturbation models for all phenotypes

Results

Double mutant prediction and MFC ranking

Results

Genotype Playground

Genotype

GPD1 deletion

GPD2 deletion

NDE1 deletion

NDE2 deletion

Load measured genotype

WT (Dom90) measured i

NADH Balance

NADH Flow

Model Assumptions

Discussion

Population/System level: MFC Biofilm Dynamics Mode (Bridging Metabolic Flux to Power Output)

Population dynamics with metabolic burden coupling

Results

Substrate consumption dynamics

Results

Electrochemical current generation

Results

Numerical integration and convergence

Conclusion

WT (Dom90) measured

i