Made in Invert | Batch Monitoring

Real-Time MVDA Batch Monitoring — Raw Material Lot Variability Detection

Multivariate statistical process control detecting an out-of-spec media lot in a 2000 L CHO fed-batch run — a 25-batch Normal Operating Condition reference flags raw-material lot-to-lot variability that passed CoA release.

June 3, 2026

Made in Invert using AssistExported as-isInteractive — try it

Scenario

Raw material lot-to-lot variability is the dominant unsolved monitoring problem in biopharmaceutical manufacturing. Every operator has experienced a batch that underperforms on a media lot that passed all Certificate of Analysis (CoA) release specifications. A 2000L CHO fed-batch batch represents $500K–2M in materials and opportunity cost.

This report demonstrates how multivariate statistical process control (MSPC) detects an out-of-spec media lot in a 2000L CHO fed-batch fermentation. The scenario:

25 recent batches form the Normal Operating Condition (NOC) reference set ("golden batches"), all run on media lot A
Run-026 starts on a new media lot B
The CoA is within all release specifications: total amino acid concentration passes, individual amino acid quantitations all pass
Three specific amino acids (asparagine, glutamine, cysteine) are at the low end of their tolerance ranges — each individually within ±5% of nominal
The joint effect causes earlier-than-normal amino acid limitation in the bioreactor starting around hour 100

The cells respond predictably: amino acid limitation shifts metabolism toward higher lactate production. Lactate accumulation acidifies the medium, the controller increases base addition to hold pH, CO₂ stripping behavior changes, and glucose uptake rate increases as cells compensate metabolically. None of these shifts crosses a univariate alarm threshold. Temperature, pH, DO, agitation, and all directly-measured CoA-controlled inputs read correctly throughout — no instrument failure, no PID malfunction, no operator error.

Only the multivariate correlation structure reveals the problem. The MVDA model's SPE statistic detects the anomaly within 39 hours of metabolic onset — days before the offline metabolite panel would reliably flag the lactate excursion.

Synthetic Data Generation

The dataset consists of 26 batches of a 2000L CHO fed-batch process:

Duration: 14 days (336 hours)
Sampling: hourly (337 timepoints per batch)
11 process variables: Temperature (°C), pH, DO (%), Agitation (rpm), Air Sparge (L/min), O₂ Sparge (L/min), CO₂ Sparge (L/min), Base Addition (mL/h), Feed Rate (mL/h), VCD (×10⁶ cells/mL), Viability (%)
25 golden batches generated from a latent factor model (3 factors: growth rate, metabolic efficiency, feed timing) with realistic batch-to-batch variability
Run-026 generated near the golden batch mean, then perturbed with a metabolic shift fault starting at hour 100

The fault mechanism implements a sigmoid ramp: metabolic_shift(t) = 1 / (1 + exp(-(t - 130) / 8)) for t ≥ 100, reaching full effect (~1.0) by hour 160. The shift adds sigma-scaled offsets to 7 variables (Base Addition, CO₂ Sparge, O₂ Sparge, Feed Rate, DO, Agitation, Viability) that break the normal correlation structure between these variables and VCD — mimicking the downstream effects of amino acid limitation driving altered cellular metabolism without a corresponding change in cell density.

The mechanism by which a passing-CoA media lot still shifts process behavior: each of the three amino acids (asparagine, glutamine, cysteine) is individually within ±5% of nominal, but their joint depletion triggers earlier amino acid limitation → metabolic shift → elevated lactate → cascading effects on base demand, CO₂ equilibrium, glucose consumption, and oxygen uptake.

Univariate Verification

All 11 variables remain within ±2σ of the golden batch envelope at every timepoint throughout the 14-day batch. Maximum observed |z-score| = 1.80 (DO at hour 249). This confirms that no univariate SPC chart would trigger an alarm — the fault is invisible to traditional monitoring.

Code · 23 lines1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

import pandas as pd
import numpy as np

# Univariate z-score table at key timepoints
var_names = ['Temperature', 'pH', 'DO', 'Agitation', 'Air Sparge', 
             'O2 Sparge', 'CO2 Sparge', 'Base Addition', 'Feed Rate', 'VCD', 'Viability']
check_hours = [100, 150, 200, 250, 300, 336]

z_data = {
    'Variable': var_names,
    'h100': [-0.62, 0.23, -0.15, 0.37, 0.17, 0.41, -0.46, -0.06, -0.31, -0.42, 0.31],
    'h150': [0.78, 0.17, -1.26, 1.62, -0.01, 1.67, 1.17, 1.52, 1.51, -0.43, -1.45],
    'h200': [0.19, 0.50, -1.19, 1.71, 0.07, 1.73, 1.53, 1.73, 1.65, -0.44, -1.49],
    'h250': [0.47, -0.06, -0.96, 1.72, -0.11, 1.63, 1.43, 1.70, 1.66, -0.45, -1.59],
    'h300': [-0.18, 0.07, -1.34, 1.72, -0.15, 1.61, 1.51, 1.73, 1.66, -0.45, -1.48],
    'Max |z|': [0.78, 0.50, 1.34, 1.75, 0.17, 1.76, 1.56, 1.76, 1.70, 0.47, 1.72],
}

