"
Team:Groningen/Modelling/Arsenic.js
From 2009.igem.org
arsenicModelMaxIter = 1000; function arsenicModelConstants() { // Prevents other scripts from overriding them.
return {
// Extra-cellular
GlpFTfactor: 2, // relative amount of GlpF
hasGlpFPlasmid: false, // set to true for experiments with GlpF on a plasmid
// Intra-cellular
ars1T: 1.6605e-9,
ars2T: 0,
pro: 0, // TODO: rename to proR
proK: 0, // Promoters in front of ArsR fusion from Kostal2004 (handled completely analogously to MBPArsR)
proM: 0, // Promoters in front of our ArsR fusion (with MBP)
proF: 0, // Promoters in front of fMT
// Reaction constants
K1d: 6e-6, // ArsR-As
KKd: 6e-6, // ArsR-As (ArsR fusion from Kostal2004)
KMd: 6e-6, // MBPArsR-As
KFd: 6e-6, // fMT-As
nf: 3, // Hill coefficient for fMT + As
K3d2: Math.pow(0.33e-6,2), // ArsR-ars
v5: 3.1863e-6, // import
K5: 27.718e-6, // import
k8: 1e2, // export (unknown)
K7: 1e-4, // export (unknown)
tauB: 1800, // half-life of ArsB
tauR: 6, // half-life of ArsR
tauG: 6, // half-life of GV
tauK: 960, // half-life of ArsR fusion of Kostal2004
tauM: 6, // half-life of MBP-ArsR
tauF: 6, // half-life of fMT
beta1: 100, // production rate of ArsR behind ars1 // TODO: Rename them all!
beta3: 100, // production rate of ArsR behind pro
beta4: 95, // production rate of ArsB behind ars1
beta5: 6.6, // production rate of GV behind ars2
betaK: 15.4, // production rate of ArsR fusion of Kostal2004 behind proK
betaM: 26.6, // production rate of MBP-ArsR behind proM
betaF: 200, // production rate of fMT behind proF
// Volumes
Vs: 1.1e-3,
Vc: 0.0073e-3
};
}
// Returns an array with all the necessary variables initialized as if a solution with cells in equilibrium suddenly gets a shot of arsenic function arsenicModelInitialization(c,AsT) {
if (c.AsT!=undefined && AsT==undefined) AsT = c.AsT; // First determine equilibrium without arsenic var x = arsenicModelEquilibrium(c,0);
// Return equilibrium with added arsenic x.AsexT = AsT/c.Vs; return x;
}
// Computes gradient of the variables in the arsenic model function arsenicModelGradient(c,x) {
// Solve 0 = As(III)in (1 + ArsR/KRd + MBPArsRT/(KMd + As(III)in)
// + nf fMTT As(III)in^(nf-1)/(KFdnf + As(III)in^nf) - As(III)inT
// and 0 = ArsR (1 + Asin/KRd + 2 ArsR ars / KAd²) - ArsRT
var arsT = c.ars1T+c.ars2T, KFdnf = Math.pow(c.KFd,c.nf);
var Asinlow = 0, Asinhigh = x.AsinT, Asintol = x.AsinT*1e-6, Asinlownf, Asinhighnf;
var ArsRlow = 0, ArsRhigh = x.ArsRT, ArsRtol = x.ArsRT*1e-6;
var arslow = arsT/(1+Math.pow(ArsRhigh,2)/c.K3d2), arshigh = arsT, arstol = arsT*1e-6;
for(var iter=1; iter<=arsenicModelMaxIter &&
(Asinhigh-Asinlow>Asintol || ArsRhigh-ArsRlow>ArsRtol || arshigh-arslow>arstol); iter++) {
Asinlownf = Math.pow(Asinlow,c.nf-1); // It is not safe to divide out Asinlow/high, so we do c.nf-1 here
Asinhighnf = Math.pow(Asinhigh,c.nf-1);
Asinlow = x.AsinT/(1+ArsRhigh/c.K1d+x.MBPArsRT/(c.KMd+Asinlow)+x.KostalArsRT/(c.KKd+Asinlow)
+c.nf*x.fMTT*Asinhighnf/(KFdnf+Asinlownf*Asinlow));
Asinhigh = x.AsinT/(1+ArsRlow/c.K1d+x.MBPArsRT/(c.KMd+Asinhigh)+x.KostalArsRT/(c.KKd+Asinhigh)
+c.nf*x.fMTT*Asinlownf/(KFdnf+Asinhighnf*Asinhigh));
ArsRlow = x.ArsRT/(1+Asinhigh/c.K1d+2*arshigh*ArsRhigh/c.K3d2);
ArsRhigh = x.ArsRT/(1+Asinlow/c.K1d+2*arslow*ArsRlow/c.K3d2);
arslow = arsT*c.K3d2/(c.K3d2+Math.pow(ArsRhigh,2));
arshigh = arsT*c.K3d2/(c.K3d2+Math.pow(ArsRlow,2));
}
//if (iter>15) alert(iter);
var Asin = Asinlow + 0.5*(Asinhigh-Asinlow);
var ArsR = ArsRlow + 0.5*(ArsRhigh-ArsRlow);
//var ars = arslow + 0.5*(arshigh-arslow);
/*alert('err='+[Asin*(1+ArsR/c.K1d+x.MBPArsRT/(c.KMd+Asin)+x.KostalArsRT/(c.KKd+Asin))
+ c.nf*x.fMTT*Math.pow(Asin,c.nf)/(KFdnf+Math.pow(Asin,c.nf)) - x.AsinT,
ArsR*(1+Asin/c.K1d+2*(c.ars1T+c.ars2T)*ArsR/c.K3d2) - x.ArsRT,
ars*(c.K3d2+Math.pow(ArsRhigh,2))/c.K3d2 - arsT]);*/
/*alert('gradient, [Asin,ArsR,arsF]='+[Asin,ArsR,ars/arsT]
+', [AsinT-Asin,ArsRT-(ArsR+2*(arsT-ars))]='+[x.AsinT-Asin,x.ArsRT-(ArsR+2*(arsT-ars))]);*/
var KostalArsR = x.KostalArsRT*c.KKd/(c.KKd+Asin); var MBPArsR = x.MBPArsRT*c.KMd/(c.KMd+Asin); var fMT = x.fMTT*c.KFd/(c.KFd+Asin); var ArsB = x.ArsBT*c.K7/(c.K7+Asin); var ArsBAs = x.ArsBT*Asin/(c.K7+Asin); var arsFraction = c.K3d2/(c.K3d2+Math.pow(ArsR,2)); var ars1 = c.ars1T*arsFraction; var ars2 = c.ars2T*arsFraction; var dAsinTdt = (c.hasGlpFPlasmid?c.GlpFTfactor:1)*c.v5*x.AsexT / (c.K5+x.AsexT) - c.k8*ArsBAs; var dArsRT = c.beta1*ars1 + c.beta3*c.pro - (Math.LN2/c.tauR)*ArsR;
return {
// Extra-cellular
AsexT: -(c.Vc/c.Vs)*dAsinTdt,
// Intra-cellular
ArsBT: c.beta4*ars1 - (Math.LN2/c.tauB)*ArsB,
AsinT: dAsinTdt,
ArsRT: dArsRT,
KostalArsRT: c.betaK*c.proK - (Math.LN2/c.tauK)*KostalArsR,
MBPArsRT: c.betaM*c.proM - (Math.LN2/c.tauM)*MBPArsR,
fMTT: c.betaF*c.proF - (Math.LN2/c.tauF)*fMT,
GV: c.beta5*ars2 - (Math.LN2/c.tauG)*x.GV,
// Subcomponents (selected)
'_ArsR': ArsR,
'_Asin': Asin,
'_arsF': arsFraction
};
}
// Computes the equilibrium from constants function arsenicModelEquilibrium(c,AsT) {
if (c.AsT!=undefined && AsT==undefined) AsT = c.AsT;
// Solve 0 = (β1 ars1T + β3 pro) (τR/ln(2)) K3d² - K3d² ArsR + β3 (τR/ln(2)) pro ArsR² - ArsR³
var fArsR = function(ArsR) { return (c.beta1*c.ars1T + c.beta3*c.pro)*(c.tauR/Math.LN2)*c.K3d2 - c.K3d2*ArsR + c.beta3*(c.tauR/Math.LN2)*c.pro*Math.pow(ArsR,2) - Math.pow(ArsR,3); }
var ArsRlow = 0, ArsRhigh = (c.beta1*c.ars1T + c.beta3*c.pro)*(c.tauR/Math.LN2);
var ArsR = findSingleZero(fArsR,ArsRlow,ArsRhigh);
// Compute ArsB and GV var arsFraction = c.K3d2/(c.K3d2+Math.pow(ArsR,2)); var ArsB = c.beta4*(c.tauB/Math.LN2)*c.ars1T*arsFraction; var KostalArsR = c.betaK*(c.tauK/Math.LN2)*c.proK; var MBPArsR = c.betaM*(c.tauM/Math.LN2)*c.proM; var fMT = c.betaF*(c.tauF/Math.LN2)*c.proF; var GV = c.beta5*(c.tauB/Math.LN2)*c.ars2T*arsFraction;
// Determine intra-cellular concentration of As(III)
if (c.K7*(c.hasGlpFPlasmid?c.GlpFTfactor:1)*c.v5==0) {
var Asin = 0, AsinT = 0;
} else if (c.k8*ArsB==0) {
var AsinT = AsT/c.Vc;
var fAsin = function(Asin) {
return Asin*(1+ArsR/c.K1d+MBPArsR/c.KMd+KostalArsR/c.KKd) + c.nf*fMT*Math.pow(Asin/c.KFd,c.nf) - AsinT;
};
var Asinlow = 0, Asinhigh = AsinT;
var Asin = findSingleZero(fAsin,Asinlow,Asinhigh);
} else {
// By solving 0 = Vs K5 As(III)in / (K7 v5/(k8 ArsB) - As(III)in)
// + Vc As(III)in (1 + ArsR/KRd + MBPArsR/KMd + fMTT As(III)in^(nf-1)/KFd^nf)
// - As(III)T
var fAsin = function(Asin) {
return c.Vs*c.K5*Asin/(c.K7*(c.hasGlpFPlasmid?c.GlpFTfactor:1)*c.v5/(c.k8*ArsB) - Asin)
+ c.Vc*Asin*(1+ArsR/c.K1d+MBPArsR/c.KMd+KostalArsR/c.KKd)
+ c.Vc*c.nf*fMT*Math.pow(Asin/c.KFd,c.nf)
- AsT;
};
var Asinlow = 0, Asinhigh = Math.min(AsT/c.Vc,c.K7*(c.hasGlpFPlasmid?c.GlpFTfactor:1)*c.v5/(c.k8*ArsB));
var Asin = findSingleZero(fAsin,Asinlow,Asinhigh);
var AsinT = Asin*(1+ArsR/c.K1d+MBPArsR/c.KMd+KostalArsR/c.KKd) + fMT*Math.pow(Asin/c.KFd,c.nf);
}
var ArsRT = ArsR*(1.0+Asin/c.K1d+2*ArsR*arsFraction*(c.ars1T+c.ars2T)/c.K3d2);
/*alert('equilibrium, [Asin,ArsR,arsF]='+[Asin,ArsR,arsFraction]
+', [AsinT-Asin,ArsRT-(ArsR+2*(arsT-ars))]='+[AsinT-Asin,ArsRT-(ArsR+2*(c.ars1T+c.ars2T)*(1-arsFraction))]);*/
return {
// Extra-cellular
'AsexT': (AsT-c.Vc*AsinT)/c.Vs,
// Intra-cellular
'ArsBT': ArsB*(1.0+Asin/c.K7),
'AsinT': AsinT,
'ArsRT': ArsRT,
'KostalArsRT': KostalArsR*(1.0+Asin/c.KKd),
'MBPArsRT': MBPArsR*(1.0+Asin/c.KMd),
'fMTT': fMT*(1.0+Math.pow(Asin/c.KFd,c.nf)),
'GV': GV,
// Subcomponents (selected)
'_ArsR': ArsR,
'_Asin': Asin,
'_arsF': arsFraction
};
}
// Assuming there is exactly one zero of f(x) for x in [low,high] this function will find it function findSingleZero(f,low,high,tol) {
if (tol===undefined) tol = 1e-6;
var ylow = f(low);
var yhigh = f(high);
var x = (low+high)/2.0, y;
for(var iter=1;
iter<=arsenicModelMaxIter && (high-low)/x>tol;
iter++)
{
y = f(x);
if (y*ylow>0) { // Checking if y and ylow have the same sign
low = x;
ylow = y;
} else {
high = x;
yhigh = y;
}
x = high - yhigh*(high-low)/(yhigh-ylow); // Secant method (roughly)
if (isNaN(x)) x = low+0.5*(high-low);
x = Math.min(Math.max(x,low+0.1*(high-low)),low+0.9*(high-low)); // Ensures convergence
}
return x;
}
