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

4 views (last 30 days)
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
  2 Comments
Bjorn Gustavsson
Bjorn Gustavsson on 28 Sep 2022
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.

Sign in to comment.

Answers (3)

Davide Masiello
Davide Masiello on 20 Sep 2022
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
  2 Comments
Bjorn Gustavsson
Bjorn Gustavsson on 20 Sep 2022
@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...

Sign in to comment.


David Goodmanson
David Goodmanson on 20 Sep 2022
Edited: David Goodmanson on 20 Sep 2022
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.

Bjorn Gustavsson
Bjorn Gustavsson on 20 Sep 2022
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

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!