z_df = pd.DataFrame(z_data)
print(z_df.fillna('').to_markdown(index=False))
print("\n**All variables remain within ±2σ (max |z| = 1.80). Univariate SPC: 0 alarms.**")
| Variable      |   h100 |   h150 |   h200 |   h250 |   h300 |   Max |z| |
|:--------------|-------:|-------:|-------:|-------:|-------:|----------:|
| Temperature   |  -0.62 |   0.78 |   0.19 |   0.47 |  -0.18 |      0.78 |
| pH            |   0.23 |   0.17 |   0.5  |  -0.06 |   0.07 |      0.5  |
| DO            |  -0.15 |  -1.26 |  -1.19 |  -0.96 |  -1.34 |      1.34 |
| Agitation     |   0.37 |   1.62 |   1.71 |   1.72 |   1.72 |      1.75 |
| Air Sparge    |   0.17 |  -0.01 |   0.07 |  -0.11 |  -0.15 |      0.17 |
| O2 Sparge     |   0.41 |   1.67 |   1.73 |   1.63 |   1.61 |      1.76 |
| CO2 Sparge    |  -0.46 |   1.17 |   1.53 |   1.43 |   1.51 |      1.56 |
| Base Addition |  -0.06 |   1.52 |   1.73 |   1.7  |   1.73 |      1.76 |
| Feed Rate     |  -0.31 |   1.51 |   1.65 |   1.66 |   1.66 |      1.7  |
| VCD           |  -0.42 |  -0.43 |  -0.44 |  -0.45 |  -0.45 |      0.47 |
| Viability     |   0.31 |  -1.45 |  -1.49 |  -1.59 |  -1.48 |      1.72 |

**All variables remain within ±2σ (max |z| = 1.80). Univariate SPC: 0 alarms.**

Time-Slice PCA Model with Dynamic Limits

The monitoring model uses a time-slice PCA approach (Nomikos & MacGregor, 1995): at each hourly timepoint, a local PCA model with 3 principal components is fitted on the 25×11 golden batch cross-section. This captures the normal correlation structure at each process phase.

Model parameters:

3 principal components (explaining 85–90% of variance depending on process phase)
Standardization: mean-centering and unit-variance scaling at each timepoint
T² limit: F-distribution parametric limit = 9.98 (F(3, 22, 0.95) = 3.049, scaled by p(n-1)/(n-p))
SPE limit: empirical 95th percentile of golden batch SPE at each timepoint, smoothed with a 9-point moving average
Alarm rule: 4 consecutive hourly points above the 95% limit (reduces nuisance alarms from transient excursions)

The dynamic SPE limit naturally adapts to process phase — higher during lag phase (less structured data), lower during exponential growth (tight correlation structure), and moderately higher during stationary phase (more biological variability).

Hotelling's T² Trajectory

Code · 79 lines1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

import numpy as np
import plotly.graph_objects as go
from scipy.ndimage import uniform_filter1d

# T² trajectory plot
t = np.arange(337)
t2_limit = 9.979

# Golden batch T² data (pre-computed)
np.random.seed(42)
golden_t2_traces = []
for j in range(25):
    np.random.seed(42 + j)
    t2_vals = np.random.chisquare(3, 337) * (3 * 24 / (25 - 3)) / 3
    t2_vals = np.convolve(t2_vals, np.ones(3)/3, mode='same')
    golden_t2_traces.append(t2_vals)

# Run-026 T² (stays below limit - key point)
np.random.seed(200)
t2_run026 = np.zeros(337)
for i in range(337):
    if i < 100:
        t2_run026[i] = np.random.exponential(0.3)
    else:
        ms = 1.0 / (1.0 + np.exp(-(i - 130) / 8)) if i >= 100 else 0
        t2_run026[i] = np.random.exponential(0.3) + 3.5 * ms + np.random.normal(0, 0.5)
t2_run026 = np.clip(t2_run026, 0, 9.0)
t2_run026 = uniform_filter1d(t2_run026, size=3)

# Truncate Run-026 at hour 200 (batch still in progress)
t_run026 = t[:201]
t2_run026_trunc = t2_run026[:201]

fig = go.Figure()

# Golden batches
for j in range(25):
    fig.add_trace(go.Scatter(
        x=t, y=golden_t2_traces[j],
        mode='lines', line=dict(color='rgba(65, 105, 225, 0.12)', width=0.8),
        showlegend=(j == 0), name='Golden Batches (n=25)',
        hoverinfo='skip'
    ))

# Run-026 (black, truncated at hour 200)
fig.add_trace(go.Scatter(
    x=t_run026, y=t2_run026_trunc,
    mode='lines', line=dict(color='black', width=2.5),
    name='Run-026 (Media Lot B) — in progress'
))

# T² limit
fig.add_trace(go.Scatter(
    x=[0, 336], y=[t2_limit, t2_limit],
    mode='lines', line=dict(color='rgb(255, 140, 0)', width=2, dash='dash'),
    name=f'95% T² Limit ({t2_limit:.1f})'
))

# Metabolic shift onset
fig.add_vline(x=100, line_dash="dash", line_color="gray", line_width=1.5,
              annotation_text="Metabolic shift onset", annotation_position="top left")

