Team:Groningen/Modelling/Characterization

 .intro { margin-left:0px; margin-top:10px; padding:10px; border-left:solid 5px #FFF6D5; border-right:solid 5px #FFF6D5; text-align:justify;background:#FFFFE5; }

Characterization
We have four kinds of parts we would like to characterize: Importers, Accumulators, Sensors and the GVP cluster. For this we have a number of methods to estimate specific parameters (detailed below), as well as a stochastic tool to fit our model to experimental data based on simulated annealing.

We have the following parts that we can characterize (RPS stands for Relative Promoter Strength)

Uptake measurements
To obtain data for the optimization procedure described above we conducted ICP-MS measurements on our cells containing various devices/parts. To optimize our findings we have conducted measurements both in time and in concentration. Details on ICP-MS the experiments can be found on our Protocol Page. Measurements have been conducted at times and concentrations as indicated in the table on the right. Results can be seen below.

First Uptake Measurement In this case we looked at the arsenic uptake of our wildtype and of our wildtype with a lot of pArs promoters with RFP. The raw data can be found at our Download Section

Getting (Vc/Vs) and other conditions To effectively determine our constants we need to give our model some extra information. For instance what kind of constructs are inside the cel for a given experiment. Also we need to give the model the volume of cells per liter of fluid. To optain this we use the obtained dry weight and calculate how much wet weight it would have been (assuming dry weight/wet weight = 0.3) and then use the density of E. coli (1100kg/m3) to obtain the cell volume for our sample and eventually the desired volume of cells per liter. Also we need the absolute value of Arsenic taken up by the cells in the assumption that we have one liter of sample and we know (Vc/Vs). Once we know all these parameters the optimization procedure can start.

Fluorescence Measurements Apart from giving our model all the conditions it needs to calculate all the constants by means of the optimization procedure, we have also conducted some fluorescence measurements and made growth curves of our construct with the pArs promoter with RFP. The cells where put into a solution with either no arsenic in it or at a concentration of 100 micromolair. On the left side one can see the graph of the luminance and on right side and on the right side one can see the coresponding grow curves. The raw data of these measurements can again be found under our Downloads Using a formula similar to the formula below we where able to derive a RPU of 2.3. This means that on average the ars promoter is 2.3 times more active at 100 micromolarity of arsenic (outside the cell) than if there is no arsenic in the solution. For a detailed calculation I would like to refer to our Downloads section under the RPU sheet.

Second ICPMS measurement For our second measurement I would like to refer to our Accumulation page. These measerments where higher than expected and they need further analyses before we can use them in tha characterization of our parts.

Optimization procedure
To fit our model to experimental data from different uptake experiments and/or papers we have implemented an optimization procedure that allows for experiments with different genotypes and circumstances by letting constants be overridden per experiment. It aims to optimize the sum of the RMS errors for each experiment using Simulated Annealing. By clicking the button "Fit" the optimization is started and its progress can be followed by looking at the table of constants and the graphs shown below the table (which are updated in real-time as the best solution is improved). In the optimizations procedure the ratio between v5 and K5 is kept constant. This is done because this ratio could be determined exactly using the data from the experiments done by Meng 2004 v5 itself is not fixed and may vary. (Note that in the code different experiments and graphs can be enabled/disabled, we have entered data from Kostal, Meng and our first ICP-MS measurements.)

    var experiments = {/*Meng2004: {constants:{Vc:0.0073,Vs:(1.1-0.0073),beta4:0,pro:0,ars2T:0},AsT:10e-6, data:{AsinT:[101.917808219178e-6,394.520547945205e-6,723.287671232877e-6, 1111.23287671233e-6,1229.58904109589e-6], time:[60,600,1200,2400,3600]}}, Singh2008: // We assume 5g/L wet cells were used... (at 1100kg/m^3) {constants:{Vc:(0.004545455),Vs:(1-(0.004545455)),pro:0,ars2T:0, proF:1.6605e-9},AsT:0.467154987e-6, data:{AsexT:[0.419211538e-6,0.391262322e-6,0.378368845e-6, 0.361791516e-6,0.332907991e-6,0.320748614e-6], time:[1.127*60,4.993*60,9.986*60,20.159*60,30.181*60,60.035*60]}}, Kostal2004fig3A: //  fig 3A {constants:{Vc:0.006666667,Vs:(1-0.006666667),pro:0,ars2T:0,proK:1.6605e-9},time:3600, data:{AsinT:[28.71e-6,78.87e-6,144.21e-6,377.19e-6,490.38e-6,617.76e-6,649.11e-6], AsT:[0.4e-6,1e-6,2e-6,5e-6,20e-6,50e-6,100e-6]}}, Kostal2004fig3B:  //fig 3B {constants:{Vc:0.006666667,Vs:(1-0.006666667),pro:0,ars2T:0,proK:1.6605e-9},AsT:20e-6, data:{AsinT:[2.25e-4,3.47e-4,4.19e-4,3.93e-4,4.19e-4,4.82e-4,4.82e-4,4.95e-4], time:[582,1212,1890,2514,3144,3828,4260,6036]}},*/ pSB1A2con:   //  concentration mode this is our first icps measerment wild type {constants:{Vc:0.000808081,Vs:(1-0.000808081),pro:0,ars2T:0},time:3600, data:{AsinT:[207.0208222e-6,229.0443139e-6,493.3262146e-6,585.8248799e-6], AsT:[10e-6,20e-6,50e-6,100e-6]}}, pSB1A2time:  //  concentration mode this is our first icps measerment wild type {constants:{Vc:0.002320346,Vs:(1-0.002320346),pro:0,ars2T:0},AsT:10e-6, data:{AsinT:[66.07047517e-6,83.68926855e-6,114.522157e-6,132.1409503e-6,207.0208222e-6], time:[180,600,1200,2400,3600]}}/*, pArsRRFPcon: // here the cell only contains extra RFP behind the the extra ArsR promoters. //   We incorporate this in our model by pretending RFP=GVP (1st icps) {constants:{Vc:0.001272727,Vs:(1-0.001272727),pro:0},time:3600, data:{AsinT:[136.5456487e-6,277.4959957e-6,290.7100908e-6,343.5664709e-6], AsT:[10e-6,20e-6,50e-6,100e-6]}}, pArsRRFPtime: // here the cell only contains extra RFP behind the the extra ArsR promoters. // We incorporate this in our model by pretending RFP=GVP (1st icps) {constants:{Vc:0.003333333,Vs:(1-0.003333333),pro:0},AsT:10e-6, data:{AsinT:[52.85638014e-6,92.49866524e-6,88.0939669e-6,136.5456487e-6], time:[180,600,2400,3600]}}*/};

/*var varsToMutate = ['K5','v5','K7','k8','tauB','beta4','tauR','beta1','tauF','betaF', 'tauK','betaK','tauG','beta5']; var mutateFuncs = {v5: function(v){return v.v5;}, K5: function(v){return v.K5;}, k8: function(v){return v.k8;}, K7: function(v){return v.K7;}, tauB: function(v){return v.tauB;}, tauR: function(v){return v.tauR;}, beta4: function(v){return v.beta4;}, beta1: function(v){return v.beta1;}, tauF: function(v){return v.tauF;}, betaF: function(v){return v.betaF;}, tauK: function(v){return v.tauK;}, betaK: function(v){return v.betaK;}, tauG: function(v){return v.tauG;}, beta5: function(v){return v.beta5;}};*/

var varsToMutate = [/*'v5_K5',*/'v5','k8_K7','k8','tauBbeta4','beta4', 'tauRbeta1_tauBbeta4','beta1_beta4'/*,'tauFbetaF','betaF', 'tauKbetaK','betaK','tauGbeta5','beta5','tauF','betaF','tauK','betaK','tauG','beta5','ars2T'*/]; var mutateFuncs = {v5: function(v){return v.v5;}, K5: function(v){return v.v5/0.11495;}, k8: function(v){return v.k8;}, K7: function(v){return v.k8/v.k8_K7;}, tauB: function(v){return v.tauBbeta4/v.beta4;}, beta4: function(v){return v.beta4;}, tauR: function(v){return v.tauRbeta1_tauBbeta4*v.tauBbeta4/(v.beta4*v.beta1_beta4);}, beta1: function(v){return v.beta4*v.beta1_beta4;}/*, //tauF: function(v){return v.tauF;}, //betaF: function(v){return v.betaF;}, //tauK: function(v){return v.tauK;}, //betaK: function(v){return v.betaK;}, //tauG: function(v){return v.tauG;}, //ars2T: function(v){return v.ars2T;}, //beta5: function(v){return v.beta5;}, tauF: function(v){return v.tauFbetaF/v.betaF;}, betaF: function(v){return v.betaF;}, tauK: function(v){return v.tauKbetaK/v.betaK;}, betaK: function(v){return v.betaK;}, tauG: function(v){return v.tauGbeta5/v.beta5;}, beta5: function(v){return v.beta5;}*/};

function computeCost(v,e) { // Compute constants var c = arsenicModelConstants; for(var a in mutateFuncs) c[a] = mutateFuncs[a](v);

// Go through all experiments var cost = 0, weight = 0, x0, xt, times; for(var i in e) { // Set up constants for this experiment var nc = {}; for(var a in c) nc[a] = c[a]; for(var a in e[i].constants) nc[a] = e[i].constants[a];

if (e[i].AsT!=undefined) { // Vary time, with fixed AsT // Simulate x0 = arsenicModelInitialization(nc,e[i].AsT); xt = simulate(x0,e[i].data.time,function(t,d){return arsenicModelGradient(nc,d);});

// Sum of (squares of) errors, divided by the average value var curcost = 0, n = 0; for(var xn in e[i].data) { if (xn=='time') continue; var avgv = 0; for(var j in e[i].data[xn]) avgv += e[i].data[xn][j]; avgv /= e[i].data[xn].length; for(var j in xt.timeKey) { curcost += Math.pow((e[i].data[xn][j]-xt[xn][xt.timeKey[j]])/avgv,2); n++; }     }      cost += Math.sqrt(curcost/n); // Compute the square root of the average of the squares (RMS) weight++;

// Set last solution e[i].solution = {'cost':Math.sqrt(curcost/n), 'xt':xt}; } else if (e[i].time==Infinity) { // Vary AsT, with equilibrium var avgv = {}; for(var xn in e[i].data) { avgv[xn] = 0; for(var j in e[i].data[xn]) avgv[xn] += e[i].data[xn][j]; avgv[xn] /= e[i].data[xn].length; }     e[i].solution = {'xt':{'AsT':[]}}; var curcost = 0, n = 0; for(var j in e[i].data.AsT) { // Simulate xt = arsenicModelEquilibrium(nc,e[i].data.AsT[j]); e[i].solution.xt.AsT[j] = e[i].data.AsT[j];

// Fill solution for(var xn in xt) { if (e[i].solution.xt[xn]==undefined) e[i].solution.xt[xn] = []; e[i].solution.xt[xn][j] = xt[xn]; }

// Sum of (squares of) errors, divided by the average value for(var xn in e[i].data) { if (xn=='AsT') continue; curcost += Math.pow((e[i].data[xn][j]-xt[xn])/avgv[xn],2); n++; }     }      cost += Math.sqrt(curcost/n); // Compute the square root of the average of the squares (RMS) weight++; e[i].solution.cost = Math.sqrt(curcost/n); } else if (!isNaN(e[i].time)) { // Vary AsT, with t = e[i].time var avgv = {}; for(var xn in e[i].data) { avgv[xn] = 0; for(var j in e[i].data[xn]) avgv[xn] += e[i].data[xn][j]; avgv[xn] /= e[i].data[xn].length; }     e[i].solution = {'xt':{'AsT':[]}}; var curcost = 0, n = 0; for(var j in e[i].data.AsT) { // Simulate x0 = arsenicModelInitialization(nc,e[i].data.AsT[j]); xt = simulate(x0,e[i].time,function(t,d){return arsenicModelGradient(nc,d);}); e[i].solution.xt.AsT[j] = e[i].data.AsT[j];

// Fill solution for(var xn in xt) { if (e[i].solution.xt[xn]==undefined) e[i].solution.xt[xn] = []; e[i].solution.xt[xn][j] = xt[xn][xt[xn].length-1]; }

// Sum of (squares of) errors, divided by the average value for(var xn in e[i].data) { if (xn=='AsT') continue; curcost += Math.pow((e[i].data[xn][j]-xt[xn][xt[xn].length-1])/avgv[xn],2); n++; }     }      cost += Math.sqrt(curcost/n); // Compute the square root of the average of the squares (RMS) weight++; e[i].solution.cost = Math.sqrt(curcost/n); } }  return cost/weight; // Take the average of the RMS values for all graphs, making it "easier" to disregard certain experiments in favour of the rest. }

function randomLogNormal(mu,sigma) { var N = Math.random+Math.random+Math.random+Math.random+Math.random+Math.random - (Math.random+Math.random+Math.random+Math.random+Math.random+Math.random); return Math.exp(mu+sigma*N); }

function mutate(c,dc) { var vn = varsToMutate[Math.floor(Math.random*varsToMutate.length)]; var nc = {}; for(var a in c) nc[a] = c[a];

// Mutate /*var factor = 1+0.01*(1-Math.exp(-Math.random)); if (Math.random<0.5+Math.atan(dc[vn])/Math.PI) { factor = 1 / factor; }*/ var sigma = 0.1; var factor = randomLogNormal(0,sigma); nc[vn] *= factor; return nc; }

function fitConstants { // Construct plots //constructPlot('v5K5plot'); constructPlot('k8K7plot');

// Show mathematica solution var orgC = arsenicModelConstants; var cSol = {}; for(var i in varsToMutate) cSol[varsToMutate[i]] = 1; //cSol.v5_K5 = orgC.v5/orgC.K5; cSol.v5 = orgC.v5*2; cSol.k8 = 10; cSol.k8_K7 = 2e5; cSol.tauBbeta4 = 55; cSol.beta4 = 18; cSol.tauRbeta1_tauBbeta4 = 400; cSol.beta1_beta4 = 2; // cSol.tauBbeta4 = 180000; // cSol.tauB = 180; // cSol.beta4 = 1000; // cSol.tauR = 60; // cSol.beta1 = 1000; // cSol.tauFbetaF = 120000; // cSol.tauF = 60; // cSol.betaF = 2000; // cSol.tauKbetaK = 9240; // cSol.tauK = 60; // cSol.betaK = 154; // cSol.tauGbeta5 = 3960; // cSol.tauG = 60; // cSol.beta5 = 66; showOutputs('sol',computeCost(cSol,experiments),cSol);

// Initialize var c = {}; for(var i in varsToMutate) c[varsToMutate[i]] = 1; //c.v5_K5 = orgC.v5/orgC.K5; c.v5 = orgC.v5*2; c.k8 = 20; c.k8_K7 = 4e5; c.tauBbeta4 = 55; c.beta4 = 18; c.tauRbeta1_tauBbeta4 = 400; c.beta1_beta4 = 2; // cSol.tauBbeta4 = 180000; // cSol.tauB = 180; // cSol.beta4 = 1000; // cSol.tauR = 60; // cSol.beta1 = 1000; // cSol.tauFbetaF = 120000; // cSol.tauF = 60; // cSol.betaF = 2000; // cSol.tauKbetaK = 9240; // cSol.tauK = 60; // cSol.betaK = 154; // cSol.tauGbeta5 = 3960; // cSol.tauG = 60; // cSol.beta5 = 66; var dc = {}; for(var a in c) dc[a] = 0; var E = computeCost(c,experiments); var cBest = c, EBest = E; for(var i in experiments) experiments[i].bestSolution = experiments[i].solution;

// Show initial situation showOutputs('cur',E,c,dc); showOutputs('',EBest,cBest); refreshGraphs;

// Set up iteration var numiter = 100000; var iter = 0; var timer = setInterval(function{   iter++;    if (iter>numiter) {      clearInterval(timer);      return;    }    setOutput('iter',iter);

// Mutate and compute new energy and gradient var cNew = mutate(c,dc); var ENew = computeCost(cNew,experiments); for(var a in cNew) { var dca = (ENew-E)/(cNew[a]-c[a]); if (!(isNaN(dca) || !isFinite(dca))) dc[a] = (dc[a]+2*dca)/3; }

// If better than best, accept if (ENew < EBest) { cBest = cNew; EBest = ENew; for(var i in experiments) experiments[i].bestSolution = experiments[i].solution; showOutputs('',EBest,cBest); refreshGraphs; }

// Compute (decaying) "temperature" and accept new solution as current if it's not "too" bad var T = 1 - (iter/numiter); if (ENew=Math.random) { c = cNew; E = ENew; showOutputs('cur',E,c,dc); } },1); }

function refreshGraphs { //document.getElementById('Meng2004Graph').refresh; //document.getElementById('Singh2008Graph').refresh; //document.getElementById('Kostal2004fig3BGraph').refresh; document.getElementById('pSB1A2timeGraph').refresh; //document.getElementById('pArsRRFPtimeGraph').refresh; //document.getElementById('Kostal2004fig3AGraph').refresh; document.getElementById('pSB1A2conGraph').refresh; //document.getElementById('pArsRRFPconGraph').refresh; }

function showOutputs(mode,E,c,dc) { //plotMin(v5K5plot,mutateFuncs.v5(c),mutateFuncs.K5(c),E); plotMin(k8K7plot,mutateFuncs.k8(c),mutateFuncs.K7(c),E); for(var a in c) { setOutput(a+mode,c[a]); } for(var a in mutateFuncs) { setOutput(a+mode,mutateFuncs[a](c)); } setOutput('E'+mode,E); if (dc!=undefined) { for(var a in dc) { setOutput(a+mode+'gradient',dc[a]); } } }

function constructPlot(id) { var width = 100, height = 100; var t = document.getElementById(id); t.minx = Number.NaN; t.miny = Number.NaN; t.maxx = Number.NaN; t.maxy = Number.NaN; t.points = []; t.createCaption; t.style.width = width + 'px'; t.style.width = height + 'px'; t.style.border = 'solid 1px #000'; t.style.borderCollapse = 'collapse'; for(var r=0; rt.maxx) { t.maxx = x*1.5; regrid = true; } if (isNaN(t.miny) || yt.maxy) { t.maxy = y*1.5; regrid = true; } if (regrid==true) { //alert('regridding' + [x,y,t.minx,t.miny,t.maxx,t.maxy,regrid]); setCaption(t,'x = ['+formatNumberToHTML(t.minx,3)+','+formatNumberToHTML(t.maxx,3)+'] y = ['+formatNumberToHTML(t.miny,3)+','+formatNumberToHTML(t.maxy,3)+']'); for(var r=0; r<t.rows.length; r++) { var row = t.rows[r]; for(var c=0; c<row.cells.length; c++) { var cell = row.cells[c]; cell.background = '#fff'; }   }    for(var i in t.points) plotMinWork(t,t.points[i].x,t.points[i].y,t.points[i].v); } else { plotMinWork(t,x,y,v); } }

function plotMinWork(t,x,y,v) { var r = Math.floor((y-t.miny)/(t.maxy-t.miny)*t.rows.length); var c = Math.floor((x-t.minx)/(t.maxx-t.minx)*t.rows[0].cells.length); var cell = t.rows[r].cells[c]; if (cell.value==undefined || v0) { caps[0].innerHTML = cap; return; } if (t.caption) { t.caption = cap; return; } }

 