ode45 won't produce te or xe

I'm trying to run ode45 with events on a function (cancerEqsVaryBurden). The events script is critCells, and the main script is CancerModelFinite. I'm having trouble in the innermost for loop with getting ode45 to give xe or te. In my experience this has been caused by running through tspan without triggering an event, but that shouldn't be the case here. I can increase tspan as much as I want and it still throws the error.

Attached are the relevant scripts plus others to make the main code work.

Thanks for helping!

Answer 1

[Star Strider]

Your code didn’t post, but in your previous post, your ODE didn’t do anything that might trigger an ‘event’, like a zero-crossing.

When you plot the output of your ODE integration, does it look as though it should have triggered an event?

[Katherine]

Whoops! Fixing that now. There's a lot of new code, and it should trigger the event, since that particular equation models cancer cells growing unrestrained. I'm plotting it now to check.

[Star Strider]

I’ll wait for you to troubleshoot. I’m here if you’re still having problems.

[Katherine]

It doesn't appear to respect the option condition at all (the condition is whether number of cells - fatality limit = 0 (ie whether there are enough cells to kill the patient).). The following piece of code is the one where the issue arises:

x0 = afterKill(end,:);
tspan = [0 50000];
options = odeset('Events', @critCells);

f = @(t, x) cancerEqsVaryBurden(t, x);
[t,x,te,xe,ie] = ode45(f, tspan, x0, options);

te and xe are empty vectors at the end of this, but x has satisfied the option. This only happens when dch = 400, however.

[Star Strider]

I’ll look and see if I can figure it out...

You have 4 script files. Which one should I run as the primary?

[Katherine]

CancerModelFinite is the primary, and the the other four are either functions or option conditions :)

[Star Strider]

For loop1 = 3 and loop2 = 1, it works fine for me, producing:

ivtime =

     282.7709e+000

ivx =

       9.9998e+009   230.2515e+003

