Team:Groningen/Modelling/Arsenic.js
From 2009.igem.org
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 (unknown) 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: 1, // export (unknown) 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(nc);
// 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); } alert(iter);*/ // Asin = x.AsinT*c.K1d/(c.K1d+ArsR) // Asin = x.AsinT/(1+ArsR/c.K1d) // Asin = x.AsinT/(1+x.ArsRT/(c.K1d+Asin)) // Asin = x.AsinT*(c.K1d+Asin)/(c.K1d+Asin+x.ArsRT) // (c.K1d+Asin+x.ArsRT)*Asin = x.AsinT*(c.K1d+Asin) // c.K1d*Asin+Asin*Asin+x.ArsRT*Asin = x.AsinT*(c.K1d+Asin) // Asin*Asin + (c.K1d + x.ArsRT - x.AsinT)*Asin - x.AsinT*c.K1d = 0; var bAsin = c.K1d + x.ArsRT - x.AsinT; var cAsin = x.AsinT*c.K1d; var Asin = (-bAsin + Math.sqrt(Math.pow(bAsin,2)+4*cAsin))/2; var ArsR = x.ArsRT-(x.AsinT-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, 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 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 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,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;
}