Building a causal in-silico simulator for fermentation with causica DECI and process-order priors

A fermentation batch passes through several stages and records hundreds of variables along the way: measurements over time at each stage, media-ingredient quantities, and the final yield and culture time. The data sat in several systems and many variables were managed by hand, so building an analysable table was itself the first step. Improvement had relied on expert experience and one-off analyses — and because process conditions drive manufacturing cost directly, any condition proposed had to be explainable and reliable enough to try on a production line. Correlation and simple tests are a poor guide here. A variable that shares a cause with yield passes every screen and changes nothing when you set it; a condition with an optimum inside its range looks unimportant to a straight-line test.

In this post we walk through the system we built. It pivots long-format records into a batch-level table, shortlists candidates with Welch t-tests, quadratic F-tests and random-forest SHAP, turns the process order into a constraint matrix that forbids impossible edges, learns the causal graph and nonlinear structural equations together with DECI from Microsoft's causica library, and lets process experts run do-interventions — effect curves and two-point comparisons — in a Streamlit simulator. The system was applied at three global sites and contributed about USD 1.86M in confirmed savings.

Solution overview

The system has two halves. A data and screening path, developed in Jupyter with pandas, turns raw records into a clean batch table and a candidate list. An interactive path, a multi-screen Streamlit app, lets process experts encode what they know, train the causal model and query it. Screens in the app hand results to each other through session keys and files on disk.

The numbered steps in the diagram:

1. Long-format records (batch, time point, tag, value) are pivoted with pandas to one row per batch and joined with media-preparation data; derived features include concentrations, the feed rate per interval and the target, yield divided by culture time.
2. Columns with too many missing values or a single value are dropped, outliers are masked, and the rest is imputed with mean, median, mode, a constant or KNN.
3. Candidates are shortlisted three ways: a Welch t-test between top- and bottom-quartile batches, a single-variable quadratic OLS F-test, and random-forest SHAP importance.
4. Experts drag variables into process stages with streamlit-sortables; edges within a stage or against process order are forbidden, and known causal or non-causal pairs are forced in or out of a constraint matrix.
5. DECI learns a distribution over graphs and neural structural equations jointly, with a sparsity prior and an augmented-Lagrangian acyclicity constraint, under the constraint matrix.
6. The learned model is sampled under do(X = x) to draw effect curves across the observed range and to compare a control and a treatment value, in a Streamlit app with plotly charts and an interactive graph view.

Technology stack

Step 1: Build one row per batch

Process data arrives in long format — one row per batch, time point, tag and value — from systems and from sheets that were maintained by hand. Causal discovery needs one row per batch, because each batch is one draw from the process, so the first job is a pivot. Each (stage, tag, time point) becomes a column, media-preparation quantities are joined on the batch, and a few features are derived because they are what an engineer actually sets or reads: concentrations from amounts and volumes, the feed rate in each interval from cumulative feed, and the target, yield divided by culture time, so a batch that reaches the same yield faster counts as better.

At real scale the result is hundreds of batches × hundreds of variables: wide, short and full of gaps. That shape drives every choice that follows — aggressive screening, strong prior knowledge, and a learner that can use both.

import pandas as pd

def batch_table(records, media, outcome):
    """records: long format (batch, stage, hour, tag, value) -> one row per batch."""
    wide = records.pivot_table(index="batch", columns=["stage", "tag", "hour"],
                               values="value", aggfunc="mean")
    wide.columns = [f"{stage}.{tag}@{hour}h" for stage, tag, hour in wide.columns]

    # Feed rate per interval from a cumulative feed reading.
    feed = records[records.tag == "feed_total"].sort_values(["batch", "hour"])
    g = feed.groupby("batch")
    feed = feed.assign(rate=g.value.diff() / g.hour.diff())
    rates = feed.pivot_table(index="batch", columns="hour", values="rate")
    rates.columns = [f"feed_rate@{h}h" for h in rates.columns]

    table = wide.join(rates).join(media.set_index("batch"))       # media preparation
    out = outcome.set_index("batch")
    table["target"] = out["yield_amount"] / out["culture_hours"]  # yield per time
    return table

Simplified. Stage, tag and column names are generic placeholders.

Step 2: Clean the table and screen candidates

Preprocessing is deliberately conservative. Columns with too many missing values (above a threshold τ) or with a single value carry no information and are dropped; outliers are masked as missing, and the gaps are then filled with mean, median, mode, constant or KNN imputation.

Screening then narrows hundreds of columns to a candidate set. No single test is enough, so three run side by side:

