Team:Groningen/Modelling/Arsenic.js

From 2009.igem.org

(Difference between revisions)

Latest revision as of 08:57, 18 October 2009

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;

}