three times, just after the cancerEqsVaryBurden statement (Line #23 in CancerModelFinite in my Editor), until it gets to:

endT(loop1, loop2) = finalt;

where it throws an error because:

endT = zeros(3, 501);

was initialised as having three columns.

You can either add a dimension to endT here or reduce one in the preallocation statement. Tell me, I’ll make the appropriate changes, and keep working on this.

None of the subsequent events statements triggers, at least to this point.

[Katherine]

Firstly, thank you so much for helping me with this!

I don't quite understand why the error is being thrown, however. What is the issue with endT?

Also, my experiments told me that te and xe were not being assigned values by the last ode45 , so finalt had no value.

[Star Strider]

My pleasure! I’m interested specifically in biomedical simulations and modeling. Besides, this is obviously research and not homework, so there’s a good chance I’ll learn something. Oncology isn’t my speciality.

The finalt line works fine for me. I was printing it to the Command Window to be sure, and finally did CTRL+C to stop it at iteration 423 (so I could put the semicolon back). Some of the events are triggering. I put te1 = te after line #42 and |te2 = te after line #59. With that trace, te1 isn’t triggering, te2 is.

P.S. — I’m now getting:

Subscripted assignment dimension mismatch.

Error in CancerModelFinite (line 64)
        endT(loop1, loop2) = finalt;

for loop1 = 3, loop2 = 1. Do you want to add a third row to the assignment, or delete a row from the preallocation?

[Katherine]

That's really cool! If you don't mind, what is your specialty?

Hmm, that's strange- it hasn't worked for me. I would probably add a row to the assignment, if it doesn't matter much either way.

Thanks!

[Star Strider]

Internal Medicine (Board-Certified). Fellowship in Endocrinology, but decided to go back to grad school and get a M.S. in Biomedical Engineering rather than take my Endocrine boards. No regrets!

When I looked further, it turned out the actual dimensions weren’t the problem, so I didn’t change that line.

The problem with the Subscripted assignment dimension mismatch at loop1 = 3 and loop2 = 1 is that the second te (from cancerEqsVaryBurden) is empty. I’m accumulating te1 (first te) and te2 (second te). The te1 vector has only one entry, but te2 has 1002 when it crashes. The a and maxidx vectors have 200 elements. The value for afterKillt = 508.9397.

Also, running R2014a here on a Win 8 64 laptop. There could be version differences if you’re running a different release or a Mac.

[Katherine]

I'm actually still a student, but I hope to go into medicine eventually. This is an independent project that I'm doing, since I find this really interesting.

That's funny. I see no reason for the second te to be empty, since cancerEqsVaryBurden should trigger the event. It's only a problem for loop1 = 3, not for the others, so I don't think it's a problem with cancerEqsVaryBurden. It's also the same syntax as

I'm running R2013a on a Mac, but I'm not sure how much difference that makes.

[Star Strider]

I take it by ‘student’ you mean undergraduate. Consider an M.D.-Ph.D. program, since I consider what you’re doing to be impressively advanced. You might have an academic career if you want it.

I added these lines after the cancerEqsVaryBurden:

if ~isempty(te2)
    te2v = [te2v; te2 loop1 loop2];
else
    te = 0;
end

to get around the problem (I’m tracking the second te here), and got this plot:

Is that what you were expecting? Should there be two curves or three (or did one overplot another)? After that run, the indices loop1 = 3, loop2 = 501, te1v has 44 rows, and te2 has 1468 rows. Only the second te is ever added to finalt.

[Katherine]

I am an undergraduate! I have a little ways to go before I need to make decisions about what I'm going to do after, but I'm definitely considering an MD-PhD program. For a while, before I knew that existed, I couldn't decide if I wanted to go into science or medicine, but then I realized I could do both!

That is what I'm expecting- there are three curves, but, since the x-axis is a function of loop1, the shortest curve (loop1 = 40) doesn't have as many values as the other two (loop1 = 200 and loop1 = 400). That's fixable easily enough by evaluating at more values of loop1.

When I paste that in, I get the error Undefined function or variable 'te2v'. What exactly does te2v do?

[Star Strider]

I was considering an M.D.-Ph.D. in Pharmacology (B.S. Chemistry major, so I was all set), but my (successful) interview with the department chairman was in his lab, where he was decapitating white mice to do diaphragm preparations. I had a white mouse as a pet when I was a kid, got slightly sick, and my Pharmacology career ended there. Several years later, I was doing my Endocrine fellowship at UCSF, bought an early Apple ][ to do statistics that I taught myself and programmed myself (learned FORTRAN in college, so BASIC was easy), decided that I needed to learn a lot more about computers and statistics, and BME was my best option. Besides, we had one biostatistics lecture in med school ( p < 0.05 is good ), and I quickly realised that my colleagues and I were hopelessly ignorant in quantitative stuff, and for my part, I set out to remedy that. My BME concentration was stats, signal processing, control, and instrumentation. Still learning...

Sorry about tev2. I put:

te2v = [];

just before loop1, but forgot to mention it. The tev2 matrix just keeps track of the second te so I would know when it was executing. (I did the same thing for tev1.) Those lines aren’t absolutely necessary, but setting te = 0 if it returns as empty circumvents the ‘Subscripted assignment dimension mismatch’ problem. That’s the only way I could get it to run all the way through. I hope setting it to zero in finalt didn’t break your code.

BTW, did we solve your problem or are we still working on it? I’m good for another hour or so tonight. Otherwise, I’ll look again tomorrow morning.

[Katherine]

I'm not sure exactly what I want to do eventually, but I do know that I particularly like oncology, genetics, and neuroscience. I'll have to decide eventually, but I like having options! I'm not sure what would give me the best balance of research and application, though.

I fixed that, and I'm currently running the script with a larger range of loop1. A test case with a smaller range worked fine, so I think it's good for now! Setting te = 0 shouldn't be too much of a problem, if the test case is anything to go by.

I think this part of it is done! I have a couple more aspects I'd like to model (such as multiple doses), but I'll probably get to those tomorrow.

Thanks for all your help!

[Star Strider]

My pleasure!

There’s no ‘perfect’ solution. Every choice has trade-offs, but if you want to most closely merge research and application, genetics might be the best you can get. Talk to professionals in your potential fields-of-interest to get an idea of whether it would be a good fit for you.

If the other aspects are part of the same problem, I’ll keep this open.

[Katherine]

I'll definitely consider that, since I do want practical application. Genetics has been a big thing for me for a while, and I think I'd be happy doing that.

Yes, there are a few more aspects :) I think keeping this open would be good!

