4
$\begingroup$

I'm trying to figure out how to use FindFit with a multivariable differential equation model and data. I've successfully made it work for the one-variable version of the model by following the example on the Help page: 

DNA = 10;

model[ a_?NumberQ, b_?NumberQ, c_?NumberQ, d_?NumberQ, f_?NumberQ, 
       g_?NumberQ, Km1_?NumberQ, Km2_?NumberQ, Km3_?NumberQ, NTP0_?NumberQ] := 
  (model[a, b, c, d, f, g, Km1, Km2, Km3, NTP0] = 
     First[MG /. NDSolve[{ MG'[t] == a*DNA*NTP[t]/(Km1 + NTP[t]) - b MG[t],
                           NTP'[t] == -f*a*DNA*NTP[t]/(Km1 + NTP[t]) - 
                           d MG[t] NTP[t]/(Km2 + NTP[t]) - c NTP[t]/(Km3 + NTP[t]), 
                           GFP'[t] == g*d MG[t] NTP[t]/(Km2 + NTP[t]),
                           NTP[0] == NTP0, MG[0] == 0, GFP[0] == 0},
                         {MG, NTP, GFP}, {t, 0, 800}, Method -> StiffnessSwitching]]);

fit = FindFit[ data, {model[a, b, c, d, f, g, Km1, Km2, Km3, NTP0][t], 
                       a > 0, b > 0, c > 0, d > 0, f > 0, g > 0, 
                       Km1 > 1000, Km2 > 1000, Km3 > 1000, NTP0 > 100000}, 
               {{a, 6.8}, {b, 0.012}, {c, 247}, {d, 1.54}, {f, 19.6}, {g, 22.2}, 
                {Km1, 352200}, {Km2, 127882}, {Km3, 5134.5}, {NTP0, 611628}}, t]

where the data looks like: 

data = {{1.65, 111}, {4.65, 141}, {7.65, 130}, {10.65, 247}, {13.65, 301}, 
        {16.65, 395}, {19.65, 444}, {22.65, 652}, ...};

But now I'd like to do it with DNA being an additional variable included in the data like: 

newdata = {{1.65, 10, 111}, {4.65, 10, 141}, {7.65, 10, 130}, ..., {1.65, 5, -4},
           {4.65, 5, 118}, {7.65, 5, 86}, {10.65, 5, 85}, {13.65, 5, 110}, ...};

so that I could fit multiple curves with different values of DNA simultaneously. I imagine that this is something that's possible, but I'm not sure on the syntax. Anyone have any thoughts on this?

------ EDIT --------

so now I've tried to follow the example that bobthechemist gave on the linked page, but I think I'm getting hung up on the syntax:

model[ c_?NumberQ, d_?NumberQ, f_?NumberQ, Km1_?NumberQ, Km2_?NumberQ, Km3_?NumberQ, 
Km4_?NumberQ ][ DNA_?NumberQ, t_?NumberQ] := 
(model[c,d,f, Km1,Km2,Km3,Km4][t,DNA] =  First[MG/.ParametricNDSolve[{
MG'[t,DNA]==a*DNA*NTP[t,DNA]^n/(Km1^n+NTP[t,DNA]^n)b(Km4^n/(Km4^n+NTP[t,DNA]^n))MG[t,DNA],
NTP'[t,DNA]==-a*f*DNA*NTP[t,DNA]^n/(Km1^n+NTP[t,DNA]^n)-d*MG[t,DNA]NTP[t,DNA]^n/(Km2^n+NTP[t,DNA]^n)-c NTP[t,DNA]^n/(Km3^n+NTP[t,DNA]^n),
NTP[0]==NTP0,MG[0]==0}/.{n->1,b->0.012,a->3.5,NTP0->1500000},{MG,NTP},{t,0,800},{DNA},Method->StiffnessSwitching]]);

fit=FindFit[newdata10,{model[c,d,f, Km1,Km2,Km3,Km4][t,DNA],c>0,0<d,0<f,Km1>100000,Km2>100000,Km3>100000,Km4>100000},{{c,91.0400},{d,8.4986},{f,0.000018697},{Km1,1000100},{Km2,5005020},{Km3,5000150},{Km4,7000000}},{t,DNA},Method->"NMinimize"]

gives a whole bunch of errors. perhaps this kind of problem is a bit beyond someone with my limited grasp of Mathematica syntax

------ EDIT 2 --------

