<<xlr8r.m;

Suppose that we have a time-series of values for two species, which looks like this:

data = {a->{1,1,1,2,2,2,2,3,3,3,3,3,4,4,4,4,4,5,5,5,5,6,6,6,7,7,8,8,9,10,11,11,12,13,14,15,16,16,17,18,18,19,19,20,20,20,20,21,21,21,21,21,21,21,21,20,20,20,20,20,19,19,19,19,18,18,18,18,17,17,17,17,16,16,16,16,15,15,15,15,15,14,14,14,14,13,13,13,13,13,12,12,12,12,12,12,11,11,11,11,11,10,10,10,10,10,9,9,9,9,9,8,8,8,8,7,7,7,7,6,6,6,6,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,8,8,8,9,9,10,10,11,11,12,12,13,13,14,14,15,15,15,16,16,16,16,16,16,17,17,17,17,17,16,16,16,16,16,16,16,16,16,15,15,15,15,15,15,14,14,14,14,14,13,13,13,13,13,13,12,12,12,12,12,12,11,11,11,11,11,10,10,10,10,10,10,9,9,9,9,9,8,8,8,8,7,7,7,7,6,6,6,6,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,8,8,8,9,9,10,10,11,11,12,12,13,13,14,14,15,15,15,16,16,16,16,16,16,17,17,17,17,17,16,16,16,16,16,16,16,16,16,15,15,15,15,15,15,14,14,14,14,14,13,13,13,13,13,13,12,12,12,12,12,11,11,11,11,11,11,10,10,10,10,10,9,9,9,9,9,8,8,8,8,8,7,7,7,7,6,6,6,6,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,8,8,8,9,9,9,10,10,11,11,12,13,13,14,14,14,15,15,15,16,16,16,16,16,16,17,17,17,17,16,16,16,16,16,16,16,16,16,15,15,15,15,15,15,14,14,14,14,14,14,13,13,13,13,13,12,12,12,12,12,12,11,11,11,11,11,11,10,10,10,10,10,9,9,9,9,9,8,8,8,8,7,7,7,7,7,6,6,6,6,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,8,8,8,9,9,10,10,11,11,12,12,13,13,14,14,15,15,15,16,16,16,16,16,16,16,17,17,17,17,16,16,16,16,16,16,16,16,16,15,15,15,15,15,15,14,14,14,14,14,14,13,13,13,13,13,12,12,12,12,12,12,11,11,11,11,11,11,10,10,10,10,10,9,9,9,9,9,8,8,8,8,7,7,7,7,6,6,6,6,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,8,8,8,9,9,10,10,11,11,12,12,13,13,14,14,15,15,15,16,16,16,16,16,16,16,17,17,17,17,16,16,16,16,16,16,16,16,16,15,15,15,15,15,15,14,14,14,14,14,13,13,13,13,13,13,12,12,12,12,12,12,11,11,11,11,11,10,10,10,10,10,10,9,9,9,9,9,8,8,8,8,7,7,7,7,6,6,6,6,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6},b->{2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3}};

Suppose further that we're pretty sure they were generated by a model something like this:

model = {
{a↦a,hill[1,9,4,0,15]},
{a↦b,hill[1,8,4,0,1.5]},
{a + b -> b, .2},
{a -> ∅, γ2},
{b -> ∅, γ3},
{∅ -> a, 3},
{∅ -> b, .05}
};

But we don't know values for γ2 or γ3.  What can we do for such a system?

{eqns, vars} = interpret[model];

To give some idea of the variability in such a model, suppose that we fix γ2 at 0.2, and vary γ3 from .04 to .25, a relatively small-seeming range. We can use this to generate a 3D plot with values for γ3 on the x-axis, time on the y-axis, and the value of the predicted "a" concentration on the z.  That looks like this:

Plot3D[
Quiet[NDSolve[Join[eqns /. {γ2->.02,γ3->x} ,{a[0]==1,b[0]==2} ],vars,{t,0,100},MaxSteps->100000]][[1]][[2,2]][y],{x,.04,.25},{y,0,100}]

Another thing you might like to do is pick some values for both γ2 and γ3 and plot predicted values of "a" versus the data. That might look like this:

plotRange[myg2_,myg3_,low_,high_]:=
Module[{myode,myinterpolatingFuncs, d,ymin,ymax},
myode = Join[eqns /. {γ2->myg2,γ3->myg3} ,{a[0]==1,b[0]==2} ];
myinterpolatingFuncs=Quiet[NDSolve[myode,vars,{t,0,100},MaxSteps->100000]][[1]];
d=(a/.data)[[Floor[low*10]+1;;Floor[high*10]+1]];
ymin = Min[d]-.2;
ymax = Max[d]+.2;
Show[
ListPlot[Thread[{Range[low,high,.1], d}],PlotRange->{{low,high},{ymin,ymax}}], 
Plot[(a/.myinterpolatingFuncs)[t],{t,low,high}]]
];


plotRange[.1,.2,0,100]
plotRange[.3,.1,0,100]

Clearly, there is a lot of variability in the model, which depends on the values of these parameters. What we would like is to seek parameter values that produce predictions which match the data. To do this , we create a function that takes parameter values, compares the model with these parameters to the data, and returns the difference. This "objective function" acts as a score for proposed parameter values.

evaluate[proposedγ2_?NumberQ, proposedγ3_?NumberQ]:=
Module[{interpolatingFuncs,predictedA,predictedB,errorA,errorB,ode},
ode = Join[eqns /. {γ2->proposedγ2,γ3->proposedγ3} ,{a[0]==1,b[0]==2} ];
interpolatingFuncs=Quiet[NDSolve[ode,vars,{t,0,100},MaxSteps->100000]][[1]];
predictedA = Table[Quiet[(a/.interpolatingFuncs)[iii]],{iii,0,100,.1}];
predictedB = Table[Quiet[(b/.interpolatingFuncs)[iii]],{iii,0,100,.1}];
errorA =Norm[(predictedA- (a/.data)), 2];
errorB = Norm[(predictedB- (b/.data)), 2];
errorA + errorB 
];

For problems with small dimension, we can plot this objective function. Here it looks like this. Smaller values are better matches with the data. But where is the lowest point?

Plot3D[evaluate[i,j],{i,0,.1},{j,0,.3}]

Fixing each of the parameters in turn, we can plot the behavior of the objective function with respect to the other. But we want to optimize them together.

Plot[evaluate[.02,j],{j,0,.2}]
Plot[evaluate[i,.15],{i,0,.2}]

For this, we turn to simulated annealing. This takes our objective function, and an idea of the valid range for each variable, and attempts to find a globally optimal value for the parameters. The return type of LamDelosmeSA is 
	{{a,b},{c,d},e}
Where:
  a is the value of the objective function at b
  b is the set of variable assignments which were current when the optimization routine ended
  c is the best value of the objective function encoutered
  d is the set of variable assignments associated with objective value c
  e is the number of iterations used

SetDirectory["/shared/uci/courses/284_compbio_c/project/"]; (* you should change this to whatever directory on your computer has the file LamDelosmeSA.m *)
<<LamDelosmeSA.m

rep[s_,e_,i_]:=Print["State<",s,"> has energy <",e,"> at iteration<", i, ">"];
timeres = Timing[LamDelosmeSA[evaluate[g2,g3],{{g2,0,1},{g3,0,1}},{g2≥0,g2≤1,g3≥0,g3≤1},ReportFunc->rep,NumForgetMoves->100,NumInitMoves->100,MaximumFrozenIterations->1]]

results= timeres[[2]]

Having run the optimization algorithm, we can compare the model of "a" given the resulting parameters with the data:

plotRange[g2/.results[[2,2]],g3/.results[[2,2]],0,100]