[Star Strider]

Everything else is rather dichotomous between clinical and basic science activities. It seems to me that with clinical genetics, you’d likely spend a relatively large proportion of your time sequencing and discovering significant new mutations, most of them publishable. Nearly all of the sequencing, gene amplification, and analysis are computerised.

Tomorrow, then...

[Katherine]

Okay! This new main script replaces the old one, and I'm trying to model chemotherapy with multiple rounds of treatment. The picture is supposed to resemble the first colored picture in the paper, but it doesn't. I'm using a series of nested for loops, varying the number of doses and the time between doses. Do you think this is a problem with my keeping track of time or number of cells?

Thank you :)

[Star Strider]

The new main script, CancerModelFiniteRest, replaces CancerModelFinite, correct? Thanks for the paper — I’m growing my library.

By ‘the first colored picture in the paper’, I assume you’re referring to Figure 8A.

This is going to take me a while to get through. It looks like the ODE system cancerEqs are equations (1) and (2) in the paper. I can’t find the cancerEqsVaryBurden ODE system in this paper though. Do they belong in your code? (I may be missing something.) Are they in another paper referenced here but not stated? (I haven’t seen it referred to.) I understand that this is research, and if you’re introducing new equations to test the hypotheses of this paper and are therefore not reproducing the results in this paper, that could certainly be a reason.

We could still be dealing with differences between solvers, since again they don’t mention the necessary characteristics of the ones they are using. I can’t estimate what may be due to solver differences, but I trust the MATLAB solvers.

In order to figure out what the problem may be, I’m going to have to read the paper in some detail and compare it to your code. That could take a while.

Running your code now to see how it goes.

[Katherine]

Yes- the new script replaces CancerModelFinite.

Oops- I meant 7 (although 7-9 as a whole represent one type of therapy).

cancerEqsVaryBurden is what the equations for cell grown would be in the absence of therapy. I simply deleted the therapy term from the given equations, and used the result for the growth period before therapy is begun and between doses. I've used it successfully for earlier simulations, and it hasn't given me issues so far.

Great! I'm really glad you're helping, since I don't think I would have gotten through this without you.

[Star Strider]

Apologise for the delay. I was running your code (while watching Mr Selfridge). It took an hour, so I saved what I hope are the relevant variables and will plot them tomorrow. On the chance that I saved the wrong ones, which ones should I save and plot? I’m familiar with pcolor and the other surface plots. Consider also contourf. (I haven’t read the paper in full yet, so I haven’t figured this out.)

Your approach is intriguing! You’re running one simulation as untreated (control) and the other as treated (intervention). This is actually fun, and I’m learning about an aspect of biomedical simulation I’d definitely not have otherwise encountered.

I’ll finish reading the paper and do the plots tomorrow.

My pleasure!

[Katherine]

It's completely fine!

The relevant variables are endT, and, for the axes, rest (in days), and Ncyc. The latter two are the axes, while endT is what is being plotted. Fig. 8 apparently has three independent variables, and I'm not quite sure how they're plotted on one graph, but max survival time (for various C0 (lambdaCt in my code), Ncyc, and trest (rest in my code)) becomes relevant there.

I'm glad I could show this to you, then :)

Thanks!

[Star Strider]

Thanks for the clarification on the variables. Figures 7 and 8 actually look to me to be filled contour plots with two independent variables and the dependent variable presented as grayscale or spectral variation. The dependent variable in that situation has to be a matrix defined at every combination of the x and y vectors. (No worries! MATLAB makes this easy.)

When you designed your simulation experiment, did you do so with a statistical design in mind, so that you could demonstrate that the difference between the treated and treated cancers was due to your intervention and not due to a random process? Also, look through PubMed to see if anyone else has done the sort of paired simulation study you’re doing. Seems to me that not even the authors of your paper did, though I’m not familiar with that literature.

[Katherine]

I have a small bit of code that I think would produce the type of plot in Figure 7:

hourRest = restArr * 24;

pcolor(hourRest, NcycArr, endT);
colormap('gray');
colorbar;
shading interp;
xlabel('trest (hours)');
ylabel('Ncyc');
title('survival time with various rest periods and numbers of cycles');

