Reproducing results of paper gives inaccurate graph

I have a project that involves trying to reproduce the results of a paper. When I run CancerModel.m, I get a picture that looks like the one in the paper, but the values on the axes are off. I'm plotting cell kill rate vs. time of death, and the paper produces a much higher lifespan and a much lower optimal cell kill rate. I've checked my equation for accuracy, and they appear to be consistent. I used a value of beta (rate of cancer cell growth retardation) that corresponds to cells/day instead of cells/hour, to be consistent with the first graph in the paper, and the

if i > 4
  endT(1, i+1) = te;
end

is to prevent errors (I don't know why they occur). If you remove that code, te (the time when the patient dies for each lambdaCt (cancer cell kill rate)) has no value for lambdaCt = 0 to 4. Also, complex numbers are appearing in x and te for no apparent reason.

Here are my three files:

CancerModel.m:

endT = zeros(1, 41);
  allLambdaCt = [0:40]
  for i = [0:40]
  lambdaCt = i;
  x0 = [1 0];
  tspan = [0 19200];
  xtot = x(1) + x(2);
  options = odeset('Events', @critCells);
  f = @(t, x) cancerEqs(t, x, lambdaCt);
  [t,x,te,xe,ie] = ode45(f, tspan, x0, options);
  x
  te
  if i > 4
      endT(1, i+1) = te;
  end
end
plot(allLambdaCt, endT);

cancerEqs.m:

function xprime = cancerEqs(t, x, lambdaCt)
%cancerEqs: Computes the derivatives involved in solving the system of equations modeling chemotherapy equations.
beta = 59.28 * 10^-4;
Ninf = 2 * 10^12;
tau1 = 10^-6;
tau2 = 10^-6;

xprime = [-beta * log((x(1) + x(2)) / Ninf) * (x(1) - lambdaCt * x(1) + tau2 * x(2) - tau1 * x(1)); -beta * log((x(1) + x(2)) / Ninf) * (x(2) + tau1 * x(1) - tau2 * x(2))];

critCells.m

function [value, isterminal, direction] = critCells(t,x);
% The event occurs when the number of cells increases through the critical level.
value = x(1) + x(2) - 5 * 10^11;
isterminal = 1;
direction = 0;

The paper is attached. I'm using the same equations and parameters.

What exactly is the problem?

EDIT: Using ode15s instead of ode45 produces a better survival time, but doesn't improve the optimal cell kill.

Answer 1

[Star Strider]

I didn’t see anything in your code that differed in any way from what the paper described. I didn’t read the paper as exhaustively as you must have, but I did read through it looking for a ‘Materials and Methods’ section, and didn’t find even the briefest mention of one. As you have seen with ode45 and ode15s, different solvers can give different results. The paper didn’t specify what solver they used. (It may have been the relatively popular and free Berkeley Madonna from U.C. Berkeley, but that’s simply a guess on my part. That is a fixed-step solver as best I can determine, but I haven’t examined it in detail because I have the MATLAB solvers.)

As for your getting complex results for x, the only possibility I can think of that could cause that is for (x(1) + x(2)) to become negative at some point. (I caution against using ‘i’ and ‘j’ as loop indices, because MATLAB uses them for its imaginary operators, though I don’t believe that’s causing those problems.)

I suggest you also look up the papers in the references by Gaffney, and Piccart-Gebhart that may help elucidate the differences you found. You can also write the authors.

That a paper is published in a well-respected, refereed journal does not guarantee that it is possible to reproduce the results, especially if there is no mention of methods.

[Katherine]

Thank you for your response! The solution does vary wildly among MATLAB's ode solvers, and I can imagine it would vary with different programs as well. I've looked at some other papers, and the difference appears to be because of the solvers. I suppose I've shown the general point of the paper, though, and that should be enough.

The complex results are now all in the form a + 0i, so they're not causing problems.

Thanks!

[Star Strider]

My pleasure!

If the results have the form ‘a + 0i’, the imaginary component may be on the order of eps. Since you mentioned that this occurs with x, just use real(x) to eliminate them. It’s highly unlikely you’re deleting any meaningful information.

[Katherine]

That worked, thanks!

[Star Strider]

My pleasure!

Answer 2

[Sean de Wolski]

Have you tried using ode15s rather than ode45? The presence of log in your function indicates the function could be stiff (as are some of the plots in the paper).

[Katherine]

That gives a much better maximum survival time, but the optimal cell kill is still too high (about 5 vs .9).

(Also, thank for the quick response!)