Bio-Manufacturing Unit Operations

This document provides comprehensive documentation for all bio-manufacturing unit operations available in the difflow_bio module.

---

Bioreactors

ContinuousBioreactor (Chemostat)

Location: difflow_bio/units/bioreactors.py

Class: ContinuousBioreactor

Description: Models a continuous stirred-tank bioreactor (chemostat) at steady state. The reactor maintains constant volume with continuous feed and withdrawal, allowing for steady-state cell cultivation.

Process Role

Continuous bioreactors are used for:

• Large-scale production of microbial products

• Maintaining cells in exponential growth phase

• Steady-state operation for consistent product quality

• High-volume, low-value products (ethanol, organic acids)

Parameters

@dataclass
class BioreactorParams:
    V: float               # Reactor volume (L)
    Y_xs: float            # Biomass yield on substrate (g cell/g substrate)
    kinetic_fn: Callable   # Growth kinetics mu(S, params) or mu(S, X, params)
    kinetic_params: dict   # Passed to kinetic_fn, e.g. {"mu_max": 0.4, "K_s": 0.5}
    k_d: float = 0.0       # Cell death rate (1/h)
    m_s: float = 0.0       # Maintenance coefficient (g substrate/g cell/h)
    alpha: float = 0.0     # Growth-associated product formation (g/g)
    beta: float = 0.0      # Non-growth-associated product formation (g/g/h)
    species_order: list = ["cells", "substrate", "product"]

The growth law is the kinetic_fn: monod_kinetics reads mu_max and K_s from kinetic_params, and the inhibition models (substrate_inhibition_kinetics, …) read their own constants from the same dict.

Inputs

Parameter	Type	Units	Description
inlet	Stream	-	Feed stream; its substrate entry is the feed substrate concentration \(S_f\)
D	float	1/h	Dilution rate (F/V)
F	float	L/h	Volumetric feed rate, used when D is None (\(D = F/V\))

Outputs

Parameter	Type	Units	Description
outlet	Stream	-	Outlet stream
info['X']	float	g/L	Cell concentration
info['S']	float	g/L	Substrate concentration
info['P']	float	g/L	Product concentration
info['mu']	float	1/h	Specific growth rate
info['productivity']	float	g/L/h	Volumetric productivity

Governing Equations

Cell Balance (steady-state):

\[\frac{dX}{dt} = 0 = (\mu - D - k_d) X\]

At steady state: \(\mu = D + k_d\) (washout condition: \(D < \mu_{max}\))

Substrate Balance:

\[\frac{dS}{dt} = 0 = D(S_f - S) - \frac{\mu X}{Y_{X/S}} - m_s X\]

Solving for steady-state substrate:

\[S = \frac{K_s D}{(\mu_{max} - D)}\]

Cell Concentration:

\[X = Y_{X/S}(S_f - S) - m_s \frac{X}{D}\]

Simplified (negligible maintenance):

\[X = Y_{X/S}(S_f - S)\]

Product Formation (Luedeking-Piret):

\[\frac{dP}{dt} = (\alpha \mu + \beta) X - D \cdot P\]

At steady state:

\[P = \frac{(\alpha \mu + \beta) X}{D}\]

Where:

• \(\alpha\): Growth-associated product formation coefficient

• \(\beta\): Non-growth-associated product formation rate

Productivity:

(1)\[P_r = D \cdot X\]

(2)\[P_r = D \cdot P\]

Critical Dilution Rate

Washout occurs when \(D > \mu_{max}\). The critical dilution rate:

\[D_{crit} = \mu_{max} \frac{S_f}{K_s + S_f}\]

Optimal Dilution Rate (maximum productivity):

\[D_{opt} = \mu_{max} \left(1 - \sqrt{\frac{K_s}{K_s + S_f}}\right)\]

Example Usage

from difflow.streams import make_stream
from difflow_bio.units.bioreactors import (
    BioreactorParams, ContinuousBioreactor, monod_kinetics)

params = BioreactorParams(
    V=10.0,                                       # L
    Y_xs=0.5,                                     # g/g
    kinetic_fn=monod_kinetics,
    kinetic_params={"mu_max": 0.4, "K_s": 0.5},   # 1/h, g/L
    alpha=0.05,         # growth-associated
    beta=0.01,          # non-growth-associated
    species_order=["cells", "substrate", "product"],
)

bioreactor = ContinuousBioreactor(params)
# S_f = 50 g/L is the substrate entry of the feed, not a call argument
feed = make_stream({'substrate': 50.0, 'cells': 0.0, 'product': 0.0}, T=310.0, P=101325.0)

outlet, info = bioreactor(feed, D=0.2)
print(f"Cell concentration: {info['X']:.2f} g/L")
print(f"Productivity: {info['productivity']:.3f} g/L/h")

---

FedBatchBioreactor

Location: difflow_bio/units/bioreactors.py

Class: FedBatchBioreactor

Description: Models a fed-batch bioreactor where substrate is added during the batch to control growth rate and avoid substrate inhibition.

Process Role

Fed-batch bioreactors are used for:

• High-value products (therapeutic proteins, antibodies)

• Avoiding substrate/product inhibition

• Maximizing final product titer

• Most industrial mAb production

Parameters