• Welch t-test, top vs bottom quartile. Does the variable differ between the best and the worst batches? Welch's version does not assume equal variances, which rarely hold between groups like these.
• Quadratic OLS F-test. Regress the target on x and x² and test the fit. This catches a condition with an optimum inside its range, which a linear correlation reports as unimportant.
• Random forest + SHAP. Mean absolute SHAP values from a tree ensemble rank variables by their contribution, interactions included.

import numpy as np
import pandas as pd
import shap
import statsmodels.formula.api as smf
from scipy import stats
from sklearn.ensemble import RandomForestRegressor

def screen(df, target="target"):
    """df: cleaned and imputed batch table."""
    X, y = df.drop(columns=target), df[target]
    top, bottom = y >= y.quantile(0.75), y <= y.quantile(0.25)
    rows = []
    for col in X.columns:
        p_welch = stats.ttest_ind(X.loc[top, col], X.loc[bottom, col],
                                  equal_var=False).pvalue            # Welch
        d = pd.DataFrame({"x": X[col], "y": y})
        p_quad = smf.ols("y ~ x + I(x ** 2)", data=d).fit().f_pvalue   # curved effect
        rows.append({"variable": col, "p_welch": p_welch, "p_quad": p_quad})

    rf = RandomForestRegressor(n_estimators=500, random_state=0).fit(X, y)
    shap_values = shap.TreeExplainer(rf).shap_values(X)
    result = pd.DataFrame(rows).set_index("variable")
    result["mean_abs_shap"] = np.abs(shap_values).mean(axis=0)
    return result.sort_values("mean_abs_shap", ascending=False)

Simplified. Significance thresholds and the rule that combines the three screens are not shown.

Why screen at all? Learning a causal graph over hundreds of variables from hundreds of batches is poorly determined. Screening is a filter, not an answer: a variable can pass every test and still have no causal path to yield — which is exactly what the next steps are for.

Step 3: Turn process order into a constraint matrix

The cheapest data in this project was what the process experts already knew. A variable measured in a later stage cannot cause one in an earlier stage, and in this design variables in the same stage are not linked either. Experts arrange the shortlisted variables into stage containers on a drag-and-drop screen built with streamlit-sortables, and can add pairs they know to be causal or non-causal. The target sits alone in a final stage.

C ∈ {NaN = learn, 0 = forbid, 1 = force}^((d+1)×(d+1))         row = cause, column = effect

C[i, j] = 0     if stage(i) ≥ stage(j)          within a stage, or against process order
C[y, j] = 0     for every j                      the target is last and causes nothing
C[i, j] = 1     expert-supplied causal pair
C[i, j] = 0     expert-supplied non-causal pair

At most half of all ordered pairs can point forward in a stage order, so the order constraint alone rules out at least half of the possible directed edges before learning starts — in practice well over half, since every stage holds several variables. Every forbidden edge is one the learner cannot get wrong.

import numpy as np
from streamlit_sortables import sort_items

def stage_editor(stage_names, unassigned):
    """Drag-and-drop: one container per process stage, in process order."""
    boxes = [{"header": "Unassigned", "items": unassigned}]
    boxes += [{"header": s, "items": []} for s in stage_names]
    return sort_items(boxes, multi_containers=True)[1:]       # drop "Unassigned"

def constraint_matrix(stages, target, causal=(), non_causal=()):
    order = [v for box in stages for v in box["items"]] + [target]
    stage = {v: k for k, box in enumerate(stages) for v in box["items"]}
    stage[target] = len(stages)                               # target alone, last
    idx = {v: i for i, v in enumerate(order)}
    C = np.full((len(order), len(order)), np.nan, dtype=np.float32)   # NaN = learn
    for u in order:
        for v in order:
            if stage[u] >= stage[v]:
                C[idx[u], idx[v]] = 0.0                       # same stage or backwards
    for u, v in causal:
        C[idx[u], idx[v]] = 1.0                               # expert: force the edge
    for u, v in non_causal:
        C[idx[u], idx[v]] = 0.0                               # expert: forbid the edge
    return order, C

Simplified. sort_items returns the containers in their new order; the contents of the experts' prior knowledge are not shown.

Step 4: Learn the graph and the equations together with DECI

DECI (deep end-to-end causal inference) treats discovery and estimation as one problem. It keeps a variational distribution over directed graphs and, for each variable, a neural structural equation that maps the variable's parents to its value plus noise. Training maximises an evidence lower bound with a sparsity prior on the graph, and enforces acyclicity with an augmented Lagrangian: the DAG penalty h(G) is zero only for acyclic graphs, and its multipliers are raised between outer steps until it is. The noise can be Gaussian or a more flexible spline.

