Estimating parameters with ODE system with changing variables

Hello!

I am working on a ODE system for a disease dynamics that includes, at a certain point, the deployment of a vaccine.

The parameters I am trying to estimate are three infection rates and the deployment of the vaccine. In order to be more biologically accurate, the system starts with just one infected element.

The output data to be compared with the field data is x(5)/(x(4)+x(5))

My code is:

File 1

global coff fin

time = [2007; 2008; 2009; 2010; 2011]; 
nip_r = [.53;.55;.45;.14;.13]; % field data

coff=2007; %2007 is the first year of vaccine deployment
fin=2020;
b0=[0.391;0.391;0.535;0.5]; 

[b_fit,resnorm,exitflag,jacobian,message] = lsqcurvefit(@V_ODE_fit_gg,b0,time,nip_r,[],[],options)

File 2

function M = V_ODE_fit_gg(b, time, x)
global a1 a2 a3 a4 a5 gm p1 p2 coff fin
a1=4;
a2=1;
a3=1;
a4=4;
a5=1;
gm=0.025;
p1=1/30;
p2=0.3;
t0=[0:coff];
t1=[coff:fin];
x0 = [100,400,0,400,0,100,100,1,0];

[~,Dv] = ode15s(@DifEq1,t0,x0);

function s=DifEq1(~, x)
    xp=ones(9,1);
    xp(1)=a4*x(6) -a1*x(1);
    xp(2)=(a1*x(1))-x(2)*(b(1)*(x(8)/(x(7)+x(8)+x(9)))+a2);

      xp(3)=x(2)*(b(1)*(x(8)/(x(7)+x(8)+x(9))))-a2*x(3);
      xp(4)=(a2*x(2))-x(4)*(b(2)*(x(8)/(x(7)+x(8)+x(9)))+a3);
      xp(5)=a2*x(3)+(b(2)*(x(8)/(x(7)+x(8)+x(9)))*x(4))-a3*x(5);
      xp(6)=a3*(x(4)+x(5))-a4*x(6);
      xp(7)= (a5*(x(7)+x(8)+x(9)))x(7(x(5)/(x(7)+x(8)+x(9)))*b(3))-a5*x(7);
      xp(8)=x(7)*b(3)*(x(5)/(x(7)+x(8)+x(9)))-gm*x(8)-a5*x(8);
      xp(9)=gm*x(8)-a5*x(9);

      s=xp;
      end

D=[Dv(end,1) Dv(end,2) Dv(end,3) Dv(end,4) Dv(end,5) Dv(end,6) Dv(end,7) Dv(end,8) Dv(end,9) 0];

[~,Mv]= ode15s(@DifEq2,t1,D);

      function s=DifEq2(~, x)
      xp=ones(10,1);

      xp(1)=a4*x(6) -a1*x(1); 

      xp(2)=(a1*x(1))-x(2)*(b(1)*(x(8)/(x(7)+x(8)+x(9)))+a2);

      xp(3)=x(2)*(b(1)*(x(8)/(x(7)+x(8)+x(9))))-a2*x(3);
      xp(4)=(a2*x(2))-x(4)*(b(2)*(x(8)/(x(7)+x(8)+x(9)))+a3);
      xp(5)=a2*x(3)+(b(2)*(x(8)/(x(7)+x(8)+x(9)))*x(4))-a3*x(5);
      xp(6)=a3*(x(4)+x(5))-a4*x(6);
      xp(7)= (a5*(x(7)+x(8)+x(9)+x(10)))-(x(7)*(x(5)/(x(7)+x(8)+x(9)+x(10)))*b(3))-b(4)*x(7)-a5*x(7);
      xp(8)=b(3)*(x(5)/(x(7)+x(8)+x(9)+x(10)))*(x(7)+b(4)*p1*x(10))-gm*x(8)-a5*x(8);
      xp(9)=gm*x(8)+b(4)*p2*x(10)-a5*x(9);
      xp(10)=b(4)*x(7)-b(4)*p1*b(3)*(x(5)/(x(7)+x(8)+x(9)+x(10)))*x(10)-b(4)*p2*x(10)-a5*x(10);

      s=xp;
      end

      I = 0*time;
      for i=1:length(time)
      [C I(i)] = min(abs(t1-time(i)));
      end

      g=Mv(I,5)/(Mv(I,4)+Mv(I,5));
      M=g(:,4);

end

Running lsqcurvefit, the results obtained are the sames as the initial estimates and the plot is a flat line. Could someone help me with this?

tnks

[Star Strider]

Three observations:

1. Avoid global calls, since global variables cause more problems than they solve. Pass extra parameters (using anonymous functions where necessary), instead.
2. You don’t pass your time variable to your ode15s call. You pass it ‘t0’ instead, although I can’t tell if that is important.
3. The I = 0*time; assignment sets ‘I’ to be a vector of zeros. I’m surprised that did not throw an error in your ‘g’ assignment.

[André_88]

Thanks for the suggestions!

I will try to reduce on the global calls

Using t0 does not create any problems on a similar model I have ran but I will try to use this suggestion.

Concerning point 3, do you suggest any modification?

Cheers

[Star Strider]

My pleasure.

It’s difficult for me to follow the structure of your code. However, one reason parameters don’t change in a regression is that the initial estimates are too far from the correct ones, and the parameter ‘surface’ is flat, so the gradient descent algorithm can’t find a direction to the global minimum. Experiment with different initial estimates. You have five data pairs and four parameters, making this a challenge.

When using ODEs as part of a regression objective function, I always test them first with realistic parameters to be certain they produce the sort of output that makes sense in the context of the problem, and that is appropriate for the regression I want to use them in. Do they produce the appropriate output with your test parameters? Does ‘g’ look as it should?

The problem with globals is that other code can change them, and you can lose track of those changes. If you pass those variables as arguments instead, their values are fixed (unless you change them in your function). You always know what they are, nevertheless.

W.R.T. #3., I misread your code. It’s an appropriate preallocation, just not a structure I’m used to seeing.

[André_88]

I will try to change the initial estimates then, to see if it makes any change in the results and in the plot.

I've already tried to use the ODE system with realistic parameters and the results do indeed make sense. 'g' makes sense indeed.

[Star Strider]

Good.

If your ODE produces reliable results, experiment with the initial estimates to your regression. If you have the Global Optimization Toolbox, consider using genetic algorithms if you can’t easily estimate a reliable set of parameters.