@dataclass
class FedBatchParams:
    V0: float              # Initial volume (L)
    Y_xs: float            # Yield coefficient (g cells / g substrate)
    kinetic_fn: Callable   # Growth kinetics function
    kinetic_params: dict   # Parameters for kinetic function
    k_d: float = 0.0       # Death rate constant (1/h)
    m_s: float = 0.0       # Maintenance coefficient (g/g/h)
    K_m: float = 0.05      # Half-saturation of maintenance uptake (g/L)
    alpha: float = 0.0     # Growth-associated product formation (g/g)
    beta: float = 0.0      # Non-growth-associated product formation (g/g/h)
    species_order: list = ["cells", "substrate", "product"]

Integration: Uses diffrax (adaptive Tsit5 by default) or RK4 fallback

Inputs

Parameter	Type	Units	Description
inlet	Stream	-	Feed stream composition
feed_rate	Callable	L/h	Feed rate as function of time
X_0	float	g/L	Initial cell concentration
S_0	float	g/L	Initial substrate concentration

Governing Equations

Volume Change:

\[\frac{dV}{dt} = F(t)\]

Cell Balance:

\[\frac{d(VX)}{dt} = V(\mu - k_d)X - V m_s Y_{X/S} (1 - \phi) X\]

Substrate Balance:

\[\frac{d(VS)}{dt} = F \cdot S_f - V\left(\frac{\mu X}{Y_{X/S}} + m_s \phi X\right), \qquad \phi = \frac{S}{K_m + S}\]

Maintenance under starvation: maintenance draws on substrate while there is any (\(\phi \approx 1\) for \(S \gg K_m\)) and on biomass when it runs out, at the endogenous decay rate \(b = m_s Y_{X/S}\) (Pirt, 1965). Without the \(\phi\) factor, uptake would continue at \(m_s X\) after the substrate is gone and drive \(S\) negative while product kept accumulating. info["S_min"] reports the lowest substrate concentration reached, so a starved run is visible.

Product Balance:

\[\frac{d(VP)}{dt} = V(\alpha \mu + \beta)X\]

Common Feed Strategies:

1. Constant Feed: \(F(t) = F_0\)

2. Exponential Feed: \(F(t) = F_0 e^{\mu_{set} t}\) to maintain constant \(\mu\)

3. Feedback Control: Adjust F based on measured substrate or DO

Example Usage

from difflow_bio.units.bioreactors import FedBatchBioreactor, FedBatchParams, monod_kinetics
import jax.numpy as jnp

params = FedBatchParams(
    V0=5.0,            # Initial volume (L)
    Y_xs=0.5,          # Yield coefficient
    kinetic_fn=monod_kinetics,
    kinetic_params={'mu_max': 0.4, 'K_s': 0.5},
    k_d=0.01,          # Death rate
    alpha=0.1,         # Growth-associated product formation
    beta=0.02          # Non-growth-associated product formation
)

bioreactor = FedBatchBioreactor(params)

# Exponential feed to maintain mu = 0.1/h
def exponential_feed(t):
    return 0.1 * jnp.exp(0.1 * t)

# Run simulation
outlet, info = bioreactor(
    X0=0.5,              # Initial cell concentration (g/L)
    S0=20.0,             # Initial substrate (g/L)
    P0=0.0,              # Initial product (g/L)
    t_final=72.0,        # hours
    feed_rate_fn=exponential_feed,
    S_feed=200.0,        # Feed substrate concentration (g/L)
)
print(f"Final cell concentration: {info['X_final']:.2f} g/L")
print(f"Final product: {info['P_final']:.2f} g/L")

---

Kinetic Models

Location: difflow_bio/units/bioreactors.py

Monod Kinetics

\[\mu = \mu_{max} \frac{S}{K_s + S}\]

from difflow_bio.units.bioreactors import monod_kinetics

mu = monod_kinetics(5.0, {"mu_max": 0.4, "K_s": 0.5})

Substrate Inhibition (Andrews/Haldane)

\[\mu = \mu_{max} \frac{S}{K_s + S + S^2/K_i}\]

from difflow_bio.units.bioreactors import substrate_inhibition_kinetics

mu = substrate_inhibition_kinetics(5.0, {"mu_max": 0.4, "K_s": 0.5, "K_i": 50.0})

Product Inhibition

\[\mu = \mu_{max} \frac{S}{K_s + S} \left(1 - \frac{P}{P_{max}}\right)^n\]

from difflow_bio.units.bioreactors import product_inhibition_kinetics

mu = product_inhibition_kinetics(
    5.0, 10.0, {"mu_max": 0.4, "K_s": 0.5, "P_max": 100.0, "n": 1})

Contois Kinetics (High Cell Density)

\[\mu = \mu_{max} \frac{S}{K_{sX} X + S}\]

from difflow_bio.units.bioreactors import contois_kinetics

mu = contois_kinetics(5.0, 50.0, {"mu_max": 0.4, "K_s": 0.1})   # K_s is K_sX here

Utility Functions

from difflow_bio.units.bioreactors import (
    dilution_rate,
    residence_time,
    optimal_dilution_rate,
)

D = dilution_rate(F=100.0, V=500.0)  # 0.2 1/h
tau = residence_time(D)  # 5 h
D_opt = optimal_dilution_rate({"mu_max": 0.4, "K_s": 0.5, "S_f": 50.0})

The washout dilution rate is a method of the reactor, ContinuousBioreactor.washout_dilution_rate().