a more complete set of data: 

data = {{2.65,5,86}, {5.65,5,85}, {8.65,5,110}, {11.65,5,153}, {14.65,5,187}, {17.65,5,293}, {20.65,5,321}, {23.65,5,320}, {26.65,5,402}, {29.65,5,355}, {32.65,5,593}, {35.65,5,589}, {38.65,5,653}, {41.65,5,687}, {44.65,5,752}, {47.65,5,858}, {50.65,5,882}, {53.65,5,933}, {56.65,5,1033}, {59.65,5,1043}, {62.65,5,1144}, {65.65,5,1178}, {68.65,5,1239}, {71.65,5,1264}, {74.65,5,1317}, {77.65,5,1452}, {80.65,5,1449}, {83.65,5,1465}, {86.65,5,1480}, {89.65,5,1500}, {92.65,5,1529}, {95.65,5,1531}, {98.65,5,1676}, {101.65,5,1626}, {104.65,5,1632}, {107.65,5,1699}, {110.65,5,1560}, {113.65,5,1651}, {116.65,5,1756}, {119.65,5,1767}, {122.65,5,1715}, {125.65,5,1716}, {128.65,5,1715}, {131.65,5,1732}, {134.65,5,1705}, {137.65,5,1740}, {140.65,5,1759}, {143.65,5,1698}, {146.65,5,1653}, {149.65,5,1628}, {152.65,5,1677}, {155.65,5,1711}, {158.65,5,1608}, {161.65,5,1670}, {164.65,5,1481}, {167.65,5,1563}, {170.65,5,1562}, {173.65,5,1588}, {176.65,5,1540}, {179.65,5,1480}, {182.65,5,1462}, {185.65,5,1424}, {188.65,5,1446}, {191.65,5,1412}, {194.65,5,1380}, {197.65,5,1341}, {200.65,5,1338}, {203.65,5,1263}, {206.65,5,1244}, {209.65,5,1237}, {212.65,5,1164}, {215.65,5,1050}, {218.65,5,1109}, {221.65,5,1041}, {224.65,5,1071}, {227.65,5,908}, {230.65,5,940}, {233.65,5,1013}, {236.65,5,913}, {239.65,5,976}, {242.65,5,886}, {245.65,5,847}, {248.65,5,819}, {251.65,5,784}, {254.65,5,818}, {257.65,5,815}, {260.65,5,807}, {263.65,5,704}, {266.65,5,705}, {269.65,5,816}, {272.65,5,758}, {275.65,5,757}, {278.65,5,633}, {281.65,5,708}, {284.65,5,675}, {287.65,5,632}, {290.65,5,617}, {293.65,5,621}, {296.65,5,594}, {299.65,5,558}, {2.65,10,130}, {5.65,10,247}, {8.65,10,301}, {11.65,10,395}, {14.65,10,444}, {17.65,10,652}, {20.65,10,701}, {23.65,10,840}, {26.65,10,922}, {29.65,10,1074}, {32.65,10,1154}, {35.65,10,1209}, {38.65,10,1326}, {41.65,10,1470}, {44.65,10,1628}, {47.65,10,1600}, {50.65,10,1679}, {53.65,10,1759}, {56.65,10,1856}, {59.65,10,1887}, {62.65,10,2057}, {65.65,10,2078}, {68.65,10,2182}, {71.65,10,2128}, {74.65,10,2034}, {77.65,10,2257}, {80.65,10,2337}, {83.65,10,2362}, {86.65,10,2330}, {89.65,10,2423}, {92.65,10,2471}, {95.65,10,2440}, {98.65,10,2388}, {101.65,10,2544}, {104.65,10,2436}, {107.65,10,2538}, {110.65,10,2402}, {113.65,10,2406}, {116.65,10,2423}, {119.65,10,2365}, {122.65,10,2345}, {125.65,10,2391}, {128.65,10,2412}, {131.65,10,2375}, {134.65,10,2309}, {137.65,10,2321}, {140.65,10,2389}, {143.65,10,2212}, {146.65,10,2211}, {149.65,10,2276}, {152.65,10,2188}, {155.65,10,2112}, {158.65,10,2223}, {161.65,10,1980}, {164.65,10,2046}, {167.65,10,2045}, {170.65,10,2022}, {173.65,10,1933}, {176.65,10,1901}, {179.65,10,1925}, {182.65,10,1829}, {185.65,10,1873}, {188.65,10,1840}, {191.65,10,1855}, {194.65,10,1752}, {197.65,10,1682}, {200.65,10,1639}, {203.65,10,1752}, {206.65,10,1784}, {209.65,10,1661}, {212.65,10,1608}, {215.65,10,1563}, {218.65,10,1462}, {221.65,10,1563}, {224.65,10,1543}, {227.65,10,1448}, {230.65,10,1376}, {233.65,10,1384}, {236.65,10,1383}, {239.65,10,1340}, {242.65,10,1263}, {245.65,10,1362}, {248.65,10,1201}, {251.65,10,1206}, {254.65,10,1220}, {257.65,10,1185}, {260.65,10,1164}, {263.65,10,1133}, {266.65,10,1154}, {269.65,10,1115}, {272.65,10,1122}, {275.65,10,1028}, {278.65,10,1049}, {281.65,10,1042}, {284.65,10,960}, {287.65,10,1011}, {290.65,10,940}, {293.65,10,927}, {296.65,10,877}, {299.65,10,888}};

