Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Fitting coupled differential equations to experimental data

  • To: mathgroup at smc.vnet.net
  • Subject: [mg83852] Re: Fitting coupled differential equations to experimental data
  • From: Szabolcs Horvát <szhorvat at gmail.com>
  • Date: Sun, 2 Dec 2007 04:14:03 -0500 (EST)
  • References: <firedg$rq6$1@smc.vnet.net>

Ausman, Kevin wrote:
> I have been beating my head against a problem for several weeks now, so
> I'm at the point of begging for help. Please help!

I don't have your data, so I will use a simplified problem with this data:

data = Table[{x, Exp[2 x]}, {x, 0, 1, .1}]

> 
> I have a set of coupled differential rate equations with several
> parameters (this is a chemical system, and these variable parameters are
> rate constants and initial concentrations):
> 
> {DHR'[t] == -kadD*DHR[t]*C60[t] + kdeD*DHR1[t],
>  DHR1'[t] == kadD*DHR[t]*C60[t] - kdeD*DHR1[t] - kox*DHR1[t],
>  RH1'[t] == kox*DHR1[t] - kdeR*RH1[t] + kadR*RH[t]*C60[t],
>  RH'[t] == kdeR*RH1[t] - kadR*RH[t]*C60[t],
>  C60'[t] == -kadD*DHR[t]*C60[t] + kdeD*DHR1[t] + kdeR*RH1[t] -
>    kadR*RH[t]*C60[t], Fluorescence[t] == ScaleFactor*(RH[t] + =
> RH1[t]),
>   DHR[0] == initDHR, DHR1[0] == 0, RH[0] == initRH, RH1[0] =
> == 0,
>  C60[0] == initC60}
> 
> My variables in this case are:
> 
> variables={DHR, DHR1, RH, RH1, C60, Fluorescence}
> 
> And my parameters are:
> 
> {initC60,initRH,initDHR,ScaleFactor,kadD,kdeD,kox,kadR,kdeR}