---

Centrifugation

Centrifuge

Location: difflow_bio/units/centrifuge.py

Class: Centrifuge

Description: Models centrifugal separation of cells/particles from liquid based on density difference. Uses Sigma factor theory for scale-up.

Process Role

Centrifugation is used for:

• Cell harvest (primary recovery)

• Cell debris removal after homogenization

• Precipitate collection

• Clarification before chromatography

Parameters

@dataclass
class CentrifugeParams:
    sigma: float                    # Sigma factor (equivalent settling area, m²)
    efficiency: float = 0.7         # Separation efficiency (0-1)
    species_order: list = None
    cell_species: str = "cells"
    lysis_threshold_g: float = None # RCF above which cells lyse (x g)
    lysis_coefficient: float = 0.0  # Lysed fraction per unit RCF above threshold
    lysate_species: str = "lysate"

The particle and fluid properties are call arguments, not parameters: d_particle (m), rho_particle and rho_fluid (kg/m³), viscosity (Pa·s), plus concentrate_fraction and rcf.

Inputs

Parameter	Type	Units	Description
inlet	Stream	-	Feed with cells and broth
Q	float	m³/s	Volumetric throughput

Outputs

Parameter	Type	Units	Description
concentrate	Stream	-	Concentrated cells/solids
supernatant	Stream	-	Clarified liquid
info['cell_recovery']	float	-	Fraction of the cells sent to the concentrate
info['separation_efficiency']	float	-	Sigma-theory capture efficiency
info['stokes_velocity']	float	m/s	Settling velocity of the mean particle
info['critical_diameter']	float	m	Particle size captured at 50% efficiency at this throughput
info['Q_critical']	float	m³/s	Sigma-theory critical throughput, \(2 v_s \Sigma\)
info['concentration_factor']	float	-	Cell concentration, concentrate over feed
info['lysis_fraction'], info['lysate_released']	float	-	Shear lysis (zero unless lysis_threshold_g is set)

Governing Equations

Stokes Settling Velocity (single particle, laminar flow):

\[v_s = \frac{d_p^2 (\rho_p - \rho_f) g}{18 \mu}\]

Where:

• \(d_p\): Particle diameter (m)

• \(\rho_p\): Particle density (kg/m³)

• \(\rho_f\): Fluid density (kg/m³)

• \(\mu\): Fluid viscosity (Pa·s)

• \(g\): Gravitational acceleration (9.81 m/s²)

Centrifugal Enhancement:

In a centrifuge, \(g\) is replaced by centrifugal acceleration:

\[a_c = \omega^2 r\]

Sigma Factor (equivalent settling area):

The Sigma factor allows scale-up between different centrifuge types:

\[\Sigma = \frac{Q}{2 v_s}\]

For a given separation: \(Q_1/\Sigma_1 = Q_2/\Sigma_2\)

Critical Particle Diameter (100% capture):

\[d_{crit} = \sqrt{\frac{18 \mu Q}{(\rho_p - \rho_f) \omega^2 \Sigma}}\]

Separation Efficiency (particle size distribution):

\[E = 1 - \exp\left(-\frac{d^2}{d_{crit}^2}\right)\]

G-Force:

\[G = \frac{\omega^2 r}{g} = \frac{(2\pi N)^2 r}{g}\]

Where N is rotational speed (rev/s).

Example Usage

from difflow.streams import make_stream
from difflow_bio.units.centrifuge import Centrifuge, CentrifugeParams

params = CentrifugeParams(sigma=5000.0, efficiency=0.95)   # m², -

centrifuge = Centrifuge(params)
feed = make_stream({'cells': 50.0, 'broth': 950.0}, T=298.0, P=101325.0)

# particle and fluid properties are call arguments
concentrate, supernatant, info = centrifuge(
    feed, Q=0.001,              # 1 L/s = 3.6 m³/h
    d_particle=5e-6,            # 5 μm cells
    rho_particle=1050.0, rho_fluid=1000.0, viscosity=0.001,
)
print(f"Cell recovery: {info['cell_recovery']:.2%}")

---

DiscStackCentrifuge

Location: difflow_bio/units/centrifuge.py

Class: DiscStackCentrifuge

Description: Specialized model for disc-stack centrifuges, the most common type in bioprocessing for cell harvest.

Parameters

@dataclass
class DiscStackParams:
    n_discs: int           # Number of discs
    r_inner: float         # Inner disc radius (m)
    r_outer: float         # Outer disc radius (m)
    cone_angle: float      # Disc cone angle (radians)
    rpm: float             # Rotational speed (rev/min)

Sigma Factor Calculation

\[\Sigma = \frac{2\pi n \omega^2 (r_o^3 - r_i^3)}{3g \tan\theta}\]

Where:

• \(n\): Number of discs

• \(\omega\): Angular velocity (rad/s)

• \(r_o, r_i\): Outer and inner disc radii

• \(\theta\): Half cone angle

Example Usage

from difflow_bio.units.centrifuge import DiscStackCentrifuge, DiscStackParams

params = DiscStackParams(
    n_discs=100,
    r_inner=0.05,    # 5 cm
    r_outer=0.15,    # 15 cm
    half_angle=0.785, # 45 degrees
    rpm=7000
)

centrifuge = DiscStackCentrifuge(params)
# Sigma is calculated internally based on geometry

