Team:Groningen/Modelling/Arsenic.js

From 2009.igem.org

(Difference between revisions)

Revision as of 18:09, 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 = 0, Asinhigh = x.AsinT;
 var ArsRlow = 0/*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];

}