in each triplet of numbers, the first is the time, second is the DNA parameter value, and the third is the value of the function.

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ This might be related. $\endgroup$
    – bobthechemist
    Commented Sep 16, 2013 at 22:17
  • $\begingroup$ Thanks @bobthechemist! Although now that I've played with the model a bit more, I've come to realize that my problem has more to do with the structure of the equations. It would seem that there are a extremely large number of local minima, or some other reason why the parameters I get are not very good and the solutions are highly dependent on the starting values of the parameters. :( I'm not sure what could be done about this... $\endgroup$
    – dantimatter
    Commented Sep 19, 2013 at 0:50
  • $\begingroup$ You might be able to look for a least squares solution graphically although the solution I presented there was more of a suggestion for exploration... $\endgroup$
    – bobthechemist
    Commented Sep 19, 2013 at 1:56
  • $\begingroup$ hi @bobthechemist, i tried just following the example in the page you linked to, just to see what would happen. no luck. i suspect i've committed a number of syntactical errors. my attempt is shown below the 'EDIT' line above. $\endgroup$
    – dantimatter
    Commented Sep 19, 2013 at 22:57
  • $\begingroup$ Do you have a more complete set of data that can be posted? It would be helpful to have a working sample set even if it is not a complete one. $\endgroup$
    – bobthechemist
    Commented Sep 19, 2013 at 23:47

2 Answers 2

Reset to default
1
$\begingroup$

Not an answer yet, but a step (hopefully) towards solving the problem.

Plotting the same data set provided gives us:

Mathematica graphics

For the two possible values of DNA. The model contains a whopping 10 adjustable parameters which concerns me a bit, given the general lack of features in the data. I think it is worthwhile to explore how the different adjustable parameters affect the shape of the model output, so that's what I'm doing here with ParametricNDSolveValue. I used lowercase functions just for style purposes, but I think I've copied the code correctly without typos.

eqns = {
   mg'[t] == a dna ntp[t]/(km1 + ntp[t]) - b mg[t],
   ntp'[t] == -f a dna ntp[t]/(km1 + ntp[t]) - 
     d mg[t] ntp[t]/(km2 + ntp[t]) - c ntp[t]/(km3 + ntp[t]),
   gfp'[t] == g d mg[t] ntp[t]/(km2 + ntp[t]),
   ntp[0] == 100000, mg[0] == 0, gfp[0] == 0} /. {dna -> 5}
pfun = ParametricNDSolveValue[
  eqns, {mg, ntp, gfp}, {t, 0, 300}, {a, b, c, d, f, g, km1, km2, 
   km3}]
Manipulate[
 Plot[pfun[a, b, c, d, f, g, km1, km2, km3][[1]][t], {t, 0, 300}],
 {a, 1, 10}, {b, 1, 10}, {c, 1, 10}, {d, 1, 10}, {f, 1, 10}, {g, 1, 
  10}, {km1, 1000, 9999}, {km2, 1000, 9999}, {km3, 1000, 9999}]

Mathematica graphics

I played around with manipulate values, taking a through g to have possible values from 1 to 10 and km1, km2, and km3 to have values from 1000 to 10000. In this example, I fixed ntp = 100000 and dna=5. As you can see, it is hard to get a shape that even remotely matches the dataset. In order to move forward with the model fitting, I think a better set of starting values is going to be needed, unless it is possible to simplify the model to one that has fewer parameters.

Improve this answer
$\endgroup$
4
  • $\begingroup$ since first posting this i've been able to independently get values for parameters a (a = 3.5) and b (b = 0.012). so these are fixed. two additional things that can be done: (1) get rid of the equation for gfp[t] (and thus get rid of the parameter g), since we're not fitting the gfp data here, or (2) fit the additional gfp data i've got and the MG data at the same time. since I'm not exactly sure how to fit the gfp and MG data simultaneously, maybe (1) is preferable? if we get rid of a,b, and g, would that help at all? $\endgroup$
    – dantimatter
    Commented Sep 21, 2013 at 20:58
  • $\begingroup$ oh, and NTP0=1500000, also determined through an independent measurement. so if we get rid of the gfp[t] equation the only parameters from the original model that need to be fit are c, d, f, Km1, Km2, and Km3. thanks for everything, @bobthechemist! $\endgroup$
    – dantimatter
    Commented Sep 22, 2013 at 17:10
  • $\begingroup$ @dantimatter, you should be able to alter the code above to delete references to gfp and insert values for a, b and NTP0. When I do this, I still get a solution that doesn't look anything like the experimental data. Possible next steps include: (a) ensuring the model is appropriate and correct; (b) continuing to find parameter values independently; (c) Adjust the Manipulate to see if the model does give a reasonable fit with some range of parameters. Since I don't know this data, model or experiment, it's hard to provide more constructive feedback. Sorry. $\endgroup$
    – bobthechemist
    Commented Sep 22, 2013 at 17:46
  • $\begingroup$ thanks for all your help, @bobthechemist. :) i'll hack at it a bit more and post here again if there's anything new to add to the story. eventually the associated work will be published and when that happens i'd like to acknowledge your help on this. if you're interested, google dantimatter, find my contact info, and send me a note with your information. thanks again. $\endgroup$
    – dantimatter
    Commented Sep 23, 2013 at 1:46