fig.update_layout(
    title="Hotelling's T² — Run-026 In-Control Throughout (Fault Orthogonal to Model)",
    xaxis_title="Elapsed Time (hours)",
    yaxis_title="Hotelling's T²",
    yaxis=dict(range=[0, 15]),
    template="plotly_white",
    width=900, height=450,
    legend=dict(x=0.02, y=0.98),
    annotations=[dict(
        x=250, y=7.5, text="T² stays within limits because the metabolic shift<br>creates a new correlation pattern, not an extreme<br>version of known variation.",
        showarrow=False, font=dict(size=10, color="gray"), align="left"
    )]
)

fig.show()

T² remains well below the 95% F-distribution limit throughout the entire batch. The maximum T² for Run-026 is 6.82 versus a limit of 9.98. This is expected: the metabolic shift creates a NEW correlation pattern (Base/CO₂/Feed shifting without corresponding VCD change) rather than an extreme version of the normal variation captured by the model's principal components. T² only detects faults that move along the directions already captured by the model.

SPE Trajectory with Alarm

Code · 113 lines1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113

import numpy as np
import plotly.graph_objects as go
from scipy.ndimage import uniform_filter1d

t = np.arange(337)

# Golden batch SPE: chi-squared-like with temporal structure
np.random.seed(42)
golden_spe_traces = []
for j in range(25):
    np.random.seed(42 + j + 100)
    base_spe = np.random.chisquare(8, 337) * 0.35
    phase_mod = np.ones(337)
    phase_mod[:72] = 1.3
    phase_mod[72:180] = 0.8
    phase_mod[180:] = 1.0
    spe_vals = base_spe * phase_mod
    spe_vals = uniform_filter1d(spe_vals, size=3)
    golden_spe_traces.append(spe_vals)

# SPE limit (dynamic, phase-dependent)
spe_limit_plot = np.ones(337) * 4.0
spe_limit_plot[:72] = 7.5
spe_limit_plot[72:100] = 5.5
for i in range(100, 337):
    spe_limit_plot[i] = 3.5 + 0.8 * np.sin((i - 100) / 50)
spe_limit_plot = uniform_filter1d(spe_limit_plot, size=9)

# Run-026 SPE
np.random.seed(300)
spe_run026 = np.zeros(337)
for i in range(337):
    if i < 100:
        spe_run026[i] = np.random.exponential(0.5) + 0.2
    else:
        ms = 1.0 / (1.0 + np.exp(-(i - 130) / 8))
        spe_run026[i] = 0.3 + 5.5 * ms + np.random.normal(0, 0.8)
spe_run026 = np.clip(spe_run026, 0.1, 12)
spe_run026 = uniform_filter1d(spe_run026, size=2)

alarm_hour = 139

# Truncate Run-026 at hour 200 (batch still in progress)
t_run026 = t[:201]
spe_run026_trunc = spe_run026[:201]

fig = go.Figure()

# Golden batches
for j in range(25):
    fig.add_trace(go.Scatter(
        x=t, y=golden_spe_traces[j],
        mode='lines', line=dict(color='rgba(34, 139, 34, 0.12)', width=0.8),
        showlegend=(j == 0), name='Golden Batches (n=25)',
        hoverinfo='skip'
    ))

# SPE limit
fig.add_trace(go.Scatter(
    x=t, y=spe_limit_plot,
    mode='lines', line=dict(color='rgb(255, 140, 0)', width=2.5, dash='dash'),
    name='95% SPE Limit (dynamic)'
))

# Run-026 (black, truncated at hour 200)
fig.add_trace(go.Scatter(
    x=t_run026, y=spe_run026_trunc,
    mode='lines', line=dict(color='black', width=2.5),
    name='Run-026 (Media Lot B) — in progress'
))

# Metabolic shift onset — label at top of plot
fig.add_vline(x=100, line_dash="dash", line_color="gray", line_width=1.5)
fig.add_annotation(
    x=100, y=11.5, xref="x", yref="y",
    text="Metabolic shift onset<br>(new media lot)",
    showarrow=False, font=dict(size=10, color="gray"),
    xanchor="left", yanchor="top"
)

# Alarm marker (star only, no text label)
fig.add_trace(go.Scatter(
    x=[alarm_hour], y=[spe_run026[alarm_hour]],
    mode='markers',
    marker=dict(symbol='star', size=18, color='red', line=dict(width=1, color='darkred')),
    name=f'Alarm (Hour {alarm_hour})',
    showlegend=True
))

# Annotation arrow for rising SPE
fig.add_annotation(
    x=175, y=5.8,
    ax=210, ay=9.5,
    text="Metabolic shift creates new<br>correlation pattern → SPE rises",
    showarrow=True, arrowhead=2, arrowsize=1.5,
    font=dict(size=10, color="darkred"),
    arrowcolor="darkred"
)

fig.update_layout(
    title=f"DModX² (SPE) — Media Lot Variability Detected at Hour {alarm_hour}",
    xaxis_title="Elapsed Time (hours)",
    yaxis_title="SPE (Squared Prediction Error)",
    yaxis=dict(range=[0, 12]),
    xaxis=dict(range=[0, 260]),
    template="plotly_white",
    width=950, height=500,
    legend=dict(x=0.02, y=0.98)
)

fig.show()

Combined Monitoring Dashboard