Utility Functions

from difflow_bio.units.centrifuge import (
    stokes_velocity,
    critical_particle_diameter,
    disc_stack_sigma,
    tubular_bowl_sigma,
    centrifuge_scale_up,
    g_force
)

# Calculate Stokes velocity
v_s = stokes_velocity(d=5e-6, rho_p=1050, rho_f=1000, mu=0.001)

# Calculate Sigma for disc stack
sigma = disc_stack_sigma(n_discs=100, r_inner=0.05, r_outer=0.15, half_angle=0.785, rpm=7000.0)

# Scale-up calculation
Q2 = centrifuge_scale_up(sigma_1=1000.0, Q_1=1.0, sigma_2=10000.0)

# G-force
G = g_force(rpm=7000, r=0.15)  # ~7900 G

---

Membrane Filtration

Ultrafiltration

Location: difflow_bio/units/filtration.py

Class: Ultrafiltration

Description: Pressure-driven membrane separation that concentrates proteins based on molecular weight cutoff (MWCO).

Process Role

Ultrafiltration is used for:

• Protein concentration

• Buffer exchange (combined with diafiltration)

• Virus removal (as secondary barrier)

• Harvest clarification (with larger MWCO)

Parameters

@dataclass
class UltrafiltrationParams:
    membrane_area: float           # Membrane area (m²); sets the process time only
    MWCO: float = 30.0             # Molecular weight cutoff (kDa), rejected at 90%
    rejection: dict = {}           # species -> R, overrides the MWCO-derived value
    molecular_weights: dict = {}   # species -> MW (kDa), overrides the tables
    Lp: float = 50.0               # Permeability (L/m²/h/bar)
    k_mass: float = 5e-6           # Mass transfer coefficient (m/s)
    sigma: float = 1000.0          # Osmotic pressure coefficient (Pa·m³/kg)
    fouling_coefficient: float = 0.0  # Flux decline per L permeate (1/L)
    species_order: list = None

Each species’ rejection is its entry in rejection if there is one, and otherwise follows the MWCO against its molecular weight, looked up in molecular_weights, then difflow_bio.database.BIO_SPECIES_MW_KDA (mAb 150 kDa, aggregates 300, fragments 50, HCP 50, DNA 330, and small solutes; representative values), then the core species database. A species with no molecular weight anywhere passes freely (\(R = 0\)), the right default for buffer components. membrane_area does not change the split (sieving is per unit area); given feed_volume (L) in the call, the unit reports info['process_time_h'], the permeate volume over flux times area.

Inputs

Parameter	Type	Units	Description
inlet	Stream	-	Feed stream
concentration_factor	float	-	Volume concentration factor (required)
TMP	float	bar	Transmembrane pressure (default 1.0)
mode	str	-	"batch" (default) or "continuous"

Outputs

Parameter	Type	Units	Description
retentate	Stream	-	Concentrated product
permeate	Stream	-	Permeate (removed species)
info['flux']	float	LMH	Permeate flux (L/m²/h)
info['recovery']	dict	-	Fraction of each species kept in the retentate
info['concentration_factor'], info['volume_reduction'], info['retentate_volume_fraction']	float	-	Volume bookkeeping
info['fouling_factor']	float	-	Flux decline from fouling (1 when off)
info['rejection']	dict	-	Rejection used for each species: the rejection override, else the MWCO-derived value

The call returns ((retentate, permeate), info): the two streams are nested in a tuple of their own.

Governing Equations

Permeate Flux (resistance model):

\[J = \frac{TMP}{\mu (R_m + R_c + R_g)}\]

Where:

• \(J\): Permeate flux (L/m²/h or LMH)

• \(TMP\): Transmembrane pressure (Pa)

• \(\mu\): Permeate viscosity (Pa·s)

• \(R_m\): Membrane resistance (1/m)

• \(R_c\): Concentration polarization resistance

• \(R_g\): Gel layer resistance

Rejection Coefficient:

\[R_i = 1 - \frac{C_{i,permeate}}{C_{i,retentate}}\]

From the MWCO (rejection_from_mw), a logistic sieving curve in \(\ln MW\) with \(R = 0.9\) at the cutoff, the nominal rating:

\[R = \left[1 + \tfrac{1}{9}\left(\frac{MW}{MWCO}\right)^{-3}\right]^{-1}\]

so a solute five times the cutoff (a 150 kDa mAb on a 30 kDa membrane) is rejected at more than 99.9% and a small solute passes freely.

Concentration Factor:

(3)\[CF = \frac{V_{feed}}{V_{retentate}} = \frac{C_{retentate}}{C_{feed}}\]

Mass Balance:

\[V_f C_{f,i} = V_r C_{r,i} + V_p C_{p,i}\]

\[C_{r,i} = \frac{C_{f,i}}{1 - (1-R_i)(1 - 1/CF)}\]

Concentration Polarization:

\[\frac{C_w - C_p}{C_b - C_p} = \exp\left(\frac{J}{k}\right)\]

Where:

• \(C_w\): Wall concentration

• \(C_b\): Bulk concentration

• \(k\): Mass transfer coefficient

Gel Polarization Model (flux-limited):

\[J = k \ln\left(\frac{C_g}{C_b}\right)\]

Where \(C_g\) is the gel concentration (limiting).

Example Usage