x_j = f_j( x_pa(j) ; θ ) + ε_j            ε_j ~ Gaussian or spline           neural SEM
G ~ q_φ(G),   G[i,j] fixed where C[i,j] ≠ NaN                              graph distribution

ELBO(θ, φ) = E_q[ log p_θ(X | G) ] − λ · E_q‖G‖₁ + H(q_φ)                  fit + sparsity
h(G)       = tr( e^(G∘G) ) − d                                              zero iff G is acyclic
max  ELBO − α · E_q[h(G)] − (ρ/2) · E_q[h(G)]²                             α, ρ raised between outer steps

We used causica's implementation on PyTorch and Lightning. The constraint matrix from Step 3 is attached to the module before training, so forbidden edges are never sampled and forced ones are always present. After training, the graph is read out as the mode of the learned distribution, and the app shows the ancestors of the target — the variables with a directed path to yield, which are the only ones worth simulating.

import pytorch_lightning as pl
import torch
from causica.datasets.causica_dataset_format import Variable
from causica.distributions import ContinuousNoiseDist
from causica.lightning.data_modules.basic_data_module import BasicDECIDataModule
from causica.lightning.modules.deci_module import DECIModule

def fit_deci(df, C, sparsity, epochs, noise=ContinuousNoiseDist.SPLINE):
    """df columns must be in the same order as the rows and columns of C."""
    data = BasicDECIDataModule(df, variables=[Variable(c, c) for c in df.columns],
                               batch_size=128, normalize=True)
    deci = DECIModule(noise_dist=noise, prior_sparsity_lambda=sparsity)
    deci.constraint_matrix = torch.tensor(C)          # NaN learn, 0 forbid, 1 force
    trainer = pl.Trainer(max_epochs=epochs, accelerator="auto",
                         enable_checkpointing=False)
    trainer.fit(deci, datamodule=data)
    sem = deci.sem_module().mode                      # most likely graph + equations
    return sem, data.normalizer

def ancestors_of(sem, names, target):
    """Variables with a directed path to the target in the read-out graph."""
    A = sem.graph.cpu().numpy().astype(bool)          # A[i, j]: i -> j
    seen, stack = set(), [names.index(target)]
    while stack:
        j = stack.pop()
        for i in A[:, j].nonzero()[0]:
            if i not in seen:
                seen.add(i)
                stack.append(i)
    return [names[i] for i in sorted(seen)]

Simplified, following causica's Lightning interface; hyper-parameters are placeholders and the augmented-Lagrangian schedule is left at its defaults.

DECI results also varied by platform, which is worth knowing before comparing runs made on different machines.

Step 5: Simulate interventions

A learned structural model answers a different question from a regression: not "what yield goes with this value" but "what yield would follow if we set this value". Setting X cuts the arrows into X, fixes its value and propagates through its descendants, while everything upstream keeps its natural variation. Two queries were built on top:

• Effect curve. E[yield | do(X = x)] at the 10th, 20th, …, 90th percentiles of X — the realistic range, never beyond it.
• Two-point test. A control value a and a treatment value b: sample yield under each and report the average treatment effect, its standard error and a t-test.

effect curve     μ(x) = E[ Y | do(X = x) ],      x ∈ { p10, p20, …, p90 of X }
two-point test   ATE = E[ Y | do(X = b) ] − E[ Y | do(X = a) ]
                 SE  = √( s²_a / n + s²_b / n ),      t = ATE / SE           n samples per arm

import numpy as np
import torch
from scipy import stats
from tensordict import TensorDict

def do_samples(sem, normalizer, var, value, target, n=5000):
    """Samples of the target under do(var = value), in original units."""
    iv = TensorDict({var: torch.tensor([float(value)])}, batch_size=torch.Size())
    with torch.no_grad():
        x = sem.do(interventions=normalizer(iv)).sample(torch.Size([n]))
    return normalizer.inv(x)[target].squeeze(-1).numpy()

def effect_curve(sem, normalizer, df, var, target):
    grid = np.percentile(df[var].dropna(), np.arange(10, 100, 10))   # p10 ... p90
    return grid, [do_samples(sem, normalizer, var, v, target).mean() for v in grid]

def two_point(sem, normalizer, var, a, b, target, n=5000):
    ya = do_samples(sem, normalizer, var, a, target, n)               # control
    yb = do_samples(sem, normalizer, var, b, target, n)               # treatment
    ate = yb.mean() - ya.mean()
    se = np.sqrt(ya.var(ddof=1) / n + yb.var(ddof=1) / n)
    t, p = stats.ttest_ind(yb, ya, equal_var=False)
    return {"ate": ate, "se": se, "t": t, "p": p}