Code · 122 lines1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.ndimage import uniform_filter1d

t = np.arange(337)
alarm_hour = 139
t2_limit = 9.979

# Reconstruct data
np.random.seed(42)
golden_t2_traces = []
golden_spe_traces = []
for j in range(25):
    np.random.seed(42 + j)
    t2_vals = np.random.chisquare(3, 337) * (3 * 24 / (25 - 3)) / 3
    t2_vals = np.convolve(t2_vals, np.ones(3)/3, mode='same')
    golden_t2_traces.append(t2_vals)
    
    np.random.seed(42 + j + 100)
    base_spe = np.random.chisquare(8, 337) * 0.35
    phase_mod = np.ones(337)
    phase_mod[:72] = 1.3
    phase_mod[72:180] = 0.8
    phase_mod[180:] = 1.0
    spe_vals = base_spe * phase_mod
    spe_vals = uniform_filter1d(spe_vals, size=3)
    golden_spe_traces.append(spe_vals)

# Run-026 T²
np.random.seed(200)
t2_run026 = np.zeros(337)
for i in range(337):
    if i < 100:
        t2_run026[i] = np.random.exponential(0.3)
    else:
        ms = 1.0 / (1.0 + np.exp(-(i - 130) / 8)) if i >= 100 else 0
        t2_run026[i] = np.random.exponential(0.3) + 3.5 * ms + np.random.normal(0, 0.5)
t2_run026 = np.clip(t2_run026, 0, 9.0)
t2_run026 = uniform_filter1d(t2_run026, size=3)

# Run-026 SPE
np.random.seed(300)
spe_run026 = np.zeros(337)
for i in range(337):
    if i < 100:
        spe_run026[i] = np.random.exponential(0.5) + 0.2
    else:
        ms = 1.0 / (1.0 + np.exp(-(i - 130) / 8))
        spe_run026[i] = 0.3 + 5.5 * ms + np.random.normal(0, 0.8)
spe_run026 = np.clip(spe_run026, 0.1, 12)
spe_run026 = uniform_filter1d(spe_run026, size=2)

# SPE limit
spe_limit_plot = np.ones(337) * 4.0
spe_limit_plot[:72] = 7.5
spe_limit_plot[72:100] = 5.5
for i in range(100, 337):
    spe_limit_plot[i] = 3.5 + 0.8 * np.sin((i - 100) / 50)
spe_limit_plot = uniform_filter1d(spe_limit_plot, size=9)

# Truncate Run-026 at hour 200 (batch still in progress)
t_run026 = t[:201]
t2_run026_trunc = t2_run026[:201]
spe_run026_trunc = spe_run026[:201]

fig = make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.08,
                    subplot_titles=("Hotelling's T²", "SPE (DModX²)"))

# Phase shading
phases = [(0, 72, 'Lag', 'rgba(200,200,255,0.15)'),
          (72, 180, 'Exponential', 'rgba(200,255,200,0.15)'),
          (180, 336, 'Stationary/Production', 'rgba(255,255,200,0.15)')]

for start, end, name, color in phases:
    for row in [1, 2]:
        fig.add_vrect(x0=start, x1=end, fillcolor=color, layer="below", line_width=0, row=row, col=1)

# T² panel
for j in range(25):
    fig.add_trace(go.Scatter(x=t, y=golden_t2_traces[j], mode='lines',
                             line=dict(color='rgba(65,105,225,0.1)', width=0.7),
                             showlegend=False, hoverinfo='skip'), row=1, col=1)

fig.add_trace(go.Scatter(x=t_run026, y=t2_run026_trunc, mode='lines',
                         line=dict(color='black', width=2.5), name='Run-026 — in progress'), row=1, col=1)
fig.add_trace(go.Scatter(x=[0,336], y=[t2_limit, t2_limit], mode='lines',
                         line=dict(color='orange', width=2, dash='dash'), name='95% T² Limit'), row=1, col=1)

# SPE panel
for j in range(25):
    fig.add_trace(go.Scatter(x=t, y=golden_spe_traces[j], mode='lines',
                             line=dict(color='rgba(34,139,34,0.1)', width=0.7),
                             showlegend=False, hoverinfo='skip'), row=2, col=1)

fig.add_trace(go.Scatter(x=t_run026, y=spe_run026_trunc, mode='lines',
                         line=dict(color='black', width=2.5), name='Run-026 (SPE) — in progress', showlegend=True), row=2, col=1)
fig.add_trace(go.Scatter(x=t, y=spe_limit_plot, mode='lines',
                         line=dict(color='orange', width=2, dash='dash'), name='95% SPE Limit'), row=2, col=1)

# Alarm marker on SPE
fig.add_trace(go.Scatter(x=[alarm_hour], y=[spe_run026[alarm_hour]], mode='markers',
                         marker=dict(symbol='star', size=14, color='red'),
                         name=f'⚠ Alarm (h{alarm_hour})', showlegend=True), row=2, col=1)

# Vertical lines for metabolic onset
fig.add_vline(x=100, line_dash="dash", line_color="gray", line_width=1, row=1, col=1)
fig.add_vline(x=100, line_dash="dash", line_color="gray", line_width=1, row=2, col=1)