from difflow.streams import make_stream
from difflow_bio.units.filtration import Ultrafiltration, UltrafiltrationParams

params = UltrafiltrationParams(
    membrane_area=5.0,                       # m²
    MWCO=30.0,                               # kDa
    rejection={"mAb": 0.999, "HCP": 0.5},    # buffer passes freely
    Lp=100.0,                                # LMH/bar
)

uf = Ultrafiltration(params)
feed = make_stream({'mAb': 2.0, 'HCP': 0.5, 'buffer': 997.5}, T=298.0, P=101325.0)

(retentate, permeate), info = uf(feed, concentration_factor=10.0, TMP=2.0)
print(f"mAb in retentate: {retentate['F_mAb']:.3f}")
print(f"Flux: {info['flux']:.1f} LMH")

---

Diafiltration

Location: difflow_bio/units/filtration.py

Class: Diafiltration

Description: Ultrafiltration with continuous buffer addition to exchange buffer or remove small molecules while retaining proteins.

Process Role

Diafiltration is used for:

• Buffer exchange

• Salt removal (desalting)

• Small molecule impurity removal

• Formulation preparation

Parameters

@dataclass
class DiafiltrationParams:
    membrane_area: float           # Membrane area (m²); sets the process time only
    MWCO: float = 30.0             # kDa; rejection as for Ultrafiltration
    rejection: dict = {}
    molecular_weights: dict = {}
    Lp: float = 50.0
    k_mass: float = 5e-6
    sigma: float = 1000.0
    fouling_coefficient: float = 0.0

The call is df(inlet, buffer, n_diavolumes=None, TMP=1.0, feed_volume=None) and returns ((retentate, permeate), info). With n_diavolumes=None the buffer stream is the buffer fed: it enters the balance and sets \(N_{DV}\) = buffer / inlet (amounts as the volume proxy). With n_diavolumes given, only the buffer’s composition is used, scaled to \(N_{DV}\) times the inlet, and info['buffer_consumed'] is the buffer that entered, so retentate + permeate = inlet + buffer consumed for every species.

Governing Equations

Diavolume Definition:

\[N_{DV} = \frac{V_{buffer\,added}}{V_{retentate}}\]

Impurity Removal (constant volume diafiltration):

\[\frac{C}{C_0} = \exp(-N_{DV}(1-R))\]

For freely permeating species (\(R = 0\)):

\[\frac{C}{C_0} = \exp(-N_{DV})\]

Buffer wash-in of a species fed at \(C_b\) with sieving \(s = 1 - R\):

\[\frac{C}{C_b} = \frac{1 - \exp(-s N_{DV})}{s}\]

which tends to \(N_{DV}\) (everything added stays) as \(R \to 1\).

Required Diavolumes for target removal:

\[N_{DV} = -\frac{\ln(C/C_0)}{1-R}\]

Diavolumes	Removal (R=0)
1	63.2%
2	86.5%
3	95.0%
5	99.3%
7	99.9%

Example Usage

from difflow_bio.units.filtration import Diafiltration, DiafiltrationParams

params = DiafiltrationParams(MWCO=30.0, membrane_area=5.0)  # 30 kDa

df = Diafiltration(params)
feed = make_stream({'mAb': 10.0, 'salt': 150.0, 'buffer': 840.0}, T=298.0, P=101325.0)
new_buffer = make_stream({'new_buffer': 5000.0}, T=298.0, P=101325.0)  # 5 diavolumes

(retentate, permeate), info = df(feed, new_buffer)
print(f"Diavolumes: {float(info['n_diavolumes']):.1f}")
print(f"Salt removal: {1 - float(retentate['F_salt'] / feed['F_salt']):.1%}")
print(f"mAb kept: {float(retentate['F_mAb']) / 10.0:.2%}")

---

TFF (Tangential Flow Filtration)

Location: difflow_bio/units/filtration.py

Class: TFF

Description: Combines ultrafiltration and diafiltration in a single system with tangential flow to minimize fouling.

Key Features

• Recirculation loop reduces concentration polarization

• Feed flows parallel to membrane surface

• Multiple passes possible for high concentration

Process Modes

1. Concentration: UF mode to reduce volume

2. Diafiltration: Buffer exchange at constant volume

3. UF/DF Sequence: Concentrate → Diafiltrate → Final concentration

TFF.uf_df_uf(inlet, buffer, CF_initial, n_diavolumes, CF_final) returns ((product, uf1_permeate, df_permeate, uf2_permeate), info); the four streams sum to the inlet plus info['buffer_consumed'] for every species.

---

Chromatography

ProteinAChromatography

Location: difflow_bio/units/chromatography.py

Class: ProteinAChromatography

Description: Affinity chromatography using Protein A ligand for mAb capture. The primary capture step in most mAb processes.

Process Role

Protein A chromatography provides:

• High selectivity for Fc-containing antibodies

• 95% purity in single step

• Significant HCP and DNA clearance

• Recovery >90%

Parameters

@dataclass
class ProteinAParams:
    column_volume: float          # Column volume (L)
    q_max: float = 35.0           # Maximum binding capacity (g/L resin)
    K_d: float = 0.1              # Dissociation constant (g/L)
    target_species: str = "mAb"
    yield_factor: float = 0.95    # Elution yield (0-1)
    impurity_clearance: dict      # impurity -> LRV, default HCP 2, DNA 3, cells 4
    k_ads: float | None = None    # Adsorption rate (1/min), for kinetic DBC
    n_plates: float | None = None # Plate count, for elution pool volume
    elution_cv: float = 2.0       # Elution retention volume (CV)

