How to set up two ODE functions in matlab

Hi,I am trying to set up two ode functions fot matlab to solve. The ODEs does not have any analytical solutions, but I want to get an numerical solution and plot in matlab. The equations are describing tumor growth, and are called the Gyllenberg-Webb model, being as follows:

dP/dt=(b-r_o(N))*P+r_i(N)*Q

dQ/dt=r_o(N)*P-(r_i(N)+my_q)*Q

with initial conditions P(0)=P_0>0 and Q(0)>or equal to 0.

So I am wondering if this ODE system is possible to set up, because of the three variables; P,Q and r_o and r_i being functions of N aswell. Also P, Q and N are functions of time, and N=P+Q (read: total cell population=cell pop.1 + cell pop.2 in the tumor).

Or if it would be better to use this system of only P and N variables:

dP/dt=(b-r_i(N)-r_o(N))*P+r_i(N)*N

dN/dt=(b+my_q)*P-my_q*N

I was thinking of setting it up with the Eulers method, but dont know exactly how, and by using the ode45 function in matlab.

Any help would be great appreciated.

Answer 1

[Are Mjaavatten]

You should collect your unknowns P and N in a vector x, and define a function for the right-hand side of your ODE system, something like this:

function dxdt = rhs(t,x,b,r_i,r_o,my_q)
    P = x(1);
    N = x(2);
    dPdt = (b-r_i(N)-r_o(N))*P+r_i(N)*N;
    dNdt = (b+my_q)*P-my_q*N;
    dxdt = [dPdt;dNdt];
end

Save this as rhs.m.

Then create a script where you define all parameters and functions, as well as initial values. The ODE solvers do not take extra function parameters explicitly, so we define an inline function f of only t and x.

% Define parameters and functions going into your ODE system
b = 1;
r_i = @(N) sqrt(N) ;
r_o = @(N) exp(-N);
my_q = -4;

% Define an inline function f(t,x), where the extra parameters to rhs are 
% "baked in":
f = @(t,x) rhs(t,x,b,r_i,r_o,my_q);

% Initialize:
P0 = 1;
N0 = 2;
x0 = [P0;N0];

% Solve:
[t,x] = ode45(f,x0,[0,100]);
P = x(:,1);
N = x(:,2);
plot(t,P);

If your system diverges very quickly, it may be useful to start with a very short time interval, to see what happens. For some parameter sets you system may be stiff, and ode45 may be very slow. If so, try e.g. ode15s.

[Kristoffer Ree]

I am assuming to first save rhs.m as one file, then to use the ode45 to solve it in a new script? Also when I run the first rhs.m, I get error 'not enough input arguments: error in rhs (line 2)'.

And I dont understand how I should define a short time interval in the above suggested script with ode45? Also shouldn't it be plot(t,x) in the last line?

[Are Mjaavatten]

Sorry, I was a bit too quick. I have modified my answer. This should run ok. Modify to use your own parameters and functions. You may of course plot both P and N in the same plot, by

plot(t,x)

[Kristoffer Ree]

I got it to run now, so I assume it works fine. But I wonder what @(N) in front actually means in:

r_i = @(N) sqrt(N) ;

r_o = @(N) exp(-N);

is it just for good structure or does matlab interpret this in any way. And the actual functions of N, I guess you just used some easy examples here with the sqrt and exp?

Because some suggested functions for r_i and r_o are these:

r_i=c/(N+m) and r_o=k*N/(a*N+1)

where c,m and a all are constants, how could I change to these suggestions for r_i and r_o?

[Kristoffer Ree]

I tried this for dPdt, dQdt and dNdt, and I am expecting to get three curves, and that N always should equal P+Q, but this does not quite occur, and it takes long time to run, so I only ran it from t=0->2.

 function dxdt = pqn(t,x,b,r_i,r_o,my_q)
    P = x(1);
    N = x(2);
    Q = x(3);
    dPdt = (b-r_o(N))*P+r_i(N)*Q;
    dQdt = r_o(N)*P-(r_i(N)+my_q)*Q;
    dNdt = (b+my_q)*P-my_q*N;
    dxdt = [dPdt;dQdt;dNdt];
 end

 % Define parameters and functions going into your ODE system
b = 0.8;
r_i = @(N) 1/(N+2);
r_o = @(N) 2*N/(1*N+1);
my_q = 0.4;
% Define an inline function f(t,x), where the extra parameters to rhs are 
% "baked in":
f = @(t,x) pqn(t,x,b,r_i,r_o,my_q);
% Initialize:
P0 = 0.95;
Q0 = 0.05;
N0 = 1;
x0 = [P0;Q0;N0];
% Solve:
[t,x] = ode45(f,[0 2],x0);
P = x(:,1);
Q = x(:,2);
N = x(:,3);
plot(t,x)
title('\bf Gyllenberg-Webb model','FontSize', 20);
xlabel({'\bf Time'},'FontSize', 18);
ylabel({'\bf Cell population'},'FontSize', 18);

Does this seem to be right scripted despite not quite showing what I want(I think this might be due to the functions r_i and r_o or the assumed parameter values)?

[Are Mjaavatten]

The expression r_i = @(N) sqrt(N) defines an anonymous function and creates a function handle to this function. Read more in the links.

Your expression for dNdt does not equal dQdT + dPdt. If I correct this expression, cell population Q goes negative at t = 0.06 and you hit a singularity at around t = 0.375. You see this better if you use ode15s.

This is probably not what you want. You should re-examine your parameters and possibly your equations.

[Kristoffer Ree]

I don1t see whats wrong with dNdt:

dPdt+dQdt

((b-r_o)*P+r_i*Q)+(r_o*P-(r_i+my_q)*Q)

(b*P-r_o*P+r_i*Q)+(r_o*P-r_i*Q-my_q*Q)

b*P-my_q*Q, And N=P+Q, thus, Q=N-P

b*P-my_q*(N-P)

(b+my_q)*P-my_q*N

?

I tried the same script with ode15s, and got at least smooth curves, but I am not actually expecting negative values for any of the populations.

This is more like what I am expecting the curves to show (without the D curve, and f(r) and g(r) are r_o(N) and r_i(N).)

[Torsten]

N = P+Q follows if you set dN=b*P-my_q*Q.

You cannot already assume it here by inserting Q = N-P in the expression dN = b*P-my_q*Q.

Best wishes

Torsten.

[Are Mjaavatten]

The error comes from your assignment of N and Q in pqn. It should be:

    P = x(1);
    Q = x(2);
    N = x(3);

Now Q stays positive and the singularity disappears.

[Kristoffer Ree]

Yes, well spotted, and now if i choose, dNdt = b*P-my_q*Q or if I choose dNdt = (b+my_q)*P-my_q*N or simply dNdt = dPdt+dQdt I get the same plots for all.

Now all the plots starts to look more like they should.