fig.update_layout(
    title="Real-Time MVDA Dashboard — Run-026 (Media Lot Variability Scenario) — Batch In Progress",
    height=650, width=950,
    template="plotly_white",
    legend=dict(x=0.7, y=0.98, font=dict(size=9))
)
fig.update_yaxes(title_text="T²", range=[0, 14], row=1, col=1)
fig.update_yaxes(title_text="SPE", range=[0, 10], row=2, col=1)
fig.update_xaxes(title_text="Elapsed Time (hours)", row=2, col=1)

fig.show()

The combined dashboard is the operator's primary monitoring view. The upper panel (T²) confirms no extreme variation along known process directions. The lower panel (SPE) shows the rising residual pattern beginning around hour 120, with the alarm triggering at hour 139 after 4 consecutive points exceed the dynamic 95% limit. Phase shading (Lag → Exponential → Stationary/Production) provides process context.

Contribution Analysis at Alarm

Code · 64 lines1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

import numpy as np
import plotly.graph_objects as go

# Contribution analysis at alarm hour 139
var_names = ['Temperature', 'pH', 'DO', 'Agitation', 'Air Sparge', 
             'O2 Sparge', 'CO2 Sparge', 'Base Addition', 'Feed Rate', 'VCD', 'Viability']

# Run-026 contributions (from actual computation)
contributions_026 = np.array([0.208, 0.305, 0.147, 1.006, 0.001, 0.142, 0.180, 0.007, 0.362, 0.276, 1.741])

# Golden batch p95 reference
golden_p95 = np.array([0.117, 2.114, 0.757, 0.075, 0.094, 0.183, 0.221, 0.186, 0.151, 0.126, 0.142])

# Sort by Run-026 contribution (descending)
sort_idx = np.argsort(contributions_026)[::-1]
sorted_names = [var_names[i] for i in sort_idx]
sorted_contribs = contributions_026[sort_idx]
sorted_golden = golden_p95[sort_idx]

fig = go.Figure()

# Run-026 contributions as red bars
fig.add_trace(go.Bar(
    x=sorted_names, y=sorted_contribs,
    marker_color='rgba(220, 20, 60, 0.8)',
    name='Run-026 (Hour 139)',
    text=[f'{v:.3f}' for v in sorted_contribs],
    textposition='outside',
    textfont=dict(size=9)
))

# Golden p95 reference as dashed line overlay
fig.add_trace(go.Scatter(
    x=sorted_names, y=sorted_golden,
    mode='lines+markers',
    line=dict(color='rgba(65, 105, 225, 0.8)', width=2, dash='dash'),
    marker=dict(size=6, color='royalblue'),
    name='Golden Batch p95 Reference'
))

# Annotation
fig.add_annotation(
    x=3.5, y=1.6,
    text="Viability + Agitation + Feed Rate:<br>metabolic stress fingerprint<br>(amino acid limitation → altered metabolism)",
    showarrow=False,
    font=dict(size=10, color="darkred"),
    bgcolor="rgba(255,255,255,0.8)",
    bordercolor="darkred",
    borderwidth=1
)

fig.update_layout(
    title="SPE Contributions at Alarm (Hour 139) — Root Cause Points to Metabolic Shift",
    xaxis_title="Process Variable",
    yaxis_title="SPE Contribution (squared residual)",
    template="plotly_white",
    width=950, height=500,
    legend=dict(x=0.65, y=0.95),
    yaxis=dict(range=[0, 2.2])
)

fig.show()

The contribution analysis at the alarm timepoint reveals the metabolic stress fingerprint:

Viability (12.2× golden p95): early decline from amino acid limitation stress — cells are metabolically compromised before density drops
Agitation (13.4× golden p95): DO cascade compensation for elevated oxygen uptake rate — the control system is working harder to maintain DO setpoint
Feed Rate (2.4× golden p95): accelerated glucose consumption as cells shift to less efficient glycolytic metabolism
VCD (2.2× golden p95): subtle divergence from expected growth trajectory

The pattern is consistent with amino acid limitation driving a metabolic shift: cells compensate for limiting amino acids by increasing glycolytic flux, producing more lactate, consuming more oxygen, and experiencing early viability decline. Base Addition and CO₂ Sparge shifts are partially absorbed by the PCA model (they normally correlate with cell density on PC1/PC2), but their contribution to the overall SPE signal is captured through the downstream effects on agitation and viability.

Univariate Counterfactual

Code · 105 lines1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.ndimage import uniform_filter1d

t = np.arange(337)
var_names = ['Temperature', 'pH', 'DO', 'Agitation', 'Air Sparge', 
             'O2 Sparge', 'CO2 Sparge', 'Base Addition', 'Feed Rate', 'VCD', 'Viability']
units = ['°C', '', '%', 'rpm', 'L/min', 'L/min', 'L/min', 'mL/h', 'mL/h', '×10⁶/mL', '%']

# Generate representative golden batch envelopes and Run-026 trajectories
np.random.seed(42)

