SIR model with recovered individuals may lose their immunity and become reinfected with the disease. But came with a failure about integration tolerances

Hi everyone,

I am facing a very intrigating problem atm. I have this SIR model with individuals may lose their immunity and become reinfected with the disease. I tried to use a SIR simulater from matlab community but it seems it didn't work too. So I am trying to do myself. After I run the code, come with an warming:

Warning: Failure at t=4.135655e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.110223e-16) at time t.

In ode45 (line 352)

In Q10_SIR_model (line 15)

See bellow my code:

%% Q10_SIR_model
clc; close all ; clear;
% time
tspan = 0:0.01:100;
% Initial Conditions (Suceptible, Infected, Recovered)
q0 = [99; 1; 0]; 
% Parameters
beta = 0.5; gamma = 20; delta = 0.4; 
%% Model SIR
% Initial conditions as vector
p = [beta; gamma; delta];
[t, y] = ode45(@(t,q)fc_ode_SIR(t,q,p),tspan,q0);
%% Plot
figure;
subplot(2,1,1)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]'); 
legend('Suceptible','Infected', 'Recovered', 'Location','Best');
y = y(:,3);%*N;
% figure(2);
subplot(2,1,2)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]'); 
legend('Infected', 'Location','Best');

And my script:

%% fc_ode_SIR
function dq = fc_ode_SIR(~,q,p) % q = initial conditions (Suceptible, Infected, Recovered)
% p = parameter (beta, gamma, delta, initial conditions) 
beta = p(1); gamma = p(2); delta = p(3); 
S = q(1); I = q(2); R = q(3);
% SIR Model
S = beta*S*I + delta*R;
I = beta*S*I - gamma*I;
R = gamma*I - delta*I;
dq = [-S;
    I;
    R];
end

I don't know where exactaly is my mistake as I run other SIR models and it was fine. I am quite sure it's something silly but if someone could find that I'd appreciate.

Cheers

[Julian Pereira]

Final version:

%% Q10_SIR_model
clc; close all ; clear;
% time
tspan = 0:0.01:100;
% Initial Conditions (Suceptible, Infected, Recovered)
q0 = [99; 1; 0]; 
% Parameters
beta = 0.5; gamma = 20; delta = 0.4; 
%% Model SIR
% Initial conditions as vector
p = [beta; gamma; delta];
[t, y] = ode45(@(t,q)fc_ode_SIR(t,q,p),tspan,q0);
%% Plot
figure;
subplot(2,1,1)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]'); 
legend('Suceptible','Infected', 'Recovered', 'Location','Best');
y = y(:,3);%*N;
% figure(2);
subplot(2,1,2)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]'); 
legend('Infected', 'Location','Best');

And the function:

%% fc_ode_SIR
function dq = fc_ode_SIR(~,q,p) % q = initial conditions (Suceptible, Infected, Recovered)
% p = parameter (beta, mi, ni, initial conditions) 
beta = p(1); gamma = p(2); delta = p(3); 
S = q(1); I = q(2); R = q(3);
% SIR Model
dS = -beta*S*I + delta*R;
dI = beta*S*I - gamma*I;
dR = gamma*I - delta*R;
% dS = beta*S*I + delta*I;           % Note changes here
% dI = beta*S*I - gamma*I - delta*I; % and here
% dR = gamma*I;
dq = [dS;
    dI;
    dR];
end

And the output:

Thanks all for your help!

[Bjorn Gustavsson]

Doesn't seem unreasonable. Perhaps there's even some signs of what is necessary for recurring influenca outbreaks in the oscillations - however, the reality is obviously far more complicated.

Answer 1

[Davide Masiello]

I guess the problem is in this three lines

S = beta*S*I + delta*R;
I = beta*S*I - gamma*I;
R = gamma*I - delta*I;

where you have redefined the variables, which seem to depend on the variables values themselves.

If I had to guess, those expressions are rather the values of the time derivatives of the variables.

