+ Parameter Search estimates kinetic parameters for microbial growth
+ and substrate consumption models by fitting time-course
+ measurements of biomass and residual substrate.
+
+
+ Fill the experimental table with the sampling times (in hours) and
+ the measured concentrations (g/L). Use the algorithm controls to
+ tune the particle swarm optimization (PSO) strategy—adjust the
+ swarm size, cognitive and social weights, inertia factor, and
+ iteration limit to reflect the variability of your dataset. Model
+ bounds let you constrain feasible ranges for each kinetic
+ parameter before running the search.
+
+
+ After clicking Run search, the tool evaluates all
+ supported kinetic expressions, displays the PSO progress, and plots
+ the simulated curves against your data. Each card lists the
+ optimized coefficients together with the maintenance term derived
+ from Pirt (1965), so you can compare how
+ different formulations reproduce the experiment without leaving the
+ page.
+
+
+
+
+
+
+
+
+
+
diff --git a/www/package.json b/www/package.json
new file mode 100644
index 0000000..f99c1cc
--- /dev/null
+++ b/www/package.json
@@ -0,0 +1,7 @@
+{
+ "name": "biolab",
+ "type": "module",
+ "scripts": {
+ "test": "node tests/demo.mjs && node tests/synthetic-search.mjs"
+ }
+}
diff --git a/www/src/Objective.js b/www/src/Objective.js
new file mode 100644
index 0000000..514b0d9
--- /dev/null
+++ b/www/src/Objective.js
@@ -0,0 +1,131 @@
+import { RK4, RK4getvalue } from "./runge-kutta.js";
+export class Objective {
+ constructor(exper, f, res, mSize) {
+ /*
+ * res: a resolução utilizada para a resolução com RK4
+ * mSize: numero de variáveis do modelo atual
+ * exper: tabela de valores obtidos experimentalmente
+ */
+ this.exper = exper;
+ this.sampleSize = exper.length;
+ // Tamanho da amostra obtida experimentalmente (número de pontos)
+ this.f = f;
+ // Função f(t) a ser utilizada para o ajuste
+ // Vetor de tempos para resolução numérica
+ let tf = exper[this.sampleSize - 1][0];
+ this.timeArray = [];
+ for (let i = 0; i <= res; i++) {
+ this.timeArray[i] = (i * tf) / res;
+ }
+ // valor inicial é o primeiro ponto experimental
+ this.y0 = [exper[1][2], exper[1][1]];
+
+ // Se o número de variaveis na tabela (descontando a coluna tempo)
+ // for maior que o numero de variáveis do modelo,
+ // a função objetivo é calculada com o número de variáveis do modelo.
+ // Caso contrário, ela é calculada com o numero de variáveis da tabela.
+ if (exper[1].length - 1 >= mSize) {
+ this.nVar = mSize;
+ } else {
+ this.nVar = exper[1].length - 1;
+ }
+
+ const headerRow = Array.isArray(exper[0]) ? exper[0] : [];
+ const normalizedHeaders = headerRow.map((value) => {
+ if (typeof value !== "string") {
+ return "";
+ }
+ return value
+ .normalize("NFD")
+ .replace(/[\u0300-\u036f]/g, "")
+ .toLowerCase();
+ });
+
+ const dataColumns = exper[1] ? exper[1].length : 0;
+ const maxValidIndex = Math.max(1, dataColumns - 1);
+ const clampColumnIndex = (candidate) => {
+ if (!Number.isFinite(candidate) || candidate < 1 || candidate >= dataColumns) {
+ return Math.min(Math.max(1, candidate || 1), maxValidIndex);
+ }
+ return candidate;
+ };
+
+ const findColumnIndex = (keywords, fallback) => {
+ const idx = normalizedHeaders.findIndex((name, columnIndex) => {
+ if (columnIndex === 0) {
+ return false;
+ }
+ return keywords.some((keyword) => name.includes(keyword));
+ });
+ const candidate = idx > 0 ? idx : fallback;
+ return clampColumnIndex(candidate);
+ };
+
+ this.stateColumnIndex = [];
+ if (this.nVar >= 1) {
+ const fallbackCells = Math.min(2, maxValidIndex);
+ this.stateColumnIndex[0] = findColumnIndex(
+ ["celula", "celulas", "celula", "cell", "cells", "biomassa", "biomass", "x"],
+ fallbackCells,
+ );
+ }
+ if (this.nVar >= 2) {
+ this.stateColumnIndex[1] = findColumnIndex([
+ "substrato",
+ "substrate",
+ "substrat",
+ "s",
+ ], 1);
+ }
+ for (let j = 2; j < this.nVar; j++) {
+ this.stateColumnIndex[j] = clampColumnIndex(Math.min(j + 1, maxValidIndex));
+ }
+
+ // Média das colunas para normalização
+ this.mean = new Array(this.nVar).fill(0);
+ for (let i = 1; i < this.sampleSize; i++) {
+ for (let j = 0; j < this.nVar; j++) {
+ const columnIndex = this.stateColumnIndex[j];
+ if (columnIndex < exper[i].length) {
+ this.mean[j] += exper[i][columnIndex];
+ }
+ }
+ }
+ for (let j = 0; j < this.nVar; j++) {
+ this.mean[j] /= this.sampleSize - 1;
+ }
+ }
+ objective(params) {
+ // Resolve o modelo pelo método RK4
+ let sol = RK4(this.f, this.timeArray, this.y0, params);
+
+ // Pega predições pontuais do modelo
+ let pSol = [];
+ for (let i = 1; i < this.sampleSize; i++) {
+ pSol[i - 1] = RK4getvalue(
+ sol,
+ this.timeArray,
+ this.exper[i][0],
+ this.f,
+ params,
+ );
+ }
+
+ // Cálculo da função objetivo
+ let obj = 0;
+ for (let i = 2; i < this.sampleSize; i++) {
+ for (let j = 0; j < this.nVar; j++) {
+ const columnIndex = this.stateColumnIndex[j];
+ if (columnIndex >= this.exper[i].length) {
+ continue;
+ }
+ const observed = this.exper[i][columnIndex];
+ const predicted = pSol[i - 1][j];
+ const denom = this.mean[j] !== 0 ? this.mean[j] : 1;
+ obj += ((predicted - observed) / denom) ** 2;
+ }
+ }
+ // console.log(obj)
+ return obj;
+ }
+}
diff --git a/www/src/PSO.js b/www/src/PSO.js
new file mode 100644
index 0000000..a81ce14
--- /dev/null
+++ b/www/src/PSO.js
@@ -0,0 +1,93 @@
+import { rrandom } from "./rrandom.js";
+export class PSO {
+ constructor(obj, n, bounds) {
+ /*
+ obj: função objetivo a ser minimizada
+ n: numero de iterações a serem performadas
+ bounds: vetor com os mínimos e maximos dos parametros da busca
+ */
+ this.obj = obj;
+ this.n = n;
+ this.bounds = bounds;
+ this.pos = [];
+ this.vel = [];
+ this.err = [];
+ this.pos_best = [];
+ this.err_best = [];
+
+ this.dimensions = bounds.length; // número de dimensões da busca (parâmetros)
+
+ for (var i = 0; i < this.n; i++) {
+ /*
+ inicia o vetor vel com valores aleatórios entre -1 e 1
+ inicia o vetor pos com posições aleatórias dentro dos limites especificados
+ */
+ var v = [];
+ var x = [];
+ for (var j = 0; j < this.dimensions; j++) {
+ v.push(rrandom(-1, 1));
+ x.push(rrandom(bounds[j][0], bounds[j][1]));
+ }
+
+ this.vel.push(v.slice());
+ this.pos.push(x.slice());
+ this.pos_best.push(x.slice()); // inicia a melhor posição como a posiçaõ atual (aleatória)
+ this.err_best.push(Infinity);
+ }
+
+ this.err_best_g = Infinity;
+ // inicia o melhor erro global como infinito (quanto menor melhor)
+ this.pos_best_g = this.pos[0].slice();
+ // inicia a melhor posição global como a posição da particula 0
+ }
+
+ run(c1, c2, w, iteration) {
+ for (let i = 0; i < iteration; i++) {
+ this.update(c1, c2, w);
+ // console.log(i, "/", iteration, "\t", this.err_best_g);
+ }
+
+ console.log(this.pos_best_g);
+ // console.log(this.err_best_g)
+ // console.table(this.pos)
+ }
+
+ update(c1, c2, w) {
+ for (let i = 0; i < this.n; i++) {
+ for (let j = 0; j < this.dimensions; j++) {
+ let r1 = rrandom(0, 1);
+ let r2 = rrandom(0, 1);
+
+ let vel_cognitive = c1 * r1 * (this.pos_best[i][j] - this.pos[i][j]);
+ let vel_social = c2 * r2 * (this.pos_best_g[j] - this.pos[i][j]);
+ const updatedVelocity = w * this.vel[i][j] + vel_cognitive + vel_social;
+ let candidatePosition = this.pos[i][j] + updatedVelocity;
+
+ if (candidatePosition > this.bounds[j][1]) {
+ candidatePosition = this.bounds[j][1];
+ this.vel[i][j] = 0;
+ } else if (candidatePosition < this.bounds[j][0]) {
+ candidatePosition = this.bounds[j][0];
+ this.vel[i][j] = 0;
+ } else {
+ this.vel[i][j] = updatedVelocity;
+ }
+
+ this.pos[i][j] = candidatePosition;
+ }
+
+ // Update error value
+ this.err[i] = this.obj.objective(this.pos[i]);
+
+ // Update global values
+ if (this.err[i] < this.err_best[i]) {
+ this.pos_best[i] = this.pos[i].slice();
+ this.err_best[i] = this.err[i];
+ if (this.err[i] < this.err_best_g) {
+ this.pos_best_g = this.pos[i].slice();
+ this.err_best_g = this.err[i];
+ }
+ }
+ }
+ }
+}
diff --git a/www/src/conhecidos.js b/www/src/conhecidos.js
new file mode 100644
index 0000000..efd45d7
--- /dev/null
+++ b/www/src/conhecidos.js
@@ -0,0 +1,49 @@
+// MONOD, 1942
+export function monod(S, mu_max, K_S) {
+ return (mu_max * S) / (K_S + S);
+}
+
+// MOSER, 1958
+export function moser(S, mu_max, K_S, n) {
+ return (mu_max * S ** n) / (K_S + S ** n);
+}
+
+// CONTOIS, 1959
+export function contois(S, X, mu_max, K_S) {
+ return (mu_max * S) / (K_S * X + S);
+}
+
+// BERGTER, 1983
+export function bergter(S, t, mu_max, K_S, T) {
+ return ((mu_max * S) / (K_S + S)) * (1 - Math.exp(-t / T));
+}
+
+// TESSIER, 1942
+export function tessier(S, mu_max, K_S) {
+ return mu_max * (1 - Math.exp(-S / K_S));
+}
+
+// ANDREWS, 1968
+export function andrews(S, mu_max, K_S, K_I) {
+ return (mu_max * S) / (K_S + S + (S ** 2) / K_I);
+}
+
+// AIBA; SHODA; NAGATANI, 1968
+export function aiba(S, mu_max, K_S, K_I) {
+ return ((mu_max * S) / (K_S + S)) * Math.exp(-K_I * S);
+}
+
+// Morte celular
+
+// SINCLAIR; KRISTIANSEN, 1987
+export function death(K_d) {
+ return -K_d;
+}
+
+// Consumo do substrato limitante para manutenção
+
+// PIRT, 1965
+
+export function pirt(mu, Y_XS, m_S) {
+ return (1 / Y_XS) * mu + m_S;
+}
diff --git a/www/src/consulta/model.js b/www/src/consulta/model.js
new file mode 100644
index 0000000..5ecb9fa
--- /dev/null
+++ b/www/src/consulta/model.js
@@ -0,0 +1,13 @@
+function model(time, y, params) {
+ let K_S = params[0];
+ let mu_max = params[1];
+ let m_S = params[2];
+ let Y_XS = params[3];
+ let X = y[0];
+ let S = y[1];
+ let dmu = monod(S, mu_max, K_S);
+ let dX = X * dmu;
+ const qS = pirt(dmu, Y_XS, m_S);
+ let dS = -qS * X;
+ return [dX, dS];
+}
diff --git a/www/src/consulta/modelos.js b/www/src/consulta/modelos.js
new file mode 100644
index 0000000..97a09cf
--- /dev/null
+++ b/www/src/consulta/modelos.js
@@ -0,0 +1,170 @@
+// Crescimento num único substrato limitante:
+
+// MONOD, 1942
+
+class monod {
+ constructor() {
+ this.nPar = 2;
+ this.parName = ["mu_max", "K_S"]
+ this.varName = ["Substrato"]
+ this.varUsed = [1, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ this.equation = "\\mu_{X} = \\frac{\\mu_{max} S}{Ks + S}"
+ this.name = "MONOD, 1942"
+ }
+ model(variavel, parametros) {
+ return (parametros[0] * variavel[0]) / (parametros[1] + variavel[0])
+ }
+}
+
+// MOSER, 1958
+class moser {
+ constructor() {
+ this.nPar = 3;
+ this.parName = ["mu_max", "K_S", "n"]
+ this.varName = ["Substrato"]
+ this.varUsed = [1, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ this.equation = "\\mu_{X} = \\frac{\\mu_{max} S^n}{Ks + S^n}"
+ this.name = "MOSER, 1958"
+ }
+ model(variavel, parametros) {
+ // moser(S, mu_max, K_S, n) {
+ return (parametros[0] * variavel[0] ** parametros[2]) / (parametros[1] + variavel[0] ** parametros[2])
+ }
+}
+
+// CONTOIS, 1959
+class contois {
+ constructor() {
+ this.nPar = 2;
+ this.parName = ["mu_max", "K_S"]
+ this.varName = ["Substrato", "Células"]
+ this.varUsed = [1, 0, 1, 0, 0] // sub1, sub2, cel, prod, tempo
+ }
+ model(variavel, parametros) {
+ // model(S, X, mu_max, K_S) {
+ return (parametros[0] * variavel[0]) / (parametros[1] * variavel[1] + variavel[0])
+ }
+}
+
+// BERGTER, 1983
+class bergter {
+ constructor() {
+ this.nPar = 3;
+ this.parName = ["mu_max", "K_S", "T"]
+ this.varName = ["Substrato", "tempo"]
+ this.varUsed = [1, 0, 0, 0, 1] // sub1, sub2, cel, prod, tempo
+ }
+ model(variavel, parametros) {
+ // (S, t, mu_max, K_S, T) {
+ return ((parametros[0] * variavel[0]) / (parametros[1] + variavel[0])) * (1 - Math.exp(-variavel[1] / parametros[2]))
+ }
+}
+
+// Morte celular
+
+// SINCLAIR; KRISTIANSEN, 1987
+class death {
+ constructor() {
+ this.nPar = 1;
+ this.parName = ["k_d"]
+ this.varName = []
+ this.varUsed = [0, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ }
+ model(parametros) {
+ return -parametros[0]
+ }
+}
+
+// Crescimento em um único substrato limitante e inibidor
+
+// AIBA; SHODA; NAGATANI; 1968
+class asn {
+ constructor() {
+ this.nPar = 3
+ this.parName = ["mu_max", "K_S", "K_i"]
+ this.varName = ["Substrato"]
+ this.varUsed = [1, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ }
+ model(variavel, parametros) {
+ return ((parametros[0] * variavel[0]) / (parametros[1] + variavel[0])) * Math.exp(-variavel[0] / parametros[2])
+ }
+}
+
+// HALDANE, 1930
+class haldane {
+ constructor() {
+ this.nPar = 3
+ this.parName = ["mu_xa", "K_S", "K_i"]
+ this.varName = ["Substrato"]
+ this.varUsed = [1, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ }
+ model(variavel, parametros) {
+ return (parametros[0] * variavel[0]) / (1 + parametros[1] / variavel[0] + variavel[0] / parametros[2])
+ }
+}
+
+// ANDREWS, 1968
+class andrews {
+ constructor() {
+ this.nPar = 3
+ this.parName = ["mu_xa", "K_S", "K_i"]
+ this.varName = ["Substrato"]
+ this.varUsed = [1, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ }
+ model(variavel, parametros) {
+ return (parametros[0] * variavel[0]) / (variavel[0] + parametros[1] + (variavel[0] ** 2 / parametros[2]))
+ }
+}
+
+// EDWARDS, 1970
+class edwards {
+ constructor() {
+ this.nPar = 3
+ this.parName = ["mu_xa", "K_S", "K_i"]
+ this.varName = ["Substrato"]
+ this.varUsed = [1, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ }
+ model(variavel, parametros) {
+ return (parametros[0] * variavel[0]) / (variavel[0] + parametros[1] + (variavel[0] ** 2 / parametros[2]) + (1 + variavel[0] / parametros[1]))
+ }
+}
+
+// WEBB, 1963
+class webb {
+ constructor() {
+ this.nPar = 3
+ this.parName = ["mu_xa", "K_S", "K_i"]
+ this.varName = ["Substrato"]
+ this.varUsed = [1, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ }
+ model(variavel, parametros) {
+ return (parametros[0] * variavel[0] * (1 + variavel[0] / parametros[2])) / (variavel[0] + parametros[1] + (variavel[0] ** 2 / parametros[2]))
+ }
+}
+
+// TEISSIER, 1942
+class teissier {
+ constructor() {
+ this.nPar = 3
+ this.parName = ["mu_max", "K_S", "K_i"]
+ this.varName = ["Substrato"]
+ this.varUsed = [1, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ }
+ model(variavel, parametros) {
+ return parametros[0] * (Math.exp(-variavel[0] / parametros[1]) - Math.exp(-variavel[0] / parametros[2]))
+ }
+}
+
+// WU et al., 1988
+class wu {
+ constructor() {
+ this.nPar = 4
+ this.parName = ["mu_xa", "K_S", "K_i", "n"]
+ this.varName = ["Substrato"]
+ this.varUsed = [1, 0, 0, 0, 0] // sub1, sub2, cel, prod, tempo
+ }
+ model(variavel, parametros) {
+ // return mu_xa / (1 + K_S / S + (S / K_i) ** n)
+ return parametros[0] / (1 + parametros[1] / variavel[0] + (variavel[0] / parametros[2]) ** parametros[3])
+ }
+}
diff --git a/www/src/jquery.jslatex.js b/www/src/jquery.jslatex.js
new file mode 100644
index 0000000..e04d149
--- /dev/null
+++ b/www/src/jquery.jslatex.js
@@ -0,0 +1,69 @@
+/*
+ * jsLaTeX v1.2 - jQuery plugin
+ *
+ * Copyright (c) 2009 Andreas Grech
+ *
+ * Dual licensed under the MIT and GPL licenses:
+ * http://www.opensource.org/licenses/mit-license.php
+ * http://www.gnu.org/licenses/gpl.html
+ *
+ * http://knowledge-aholic.blogspot.com
+ */
+
+(function ($) {
+ var attachToImage = function () {
+ return $("").attr({
+ src: this.src
+ });
+ },
+ formats = {
+ 'gif': attachToImage,
+ 'png': attachToImage,
+ 'swf': function () {
+ return $("").attr({
+ src: this.src,
+ type: 'application/x-shockwave-flash'
+ });
+ }
+ },
+ sections = {
+ '{f}': 'format',
+ '{e}': 'equation'
+ },
+ escapes = {
+ '+': '2B',
+ '=': '3D'
+ };
+
+ $.fn.latex = function (opts) {
+ opts = $.extend({},
+ $.fn.latex.defaults, opts);
+ opts.format = formats[opts.format] ? opts.format : 'gif';
+ return this.each(function () {
+ var $this = $(this),
+ format, s, element, url = opts.url;
+ opts.equation = $.trim($this.text());
+ for (s in sections) {
+ if (sections.hasOwnProperty(s) && (format = url.indexOf(s)) >= 0) {
+ url = url.replace(s, opts[sections[s]]);
+ }
+ }
+ for (s in escapes) {
+ if (escapes.hasOwnProperty(s) && (format = url.indexOf(s)) >= 0) {
+ url = url.replace(s, '%' + escapes[s]);
+ }
+ }
+ opts.src = url;
+ element = formats[opts.format].call(opts);
+ $this.html('').append(element);
+ if (opts.callback) {
+ opts.callback.call(element);
+ }
+ });
+ };
+
+ $.fn.latex.defaults = {
+ format: 'gif',
+ url: 'https://latex.codecogs.com/{f}.latex?{e}'
+ };
+}(jQuery));
diff --git a/www/src/rrandom.js b/www/src/rrandom.js
new file mode 100644
index 0000000..d71318f
--- /dev/null
+++ b/www/src/rrandom.js
@@ -0,0 +1,4 @@
+export function rrandom(min, max) {
+ // cria um número aleatório dentro do intervalo
+ return Math.random() * (max - min) + min;
+}
diff --git a/www/src/runge-kutta.js b/www/src/runge-kutta.js
new file mode 100644
index 0000000..3ff81bd
--- /dev/null
+++ b/www/src/runge-kutta.js
@@ -0,0 +1,70 @@
+export function RK4(f, t, Y0, params) {
+ /*
+ * f: trata-se da função f(t) ser integrada
+ * t: trata-se do vetor de pontos no tempo
+ * Y0: trata-se do vetor de condições iniciais
+ * params: paremetros passados para o modelo de crescimento
+ */
+
+ let resolution = t.length; // resolução é o numero de pontos amostrados
+ let h = t[1]; // h é o tamanho do passo de integração
+
+ let y = []; // inicializa o vetor que será devolvido como eixo y
+
+ y[0] = Y0; //inicializa o vetor com as condições iniciais
+
+ for (let i = 0; i < resolution; i++) {
+ y[i + 1] = RK4step(f, t[i], y[i], h, params);
+ }
+
+ return y;
+}
+
+function RK4step(f, t, y0, h, params) {
+ // Algorítimo Runge-kutta 4
+
+ let s = [];
+ let y = [];
+ let nVar = y0.length;
+
+ const k1 = f(t, y0, params);
+ for (let i = 0; i < nVar; i++) {
+ s[i] = y0[i] + (k1[i] * h) / 2;
+ }
+ const k2 = f(t + h / 2, s, params);
+
+ for (let i = 0; i < nVar; i++) {
+ s[i] = y0[i] + (k2[i] * h) / 2;
+ }
+ const k3 = f(t + h / 2, s, params);
+
+ for (let i = 0; i < nVar; i++) {
+ s[i] = y0[i] + k3[i] * h;
+ }
+ const k4 = f(t + h, s, params);
+
+ for (let i = 0; i < nVar; i++) {
+ s[i] = y0[i] + k3[i] * h;
+ y[i] = y0[i] + (k1[i] / 6 + k2[i] / 3 + k3[i] / 3 + k4[i] / 6) * h;
+ }
+
+ return y;
+}
+
+export function RK4getvalue(sol, timearray, time, f, params) {
+ /* A partir de uma solução gerada pela função RK4,
+ * esta função permite pegar o valor de y para um ponto arbitrário de tempo,
+ * para tal, é pego o valor mais próximo do ponto desejado
+ * e calculado com um passo de Runge-kutta para o valor exato de tempo desejado
+ */
+
+ let stepsize = timearray[1];
+ // o tamanho do passo é igual o ponto 1 do vetor tempo (ponto 0 é igual a 0)
+ let h = time % stepsize;
+ // o passo dado é igual ao modulo entre o tempo arbitrário e o passo na resolução original
+ let y0 = sol[Math.floor(time / stepsize)];
+ // inicia-se no ponto da curva imediatamente anterior ao ponto desejado
+ let y = RK4step(f, time, y0, h, params);
+ // calcula-se o y com um único passo RK4 a partir do ponto anterior
+ return y;
+}
diff --git a/www/src/search.js b/www/src/search.js
new file mode 100644
index 0000000..d57c1ad
--- /dev/null
+++ b/www/src/search.js
@@ -0,0 +1,559 @@
+import { PSO } from "./PSO.js";
+import { RK4, RK4getvalue } from "./runge-kutta.js";
+import { Objective } from "./Objective.js";
+import {
+ monod,
+ moser,
+ contois,
+ bergter,
+ tessier,
+ andrews,
+ aiba,
+ pirt,
+} from "./conhecidos.js";
+import "https://cdn.plot.ly/plotly-2.29.1.min.js";
+
+const parameterCatalog = {
+ K_S: {
+ latex: "K_S",
+ unitText: "g/L",
+ unitLatex: "\\mathrm{g\\,L^{-1}}",
+ overrides: {
+ contois: {
+ unitText: "g_S/g_X",
+ unitLatex: "\\frac{\\mathrm{g}_{S}}{\\mathrm{g}_{X}}",
+ },
+ },
+ },
+ mu_max: {
+ latex: "\\mu_{max}",
+ unitText: "h⁻¹",
+ unitLatex: "\\mathrm{h^{-1}}",
+ },
+ K_I: {
+ latex: "K_I",
+ unitText: "g/L",
+ unitLatex: "\\mathrm{g\\,L^{-1}}",
+ overrides: {
+ aiba: {
+ unitText: "L/g",
+ unitLatex: "\\mathrm{L\\,g^{-1}}",
+ },
+ },
+ },
+ m_S: {
+ latex: "m_S",
+ unitText: "g_S/(g_X·h)",
+ unitLatex: "\\frac{\\mathrm{g}_{S}}{\\mathrm{g}_{X}\\,\\mathrm{h}}",
+ },
+ Y_XS: {
+ latex: "Y_{XS}",
+ unitText: "g_X/g_S",
+ unitLatex: "\\frac{\\mathrm{g}_{X}}{\\mathrm{g}_{S}}",
+ },
+ T: {
+ latex: "T",
+ unitText: "h",
+ unitLatex: "\\mathrm{h}",
+ },
+ n: {
+ latex: "n",
+ unitText: null,
+ unitLatex: null,
+ },
+};
+
+const modelParameters = {
+ aiba: [
+ { key: "K_S", bounds: [0.005, 2] },
+ { key: "mu_max", bounds: [0.05, 0.9] },
+ { key: "K_I", bounds: [0.01, 1] },
+ { key: "m_S", bounds: [0.0015, 0.05] },
+ { key: "Y_XS", bounds: [0.3, 0.7] },
+ ],
+ andrews: [
+ { key: "K_S", bounds: [0.005, 2] },
+ { key: "mu_max", bounds: [0.05, 0.9] },
+ { key: "K_I", bounds: [5, 150] },
+ { key: "m_S", bounds: [0.0015, 0.05] },
+ { key: "Y_XS", bounds: [0.3, 0.7] },
+ ],
+ bergter: [
+ { key: "K_S", bounds: [0.005, 2] },
+ { key: "mu_max", bounds: [0.05, 0.9] },
+ { key: "T", bounds: [5, 80] },
+ { key: "m_S", bounds: [0.0015, 0.05] },
+ { key: "Y_XS", bounds: [0.3, 0.7] },
+ ],
+ contois: [
+ { key: "K_S", bounds: [0.005, 2] },
+ { key: "mu_max", bounds: [0.05, 0.9] },
+ { key: "m_S", bounds: [0.0015, 0.05] },
+ { key: "Y_XS", bounds: [0.3, 0.7] },
+ ],
+ monod: [
+ { key: "K_S", bounds: [0.005, 2] },
+ { key: "mu_max", bounds: [0.05, 0.9] },
+ { key: "m_S", bounds: [0.0015, 0.05] },
+ { key: "Y_XS", bounds: [0.3, 0.7] },
+ ],
+ moser: [
+ { key: "K_S", bounds: [0.005, 2] },
+ { key: "mu_max", bounds: [0.05, 0.9] },
+ { key: "n", bounds: [0.8, 2.5] },
+ { key: "m_S", bounds: [0.0015, 0.05] },
+ { key: "Y_XS", bounds: [0.3, 0.7] },
+ ],
+ tessier: [
+ { key: "K_S", bounds: [0.005, 2] },
+ { key: "mu_max", bounds: [0.2, 0.9] },
+ { key: "m_S", bounds: [0.005, 0.05] },
+ { key: "Y_XS", bounds: [0.3, 0.7] },
+ ],
+};
+
+function getParamDisplayInfo(paramKey, modelKey) {
+ const baseInfo = parameterCatalog[paramKey];
+ if (!baseInfo) {
+ throw new Error(`Unknown parameter: ${paramKey}`);
+ }
+ const override = baseInfo.overrides?.[modelKey];
+ if (!override) {
+ return baseInfo;
+ }
+ return { ...baseInfo, ...override };
+}
+
+function getDefaultBounds(modelKey) {
+ return modelParameters[modelKey].map((param) => param.bounds.slice());
+}
+
+function getParamDetails(modelKey) {
+ return modelParameters[modelKey].map((param) => {
+ const { latex, unitText, unitLatex } = getParamDisplayInfo(
+ param.key,
+ modelKey,
+ );
+ return {
+ key: param.key,
+ latex,
+ unitText,
+ unitLatex,
+ };
+ });
+}
+
+function monodPirt(_, y, params) {
+ let K_S = params[0];
+ let mu_max = params[1];
+ //pirt:
+ let m_S = params[2];
+ let Y_XS = params[3];
+ let X = y[0];
+ let S = y[1];
+ let dmu = monod(S, mu_max, K_S);
+ let dX = X * dmu;
+ const qS = pirt(dmu, Y_XS, m_S);
+ let dS = -qS * X;
+ return [dX, dS];
+}
+
+function moserPirt(_, y, params) {
+ let K_S = params[0];
+ let mu_max = params[1];
+ let n = params[2];
+ //pirt:
+ let m_S = params[3];
+ let Y_XS = params[4];
+ let X = y[0];
+ let S = y[1];
+ let dmu = moser(S, mu_max, K_S, n);
+ let dX = X * dmu;
+ const qS = pirt(dmu, Y_XS, m_S);
+ let dS = -qS * X;
+ return [dX, dS];
+}
+
+function contoisPirt(_, y, params) {
+ let K_S = params[0];
+ let mu_max = params[1];
+ let m_S = params[2];
+ let Y_XS = params[3];
+ let X = y[0];
+ let S = y[1];
+ let dmu = contois(S, X, mu_max, K_S);
+ let dX = X * dmu;
+ const qS = pirt(dmu, Y_XS, m_S);
+ let dS = -qS * X;
+ return [dX, dS];
+}
+
+function bergterPirt(time, y, params) {
+ let K_S = params[0];
+ let mu_max = params[1];
+ let T = params[2];
+ // pirt:
+ let m_S = params[3];
+ let Y_XS = params[4];
+ let X = y[0];
+ let S = y[1];
+ let dmu = bergter(S, time, mu_max, K_S, T);
+ let dX = X * dmu;
+ const qS = pirt(dmu, Y_XS, m_S);
+ let dS = -qS * X;
+ return [dX, dS];
+}
+
+function tessierPirt(_, y, params) {
+ let K_S = params[0];
+ let mu_max = params[1];
+ // pirt:
+ let m_S = params[2];
+ let Y_XS = params[3];
+ let X = y[0];
+ let S = y[1];
+ let dmu = tessier(S, mu_max, K_S);
+ let dX = X * dmu;
+ const qS = pirt(dmu, Y_XS, m_S);
+ let dS = -qS * X;
+ return [dX, dS];
+}
+
+function andrewsPirt(_, y, params) {
+ let K_S = params[0];
+ let mu_max = params[1];
+ let K_I = params[2];
+ // pirt:
+ let m_S = params[3];
+ let Y_XS = params[4];
+ let X = y[0];
+ let S = y[1];
+ let dmu = andrews(S, mu_max, K_S, K_I);
+ let dX = X * dmu;
+ const qS = pirt(dmu, Y_XS, m_S);
+ let dS = -qS * X;
+ return [dX, dS];
+}
+
+function aibaPirt(_, y, params) {
+ let K_S = params[0];
+ let mu_max = params[1];
+ let K_I = params[2];
+ // pirt:
+ let m_S = params[3];
+ let Y_XS = params[4];
+ let X = y[0];
+ let S = y[1];
+ let dmu = aiba(S, mu_max, K_S, K_I);
+ let dX = X * dmu;
+ const qS = pirt(dmu, Y_XS, m_S);
+ let dS = -qS * X;
+ return [dX, dS];
+}
+
+export async function main(input, options = {}) {
+ const alg = options.alg || {
+ particles: 50,
+ c1: 1.49618,
+ c2: 1.49618,
+ w: 0.7298,
+ iterations: 150,
+ };
+ const boundsOpt = options.bounds || {};
+ const onProgress = options.onProgress || (() => {});
+
+ const time = [];
+ const subs = [];
+ const cels = [];
+ for (let i = 1; i < input.length; i++) {
+ time[i] = input[i][0];
+ subs[i] = input[i][1];
+ cels[i] = input[i][2];
+ }
+
+ const basePlot = [
+ {
+ x: time,
+ y: subs,
+ name: "Experimental substrate",
+ mode: "markers",
+ marker: { size: 10, color: "#4a90e2" },
+ type: "scatter",
+ },
+ {
+ x: time,
+ y: cels,
+ name: "Experimental cells",
+ mode: "markers",
+ marker: { size: 10, color: "#50e3c2" },
+ type: "scatter",
+ },
+ ];
+
+ function calculateAIC(obj, sol, params) {
+ let sse = 0;
+ let count = 0;
+
+ for (let i = 1; i < input.length; i++) {
+ const timePoint = input[i][0];
+ const prediction = RK4getvalue(
+ sol,
+ obj.timeArray,
+ timePoint,
+ obj.f,
+ params,
+ );
+
+ const cellResidual = prediction[0] - input[i][2];
+ const substrateResidual = prediction[1] - input[i][1];
+
+ sse += cellResidual ** 2 + substrateResidual ** 2;
+ count += 2;
+ }
+
+ const n = count;
+ const k = params.length;
+ const meanSquaredError = n > 0 ? sse / n : 0;
+ const safeMSE = Math.max(meanSquaredError, Number.EPSILON);
+
+ return {
+ aic: 2 * k + n * Math.log(safeMSE),
+ mse: meanSquaredError,
+ };
+ }
+
+ function runModel(cfg) {
+ const obj = new Objective(input, cfg.ode, 500, 2);
+ const optim = new PSO(obj, alg.particles, cfg.bounds);
+ optim.run(alg.c1, alg.c2, alg.w, alg.iterations);
+
+ const sol = RK4(obj.f, obj.timeArray, obj.y0, optim.pos_best_g);
+ const celsM = sol.map((row) => row[0]);
+ const subsM = sol.map((row) => row[1]);
+
+ const { aic, mse } = calculateAIC(obj, sol, optim.pos_best_g);
+
+ let objectiveText = "MSE: N/A";
+ if (Number.isFinite(mse)) {
+ const absValue = Math.abs(mse);
+ const formatted =
+ absValue !== 0 && (absValue >= 1e3 || absValue < 1e-2)
+ ? mse.toExponential(4)
+ : mse.toFixed(8);
+ objectiveText = `MSE: ${formatted}`;
+ }
+
+ Plotly.newPlot(
+ document.getElementById(cfg.plotId),
+ [
+ ...basePlot,
+ {
+ x: obj.timeArray,
+ y: subsM,
+ name: "Model fit for the substrate",
+ mode: "lines",
+ line: { color: "#4a90e2" },
+ },
+ {
+ x: obj.timeArray,
+ y: celsM,
+ name: "Model fit for the cells",
+ mode: "lines",
+ line: { color: "#50e3c2" },
+ },
+ ],
+ {
+ margin: { t: 10, b: 30 },
+ paper_bgcolor: "#f0f4f8",
+ plot_bgcolor: "#f0f4f8",
+ xaxis: {
+ title: { text: "Time (h)" },
+ },
+ yaxis: {
+ title: { text: "Concentration (g/L)" },
+ },
+ legend: {
+ orientation: "h",
+ yanchor: "top",
+ y: -0.2,
+ xanchor: "left",
+ x: 0,
+ font: { size: 10 },
+ },
+ annotations: [
+ {
+ text: objectiveText,
+ x: 1,
+ y: 1,
+ xref: "paper",
+ yref: "paper",
+ xanchor: "right",
+ yanchor: "top",
+ showarrow: false,
+ font: { size: 12, color: "#1f2933" },
+ bordercolor: "#d9e2ef",
+ borderwidth: 1,
+ borderpad: 6,
+ },
+ ],
+ },
+ );
+
+ const params = optim.pos_best_g.map((p) => p.toFixed(3));
+ const paramContainer = document.getElementById(cfg.paramDiv);
+ paramContainer.classList.add("model-params");
+ paramContainer.innerHTML = "";
+
+ const label = document.createElement("div");
+ label.className = "param-label";
+ label.textContent = "Parameters:";
+ paramContainer.appendChild(label);
+
+ const list = document.createElement("ul");
+ list.className = "param-list";
+
+ params.forEach((value, i) => {
+ const detail = cfg.paramDetails[i];
+ const item = document.createElement("li");
+ const mathSpan = document.createElement("span");
+ const latexUnit = detail.unitLatex;
+ const latexExpression = latexUnit
+ ? `${detail.latex} = ${value}\\,${latexUnit}`
+ : `${detail.latex} = ${value}`;
+
+ if (typeof katex !== "undefined") {
+ katex.render(latexExpression, mathSpan, { throwOnError: false });
+ } else {
+ mathSpan.textContent = latexUnit
+ ? `${detail.latex} = ${value} ${detail.unitText || ""}`
+ : `${detail.latex} = ${value}`;
+ }
+
+ item.appendChild(mathSpan);
+ list.appendChild(item);
+ });
+
+ paramContainer.appendChild(list);
+
+ return { title: cfg.title, aic, key: cfg.key };
+ }
+
+ const models = [
+ {
+ key: "aiba",
+ ode: aibaPirt,
+ bounds: boundsOpt.aiba || getDefaultBounds("aiba"),
+ plotId: "Plot1",
+ paramDiv: "AibaParam",
+ title: "Pirt-Aiba",
+ paramDetails: getParamDetails("aiba"),
+ container: document.getElementById("Plot1")?.closest(".box"),
+ },
+ {
+ key: "andrews",
+ ode: andrewsPirt,
+ bounds: boundsOpt.andrews || getDefaultBounds("andrews"),
+ plotId: "Plot2",
+ paramDiv: "AndrewsParam",
+ title: "Pirt-Andrews",
+ paramDetails: getParamDetails("andrews"),
+ container: document.getElementById("Plot2")?.closest(".box"),
+ },
+ {
+ key: "bergter",
+ ode: bergterPirt,
+ bounds: boundsOpt.bergter || getDefaultBounds("bergter"),
+ plotId: "Plot3",
+ paramDiv: "BergterParam",
+ title: "Pirt-Bergter",
+ paramDetails: getParamDetails("bergter"),
+ container: document.getElementById("Plot3")?.closest(".box"),
+ },
+ {
+ key: "contois",
+ ode: contoisPirt,
+ bounds: boundsOpt.contois || getDefaultBounds("contois"),
+ plotId: "Plot4",
+ paramDiv: "ContoisParam",
+ title: "Pirt-Contois",
+ paramDetails: getParamDetails("contois"),
+ container: document.getElementById("Plot4")?.closest(".box"),
+ },
+ {
+ key: "monod",
+ ode: monodPirt,
+ bounds: boundsOpt.monod || getDefaultBounds("monod"),
+ plotId: "Plot5",
+ paramDiv: "MonodParam",
+ title: "Pirt-Monod",
+ paramDetails: getParamDetails("monod"),
+ container: document.getElementById("Plot5")?.closest(".box"),
+ },
+ {
+ key: "moser",
+ ode: moserPirt,
+ bounds: boundsOpt.moser || getDefaultBounds("moser"),
+ plotId: "Plot6",
+ paramDiv: "MoserParam",
+ title: "Pirt-Moser",
+ paramDetails: getParamDetails("moser"),
+ container: document.getElementById("Plot6")?.closest(".box"),
+ },
+ {
+ key: "tessier",
+ ode: tessierPirt,
+ bounds: boundsOpt.tessier || getDefaultBounds("tessier"),
+ plotId: "Plot7",
+ paramDiv: "TessierParam",
+ title: "Pirt-Tessier",
+ paramDetails: getParamDetails("tessier"),
+ container: document.getElementById("Plot7")?.closest(".box"),
+ },
+ ];
+
+ const modelByKey = new Map(models.map((model) => [model.key, model]));
+
+ const results = [];
+ for (let idx = 0; idx < models.length; idx++) {
+ const cfg = models[idx];
+ const res = runModel(cfg);
+ results.push(res);
+ onProgress(idx + 1, models.length);
+ await new Promise((resolve) => setTimeout(resolve, 0));
+ }
+
+ results.sort((a, b) => a.aic - b.aic);
+ const comparisonBox = document.getElementById("comparison");
+ const chartParent = models[0]?.container?.parentElement;
+
+ if (comparisonBox && chartParent) {
+ results.forEach((result) => {
+ const modelCfg = modelByKey.get(result.key);
+ if (!modelCfg?.container) {
+ return;
+ }
+ chartParent.insertBefore(modelCfg.container, comparisonBox);
+ });
+ }
+
+ const compDiv = document.getElementById("comparison");
+ compDiv.innerHTML =
+ `
+