Global stability of an epidemiology

I want to get the graph of the global stability and the R_0 of the model. Below is the code and the error message please I need help on how to fix the error and obtain a desired result clear; clc; % Define the parameters of the epidemic model alpha1 = 0.0444; alpha3 = 0.2; alpha0 =0.06685; alpha2 = 0.0000548; alpha4 =0.001096; kappa = 10^6; mu = 2.249*10^-5; lambda = 0.0546; omega =0; %(0, 0.9) phi =0.012 ; rho = 0.0031; tua =0.85 ; eta =0.33; a=0.028;

% Define the initial conditions of the epidemic model I0 = 0.01; S0 = 1 ; R0 = 0; E0=0; H0=0.68; % Define the time vector tspan = [0 10]; % Define the differential equations of the epidemic model ode=@(t,y)[alpha0-alpha1*lambda*S-mu*S+alpha2*omega*E+alpha4*R; alpha1*lambda*S -... (alpha2 + mu)*E; alpha2*(1-omega)*E - (alpha3 + phi+ mu)*I; alpha3*I - (alpha4 + mu)*R;... rho*(1-tau )*I+ (a-eta)*H]; % Solve the system of ODEs using the ode45 solver [t,y] = ode45(ode, tspan, [S0, E0, I0, R0, H0]); % Plot the results figure; plot(t, y(:,1), 'r', t, y(:,2), 'b', t, y(:,3), 'g', y(:,4), 'k', y(:,5), 'm'); xlabel('Time'); ylabel('Population'); legend('Susceptible', 'Exposed', 'Infected', 'Recovered', 'Vibro cholera population'); % Calculate the basic reproduction number R0 n = alpha2*alpha0*(1-omega)*kappa*(eta-a)+rho*(1-tau); R0 = kappa*mu*(alpha2+mu)*(alpha3+phi+mu)*(eta-a)/n; % Determine the global stability of the epidemic model if R0 < 1 disp('The epidemic model is globally stable'); else disp('The epidemic model is globally unstable'); end

error Error in ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in glo (line 32) [t,y] = ode45(ode, tspan, [S0, E0, I0, R0, H0]);

[Sam Chak]

Hi @Sunday

From the theory of system stability, you need to find the Lyapunov function to prove if the disease equilibrium is globally asymptotically stable or not. The basic reproduction number R₀ is just an indicator about the spread, but doesn't mathematically guarantee the stability. You can find this article. It is open access.

https://link.springer.com/article/10.1007/s12648-023-02895-6

Answer 1

[Torsten]

clear; 
clc; 
% Define the parameters of the epidemic model 
alpha1 = 0.0444; % transmission rate 
alpha3 = 0.2; % fixed recovery rate 
alpha0 =0.06685; % fixed birth rate 
alpha2 = 0.0000548; % fixed exposed rate 
alpha4 =0.001096; % fixed rate 
kappa = 10^6; % fixed concentration rate of vibrio cholera in water bodies 
mu = 2.249*10^-5; %fixed Natural death rate 
lambda = 0.0546; 
omega =0; %(0, 0.9) fixed rate of strong immune system 
phi =0.012 ; % %fixed rate of death induced by cholera 
rho = 0.0031; %contribution of each infected persons to the number of vibrio cholera 
tau =0.85 ; %infectious class compliance to hygienic principles 
eta =0.33; %death rate induced by of vibrio cholera 
a=0.028; % birth rate of Vibrio cholera
% Define the initial conditions of the epidemic model 
I0 = 0.01; % initial number of infected individuals 
S0 = 1 ; % initial number of susceptible individuals 
R0 = 0; % initial number of recovered individuals 
E0=0; % initial number of exposed individuals 
H0=0.68; % initial number of Vibrao cholera bacteria 
% Define the time vector 
tspan = [0 10]; 
% Define the differential equations of the epidemic model 
Ode = @(t,S,E,I,R,H)[alpha0-alpha1*lambda*S-mu*S+alpha2*omega*E+alpha4*R; alpha1*lambda*S-... 
                    (alpha2 + mu)*E; alpha2*(1-omega)*E-(alpha3 + phi+ mu)*I; alpha3*I-(alpha4 + mu)*R;...
                    rho*(1-tau)*I+(a-eta)*H]; 
ode = @(t,y)Ode(t,y(1),y(2),y(3),y(4),y(5));
% Solve the system of ODEs using the ode45 solver 
[t,y] = ode45(ode, tspan, [S0, E0, I0, R0, H0]); 
% Plot the results figure; 
plot(t, y(:,1), 'r', t, y(:,2), 'b', t, y(:,3), 'g', t,y(:,4), 'k', t,y(:,5), 'm'); 
xlabel('Time'); 
ylabel('Population'); 
legend('Susceptible', 'Exposed', 'Infected', 'Recovered', 'Vibro cholera population'); 

% Calculate the basic reproduction number R0 
n = alpha2*alpha0*(1-omega)*kappa*(eta-a)+rho*(1-tau); 
R0 = kappa*mu*(alpha2+mu)*(alpha3+phi+mu)*(eta-a)/n; 
% Determine the global stability of the epidemic model 
if R0 < 1 
  disp('The epidemic model is globally stable'); 
else 
  disp('The epidemic model is globally unstable'); 
end

The epidemic model is globally stable