Team:Groningen/Modelling/Arsenic.js

From 2009.igem.org

(Difference between revisions)

Revision as of 14:42, 6 October 2009

arsenicModelMaxIter = 1000; function arsenicModelConstants() { // Prevents other scripts from overriding them.

 return {
   // Extra-cellular
   GlpFTfactor: 1, // relative amount of GlpF
   // Intra-cellular
   ars1T: 1.6605e-9,
   ars2T: 16.605e-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
   k8: 0, // export (unknown)
   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
   K7: 1, // export (unknown)
   tauB: 60,  // half-life of ArsB
   tauR: 60,  // half-life of ArsR
   tauG: 60,  // half-life of GV
   tauK: 120, // half-life of ArsR fusion of Kostal2004
   tauM: 120, // half-life of MBP-ArsR
   tauF: 120, // half-life of fMT
   beta1: 1000, // production rate of ArsR behind ars1 // TODO: Rename them all!
   beta3: 1000, // production rate of ArsR behind pro
   beta4: 1000, // production rate of ArsB behind ars1
   beta5: 1000, // production rate of GV behind ars2
   betaK: 500,  // production rate of ArsR fusion of Kostal2004 behind proK
   betaM: 500,  // production rate of MBP-ArsR behind proM
   betaF: 500,  // 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 + ArsRT/(KRd + As(III)in) + MBPArsRT/(KMd + As(III)in)
 //            + fMTT As(III)in^(nf-1)/(KFdnf + As(III)in^nf)
 var fAsin = function(Asin) {
               var Asinnf = Math.pow(Asin,c.nf);
               return Asin*(1 + x.ArsRT/(c.K1d + Asin) + x.MBPArsRT/(c.KMd + Asin) + x.KostalArsRT/(c.KKd + Asin))
                       + x.fMTT*Asinnf/(Math.pow(c.KFd,c.nf)+Asinnf) - x.AsinT;
             };
 var Asinlow = 0, Asinhigh = x.AsinT;
 var Asin = findSingleZero(fAsin,Asinlow,Asinhigh);

 var ArsR = x.ArsRT*c.K1d/(c.K1d+Asin);
 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.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)
 // 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.v5/(c.k8*ArsB) - Asin)
                       + c.Vc*Asin*(1+ArsR/c.K1d+MBPArsR/c.KMd+KostalArsR/c.KKd)
                       + c.Vc*fMT*Math.pow(Asin/c.KFd,c.nf)
                       - AsT;
             };
 var Asinlow = 0, Asinhigh = Math.min(AsT/c.Vc,c.K7*c.v5/(c.k8*ArsB));
 var Asin = findSingleZero(fAsin,Asinlow,Asinhigh);
 var Asinnf = Math.pow(Asin,c.nf);
 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),
   '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)
   x = Math.min(Math.max(x,low+0.1*(high-low)),low+0.9*(high-low)); // Ensures convergence
 }
 return x;

}