How to estimate parameters of fractional SIR epidemic model?

I want to know the code for finding parameter estimations of a fractional SIR epidemic model for a vector born disease in Matlab. I have some real data (no of infected individuals) on a specific vector born disease for eg dengue. How to simulate this model and estimate the parameters?

Answer 1

[Manikanta Aditya]

Hi @JAYARAM PRAKASH

To estimate parameters of a fractional SIR epidemic model for a vector-borne disease like dengue in MATLAB, you can follow these steps:

1. Define the Fractional SIR Model: Write the differential equations for the fractional SIR model. You can use the Caputo derivative for the fractional part.
2. Load Your Data: Import your real data of infected individuals into MATLAB.
3. Set Up the Optimization Problem: Use MATLAB's optimization functions like fminsearch or lsqcurvefit to estimate the parameters by fitting the model to your data.
4. Simulate the Model: Use the estimated parameters to simulate the model and compare it with your data.

Please find some resources, I found on modeling epidemics:

• GitHub - saucermanlab/modeling-epidemics: Computational models to simulate epidemics
• param-estimation-SIR/Matlab/SIR_Example_Main.m at master · epimath/param-estimation-SIR · GitHub

I hope this helps.

[Manikanta Aditya]

Basic outline of code, how it might look:

% Define the fractional SIR model using Caputo derivative
function dydt = fractionalSIR(t, y, beta, gamma, alpha)
    D = @(f, a) diff(f, a); % Caputo derivative
    S = y(1);
    I = y(2);
    R = y(3);
    dydt = [-beta * S * I;
            beta * S * I - gamma * I;
            gamma * I];
end

% Load your data
data = load('your_data.mat'); % Replace with your data file

% Initial guess for parameters [beta, gamma, alpha]
params0 = [0.5, 0.1, 0.9];

% Define the objective function for optimization
objective = @(params) sum((data - simulateModel(params)).^2);

% Estimate parameters using fminsearch
params_est = fminsearch(objective, params0);

% Simulate the model with estimated parameters
[t, y] = ode45(@(t, y) fractionalSIR(t, y, params_est(1), params_est(2), params_est(3)), [0, max(data.time)], [S0, I0, R0]);

% Plot the results
plot(t, y(:, 2), 'r', 'LineWidth', 2); % Infected individuals
hold on;
plot(data.time, data.infected, 'bo'); % Real data
legend('Model', 'Data');
xlabel('Time');
ylabel('Infected Individuals');
title('Fractional SIR Model Fit');

Answer 2

[Star Strider]

See for example Parameter Estimation for a System of Differential Equations

You have to define your differential equations and use them as the ‘DifEq’ function. The data you fit must match the dimensions of the output, so use the ‘C’ result for that, and select the columns of ‘Cv’ to match the columns in your data, if necessary.

Example —

function C=kinetics(theta,t)
    c0=theta(end-3:end);
    [T,Cv]=ode45(@DifEq,t,c0);
    
    function dC=DifEq(t,c)
        dcdt=zeros(4,1);
        dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
        dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
        dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
        dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
        dC=dcdt;
    end
    C=Cv;
end

t=[0.1
    0.2
    0.4
    0.6
    0.8
    1
    1.5
    2
    3
    4
    5
    6];
c=[0.902	    0.06997	0.02463	0.00218
    0.8072	    0.1353	0.0482	0.008192
    0.6757	    0.2123	0.0864	0.0289
    0.5569	    0.2789	0.1063	0.06233
    0.4297	    0.3292	0.1476	0.09756
    0.3774	    0.3457	0.1485	0.1255
    0.2149	    0.3486	0.1821	0.2526
    0.141	    0.3254	0.194	0.3401
    0.04921	    0.2445	0.1742	0.5277
    0.0178	    0.1728	0.1732	0.6323
    0.006431	0.1091	0.1137	0.7702
    0.002595	0.08301	0.08224	0.835]; 

theta0 = rand(10,1);

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

for k1 = 1:numel(theta)
    fprintf(1, '\t\tTheta(%2d) = %23.15E\n', k1, theta(k1))
end

		Theta( 1) =   7.648784502758464E-01
		Theta( 2) =   2.348642320704391E-01
		Theta( 3) =   2.089879571711175E-01
		Theta( 4) =   4.920779719013652E-01
		Theta( 5) =   6.221107591627785E-01
		Theta( 6) =   1.175188421416447E-13
		Theta( 7) =   9.028750693139709E-01
		Theta( 8) =   7.145169346891880E-02
		Theta( 9) =   2.840857166532932E-02
		Theta(10) =   4.784361317284549E-12

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

.