Let the diff eq be {f'[x] == a f[x], f[0] == 1}

> 
> I can successfully create a model function at a given set of values for
> these parameters using NDSolve:
> 
> model[kadD_?NumberQ, kdeD_?NumberQ, kox_?NumberQ, kadR_?NumberQ,
>   kdeR_?NumberQ, initDHR_?NumberQ, initRH_?NumberQ, initC60_?NumberQ,
>   ScaleFactor_?
>    NumberQ] := (model[kadD, kdeD, kox, kadR, kdeR, initDHR, initRH,
>     initC60, ScaleFactor] =
>    First[NDSolve[{DHR'[t] == -kadD*DHR[t]*C60[t] + kdeD*DHR1[t],
>       DHR1'[t] == kadD*DHR[t]*C60[t] - kdeD*DHR1[t] - kox*DHR1[t],
>       RH1'[t] == kox*DHR1[t] - kdeR*RH1[t] + kadR*RH[t]*C60[t],
>       RH'[t] == kdeR*RH1[t] - kadR*RH[t]*C60[t],
>       C60'[t] == -kadD*DHR[t]*C60[t] + kdeD*DHR1[t] + kdeR*RH1[t] -
>         kadR*RH[t]*C60[t],
>       Fluorescence[t] == ScaleFactor*(RH[t] + RH1[t]),
>       DHR[0] == initDHR, DHR1[0] == 0, RH[0] == initRH, =
> RH1[0] == 0,
>       C60[0] == initC60}, variables, {t, 0, 7200}]])


Clear[model]

model[a_?NumericQ] :=
  model[a] = First@NDSolve[{f'[x] == a f[x], f[0] == 1}, f, {x, 0, 1}]

It is better to use NumericQ instead of NumberQ, because 
NumberQ[Sqrt[2]] === False but NumericQ[Sqrt[2]] === True.  It is also 
good to avoid symbol names starting with capital letters to prevent 
collisions with built-in symbols.  E.g. ScaleFactor is already used by 
the VectorFieldPlots package.

> 
> solution =
>   model[0.018, 0.024, 0.394, 0.024, 0.018, 0.78, 0, 0.018, 90700];
> Plot[Fluorescence[t] /. solution, {t, 0, 7200}]
> 
> <Graphic>

solution = model[1]
Plot[f[t] /. solution, {t, 0, 1}]

Yes, this works as expected.

> 
> I can even treat the parameters as adjustible using the Manipulate
> operation, plotting it with the experimental data that I want to fit:
> 
> XPData :=
>  Import["C:\Documents and Settings\xxx\My \
> Documents\Briefcase\Research\Kinetic Fitting with \
> Mathematica\DHR123_7.PRN", "Table"]
> XPTime = XPData[[All, 1]];
> XPFluor = XPData[[All, 2]];
> 
> Manipulate[(solution =
>    model[kadD, kdeD, kox, kadR, kdeR, initDHR, initRH, initC60,
>     ScaleFactor];
>   Show[ListPlot[XPData],
>    Plot[Fluorescence[t] /. solution, {t, 0, 7200}]]), {{initC60, 0.1},
>    0, 1}, {{initDHR, 0.1}, 0, 1}, {{initRH, 0}, 0,
>   1}, {{ScaleFactor, 10000}, 0, 1000000}, {{kadD, 0.1}, 0,
>   1}, {{kdeD, 0.1}, 0, 1}, {{kox, 0.1}, 0, 1}, {{kadR, 0.1}, 0,
>   1}, {{kdeR, 0.1}, 0, 1}]


This works also:

Manipulate[solution = model[a];
  Show[ListPlot[data], Plot[f[t] /. solution, {t, 0, 1}]],
  {{a, 1}, 0, 3}]

> 
> What I cannot seem to do is fit this model to experimental data using
> FindFit, or by generating a least-squares error function and using
> Minimize, NMinimize, or FindMinimum.
> 
> For example:
> 
> FindFit[XPFluor, (Model2 =
>    model[kadD, kdeD, kox, kadR, kdeR, initDHR, initRH, initC60,
>     ScaleFactor];
>   Fluorescence[t] /. Model2), {{initC60, 0.018}, {initDHR,
>    1}, {initRH, 0}, {ScaleFactor, 134000}, {kadD, 0.012}, {kdeD,
>    0.118}, {kox, 1}, {kadR, 0.012}, {kdeR, 0.134}}, t]
 >
 > generates ReplaceAll::reps errors and InterpolatingFunction::dmval
 > errors.

The simplified version of this command is

FindFit[data, mdl = model[a]; f[t] /. mdl, {{a, 1.}}, t]

But there are several problems with this.  First, not that FindFit does 
not have any Hold* attributes, so 'mdl = model[a]; f[t] /. mdl' is 
evaluated before passed to FindFit[].  The result of the evaluation is 
f[t] /. model[a], and a ReplaceAll::reps error because model[a] is not a 
list of replace rules.

You can see this if you evaluate 'mdl = model[a]; f[t] /. mdl separately.

The solution is to define
fitfun[a_?NumericQ, t_] := f[t] /. model[a]
and use
FindFit[data, fitfun[a, t], {{a, 1.}}, t]

> 
> LSE[kadD_?NumberQ, kdeD_?NumberQ, kox_?NumberQ, kadR_?NumberQ,
>   kdeR_?NumberQ, initDHR_?NumberQ, initRH_?NumberQ, initC60_?NumberQ,
>   ScaleFactor_?NumberQ, XPT_?VectorQ,
>   XPF_?VectorQ] := (LSE[kadD, kdeD, kox, kadR, kdeR, initDHR, initRH,
>     initC60, ScaleFactor, XPT,
>     XPF] = (solution =
>      model[kadD, kdeD, kox, kadR, kdeR, initDHR, initRH, initC60,
>       ScaleFactor]; M = Fluorescence[XPT] /. solution; Norm[XPF - M]))
> 
> FindMinimum[
>  LSE[kadD, kdeD, kox, kadR, kdeR, initDHR, initRH, initC60,
>   ScaleFactor, XPTime,
>   XPFluor], {{initC60, 0.018}, {initDHR, 1}, {initRH,
>    0}, {ScaleFactor, 134000}, {kadD, 0.012}, {kdeD, 0.118}, {kox,
>    1}, {kadR, 0.012}, {kdeR, 0.134}}]
> 
> generates NDSolve::icfail, First::first, ReplaceAll::reps, and
> FindMinimum::nrnum errors.

I didn't take the time to figure out what went wrong here but

lse[a_?NumericQ, data_] :=
  Norm[fitfun[a, data[[All, 1]]] - data[[All, 2]]]
FindMinimum[lse[a, data], {a, 1.}]

works correctly.

Just one more comment: try not to use global variables inside functions. 
   If you don't want to write rep /. fun[x] because it is too long, then 
use With[{ff = fun[x]}, rep /. ff], but *not* ff = fun[x]; rep /. ff. 
And don't forget to Clear[] on functions before redefining them if you 
use memoization.

> 
> I have tried dozens of other formulations, and can evaluate fits
> one-at-a-time or by using Manipulate, and then evaluate LSE (or other
> formulations) one-at-a-time or by using Manipulate, but I cannot seem to
> get Mathematica to perform an optimization.
> 
> Does anyone have an idea of how to help here? Thanks!

-- 
Szabolcs


  • Prev by Date: Re: Vieta infinite product formula
  • Next by Date: Hiding number cell
  • Previous by thread: Fitting coupled differential equations to experimental data
  • Next by thread: Re: Re: optimization routine using conjugate gradient (or other) method?