It goes after the largest for loop (about line 129), and it should produce a smooth gray plot. (The attached picture a couple comments ago was a basic pcolor of endT, where red was higher survival time. My issue was that it didn't match with what I had expected, based on Figure 7.)

Not particularly- I was simply curious as to the effects of various types of treatment, mostly changing various parameters, and I'm following the paper's lead with that, since I'm by no means an expert in the field.

[Star Strider]

I discovered one significant fact thus far. Quoting:

• Given the model consists of a small number of ordinarydifferential equations, our use of an Euler method for numericalsimulations is amply sufficient; there are no stability difficultiesnor inconvenient computational overheads. (Page. 295)

The Wikipedia article on the Euler method is instructive. The MATLAB solvers are significantly more sophisticated, so those differences could explain some of the inconsistencies between your results and theirs.

Continuing...

[Star Strider]

I may not have the latest version of your code. I can’t find your hourRest variable in the CancerModelFiniteRest code version I have. (I can find NcycArr and endT.) Considering CancerModelFiniteRest took an hour to run the first time, I would prefer to let it do so at night, generate a save file, load the save file in the morning, and experiment with it then. Currently I only plan to save those variables (and some information variables I added), so if there are any other variables you want me to save, plot, and experiment with, let me know.

[Katherine]

Thank you :D

It's rest, I think, and it appears at line 28 of CancerModelFiniteRest. I think lambdaCt would be good as well.

Thanks again!

[Star Strider]

Thanks! I made the changes in my save statement. I’ll run CancerModelFiniteRest in a bit and look for the results in the morning.

The paper was interesting, but difficult to follow, especially the contour plots. It also employed some simplifying assumptions that while necessary, makes it impossible to apply directly. It would be useful in designing an experimental protocol to test the hypotheses it raises (it might be interesting to check to see what other papers have referenced this paper), but is not itself directly clinically applicable. The only significant difference I could find between your code and my interpretation of the paper were the solvers used. Nevertheless, your results and those of the paper should be qualitatively similar but by no means identical.

[Katherine]

I have quite a few plots that match their results (qualitatively, at least!), so I should be good with showing results. Maybe next time I'll pick a better paper :P

Thanks for all your help!

[Star Strider]

My pleasure!

I could only get one plot, that of endT, that I include here as a surface plot. I plotted it as a contourf plot first, but thought the surfc more interesting. Your data are more accurate than those in the paper, and likely could lead to different conclusions.

The endT array were the only data I could plot. Only NcycArr was a vector (all the others were scalars), and I had no idea what to do with any of them. I plotted this figure with the simple command:

surfc(endT)

It would be interesting to look at this and other models, compare them with actual patient survival data, estimate parameters, and see how well the various models explain the data.

This was an interesting project, but a difficult paper to analyse. I congratulate you on being able to decipher it to the extent of integrating the DEs and reproducing and comparing the results. I encourage you to continue your research.

I’ll help if I can!

[Katherine]

Hm. In that case, I'm inclined to believe what my simulation is telling me, since my plot looked essentially the same, although I think checking actual data would be a good way of verifying. However, I think this entire model is developed using flawed assumptions, so I might actually tweak it a bit. That could be an interesting next project.

Thank you! I didn't know much about this field until now, but I think it's actually really interesting; I'll definitely continue to explore this and other papers (and maybe try some original modeling!).

If I do need help I'll definitely email you!

Thanks :)

[Star Strider]

I definitely agree. They actually made a number of simplifying assumptions that while necessary for a model of this sort, can lead to erroneous conclusions, that I believe they overstated. (It would be interesting to see if any clinical papers cited this one, or any similar papers by the same authors.) My modeling experiences are in physiological control, not oncology, but I found your project interesting.

As for doing your own simulations, if your university has a medical school, I suggest you talk with some clinical oncologists and oncology researchers to see what they suggest. They might have some ideas that you could simulate in your models, as well as data you could fit to your simulations to identify parameters. (Identifying models characterized by differential equations are relatively easy to do in MATLAB. I’ve done a fair amount of it.)

I get the impression that you did all this on your own volition and on your own time, as well as keeping up with your university studies. Your work on this project is graduate-level, and is truly impressive.

My e-mail is on my File Exchange profile.

My pleasure!