Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2007

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

Search the Archive

Fitting coupled differential equations to experimental data

  • To: mathgroup at smc.vnet.net
  • Subject: [mg83812] Fitting coupled differential equations to experimental data
  • From: "Ausman, Kevin" <ausman at okstate.edu>
  • Date: Sat, 1 Dec 2007 05:43:26 -0500 (EST)

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 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}

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}]])

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>

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}]

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.

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 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!

Kevin Ausman

-------------------------
Prof. Kevin Ausman
Oklahoma State University
Department of Chemistry
342 Physical Science
Stillwater, OK 74078

Email:  ausman at okstate.edu
Office: 405-744-4330
Fax:    405-744-6007
--------------------------


  • Prev by Date: Re: Using ReadList to read a string
  • Next by Date: Re: Using ReadList to read a string
  • Previous by thread: Re: Using ReadList to read a string
  • Next by thread: Re: Fitting coupled differential equations to experimental data