Model

Dry Lab members initially generated mathematical models and computer simulations based purely on theory and reading. After the Wet Lab had collected real data, it became possible to compare theoretical expectations with their findings.

Once the Wet Lab's APET film degradation results came back (BHET, MHET and TPA measured by HPLC at 24, 48, 72 and 96 hours, across four ICCG conditions, i.e. dockerin fusion present or absent, crossed with scaffold present or absent), we fit a mechanistic model directly to that data rather than relying on the theoretical models below. The reaction network allows the enzyme to release BHET, MHET and TPA in parallel directly from the solid film (not only in strict sequence), lets BHET and MHET convert onward to the next product, and lets the film's effective release rate grow over time as hydrolysis exposes more surface area. The equations are given in full below the chart, and the code driving the chart is reproduced at the bottom of this section so it can be inspected line by line.

Solid lines are the model; dots are the measured HPLC concentrations at 24/48/72/96 h. Switch conditions with the buttons above the chart.

Measured data

Concentrations in µM, from HPLC quantification of the APET film degradation assay (5 mL reaction, 100 mU·mL⁻¹ ICCG, 6×6×0.5 mm APET film chip). “Dock+” = dockerin-fused ICCG–DoT; “Dock−” = unfused ICCG. “Scaf+” = co-incubated with 1 µg·mL⁻¹ ScafGVT.

Governing equations

Three coupled species, namely BHET (\(B\)), MHET (\(M\)) and TPA (\(T\)), are released from the solid film in parallel and also convert onward once in solution:

k_B0, k_M0, k_T0 are the rates at which BHET, MHET and TPA are released directly from the film. k₁ and k₂ are the secondary conversions once product is in solution: BHET→MHET and MHET→TPA.

The film does not present a constant surface to the enzyme; as hydrolysis proceeds, more chain ends and surface area become exposed, so all three release rates are scaled by a shared, saturating “erosion” term:

α sets how fast accessibility grows early on; β sets where it levels off. As β→0 this becomes unbounded linear growth; the data only pin down a finite β for the unfused (Dock−) conditions. See note below.

All three species start at zero (\(B(0)=M(0)=T(0)=0\)) and the seven parameters (\(k_{B0}, k_{M0}, k_{T0}, \alpha, \beta, k_1, k_2\)) were fit per condition by nonlinear least squares against the 24/48/72/96 h HPLC measurements, with each species' residuals normalised by its own maximum observed value so that BHET (numerically much smaller than MHET or TPA) is not swamped in the fit.

*α pinned at the 20 upper bound imposed during fitting for the two Dock+ conditions, and k₁ pinned at its lower bound of 0: neither curve shows any sign of levelling off within 96 h, so the data cannot distinguish “very fast, still-growing accessibility” from “unboundedly growing accessibility”; more timepoints beyond 96 h would be needed to pin this down. The Dock− conditions, by contrast, show a real, finite β because BHET visibly peaks and then declines by 96 h in those wells, which requires a saturating (not ever-growing) source term to reproduce.

Source code

The chart above is generated entirely client-side: a 4th-order Runge–Kutta integrator steps the three equations forward in small time increments.

// ---------------------------------------------------------------------
// PETosome parallel-release + erosion model, fit to APET film HPLC data
// (24/48/72/96 h). Three species: B = BHET, M = MHET, T = TPA.
//
//   dB/dt = kB0*s(t) - k1*B
//   dM/dt = kM0*s(t) + k1*B - k2*M
//   dT/dt = kT0*s(t) + k2*M
//   s(t)  = 1 + alpha*t/(1+beta*t)
//
// This uses a 4th-order Runge-Kutta stepper. 
// ---------------------------------------------------------------------

