Team:Groningen/Modelling/Arsenic.js

From 2009.igem.org

(Difference between revisions)
m
m
Line 119: Line 119:
     //var dy = df(x);
     //var dy = df(x);
     //alert("x="+x+", y="+y+", dy="+dy+"\nlow="+low+", ylow="+ylow+"\nhigh="+high+", yhigh="+yhigh);
     //alert("x="+x+", y="+y+", dy="+dy+"\nlow="+low+", ylow="+ylow+"\nhigh="+high+", yhigh="+yhigh);
 +
    // It can be shown that the division is always >=0 by multiplying by (ylow-yhigh)/(ylow-yhigh).
 +
    nx = Math.sqrt((Math.pow(high,2)*ylow-Math.pow(low,2)*yhigh)/(ylow-yhigh));
     if (y*ylow>0) { // Checking if y and ylow have the same sign
     if (y*ylow>0) { // Checking if y and ylow have the same sign
       low = x;
       low = x;
Line 126: Line 128:
       yhigh = y;
       yhigh = y;
     }
     }
-
     x = (low+high)/2.0
+
     x = (low<=nx && nx<=high) ? nx : -nx; // Select either + or - the solution
 +
    //x = (low+high)/2.0
     /*x = x - y/dy;
     /*x = x - y/dy;
     if (!(x<high) && !(x>low)) {
     if (!(x<high) && !(x>low)) {

Revision as of 16:27, 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);
 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 aAs = c.k8*ArsB*c.Vc;
 var bAs = c.v5*c.Vc*c.K7*(1+ArsR/c.K1d) + c.k8*ArsB*(c.Vs*c.K5 + c.AsT);
 var cAs = c.v5*c.K7*(1+ArsR/c.K1d)*c.AsT;
 var AsinT = 2*cAs/(bAs + Math.sqrt(Math.pow(bAs,2) - 4*aAs*cAs));
 if (isNaN(AsinT)) AsinT = 0;
 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 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);
   //var dy = df(x);
   //alert("x="+x+", y="+y+", dy="+dy+"\nlow="+low+", ylow="+ylow+"\nhigh="+high+", yhigh="+yhigh);
   // It can be shown that the division is always >=0 by multiplying by (ylow-yhigh)/(ylow-yhigh).
   nx = Math.sqrt((Math.pow(high,2)*ylow-Math.pow(low,2)*yhigh)/(ylow-yhigh));
   if (y*ylow>0) { // Checking if y and ylow have the same sign
     low = x;
     ylow = y;
   } else {
     high = x;
     yhigh = y;
   }
   x = (low<=nx && nx<=high) ? nx : -nx; // Select either + or - the solution
   //x = (low+high)/2.0
   /*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.");
   }*/
 }
 alert("x="+x+", iter="+iter);
 return x;

}

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