Update ODE variables in Matlab

Question

I am trying to replicate the results of a paper (Immune Response). In short, it looks at the concentration of pathogens, antibodies and plasma cells, and the effect the immune response has on an organ. Depending on the state of the organ, a variable which affects the concentration of plasma cells will change.

I have modeled this in Matlab. My issue is that the conditional statements for the plasma-concentration-affecting variable is being ignored. By commenting out portions of the code and looking at the graphs, regardless of what the second condition is, the output is only responding to the initial condition.

How can I update values so that they are interpreted by ode45 (or ode23) and is there a better approach to approach this problem?

My code is below.

function dx = iir_2(t, x)
dx = [0; 0; 0; 0];
a11 = 1;
a31 = 1;
a41 = 1;
a12 = 1;
a22 = 3;
a32 = 1.5;
a42 = 1;
a23 = 1;
a33 = 0.5;

b1 = -1;
b2 = 1;
b3 = 1;
b4 = -1;

u1 = 0;
u2 = 0;
u3 = 0;
u4 = 0;

tau = 0.1;

if x(4) >= 0.5
    a21x4 = 0;
else
    a21x4 = cos(pi * x(4));
end

if x(4) > 1
    x(4) = 1;
end

dx(1) = (a11 - a12 * x(3)) * x(1) + b1 * u1;
dx(2) = a21x4 * a22 * x(1) * (t - tau) * x(3) * (t - tau) - a23 * (x(2) - 2) + b2 * u2;
dx(3) = a31 * x(2) - (a32 + a33 * x(1)) * x(3) + b3 * u3;
dx(4) = a41 * x(1) - a42 * x(4) + b4 * u4;

This second script that calls the above function is..

    Case = ['Subclinical' 'Clinical' 'Chronic' 'Lethal'];
x1_0 = [1.5 2 2.57 3];
x2_0 = 2;
x3_0 = (1*x2_0)/1.5;
x4_0 = 0;

for i = 1:4
    if i == 1
        state = 'Subclinical';
    elseif i == 2
        state = 'Clinical';
    elseif i == 3
        state = 'Chronic';
    else
        state = 'Lethal';
    end

    options = odeset('RelTol', 1e-3, 'NonNegative', [1 2 3 4]);
    [t,x] = ode45(@iir_2, [0 10], [x1_0(i) x2_0 x3_0 x4_0], options); % use ode23???
    figure 
    plot(t,x);
    str = sprintf('Case Number = %d\nx(0) = %d\n%s', i, x1_0(i), state);
    title(str);
    axis([0,10,0,10])
    legend('Pathogen', 'Plasma Cell', 'Antibody', 'Organ');
end

The conditional statements are in the function.

if x(4) >= 0.5
        a21x4 = 0;
    else
        a21x4 = cos(pi * x(4));
    end

    if x(4) > 1
        x(4) = 1;
    end

Answer 1

There are two ways of doing this.

Either way you need to understand that you cannot just set x(4) to 1 and hope everything works out for the best. Matlab won't care about the preceding values of x(4) as they are all stored in memory. Furthermore, the upcoming value of x(4) is determined by dx(4) and the previously retained values of x(4) (which you cannot set).

You have two possible solutions to your problem:

A) set dx(4) = 0 when your condition is met

if x(4) > 1
    dx(4) = 0;
else
    dx(4) = a41 * x(1) - a42 * x(4) + b4 * u4;
end

This however will not result in a perfect x(4) = 1, rather there will be a small error.

B) Calling an event function, which can be learned about here and here, in which case you can interrupt for a new function call.