Inputs

Parameter	Type	Units	Description
inlet	Stream	-	Feed stream with mAb (component flows as mass)
load_volume	float or None	L	Volume of feed loaded. None (default) loads the whole inlet; the column capacity, \(DBC \cdot V_{column}\), then limits what binds and the excess breaks through
feed_volume	float	L	Total feed volume, so load_volume / feed_volume is the fraction loaded. Required with load_volume (the stream carries amounts, not a volume; a ValueError says so). With it, the capacity uses the Langmuir loading at the feed concentration, \(q_{max} C/(K_d + C)\)

Outputs

Parameter	Type	Units	Description
product	Stream	-	Elution pool (purified mAb)
waste	Stream	-	Flow-through, wash and column losses
info['yield']	float	-	Target in the product / target in the inlet (feed not loaded counts as lost)
info['purity']	float	-	Product purity
info['mass_loaded'], info['mass_bound']	float	mass	Target loaded and bound (bound is capped by capacity)
info['capacity_utilization']	float	-	Bound mass / column capacity
info['impurity_clearance']	dict	LRV	Clearance applied per impurity
info['DBC']	float	g/L	Dynamic binding capacity

Process Steps

1. Equilibration: Condition column with loading buffer

2. Load: Apply feed, mAb binds to resin

3. Wash: Remove unbound impurities

4. Elution: Release mAb with low pH buffer

5. Regeneration: Clean and re-equilibrate column

Governing Equations

Langmuir Isotherm (equilibrium binding):

\[q = \frac{q_{max} C}{K_d + C}\]

Where:

• \(q\): Bound concentration (g/L resin)

• \(q_{max}\): Maximum binding capacity

• \(C\): Solution concentration

• \(K_d\): Dissociation constant

Dynamic Binding Capacity (at 10% breakthrough):

\[DBC_{10\%} = q_{max} \cdot f(RT, C_{load}, K_d)\]

Typically 80-95% of equilibrium capacity.

Column Capacity Utilization:

\[\text{Load} = \frac{m_{mAb,loaded}}{V_{column} \cdot DBC}\]

Typical loading: 70-90% of DBC.

Log Reduction Value:

\[LRV = \log_{10}\left(\frac{C_{in}}{C_{out}}\right)\]

Example Usage

from difflow_bio.units.chromatography import ProteinAChromatography, ProteinAParams

from difflow import make_stream

params = ProteinAParams(
    column_volume=10.0,    # L
    q_max=35.0,            # g/L (typical for MabSelect)
)

protein_a = ProteinAChromatography(params)
load = make_stream({'mAb': 250.0, 'HCP': 5.0, 'DNA': 0.1}, T=298.0, P=101325.0)

# Load the whole batch: 250 g of mAb on a 10 L column (346.5 g capacity)
(eluate, waste), info = protein_a(load)
print(f"Yield: {float(info['yield']):.2%}")
print(f"Capacity used: {float(info['capacity_utilization']):.0%}")
print(f"HCP clearance: {info['impurity_clearance']['HCP']:.1f} LRV")

---

IonExchangeChromatography

Location: difflow_bio/units/chromatography.py

Class: IonExchangeChromatography

Description: Separation based on electrostatic interactions between charged proteins and charged resin.

Types

• Cation Exchange (CEX): Negatively charged resin, binds positively charged proteins

• Anion Exchange (AEX): Positively charged resin, binds negatively charged proteins

Process Role

Ion exchange is used for:

• Aggregate removal (CEX)

• Charge variant separation (CEX)

• DNA/HCP removal (AEX in flow-through mode)

• Viral clearance (AEX)

Parameters

@dataclass
class IEXParams:
    column_volume: float              # Column volume (L)
    mode: str = "bind_elute"          # 'bind_elute' or 'flow_through'
    q_max: float = 50.0               # Binding capacity (g/L)
    K_d: float = 0.5                  # Dissociation constant (g/L)
    target_species: str = "mAb"
    selectivity: dict = {}            # species -> binding selectivity (0-1)
    yield_factor: float = 0.90        # Step recovery of the target
    impurity_clearance: dict = {}     # impurity -> LRV across the step

CEX and AEX are the same unit; mode and the parameters say which one it is. load_volume (L) is optional when calling it, and None loads the whole inlet; a partial load needs feed_volume (L) as well. The column holds at most \(q \cdot V_{column}\), with \(q = q_{max} C / (K_d + C)\) when the feed concentration is known (feed_concentration, or the amount over feed_volume) and \(q = q_{max}\) otherwise. In bind-elute mode target beyond that capacity breaks through to waste; in flow-through mode the impurity binding is scaled down to fit it. info['yield'] is the target in the product over the target in the inlet, and info['capacity'] reports the capacity used.

Impurity clearance

An impurity listed in impurity_clearance reaches the product as \(10^{-\mathrm{LRV}}\) of what was loaded, in either mode; the rest goes to waste. Impurities not listed fall back to a selectivity rule: in flow-through mode a fraction selectivity binds, and in bind-elute mode \(s^2 (1 - Y)\) co-elutes with the product. That rule caps carry-over at \(1 - Y\), so a 90% step would remove at least 90% of every impurity, far more than a polishing step removes aggregate. State clearances explicitly when purity matters.

