Team:Groningen/Modelling/Arsenic.js

From 2009.igem.org

(Difference between revisions)
m (Debugging.)
m
Line 91: Line 91:
   var b = c.v5*c.Vc*c.K7*(1+ArsR/c.K1d) + c.k8*ArsB*(c.Vs*c.K5 + c.AsT);
   var b = c.v5*c.Vc*c.K7*(1+ArsR/c.K1d) + c.k8*ArsB*(c.Vs*c.K5 + c.AsT);
   var c = c.v5*c.K7*(1+ArsR/c.K1d)*c.AsT;
   var c = c.v5*c.K7*(1+ArsR/c.K1d)*c.AsT;
 +
  alert("a="+a+", b="+b+", c="+c);
   var AsinT = 2*c/(b + Math.sqrt(Math.pow(b,2) - 4*a*c));
   var AsinT = 2*c/(b + Math.sqrt(Math.pow(b,2) - 4*a*c));
   var Asin = AsinT*c.K1d/(c.K1d+ArsR);
   var Asin = AsinT*c.K1d/(c.K1d+ArsR);
Line 124: Line 125:
     if (!(x<high) && !(x>low)) {
     if (!(x<high) && !(x>low)) {
       x = (low+high)/2.0;
       x = (low+high)/2.0;
 +
      alert("Taking middle.");
     } else if (!(x<high)) {
     } else if (!(x<high)) {
       x = low+0.9*(high-low);
       x = low+0.9*(high-low);
 +
      alert("Taking high.");
     } else if (!(x>low)) {
     } else if (!(x>low)) {
       x = low+0.1*(high-low);
       x = low+0.1*(high-low);
 +
      alert("Taking low.");
     }
     }
   }
   }
   return x;
   return x;
}
}

Revision as of 14:16, 18 September 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,
   // Reaction constants
   k8: 0, // export
   K1d: 6e-6, // ArsR-As
   K3d2: Math.pow(0.33e-6,2), // ArsR-ars
   v5: 3.1863e-6, // import
   K5: 27.718e-6, // import
   K7: 0, // export
   tauB: 60, // half-life of ArsB
   tauR: 60, // half-life of ArsR
   tauG: 60, // half-life of GV
   beta1: 1000, // production rate of ArsR behind ars1
   beta3: 1000, // production rate of ArsR behind pro
   beta4: 1000, // production rate of ArsB behind ars1
   beta5: 1000, // production rate of GV behind ars2
   // Volumes
   Vs: 1.1e-3,
   Vc: 0.0073e-3,
   // Amount of arsenic
   AsT: 10e-9
 };

}

// 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) { 

 // First determine equilibrium without arsenic
 var nc = {};
 for(var a in c) nc[a] = c[a];
 nc.AsT = 0;
 var x = arsenicModelEquilibrium(c);

 // Return equilibrium with added arsenic
 x.AsexT = c.AsT/c.Vs;
 return x;

}

// Computes gradient of the variables in the arsenic model function arsenicModelGradient(c,x) { 

 var Asin = x.AsinT/2.0;
 var ArsR = x.ArsRT*c.K1d/(c.K1d+Asin);
 var pAsin = Number.POSITIVE_INFINITY, pArsR = Number.POSITIVE_INFINITY;
 for(var iter=1;
     iter<=arsenicModelMaxIter && (Math.abs(Asin-pAsin)/Asin>1e-6 || Math.abs(ArsR-pArsR)/ArsR>1e-6);
     iter++)
 {
   pAsin = Asin;
   pArsR = ArsR;
   Asin = x.AsinT*c.K1d/(c.K1d+ArsR);
   ArsR = x.ArsRT*c.K1d/(c.K1d+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;
 return {
   // Extra-cellular
   AsexT: -(c.Vc/c.Vs)*dAsinTdt,
   // Intra-cellular
   ArsBT: c.beta4*ars1 - (Math.LN2/c.tauB)*ArsB,
   AsinT: dAsinTdt,
   ArsRT: c.beta1*ars1 + c.beta3*c.pro - (Math.LN2/c.tauR)*ArsR,
   GV: c.beta5*ars2 - (Math.LN2/c.tauG)*x.GV
 };

}

// Computes the equilibrium from constants function arsenicModelEquilibrium(c) { 

 // 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 dfArsR = function(ArsR) { return -c.K3d2 + 2*c.beta3*(c.tauR/Math.LN2)*c.pro*ArsR - 3*Math.pow(ArsR,2); }
 var ArsRlow = 0, ArsRhigh = (c.beta1*c.ars1T + c.beta3*c.pro)*(c.tauR/Math.LN2)*c.K3d2;
 var ArsR = findSingleZero(fArsR,dfArsR,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 GV = c.beta5*(c.tauB/Math.LN2)*c.ars2T*arsFraction;

 // Determine intra-cellular concentration of As(III)
 var a = c.k8*ArsB*c.Vc;
 var b = c.v5*c.Vc*c.K7*(1+ArsR/c.K1d) + c.k8*ArsB*(c.Vs*c.K5 + c.AsT);
 var c = c.v5*c.K7*(1+ArsR/c.K1d)*c.AsT;
 alert("a="+a+", b="+b+", c="+c);
 var AsinT = 2*c/(b + Math.sqrt(Math.pow(b,2) - 4*a*c));
 var Asin = AsinT*c.K1d/(c.K1d+ArsR);

 return {
   // Extra-cellular
   'AsexT': (c.AsT-c.Vc*AsinT)/c.Vs,
   // Intra-cellular
   'ArsBT': ArsB*(1.0+Asin/c.K7),
   'AsinT': AsinT,
   'ArsRT': ArsR*(1.0+Asin/c.K1d),
   'GV': GV
 };

}

// Assuming there is exactly one zero of f(x) for x in [low,high] this function will find it function findSingleZero(f,df,low,high,tol) { 

 if (tol===undefined) tol = 1e-6;
 var x = (low+high)/2.0;
 for(var iter=1;
     iter<=arsenicModelMaxIter && (high-low)/x>tol;
     iter++)
 {
   var y = f(x);
   var dy = df(x);
   alert("x="+x+" ["+low+","+high+"], y="+y+", dy="+dy);
   if (y>0) {
     low = x;
   } else {
     high = x;
   }
   x = x - y/dy;
   if (!(x<high) && !(x>low)) {
     x = (low+high)/2.0;
     alert("Taking middle.");
   } else if (!(x<high)) {
     x = low+0.9*(high-low);
     alert("Taking high.");
   } else if (!(x>low)) {
     x = low+0.1*(high-low);
     alert("Taking low.");
   }
 }
 return x;

}

Retrieved from "http://2009.igem.org/Team:Groningen/Modelling/Arsenic.js"