Team:Groningen/Modelling/Arsenic.js
From 2009.igem.org
m |
m |
||
Line 67: | Line 67: | ||
ArsR*(1 + Asin/c.K1d) + 2*(c.ars1T+c.ars2T)*ArsR2/(c.K3d2+ArsR2) - x.ArsRT]; | ArsR*(1 + Asin/c.K1d) + 2*(c.ars1T+c.ars2T)*ArsR2/(c.K3d2+ArsR2) - x.ArsRT]; | ||
}; | }; | ||
- | var Asinlow = | + | var Asinlow = x.AsinT/(1+x.ArsRT/c.K1d+x.MBPArsRT/c.KMd+x.KostalArsRT/c.KKd |
+ | +c.nf*x.fMTT*Math.pow(x.AsinT/c.KFd,c.nf)), Asinhigh = x.AsinT; | ||
var ArsRlow = x.ArsRT/(1+x.AsinT/c.K1d+2*(c.ars1T+c.ars2T)*x.ArsRT/c.K3d2), ArsRhigh = x.ArsRT; | var ArsRlow = x.ArsRT/(1+x.AsinT/c.K1d+2*(c.ars1T+c.ars2T)*x.ArsRT/c.K3d2), ArsRhigh = x.ArsRT; | ||
var AsinArsR = findSingleZero2D(fAsin,[Asinlow,ArsRlow],[Asinhigh,ArsRhigh]); | var AsinArsR = findSingleZero2D(fAsin,[Asinlow,ArsRlow],[Asinhigh,ArsRhigh]); |
Revision as of 18:22, 15 October 2009
arsenicModelMaxIter = 1000; function arsenicModelConstants() { // Prevents other scripts from overriding them.
return { // Extra-cellular GlpFTfactor: 1, // relative amount of GlpF hasGlpFPlasmid: false, // set to true for experiments with GlpF on a plasmid // Intra-cellular ars1T: 1.6605e-9, ars2T: 166.05e-9, pro: 16.605e-9, // 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: 1000, //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 // = ArsR (1 + Asin/KRd + 2 ArsR arsT KAd²/(KAd²+ArsR²) / KAd²) - ArsRT // = ArsR (1 + Asin/KRd) + 2 arsT / (1+KAd²/ArsR²)) - ArsRT var fAsin = function(AsinArsR) { var Asin = AsinArsR[0], ArsR = AsinArsR[1]; var Asinnf = Math.pow(Asin,c.nf), ArsR2 = Math.pow(ArsR,2); return [Asin*(1 + ArsR/c.K1d + x.MBPArsRT/(c.KMd + Asin) + x.KostalArsRT/(c.KKd + Asin)) + c.nf*x.fMTT*Asinnf/(Math.pow(c.KFd,c.nf)+Asinnf) - x.AsinT, ArsR*(1 + Asin/c.K1d) + 2*(c.ars1T+c.ars2T)*ArsR2/(c.K3d2+ArsR2) - x.ArsRT]; }; var Asinlow = x.AsinT/(1+x.ArsRT/c.K1d+x.MBPArsRT/c.KMd+x.KostalArsRT/c.KKd +c.nf*x.fMTT*Math.pow(x.AsinT/c.KFd,c.nf)), Asinhigh = x.AsinT; var ArsRlow = x.ArsRT/(1+x.AsinT/c.K1d+2*(c.ars1T+c.ars2T)*x.ArsRT/c.K3d2), ArsRhigh = x.ArsRT; var AsinArsR = findSingleZero2D(fAsin,[Asinlow,ArsRlow],[Asinhigh,ArsRhigh]); var Asin = AsinArsR[0]; var ArsR = AsinArsR[1];
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 };
}
// 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); }
return { // Extra-cellular 'AsexT': (AsT-c.Vc*AsinT)/c.Vs, // Intra-cellular 'ArsBT': ArsB*(1.0+Asin/c.K7), 'AsinT': AsinT, 'ArsRT': ArsR*(1.0+Asin/c.K1d+ArsR*arsFraction*(c.ars1T+c.ars2T)/c.K3d2), '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 };
}
// 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;
}
// Assuming there is exactly one zero ([0,0]) of f(x) for x in [low,high] this function will find it // Also, f(x) should be monotone in all variables (for all components of its output). // And f(x) should be finite for all values in the given domain. function findSingleZero2D(f,low,high,tol) {
if (tol===undefined) tol = 1e-6; var busy; //alert('bla'); alert('f(x0)='+f([(low[0]+high[0])/2,(low[1]+high[1])/2])); do { busy = false; //alert([low,high]); if (2*(high[0]-low[0])/(Math.abs(low[0])+Math.abs(high[0]))>tol) { var x0l = findSingleZero(function(x0){return f([x0,low[1]])[0];},low[0],high[0],tol); var x0h = findSingleZero(function(x0){return f([x0,high[1]])[0];},low[0],high[0],tol); low[0] = Math.min(x0l,x0h); high[0] = Math.max(x0l,x0h); busy = true; } //alert([low,high]); if (2*(high[1]-low[1])/(Math.abs(low[1])+Math.abs(high[1]))>tol) { var x1l = findSingleZero(function(x1){return f([low[0],x1])[1];},low[1],high[1],tol); var x1h = findSingleZero(function(x1){return f([high[0],x1])[1];},low[1],high[1],tol); low[1] = Math.min(x1l,x1h); high[1] = Math.max(x1l,x1h); busy = true; } } while(busy); alert('f(x)='+f([(low[0]+high[0])/2,(low[1]+high[1])/2])); return [(low[0]+high[0])/2,(low[1]+high[1])/2];
}