TYPICAL_CEX_CLEARANCE and TYPICAL_AEX_CLEARANCE in the same module are the representative values the packaged trains use (CEX: HCP 0.5, DNA 1.0, aggregates 0.7; AEX: HCP 1.5, DNA 3.0, aggregates 0.1). They are not measured clearances. They are chosen so that a harvest at about 5% HCP and 3% aggregate ends near reported end-of-process levels: about 10 ppm HCP against a <100 ppm target, and aggregate below the usual 1% target. The AEX HCP value matches one reported flow-through step, 530 to 15 ppm (Liu et al., mAbs 2:480, 2010, doi:10.4161/mabs.2.5.12645). Replace them with your own process data. With them, mAbDSPTrain reports about 99.4% purity, 6 ppm HCP and 0.6% aggregate from that harvest, and returns hcp_ppm and aggregate_fraction alongside purity.

Operating Modes

Bind-Elute Mode:

1. Product binds to resin

2. Impurities wash through or elute at different salt concentrations

3. Product eluted with salt gradient

Flow-Through Mode:

1. Product passes through column

2. Impurities bind to resin

3. Common for AEX after Protein A

Governing Equations

Selectivity:

\[\alpha_{ij} = \frac{K_i}{K_j}\]

Resolution:

\[R_s = \frac{t_{R,2} - t_{R,1}}{0.5(w_1 + w_2)} = \frac{\sqrt{N}}{4} \cdot \frac{\alpha - 1}{\alpha} \cdot \frac{k'}{1 + k'}\]

Where:

• \(N\): Number of theoretical plates

• \(\alpha\): Selectivity