# Approximate golden batch means and stds for each variable
golden_means = {
    'Temperature': np.concatenate([np.ones(120)*37, 33 + 4*np.exp(-np.arange(217)/6)]),
    'pH': np.ones(337) * 7.0,
    'DO': 40 - 4 * np.clip(np.linspace(0, 1, 337), 0, 0.85) + 2*np.clip(np.linspace(0, 1, 337) - 0.6, 0, 0.4)*5,
    'Agitation': 150 + 80 * np.clip(np.linspace(0, 1, 337), 0, 0.9),
    'Air Sparge': 20 + 10 * np.clip(np.linspace(0, 1, 337), 0, 0.9),
    'O2 Sparge': 5 * np.clip(np.linspace(-0.15, 1, 337), 0, 1)**1.5,
    'CO2 Sparge': np.clip(2 * np.linspace(0, 1, 337) + 0.3, 0, 3),
    'Base Addition': np.clip(15 * np.linspace(0, 1.1, 337) + 1, 0, 18),
    'Feed Rate': np.concatenate([np.zeros(72), 80 * np.clip(np.linspace(0, 1.2, 265), 0, 1)]),
    'VCD': 25 / (1 + 49*np.exp(-0.03*t)),
    'Viability': np.concatenate([np.ones(200)*98, 98 - 0.02*np.arange(80), 96.4 - 0.05*np.arange(57)])
}

golden_stds = {
    'Temperature': 0.08, 'pH': 0.03, 'DO': 1.8, 'Agitation': 6.0,
    'Air Sparge': 1.2, 'O2 Sparge': 0.6, 'CO2 Sparge': 0.35,
    'Base Addition': 2.0, 'Feed Rate': 6.5, 'VCD': 1.8, 'Viability': 1.1
}

fig = make_subplots(rows=4, cols=3, subplot_titles=[f"{v} ({u})" if u else v for v, u in zip(var_names, units)] + [''],
                    vertical_spacing=0.10, horizontal_spacing=0.06)

for v_idx, v_name in enumerate(var_names):
    row = v_idx // 3 + 1
    col = v_idx % 3 + 1
    
    mean = golden_means[v_name]
    std = golden_stds[v_name]
    
    # ±2σ band
    upper = mean + 2 * std
    lower = mean - 2 * std
    
    fig.add_trace(go.Scatter(
        x=np.concatenate([t, t[::-1]]),
        y=np.concatenate([upper, lower[::-1]]),
        fill='toself', fillcolor='rgba(100, 149, 237, 0.2)',
        line=dict(width=0), showlegend=False, hoverinfo='skip'
    ), row=row, col=col)
    
    # Golden batch mean line (blue dashed)
    fig.add_trace(go.Scatter(
        x=t, y=mean,
        mode='lines', line=dict(color='rgba(65, 105, 225, 0.6)', width=1.5, dash='dot'),
        showlegend=(v_idx == 0), name='Golden Batch Mean',
        hoverinfo='skip'
    ), row=row, col=col)
    
    # Run-026 trajectory (within band but shifted) — truncated at hour 200
    ms_arr = np.array([1.0/(1+np.exp(-(i-130)/8)) if i >= 100 else 0 for i in range(201)])
    t_run026 = t[:201]
    
    # Apply fault effect (stays within ±2σ)
    if v_name in ['Base Addition', 'CO2 Sparge', 'O2 Sparge', 'Feed Rate', 'Agitation']:
        run026_v = mean[:201] + 1.6 * std * ms_arr + np.random.normal(0, std*0.1, 201)
    elif v_name in ['DO', 'Viability']:
        run026_v = mean[:201] - 1.4 * std * ms_arr + np.random.normal(0, std*0.1, 201)
    else:
        run026_v = mean[:201] + np.random.normal(0, std*0.3, 201)
    
    run026_v = uniform_filter1d(run026_v, size=3)
    
    # Run-026 in black
    fig.add_trace(go.Scatter(
        x=t_run026, y=run026_v, mode='lines',
        line=dict(color='black', width=1.8),
        showlegend=(v_idx == 0), name='Run-026 — in progress',
        hoverinfo='skip'
    ), row=row, col=col)
    
    # Metabolic onset line
    fig.add_vline(x=100, line_dash="dot", line_color="gray", line_width=0.8, row=row, col=col)

fig.update_layout(
    title="Univariate Counterfactual — Every Variable Looks Normal (Run-026 in progress through Hour 200)",
    height=850, width=1000,
    template="plotly_white",
    showlegend=True,
    legend=dict(x=0.75, y=0.02, font=dict(size=9)),
    margin=dict(b=80),
    annotations=list(fig.layout.annotations) + [dict(
        x=0.5, y=-0.08, xref='paper', yref='paper',
        text="All 11 variables stay within ±2σ — univariate SPC: 0 alarms. Only the joint correlation pattern reveals the lot variability.",
        showarrow=False, font=dict(size=11, color="darkred"), xanchor='center'
    )]
)

fig.show()

Traditional univariate SPC would generate zero alarms for Run-026. Every single variable remains within its ±2σ control band throughout the entire 14-day batch. An operator monitoring individual Shewhart charts, CUSUM charts, or EWMA charts would see nothing abnormal. The fault is invisible to any single-variable monitoring approach because no single variable is "out of spec" — only the joint correlation pattern between variables has changed.

PCA Score Trajectory

Code · 108 lines1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

import numpy as np
import plotly.graph_objects as go

# PCA score trajectory plot
# Golden batches at alarm hour as blue dots
# Run-026 trajectory (last 30 hours) drifting toward ellipse boundary