const PETOSOME_PARAMS = {
  "Dock+, Scaf+": {
    kB0: 0.0012, kM0: 0.0355, kT0: 0.0004, alpha: 20, beta: 0, k1: 0.0, k2: 0.01204,
    obs: { B: [0.28, 12.72, 62.67, 117.19], M: [5.19, 233.51, 1093.69, 2609.80], T: [1.37, 70.01, 401.57, 1052.61] }
  },
  "Dock+, Scaf-": {
    kB0: 0.0012, kM0: 0.0347, kT0: 0.0062, alpha: 20, beta: 0, k1: 0.0, k2: 0.00690,
    obs: { B: [5.69, 18.18, 57.51, 115.24], M: [79.38, 479.11, 1450.64, 2710.69], T: [19.89, 174.92, 574.71, 1195.90] }
  },
  "Dock-, Scaf+": {
    kB0: 0.0109, kM0: 0.2128, kT0: 0.0286, alpha: 20, beta: 0.0802, k1: 0.0279, k2: 0.00666,
    obs: { B: [15.69, 42.14, 83.46, 63.54], M: [301.15, 1289.08, 2725.47, 2596.61], T: [89.03, 402.69, 878.57, 1322.48] }
  },
  "Dock-, Scaf-": {
    kB0: 0.0117, kM0: 0.2429, kT0: 0.0502, alpha: 20, beta: 0.1034, k1: 0.0190, k2: 0.00553,
    obs: { B: [17.03, 52.23, 70.94, 80.96], M: [315.04, 1481.48, 2267.99, 2675.22], T: [93.30, 452.75, 973.19, 1362.90] }
  }
};

const PETOSOME_T_MAX = 96; // hours - the model only covers the real 24-96 h data window

function petosomeDeriv(t, y, p) {
  const surf = 1 + (p.alpha * t) / (1 + p.beta * t);
  const B = y[0], M = y[1];
  const dB = p.kB0 * surf - p.k1 * B;
  const dM = p.kM0 * surf + p.k1 * B - p.k2 * M;
  const dT = p.kT0 * surf + p.k2 * M;
  return [dB, dM, dT];
}

function petosomeIntegrate(p, tMax, steps) {
  const h = tMax / steps;
  let y = [0, 0, 0];
  let t = 0;
  const out = [{ t: 0, B: 0, M: 0, T: 0 }];
  for (let i = 0; i < steps; i++) {
    const k1v = petosomeDeriv(t, y, p);
    const y2 = y.map((v, j) => v + (h / 2) * k1v[j]);
    const k2v = petosomeDeriv(t + h / 2, y2, p);
    const y3 = y.map((v, j) => v + (h / 2) * k2v[j]);
    const k3v = petosomeDeriv(t + h / 2, y3, p);
    const y4 = y.map((v, j) => v + h * k3v[j]);
    const k4v = petosomeDeriv(t + h, y4, p);
    y = y.map((v, j) => v + (h / 6) * (k1v[j] + 2 * k2v[j] + 2 * k3v[j] + k4v[j]));
    t += h;
    out.push({ t: t, B: y[0], M: y[1], T: y[2] });
  }
  return out;
}

const PETOSOME_THIS_SCRIPT = document.currentScript;