Below, I tried to code it that way and it works, but please check the math carefully.

clear,clc

tspan = [0,100];
q0 = [99; 1; 0]; 
beta = 0.5; gamma = 20; delta = 0.4; 
p = [beta; gamma; delta];
[t, y] = ode45(@(t,q)fc_ode_SIR(t,q,p),tspan,q0);
%% Plot
figure;
subplot(2,1,1)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]'); 
legend('Suceptible','Infected', 'Recovered', 'Location','Best');
y = y(:,3);%*N;
% figure(2);
subplot(2,1,2)
plot(t, y); grid on;
title(sprintf('Model SIR mod: p = [%g %g %g].', p));
xlabel('t [dia]'); 
legend('Infected', 'Location','Best');

function dq = fc_ode_SIR(~,q,p)

beta = p(1); gamma = p(2); delta = p(3); 
S = q(1); I = q(2); R = q(3);

dS = beta*S*I + delta*R;
dI = beta*S*I - gamma*I;
dR = gamma*I - delta*I;

dq = [  -dS;
        dI;
        dR
        ];
end

[Julian Pereira]

Thank you so much! Oh it was a silly mistake indeed! Basically I was using the same variable for two diferente meanings.... Thank you!!!

[Bjorn Gustavsson]

@Julian Pereira, no it is not. The loss of R due to transition back to S couldn't be proportional to the amount of I. This makes it possible, at least in principle that the loss of R (delta*I) becomes larger than R which would lead to a negatve number of R...

Answer 2

[David Goodmanson]

Hi Julian,

I have doubts about the equations, because with dq(1) = -S as you have, the equations are equivalent to

S = -beta*S*I - delta*R;
I = beta*S*I - gamma*I;
R = gamma*I - delta*I;

and it doesn't seem like recovered but reinfectible people should decrease the number of susceptible people. Should the upper right sign be positive? Also, should the lower right term be -delta*R instead of -delta*I?

Let's say whatever changes need to be done, or not, are done. The real problem is that the second line above redefines I, and then the new value, rather than the original value of I, is plugged into the third line. And the same thing is true for S redefined on line 1 and used on line 2, actually.

[Julian Pereira]

Hi David,

thanks for your answer and question! I got the model from this book: https://thalis.math.upatras.gr/~bountis/files/def-eq.pdf , page 235. Tha was the def eq:

Answer 3

[Bjorn Gustavsson]

You have an error in the SIR-equations. They should look something like this:

%% fc_ode_SIR
function dq = fc_ode_SIR(~,q,p) % q = initial conditions (Suceptible, Infected, Recovered)
% p = parameter (beta, gamma, delta, initial conditions) 
beta = p(1); gamma = p(2); delta = p(3); 
S = q(1); I = q(2); R = q(3);
% SIR Model
S = beta*S*I + delta*R;
I = beta*S*I - gamma*I;
R = gamma*I - delta*R; % NOTE THE CHANGE HERE!!!!!!
dq = [-S;
    I;
    R];
end

The source of susceptibles should either be some fraction of the recovered as I've changed it above, that is a loss of the recovered should be proportional to the number of recovered. Or some fraction of the infected should transfer directly to the susceptible:

%% fc_ode_SIR
function dq = fc_ode_SIR(~,q,p) % q = initial conditions (Suceptible, Infected, Recovered)
% p = parameter (beta, gamma, delta, initial conditions) 
beta = p(1); gamma = p(2); delta = p(3); 
S = q(1); I = q(2); R = q(3);
% SIR Model
S = beta*S*I + delta*I;           % Note changes here
I = beta*S*I - gamma*I - delta*I; % and here
R = gamma*I;
dq = [-S;
    I;
    R];
end

That was your goof.

Also worth to not is that these simple SIR-models are decidedly un-biological. The constant transition-rates means that there will be a time-independent recovery-rate - if we at time zero turn off the infection-rate the reduction in infected will be a simple exponential decay - similar to radioactive decay. This is clearly not sensible.

HTH