How do I fit coupled differential equations to experimental data?

I have a set of differential equations that I would like to fit to experimental data to get chemical kinetics parameters out of. The equations are:

d(Mzb(t))/dt = -[Mzb(t)-Mzb0]/T1b - kb*Mzb(t) + kf*Mzf(t);

d(Mzf(t))/dt = -[Mzf(t)-Mzf0]/T1f - kf*Mzf(t) + kb*Mzb(t);

The data I have is of Mzb(t) and Mzf(t) at different time points. From that data, I can glean Mzb0 and Mzf0 (the equilibrium values Mzb(t) and Mzf(t), which will be constants), and from other data I know the constants T1b and T1f. I also know that this will be a least-squares optimization to get the best fit for the kb and kf parameters. The part I don't understand is how to code all of this-- I'm a novice with the syntax and would appreciate any and all help about it.

Answer 1

[Star Strider]

Two possibilities: Monod Kinetics and Curve Fitting and Optimizing a Simulation or Ordinary Differential Equation. Those should give you the information you need.

[CP]

Star Strider, thank you for pointing me in the direction of the link you had answered earlier. I have been trying for a couple days to wrap my head around the syntax and to adapt your function file to my problem. This is what I have so far (I also noticed an error in my equations, where I used the initial values, Mzb0 and Mzf0, when I should have used the equilibrium values of Mzbeq and Mzfeq instead):

function S = Mzbt(B,t)
%Mzbt codes the system of differential equations describing the
%chemical exchange seen in the transfer of magnetization of NMR
%active nuclei, as expressed in the Bloch-McConnell equations.
%
%d(S) = -(S-Mzbeq)/T1b - kb*S + kf*X;
%d(X) = -(X-Mzfeq)/T1f - kf*X + kb*S;
% with
% Variables: x(1) = S = Mzb(t), x(2) = X = Mzf(t);
% Parameters: B(1) = Mzbeq, B(2) = Mzfeq, B(3) = kb, B(4) = kf
global T1b T1f Mzb0 Mzf0;
T1b = 1.57;
T1f = 2.5;
Mzb0 = 301638;
Mzf0 = -122060;
x0 = [Mzb0; Mzf0];
[T,Sv] = ode45(@DiffEq, t, x0);

     function dS = DiffEq(t,x);
         xdot(1) = -(x(1)-B(1))/T1b - B(3).*x(1) + B(4).*x(2);
         xdot(2) = -(x(2)-B(2))/T1f - B(4).*x(2) + B(3).*x(1);
         dS = xdot;
     end
S = Sv(:,1);
  end

I have been hit with errors when trying to implement this code, the first being that when ode45 has a function as its first argument, tspan must have two elements. When seeing that, I changed t to [0 1000], and I get the following errors:

Error using odearguments (line 90) MZBT/DIFFEQ must return a column vector

Error in ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn...

Error in Mzbt (line 21) [T,Sv] = ode45(@DiffEq,[0 1000],x0):

I appreciate any help you can offer in this-- I have been struggling (but learning), and would appreciate a push in the right direction. Thanks!

[Star Strider]

My pleasure.

First, there are several ways to be sure your ODE function returns a column vector. Probably the easiest is to add a second dimension reference to your equations:

         xdot(1,:) = -(x(1)-B(1))/T1b - B(3).*x(1) + B(4).*x(2);
         xdot(2,:) = -(x(2)-B(2))/T1f - B(4).*x(2) + B(3).*x(1);

Those will return a column vector. Another way is to have as the first assignment statement in your function:

xdot = zeros(2,1);

choose the one you like best!

Note that you can pass several parameters to your ODE function, including the initial conditions, that will produce the best fit to your data. You just have to be sure to parse them out appropriately within your objective function (here ‘Mzbt’).

Second, I discourage your using global variables. It is best to pass them as parameters in your function:

function S = Mzbt(B,t,T1b,T1f,Mzb0,Mzf0)

Then in your script file that calls lsqnonlin or nlinfit, use an anonymous function to be sure the fitting function only ‘sees’ the arguments of the function you want it to:

objfcn = @(B,t) Mzbt(B,t,T1b,T1f,Mzb0,Mzf0);

You can use that as your objective function, or include it as the appropraite argument in, for example, lsqcurvefit:

B = lsqcurvefit(@(B,t) Mzbt(B,t,T1b,T1f,Mzb0,Mzf0), ... );

The function will pick up the other arguments from your workspace, since you have already defined them.

This is not a straightforward problem, as you have discovered. I will be here to help you as much as you need.