SEIR model for ebola (2024)

%% SEIR model function

function dydt = seir_model(t, y, params, N, fixed_gamma, fixed_delta)

beta = params(1);

q = params(2);

alpha = params(3); % Now alpha is part of the parameters to be estimated

S = y(1);

E = y(2);

I = y(3);

R = y(4);

dS = -beta * S * (q * E + I) / N;

dE = beta * S * (q * E + I) / N - fixed_delta * E;

dI = fixed_delta * E - fixed_gamma * I;

dR = fixed_gamma * I;

dydt = [dS; dE; dI; dR];

end

%% GMMP scheme to solve FDEs

function y = gmmp_scheme(f, y0, t, params, N, fixed_gamma, fixed_delta)

n = length(t);

m = length(y0);

y = zeros(m, n);

y(:, 1) = y0;

% Initial values

for i = 2:n

ti = t(i);

yi = y(:, i-1);

% Fractional derivative approximation using GMMP

yi_next = yi + (ti^(params(3)-1) - (i-1)^(params(3)-1)) * f(ti, yi, params, N, fixed_gamma, fixed_delta);

y(:, i) = yi_next;

end

end

%% Objective function for MGAM

function error = objective_function(params, f, y0, t, data, N, fixed_gamma, fixed_delta)

beta = params(1);

q = params(2);

alpha = params(3); % Now alpha is part of the parameters to be estimated

y = gmmp_scheme(f, y0, t, params, N, fixed_gamma, fixed_delta);

infected_model = y(3, :);

error = sum((infected_model - data).^2);

end

%% MGAM function for parameter estimation

function estimated_params = mgam(f, y0, t, data, N, fixed_gamma, fixed_delta)

param_guess = [0.5, 0.5, 0.5]; % Guess for beta, q, and alpha

obj_func = @(params) objective_function(params, f, y0, t, data, N, fixed_gamma, fixed_delta); % Pass data here

options = optimoptions('ga', 'Display', 'iter', 'PopulationSize', 50);

[estimated_params, ~] = ga(obj_func, 3, [], [], [], [], [0, 0, 0], [1, 1, 1], [], options);

end

%% Main script

N = 18805278;

m = 1; % <-- User input value

S0 = N * m/100;

E0 = 0;

I0 = 15;

R0 = 0;

y0 = [S0; E0; I0; R0];

t0 = 0;

tf = 250;

h = 0.1;

t = t0:h:tf;

data = y0(3) * exp(0.05 * (t - t0));

%% Define fixed gamma and delta here (replace with your desired values)

fixed_gamma = 1/7;

fixed_delta = 1/12;

estimated_params = mgam(@seir_model, y0, t, data, N, fixed_gamma, fixed_delta);

beta = estimated_params(1);

q = estimated_params(2);

alpha = estimated_params(3);

fprintf('Estimated Parameters:\n');

fprintf('Beta: %.2f\n', beta);

fprintf('Alpha: %.2f\n', alpha);

fprintf('Q: %.2f\n', q);

fprintf('\n');

y = gmmp_scheme(@seir_model, y0, t, [beta, q, alpha], N, fixed_gamma, fixed_delta);

figure;

plot(t, y(3, :), 'r-', 'LineWidth', 2);

hold on;

% Set the plot limits

xlim([0 250]); % X-axis limit from 0 to 250

ylim([0 12000]); % Y-axis limit from 0 to 12000

% Optionally, add labels and title for better visualization

xlabel('Time');

ylabel('Infections');

title('SEIR Model Simulation');

hold off;