document.addEventListener("DOMContentLoaded", function () {
  const canvas = document.getElementById("petosomeChart");
  if (!canvas) return;

  let currentCondition = "Dock+, Scaf+";

  const chart = new Chart(canvas, {
    data: { datasets: [] },
    options: {
      responsive: true,
      maintainAspectRatio: false,
      parsing: false,
      scales: {
        x: { type: "linear", min: 0, max: PETOSOME_T_MAX, title: { display: true, text: "Time (hours)" } },
        y: { title: { display: true, text: "Concentration (\u00B5M)" } }
      },
      plugins: {
        legend: { display: false },
        tooltip: { callbacks: { label: c => c.dataset.label + ": " + Math.round(c.parsed.y) + " \u00B5M" } }
      }
    }
  });

  function renderChart() {
    const p = PETOSOME_PARAMS[currentCondition];
    const sol = petosomeIntegrate(p, PETOSOME_T_MAX, 200);
    const line = (key, color) => ({
      label: key, type: "line", showLine: true, pointRadius: 0, borderWidth: 2, borderColor: color,
      data: sol.map(d => ({ x: d.t, y: d[key] }))
    });
    const dots = (key, color, obsArr) => ({
      label: key + " (measured)", type: "scatter", pointRadius: 4,
      pointBackgroundColor: color, pointBorderColor: color,
      data: [24, 48, 72, 96].map((t, i) => ({ x: t, y: obsArr[i] }))
    });
    chart.data.datasets = [
      line("B", "#d85a30"), line("M", "#2a78d6"), line("T", "#1b9e70"),
      dots("B", "#d85a30", p.obs.B), dots("M", "#2a78d6", p.obs.M), dots("T", "#1b9e70", p.obs.T)
    ];
    chart.update();
  }

  const btnRow = document.getElementById("condBtns");
  Object.keys(PETOSOME_PARAMS).forEach(function (name) {
    const btn = document.createElement("button");
    btn.type = "button";
    btn.textContent = name;
    btn.className = "btn btn-sm btn-outline-secondary";
    btn.dataset.name = name;
    btn.addEventListener("click", function () {
      currentCondition = name;
      updateButtons();
      renderChart();
    });
    btnRow.appendChild(btn);
  });

  function updateButtons() {
    Array.prototype.forEach.call(btnRow.children, function (btn) {
      if (btn.dataset.name === currentCondition) {
        btn.classList.remove("btn-outline-secondary");
        btn.classList.add("btn-secondary");
      } else {
        btn.classList.remove("btn-secondary");
        btn.classList.add("btn-outline-secondary");
      }
    });
  }

  updateButtons();
  renderChart();

  // Print this entire script into the visible source-code block at the
  // bottom of the section, so the model is auditable without view-source.
  const sourceEl = document.getElementById("petosome-source");
  if (sourceEl) {
    sourceEl.textContent = PETOSOME_THIS_SCRIPT ? PETOSOME_THIS_SCRIPT.textContent.trim() : "";
  }
});

Every simulation reviewed here tracks the same chemistry; PET broken down by PETase, its intermediate cleaved by MHETase into terephthalic acid (TPA) and ethylene glycol (EG), but each encodes a different mathematical assumption about how fast that happens. This section lays out the governing equations side by side for the five models retained after review.

Saturable (Michaelis–Menten) kinetics · Simulation X · baseline

The baseline model treats both enzymes as classic single-substrate catalysts. Reaction velocity saturates as substrate becomes abundant, because enzyme active sites are the limiting resource:

V_max is the top speed, i.e., the fastest the enzyme goes if you flood it with substrate. K_m is the substrate level where it's already running at half that speed. Small K_m means the enzyme is greedy; it doesn't need much substrate to get going.

Crystalline PET is harder to attack than amorphous PET, so velocity is split by phase and the crystalline fraction is damped exponentially by a crystallinity index \(x = S_{\text{cryst}}/S\):

x is just "how crystalline is the plastic," 0 to 1. The e^-2.5x term is what makes crystalline PET so much harder to digest.

In general, we expect chemical reactions to go faster when the temperature is high. However, excessively high proteins may damage the shape of the proteins. Based on a survey of previously published books, our first conjecture was that 30°C might be optimal. Temperature was modeled through an Arrhenius term centered on a 30°C optimum, with an extra penalty past 55°C:

This is "reactions speed up when it's warmer, until the enzyme cooks." The exp(...) part is the speed-up; the piecewise case after it is where things fall apart. For X, that happens past 55°C.

MHET cleavage splits into TPA and EG at a fixed, mass-conserving ratio (molar composition of MHET): \( \Delta\text{TPA} = 0.73\,v_{\text{MHETase}}, \; \Delta\text{EG} = 0.27\,v_{\text{MHETase}} \). Time is advanced in discrete ticks (an explicit Euler-style update), not continuous integration.

Cooperative (Hill) kinetics with a sharper thermal cliff · Simulation Z · variant

Here PETase activity is modeled as cooperative rather than simple Michaelis–Menten, using a Hill coefficient \(n = 1.8\); velocity rises more sigmoidally with substrate concentration:

n is the Hill coefficient. n=1 is plain Michaelis–Menten. n>1 means the enzyme gets more effective once it's already working.

MHETase kinetics stay Michaelis–Menten (identical to X). What changes most is thermal behavior: instead of X's gentle penalty above 55°C, Z applies a much steeper denaturation cliff starting ten degrees earlier:

Same exp(...) term as X, but the cutoff moved ten degrees earlier (45°C vs 55°C) and the fall-off is steeper (−0.8 vs −0.3).

Crystallinity hindrance is also linearized rather than exponential: \( e^{-2.5x} \to 1 - 0.95x \). Stoichiometry reverts to the correct 0.73/0.27 TPA/EG split, so unlike Y, Z conserves mass.

Continuous-time ODE integration · alternate_simulation.html

Rather than stepping forward in discrete ticks, this model expresses the same PET→MHET→TPA pathway as a coupled system of ordinary differential equations and integrates it continuously with an explicit Runge–Kutta solver (scipy.integrate.solve_ivp, RK45):

k₁ and k₂ aren't measured constants here; they function as sliders. RK45 just means "solve continuously," instead of nudging forward in fixed steps like X and Z do.

Here \(k_1, k_2\) (PETase/MHETase expression levels) are the free parameters a user drags on sliders; the interface also perturbs each by a small \(\delta\) and re-integrates to show a live finite-difference sensitivity, \( \partial Y/\partial k_i \approx \big[Y(k_i+\delta)-Y(k_i)\big]/\delta \).

Gene-circuit-coupled degradation · petase_circuit.html

This is the only model where enzyme concentration is not a fixed input. In this model, enzyme concentration is the output of an upstream transcription/translation cascade, driven by promoter strength \(p\), plasmid copy number \(c\), and cell density \(d\):

w₁ and w₂ are just relative production rates; PETase's mRNA gets weight 1, MHETase's gets 0.8, so MHETase is made a bit slower per unit of the same upstream signal.

PET breakdown itself is then zero-order in substrate and first-order only in enzyme concentration; there is no \(K_m\) term at all, unlike some other models here:

"Zero-order in substrate" means there's so much PET around that adding more doesn't speed anything up; the enzyme amount is the only thing that matters here, which is why K_m drops out entirely.

Molar-mass stoichiometric pathway · intermediate.html (v7) & 2033_intermediate.html (v5)

The most detailed model in the set: it tracks an explicit intermediate (BHET) between PET and MHET, converts every flux through real molar masses, and feeds temperature and pH back into the reaction rates rather than treating them as fixed constants:

f₁, f₂, f₃ are the three step rates (PET→BHET→MHET→TPA); each runs through the same φ(pH) and θ(T) modifiers; therefore, a bad pH slows the whole chain, not just one step.

Every Δ here gets multiplied by a real molar mass; the other models in this set consider abstract "concentration units."

Molar masses (g/mol): PET 192.17, BHET 254.24, MHET 210.18, TPA 166.13, EG 62.07.

Two feedback loops make this model distinctive: reaction extent raises temperature (\(\Delta T \mathrel{+}= 1.5f_1\), an exothermic-heat proxy), and acid byproduct lowers pH, which in turn slows \(\phi(\text{pH})\); a self-limiting reaction absent from every other model above. The v5 file (2033_intermediate.html) adds one further wrinkle: a crystallinity step-modifier, \( \text{crystMod} = 0.3 \) once \([\text{PET}] < 0.1\times\text{PET}_0\), representing a recalcitrant crystalline core that resists the final stage of digestion.

At a glance

In short: X and Z form a matched pair testing how one modeling choice (saturation law, thermal response, crystallinity penalty) changes outcomes on an otherwise identical scaffold. X assumes standard Michaelis–Menten kinetics, Z assumes cooperative (Hill) kinetics, reflecting the two candidate assumptions considered before the team had real experimental data to test between them. Source simulations: SimulationX/Z (24 May 2026), alternate_simulation.html (13 May 2026), petase_circuit.html & intermediate.html/2033_intermediate.html (12 Apr 2026).

Simulation Y, an early linearized-kinetics variant, was removed from this comparison after review; it approximated saturating kinetics as first-order (valid only when substrate concentration is much smaller than K_m, which did not hold here) and its stoichiometric bookkeeping did not conserve total mass.

These five models offered five different guesses at the same reaction, but were made obsolete by the arrival of the actual lab data.