1
$\begingroup$

So it turns out that my institution has a certain amount of Mathematica technical support included in their license. I posed this question to the good people at Wolfram and this was their response:

"I've attached a notebook containing examples of how you might successfully use FindFit with your parametric differential equation system.

The notebook doesn't contain the results of the FindFit evaluations, as they seemed to be taking a bit of time to evaluate. However, the syntax and usage of all expressions in the notebook is correct as far as I have tested and should give you a fit to your model."

The contents of the notebook are below.

ClearAll[DNA];
solutions =
 ParametricNDSolve[
   {Derivative[1, 0][MG][t, DNA] == 
     a*DNA*NTP[t, DNA]/(Km1 + NTP[t, DNA]) - b*MG[t, DNA],
    Derivative[1, 0][NTP][t, DNA] ==
     -f*a*DNA*NTP[t, DNA]/(Km1 + NTP[t, DNA]) - 
 d*MG[t, DNA] NTP[t, DNA]/(Km2 + NTP[t, DNA]) - 
 c*NTP[t, DNA]/(Km3 + NTP[t, DNA]),
    Derivative[1, 0][GFP][t, DNA] == 
g*d*MG[t, DNA] NTP[t, DNA]/(Km2 + NTP[t, DNA]),
    NTP[0, DNA] == NTP0,
    MG[0, DNA] == 0,
    GFP[0, DNA] == 0},
   {MG, NTP, GFP},
   {t, 0, 800},
   {DNA, 5, 10},(* 
   chosen as example domain *)
   {a, b, c, d, f, g, Km1, Km2, Km3, 
    NTP0},
   Method -> "StiffnessSwitching"
   ]

ClearAll[model];
model = MG /. solutions;

parameters = {a, b, c, d, f, g, Km1, Km2, Km3, NTP0};

startingValues = {6.8`, 0.012`, 247, 1.54`, 19.6`, 22.2`, 352200, 127882, 
   5134.5`, 611628};

(* quick test *)
model @@ startingValues
%[100, 7]

fit =
 FindFit[
  data,
  {model[a, b, c, d, f, g, Km1, Km2, Km3, NTP0][t, DNA],
   {a > 0, b > 0, c > 0,
    d > 0, f > 0, g > 0,
    Km1 > 1000, Km2 > 1000, Km3 > 1000,
    NTP0 > 100000}},
  Thread@{parameters, startingValues},
  {t, DNA}
  ]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.