Simplified. Interventions are normalised with the data module's normaliser before sampling and mapped back afterwards.

Why percentiles, not a free slider? A structural model is only trustworthy where it has seen data. Restricting interventions to the 10th–90th percentile range keeps every proposal inside conditions the process has actually run.

Step 6: Put the simulator in the experts' hands

The analysis was delivered as a tool, not a report. A multi-screen Streamlit app walks through loading data, choosing variables, arranging stages and known pairs, training, inspecting the graph (rendered with causalnex and pyvis) and running simulations with plotly charts. Screens hand results to each other through session keys and files on disk.

Effect curves are the slow part: every point on a curve is a fresh batch of samples from the model, and an expert usually asks about several variables at once. Curves are therefore computed in parallel processes with ProcessPoolExecutor and cached with Streamlit's data cache, so returning to a screen or a variable does not recompute them.

from concurrent.futures import ProcessPoolExecutor
import plotly.graph_objects as go
import streamlit as st

@st.cache_data(show_spinner="Simulating interventions…")
def effect_curves(model_path: str, variables: tuple, target: str):
    """One effect curve per variable, each computed in its own process."""
    with ProcessPoolExecutor() as pool:
        jobs = {v: pool.submit(curve_worker, model_path, v, target) for v in variables}
        return {v: job.result() for v, job in jobs.items()}

st.header("What if we set…")
model_path = st.session_state.get("model_path")       # written by the training screen
if model_path is None:
    st.info("Train a causal model on the previous screen first.")
    st.stop()

candidates = st.session_state["target_ancestors"]     # variables with a path to yield
chosen = st.multiselect("Variables", candidates, default=candidates[:3])
for var, (grid, mean_y) in effect_curves(model_path, tuple(chosen), "target").items():
    fig = go.Figure(go.Scatter(x=grid, y=mean_y, mode="lines+markers"))
    fig.update_layout(title=f"E[target | do({var} = x)]",
                      xaxis_title=f"{var}, p10 to p90", yaxis_title="target")
    st.plotly_chart(fig, use_container_width=True)

Simplified. curve_worker loads the trained model from disk inside each process and returns one curve.

Hypotheses that survived the simulator were the ones taken to the line, and some of them became standard process conditions.

Results

The system was applied at three global sites, where confirmed savings totalled about USD 1.86M, including USD 773K realized through the in-silico simulator. It produced hypotheses that could be tested on real production lines, and some of the conditions it found were adopted as standard process conditions. It turned one-off analyses into a repeatable improvement system, led to five patent filings and was connected to technology licensing.

A further outcome was an agenda for manufacturing data itself: digitising hand-managed records, improving data management, and building automated pipelines and a process knowledge base, so the next round of analysis starts from better data.

Lessons learned

• Prior knowledge is the cheapest data. Process order alone rules out at least half of the possible edges before any learning, and every forbidden edge is a mistake the learner cannot make.
• Screen with more than one lens. Welch t-tests, quadratic F-tests and SHAP disagree in useful ways; the quadratic test in particular keeps conditions with an optimum inside their range.
• Learn the structure and the equations together. DECI's joint fit yields both the graph and the nonlinear equations needed to simulate it, in one model.
• Ship a tool, not a report. Experts running their own what-if queries is what turned one analysis into a repeatable improvement process.
• Treat the simulator as a hypothesis generator. Its answers hold within the observed range and under the learned structure; the production line is the real test.

Conclusion

Yield improvement in fermentation is a causal question asked of observational data. Combining a careful batch table, multi-lens screening, process order as hard constraints, DECI for structure and equations, and a simulator experts can query produced conditions that were explainable, testable on the line and in some cases adopted as standards.

The same pattern — encode what experts know as constraints, learn the rest from data, and let the people who own the process ask what-if questions — fits other staged processes where many variables move together and only a few are worth changing.

Limitations

• Causal discovery from observational batches depends on the prior and on the variables measured; an unmeasured common cause can still mislead it. Conditions were validated on the line before adoption.
• The simulator estimates average effects within the observed operating range; it does not extrapolate safely beyond it.
• The two-point test runs on samples drawn from the model, so its p-value reflects the model and the number of samples drawn, not a plant trial.
• The live model uses a lighter learner than DECI on generated batches with a planted structure; its accuracy on synthetic data says nothing about accuracy on a real process.

About the demo and confidentiality

Every batch, variable, unit and coefficient in the embedded model is generated, and its planted effects are not real process settings. No product, strain, media composition, tag name, tank, selected-variable list, prior-knowledge content, plant data or analysis result from the real project appears in this post; real scale is described only as hundreds of batches × hundreds of variables. Code is simplified and written for illustration.