• \(k'\): Capacity factor

---

SizeExclusionChromatography

Location: difflow_bio/units/chromatography.py

Class: SizeExclusionChromatography

Description: Separation based on molecular size. Large molecules elute first (excluded from pores), small molecules elute last.

Process Role

SEC is used for:

• Aggregate/fragment removal (polishing)

• Buffer exchange

• Molecular weight determination (analytical)

• Final formulation preparation

Parameters

@dataclass
class SECParams:
    column_volume: float   # Column volume (L)
    void_fraction: float   # Interparticle void (V_0/V_c)
    total_porosity: float  # Total porosity (V_t/V_c)
    exclusion_limit: float # MW exclusion limit (Da)
    permeation_limit: float # MW permeation limit (Da)

Governing Equations

Distribution Coefficient:

\[K_d = \frac{V_e - V_0}{V_t - V_0}\]

Where:

• \(V_e\): Elution volume

• \(V_0\): Void volume (excluded volume)

• \(V_t\): Total accessible volume

Calibration (MW vs elution volume):

\[\log(MW) = a - b \cdot K_d\]

Or: \(K_d = a' - b' \cdot \log(MW)\)

Resolution Limits:

• \(K_d = 0\): Totally excluded (MW > exclusion limit)

• \(K_d = 1\): Fully permeating (MW < permeation limit)

---

Adsorption Isotherms

Location: difflow_bio/units/chromatography.py

Langmuir Isotherm

\[q = \frac{q_{max} C}{K_d + C}\]

from difflow_bio.units.chromatography import langmuir_isotherm

q = langmuir_isotherm(C=5.0, q_max=35.0, K_d=0.5)

Linear Isotherm (dilute systems)

\[q = K \cdot C\]

from difflow_bio.units.chromatography import linear_isotherm

q = linear_isotherm(C=0.1, K=100.0)

Langmuir-Freundlich (heterogeneous sites)

\[q = \frac{q_{max} (C/K_d)^n}{1 + (C/K_d)^n}\]

from difflow_bio.units.chromatography import langmuir_freundlich_isotherm

q = langmuir_freundlich_isotherm(C=5.0, q_max=35.0, K_d=0.5, n=0.8)

Chromatography Utility Functions

from difflow_bio.units.chromatography import (
    dynamic_binding_capacity,
    column_productivity,
    resolution,
    plate_count,
    hetp
)

# Dynamic binding capacity at 10% breakthrough
DBC = dynamic_binding_capacity(q_max=35.0, C_feed=5.0, K_d=0.5, residence_time=6.0, k_ads=1.0)

# Column productivity
P = column_productivity(DBC=30.0, column_volume=1.0, cycle_time=4.0)  # g/L/h

# Resolution between peaks
Rs = resolution(t_R1=10.0, t_R2=12.0, w1=0.5, w2=0.6)

# Theoretical plates
N = plate_count(t_R=10.0, w=0.5)

# Height equivalent to theoretical plate
H = hetp(L=20.0, N=10000)  # cm

---

Summary: Typical mAb Downstream Process

Harvest (Bioreactor)
        │
        ▼
┌───────────────────┐
│   Centrifugation  │  Cell removal
└───────────────────┘
        │
        ▼
┌───────────────────┐
│ Depth Filtration  │  Clarification
└───────────────────┘
        │
        ▼
┌───────────────────┐
│    Protein A      │  Capture (>95% purity)
│  Chromatography   │  HCP: 3-4 LRV, DNA: 4-5 LRV
└───────────────────┘
        │
        ▼
┌───────────────────┐
│  Low pH Viral     │  Viral inactivation
│   Inactivation    │
└───────────────────┘
        │
        ▼
┌───────────────────┐
│  Cation Exchange  │  Aggregate removal
│  (Bind-Elute)     │  Charge variant control
└───────────────────┘
        │
        ▼
┌───────────────────┐
│  Anion Exchange   │  DNA/HCP polishing
│  (Flow-Through)   │  Viral clearance
└───────────────────┘
        │
        ▼
┌───────────────────┐
│  Virus Filtration │  Final viral clearance
└───────────────────┘
        │
        ▼
┌───────────────────┐
│    UF/DF          │  Concentration
│                   │  Formulation
└───────────────────┘
        │
        ▼
      Drug Substance

Typical Performance Targets

Step	Yield	Purity	HCP	DNA
Protein A	>95%	>95%	100-500 ppm	<10 ppb
CEX	>90%	>98%	<50 ppm	<10 ppb
AEX	>95%	>99%	<10 ppm	<1 ppb
Final	>70% overall	>99.5%	<10 ppm	<10 ppb

Cost of goods

difflow_bio.economics.cogs_breakdown estimates the annual cost of goods of a mAb process step by step, and cost per gram released. The cost follows the process rather than the batch count alone:

Item	How it is computed
Resin	Each chromatography step has its own resin, column volume, capacity and lifetime. Cycles per batch are load / (DBC x CV), at least one, and the resin is replaced every resin_lifetime_cycles cycles.
Product mass	The harvest (working_volume_L x titer_g_L) runs through the steps in order, each keeping step_yield of what it receives, so a step’s load depends on everything upstream.
Buffers	Column volumes per cycle for chromatography, litres per m2 for membranes, litres per batch otherwise, each priced from buffer_usd_L.
Media	Basal medium, fed-batch feed and seed train, per litre of working volume.
Labor	Fixed support staff, plus operator hours per batch for the upstream train and for each downstream step.
Other	QC per batch, single-use items per step, utilities, liquid waste, maintenance and depreciation on CAPEX.
Failed batches	Batches are charged when started and only batch_success_rate of them are released.

Cycles per batch are continuous, 1 + softplus(k (n - 1)) / k with n = load / (DBC x CV), so cost has a derivative with respect to titer, binding capacity, column volume and every yield. jax.grad of cost per gram with respect to titer is therefore not just -cost / titer: a higher titer also loads the columns harder.

import dataclasses
import jax
from difflow_bio.economics import load_cost_model, cogs_breakdown

process, basis = load_cost_model()          # the shipped reference
out = cogs_breakdown(process, basis)
out["cost_per_g"], out["steps"]["protein_a_capture"]["cycles"]

d_titer = jax.grad(lambda t: cogs_breakdown(
    dataclasses.replace(process, titer_g_L=t), basis)["cost_per_g"])(5.0)

Using your own data

Every price, rate and process specification is data. load_cost_model(path) reads a YAML or JSON file with a process section (bioreactor, batches and the ordered steps, each with a type of chromatography, filtration or yield) and a basis section (prices and labor rates). The shipped difflow_bio/economics/data/mab_reference.yaml is the template: copy it, replace the numbers with your own plant’s, and load the copy. Every number in it is tagged with its source: [database] (copied from the resin database), [carried over] (the value the older functions used) or [placeholder] (an order-of-magnitude assumption). Most are placeholders, and the reference result is not a benchmark.

The benchmark is a published process encoded in the same format: difflow_bio/economics/data/petrides2015_mab.yaml, the large-scale mAb process of D. Petrides, Bioprocess Design and Economics (Intelligen, 2015, section 11.6.3; an improved version is Chapter 11 of Harrison, Todd, Rudge and Petrides, Bioseparations Science and Engineering, 2nd ed., Oxford University Press, 2015). Every number in it is tagged with the page that states it, or marked derived or not stated. From the source’s stated inputs (a 15,000 L fed-batch harvest at about 2 g/L, Protein A, IEX and HIC columns with their volumes, capacities, yields, buffer volumes, prices and lifetimes, 80 batches a year), cogs_breakdown reproduces:

Quantity	difflow	Source
Product per batch / per year	19.28 kg / 1,542 kg	19.3 kg / 1,544 kg
Cycles per batch, Protein A / IEX / HIC	3.85 / 3.01 / 3.00	4 / 3 / 3 (whole cycles)
Protein A / HIC elution buffer per batch	9,658 / 2,846 L	10,002 / 2,990 kg (Table 11.15)
Resin replacement	$18.5M/yr	within Consumables, $23.6M/yr (Table 11.17)
Media	$11.2M/yr	within Raw Materials, $16.7M/yr (Table 11.17)

The source gives labor, facility-dependent, QC and miscellaneous costs only as totals, without the inputs behind them, so its total of $84/g is not reproduced; those categories are zero in the benchmark file.

ChromatographyStep.from_resin builds a step from the resin database (capacity derated to 80 % of q_max, price and lifetime), and ProcessSpec.from_dict / CostBasis.from_dict build the rest from plain dictionaries; unknown keys and unpriced buffers are reported by name.

The older estimate_* functions remain for coarse estimates. estimate_total_opex now prices each chromatography step with its own resin (Protein A capture, then cation and anion exchange) and warns that it is deprecated in favour of cogs_breakdown.