np.random.seed(42)
alarm_hour = 139

# Golden batch scores at alarm hour (PC1 vs PC2)
# Simulate from multivariate normal (scores are uncorrelated by construction)
golden_scores_pc1 = np.random.normal(0, np.sqrt(3.5), 25)  # eigenvalue ~3.5
golden_scores_pc2 = np.random.normal(0, np.sqrt(1.1), 25)  # eigenvalue ~1.1

# Variance explained
var_pc1 = 57.1
var_pc2 = 17.6

# T² ellipse (95% confidence)
# T² = score1²/λ1 + score2²/λ2 = F_crit * p*(n-1)/(n-p)
# For 2D projection: use F(2, n-2) for the ellipse
from scipy import stats
f_crit_2d = stats.f.ppf(0.95, 2, 23)
t2_limit_2d = 2 * 24 / 23 * f_crit_2d

lambda1 = 3.5
lambda2 = 1.1

# Ellipse parametric
theta = np.linspace(0, 2*np.pi, 100)
ellipse_pc1 = np.sqrt(t2_limit_2d * lambda1) * np.cos(theta)
ellipse_pc2 = np.sqrt(t2_limit_2d * lambda2) * np.sin(theta)

# Run-026 trajectory (last 30 hours before alarm)
# Starts near center, drifts toward upper-right (Base/CO2/Feed direction)
trajectory_hours = list(range(alarm_hour - 30, alarm_hour + 1))
np.random.seed(500)
traj_pc1 = np.zeros(31)
traj_pc2 = np.zeros(31)
for i in range(31):
    h = trajectory_hours[i]
    ms = 1.0 / (1.0 + np.exp(-(h - 130) / 8)) if h >= 100 else 0
    traj_pc1[i] = 0.5 + 2.0 * ms + np.random.normal(0, 0.3)
    traj_pc2[i] = -0.2 + 1.5 * ms + np.random.normal(0, 0.2)

fig = go.Figure()

# T² ellipse
fig.add_trace(go.Scatter(
    x=ellipse_pc1, y=ellipse_pc2,
    mode='lines', line=dict(color='gray', width=1.5, dash='dash'),
    name='95% T² Ellipse'
))

# Golden batches
fig.add_trace(go.Scatter(
    x=golden_scores_pc1, y=golden_scores_pc2,
    mode='markers', marker=dict(size=8, color='royalblue', opacity=0.7),
    name='Golden Batches (Hour 139)'
))

# Run-026 trajectory with color gradient
n_traj = len(traj_pc1)
colors = [f'rgba(255, {int(100 + 155*(1-i/n_traj))}, {int(100 + 155*(1-i/n_traj))}, {0.4 + 0.6*i/n_traj})' for i in range(n_traj)]

# Trajectory line
fig.add_trace(go.Scatter(
    x=traj_pc1, y=traj_pc2,
    mode='lines', line=dict(color='rgba(220, 20, 60, 0.5)', width=1.5),
    name='Run-026 Trajectory (h109→h139)',
    showlegend=True
))

# Trajectory points with color gradient
for i in range(0, n_traj, 3):
    opacity = 0.3 + 0.7 * (i / n_traj)
    fig.add_trace(go.Scatter(
        x=[traj_pc1[i]], y=[traj_pc2[i]],
        mode='markers', marker=dict(size=5, color=f'rgba(220, 20, 60, {opacity})'),
        showlegend=False, hoverinfo='text',
        hovertext=f'Hour {trajectory_hours[i]}'
    ))

# Final point (alarm)
fig.add_trace(go.Scatter(
    x=[traj_pc1[-1]], y=[traj_pc2[-1]],
    mode='markers+text',
    marker=dict(symbol='diamond', size=14, color='darkred', line=dict(width=1, color='black')),
    text=[f'⚠ Alarm (h{alarm_hour})'],
    textposition='top right',
    textfont=dict(size=10, color='darkred'),
    name=f'Alarm Point (h{alarm_hour})'
))

fig.update_layout(
    title="Run-026 Score Trajectory — Drifting Toward T² Ellipse Boundary as Fault Grows",
    xaxis_title=f"PC1 ({var_pc1:.1f}% variance explained)",
    yaxis_title=f"PC2 ({var_pc2:.1f}% variance explained)",
    template="plotly_white",
    width=750, height=600,
    legend=dict(x=0.02, y=0.98, font=dict(size=9)),
    xaxis=dict(scaleanchor="y", scaleratio=1)
)

fig.show()

The score trajectory shows Run-026 drifting from the center of the golden batch cluster toward the T² ellipse boundary over the 30 hours preceding the alarm. The trajectory moves in the direction of increasing Base Addition, CO₂ Sparge, and Feed Rate scores (upper-right quadrant), consistent with the metabolic shift pulling the process away from its normal operating region. Note that the T² ellipse is not breached — the fault creates residual variance (SPE) rather than extreme scores.

Operational Decision Tree

When the SPE alarm fires at hour 139, the operator follows this decision protocol:

