SEIR model for ebola

%% 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);

Single objective optimization:
3 Variables

Options:
CreationFcn:       @gacreationuniform
CrossoverFcn:      @crossoverscattered
SelectionFcn:      @selectionstochunif
MutationFcn:       @mutationadaptfeasible

                                  Best           Mean      Stall
Generation      Func-count        f(x)           f(x)    Generations
    1              100       1.628e+15       1.628e+15        0
    2              147       1.628e+15       1.628e+15        0
    3              194       1.628e+15       1.628e+15        0
    4              241       1.628e+15       1.628e+15        1
    5              288       1.628e+15       1.628e+15        0
    6              335       1.628e+15       1.628e+15        1
    7              382       1.628e+15       1.628e+15        2
    8              429       1.628e+15       1.628e+15        0
    9              476       1.628e+15       1.628e+15        1
   10              523       1.628e+15       1.628e+15        2
   11              570       1.628e+15       1.628e+15        0
   12              617       1.628e+15       1.628e+15        0
   13              664       1.628e+15       1.628e+15        1
   14              711       1.628e+15       1.628e+15        0
   15              758       1.628e+15       1.628e+15        1
   16              805       1.628e+15       1.628e+15        2
   17              852       1.628e+15       1.628e+15        3
   18              899       1.628e+15       1.628e+15        4
   19              946       1.628e+15       1.628e+15        5
   20              993       1.628e+15       1.628e+15        0
   21             1040       1.628e+15       1.628e+15        1
   22             1087       1.628e+15       1.628e+15        0
   23             1134       1.628e+15       1.628e+15        0
   24             1181       1.628e+15       1.628e+15        0
   25             1228       1.628e+15       1.628e+15        0
   26             1275       1.628e+15       1.628e+15        1
   27             1322       1.628e+15       1.628e+15        2
   28             1369       1.628e+15       1.628e+15        3
   29             1416       1.628e+15       1.628e+15        4
   30             1463       1.628e+15       1.628e+15        5

                                  Best           Mean      Stall
Generation      Func-count        f(x)           f(x)    Generations
   31             1510       1.628e+15       1.628e+15        0
   32             1557       1.628e+15       1.628e+15        0
   33             1604       1.628e+15       1.628e+15        1
   34             1651       1.628e+15       1.628e+15        0
   35             1698       1.628e+15       1.628e+15        1
   36             1745       1.628e+15       1.628e+15        0
   37             1792       1.628e+15       1.628e+15        1
   38             1839       1.628e+15       1.628e+15        2
   39             1886       1.628e+15       1.628e+15        0
   40             1933       1.628e+15       1.628e+15        1
   41             1980       1.628e+15       1.628e+15        0
   42             2027       1.628e+15       1.628e+15        0
   43             2074       1.628e+15       1.628e+15        1
   44             2121       1.628e+15       1.628e+15        0
   45             2168       1.628e+15       1.628e+15        0
   46             2215       1.628e+15       1.628e+15        1
   47             2262       1.628e+15       1.628e+15        0
   48             2309       1.628e+15       1.628e+15        0
   49             2356       1.628e+15       1.628e+15        1
   50             2403       1.628e+15       1.628e+15        2
   51             2450       1.628e+15       1.628e+15        3
ga stopped because the average change in the fitness value is less than options.FunctionTolerance.

beta = estimated_params(1);
q = estimated_params(2);
alpha = estimated_params(3);

fprintf('Estimated Parameters:\n');

Estimated Parameters:

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

Beta: 1.00

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

Alpha: 1.00

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

Q: 1.00

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;

this is the code which i used for plotting. but while estimating the parameters, each time I run the code i get different parameter values also the parameter values estimated are wrong. what would possibly be the problem?

Answer 1

[Kautuk Raj]

It looks like you are encountering issues with parameter estimation in your SEIR model using a genetic algorithm (GA). Here are a few suggestions which you can try:

1. Initial Guess and Bounds:

Expand the bounds to realistic ranges. For example:

lower_bounds = [0, 0, 0]; 
upper_bounds = [5, 5, 2]; 

2. Objective Function Scaling:

Normalize the error in your objective_function to improve optimization:

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

3. Optimization Options:

Adjust GA options to improve convergence:

options = optimoptions('ga', 'PopulationSize', 100, 'MaxGenerations', 200, 'FunctionTolerance', 1e-6); 

4. Stochastic Nature of GA:

Set a random seed for reproducibility:

rng(1); 

I hope these adjustments help improve the parameter estimation process in your SEIR model.