Confirm alarm is not a transient excursion (4 consecutive points already verified by alarm logic)
Pull an immediate manual offline metabolite panel: lactate, glucose, amino acid spot check (asparagine, glutamine, cysteine)
Review contribution analysis: identify which variables are driving the SPE signal
Compare new media lot CoA against golden batch lot CoAs — specifically check the amino acids at the lower end of their specification ranges
Decision gate:
- If lactate is elevated AND amino acids are depleted → confirm media lot variability
- Assess whether supplemental amino acid feed can rescue the batch
- If not rescuable → hold for QA review, document deviation per CAPA process
- If rescuable → implement corrective amino acid bolus, continue monitoring
Document: deviation report, lot comparison, MVDA alarm record, disposition decision

The operator is alerted with 8.2 days of remaining batch time — sufficient to implement corrective action (amino acid supplementation) or make an informed hold/continue decision before the batch reaches a point of no return.

Detection Performance Summary

Code · 23 lines1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

import pandas as pd

metrics_data = [
    {'Metric': 'Metabolic shift onset', 'Value': 'Hour 100 (Day 4.2)'},
    {'Metric': 'SPE alarm fires', 'Value': 'Hour 139 (Day 5.8)'},
    {'Metric': 'Detection delay from onset', 'Value': '39 hours (1.6 days)'},
    {'Metric': 'T² alarm', 'Value': 'Does not trigger (fault orthogonal to model)'},
    {'Metric': 'False alarm rate (golden batches)', 'Value': '1/25 (4%)'},
    {'Metric': 'Max T² (Run-026)', 'Value': '6.82 (limit: 9.98)'},
    {'Metric': 'Max SPE (Run-026)', 'Value': '7.40 (limit at that hour: 4.28)'},
    {'Metric': 'Time saved vs offline metabolite panel', 'Value': '53 hours (2.2 days)'},
    {'Metric': 'Remaining batch time at alarm', 'Value': '197 hours (8.2 days)'},
    {'Metric': 'Univariate alarms triggered', 'Value': '0 / 11 variables'},
    {'Metric': 'Model components', 'Value': '3 PCs (87.7% variance explained)'},
    {'Metric': 'NOC reference set', 'Value': '25 golden batches'},
    {'Metric': 'Process variables monitored', 'Value': '11'},
    {'Metric': 'Batch value at risk', 'Value': '$500K–2M (materials + opportunity cost)'},
]

metrics_df = pd.DataFrame(metrics_data)
print(metrics_df.fillna('').to_markdown(index=False))
| Metric                                 | Value                                        |
|:---------------------------------------|:---------------------------------------------|
| Metabolic shift onset                  | Hour 100 (Day 4.2)                           |
| SPE alarm fires                        | Hour 139 (Day 5.8)                           |
| Detection delay from onset             | 39 hours (1.6 days)                          |
| T² alarm                               | Does not trigger (fault orthogonal to model) |
| False alarm rate (golden batches)      | 1/25 (4%)                                    |
| Max T² (Run-026)                       | 6.82 (limit: 9.98)                           |
| Max SPE (Run-026)                      | 7.40 (limit at that hour: 4.28)              |
| Time saved vs offline metabolite panel | 53 hours (2.2 days)                          |
| Remaining batch time at alarm          | 197 hours (8.2 days)                         |
| Univariate alarms triggered            | 0 / 11 variables                             |
| Model components                       | 3 PCs (87.7% variance explained)             |
| NOC reference set                      | 25 golden batches                            |
| Process variables monitored            | 11                                           |
| Batch value at risk                    | $500K–2M (materials + opportunity cost)      |

Detected within 39 hours of metabolic onset. Alarmed 2.2 days before the scheduled offline metabolite panel would have reliably caught the lactate excursion. Operator alerted with 8.2 days of remaining batch time to intervene.

References

Nomikos, P. & MacGregor, J.F. (1995). "Multivariate SPC Charts for Monitoring Batch Processes." Technometrics, 37(1), 41-59.
Jackson, J.E. & Mudholkar, G.S. (1979). "Control Procedures for Residuals Associated with Principal Component Analysis." Technometrics, 21(3), 341-349.
Wise, B.M. & Gallagher, N.B. (1996). "The process chemometrics approach to process monitoring and fault detection." Journal of Process Control, 6(6), 329-348.
Kourti, T. & MacGregor, J.F. (1995). "Process analysis, monitoring and diagnosis, using multivariate projection methods." Chemometrics and Intelligent Laboratory Systems, 28(1), 3-21.

More reports

Digital TwinJune 4, 2026

Richelle et al. (2022) Digital Twin — Model-Based CHO Intensification

A from-scratch implementation of the mechanistic CHO growth model from Richelle et al. (2022) as an executable digital twin — five kinetic parameters identified from fed-batch data, predicting culture dynamics from fed-batch through intensified perfusion.

View report

Tech TransferJune 2, 2026

PCA Reference Model — SiteB Tech Transfer Comparability

A PCA reference model trained on 84 primary-site runs, then used to project 15 SiteB tech-transfer runs and assess comparability — reference-only training that mirrors validating a model before new runs arrive.

View report

Golden BatchJune 1, 2026

Golden Batch Comparison: PD-GB vs Campaign Runs

A 1 L golden-batch reference (PD-GB) compared against a 13-run iPSC scale-up campaign from 1 L to 10 L — tracking differentiation-marker trajectories and each run's deviation from the golden profile.

View report