Parameter estimation of SEIHRDBP model

% I want simulate "parameter estimation of SEIHRDBP Model", where
% S=susceptible, E=exposed, I=infected, H=hospitalised, R=recovered,
% D=death, B=buried, P=pathogen. I tried Matlab codes from the Matlab
% community but couldn't get right.  The parameters here are assumed values

        % Parameter value
        % Delta =  params(1);   mu    =  params(2); delta =  params(3); gamma =  params(4);  gamma1 = params(5); 
        % eta  =   params(6);   xi   =   params(7);   d    =   params(8); b    =   params(9); phip =   params(10);
        % d1   =   params(11);  b1   =   params(12); alpha =  params(13);  sigma =  params(14);
        % eta1  =  params(15); beta  =  params(16); beta1  = params(17);  beta2  = params(18);  beta3 =  params(19); 
        % betap =  params(20);
       % params=  [10, 0.5, 0.5, 0.005, 0.12, 0.12, 0.6, 0.8, 0.5, 0, 0.04, 0.04, 0.03, 0.04, 0.04, 0.006, 0.002, 0.012, 0.0012; 0.01];


        % dxdt =   zeros(8,  1);    % initializing the state variables
        % params = zeros(20, 1);    % initializing the parameters
        % 
        % dxdt(1) = params(1) - (params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) +  params(20)*x(8))*x(1)-params(2)*x(1); 
        % dxdt(2) =(params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) +  params(20)*x(8))*x(1) - (params(2) + params(3))*x(2);
        % dxdt(3) = params(3)*x(2) - (params(2) + params(4) + params(5))*x(1);
        % dxdt(4) = params(5)*x(3) - (params(6) + params(7) + params(2))*x(4);
        % dxdt(5) = params(7)*x(4) - params(2)*x(5);
        % dxdt(6) = params(4)*x(3) + params(6)*x(4) -params(8)*x(6);
        % dxdt(7) = params(8)*x(6) - params(9)*x(7);
        % dxdt(8) = params(10) +  params(14)*x(3) + params(15)*x(4) + params(11)*x(6) + params(12)*x(7) -params(13)*x(8);
    
    % This is code I tried. I basically used this help site: https://uk.mathworks.com/matlabcentral/answers/1710260-need-help-with-ode-parameter-estimation
    % The c1---c8 are the state variables, while the theta1---theta20 are the parameter values.   
    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.902	0.06997	0.02463	0.00218
    0.8072	0.1353	0.0482	0.008192 0.8072	0.1353	0.0482	0.008192
    0.6757	0.2123	0.0864	0.0289   0.6757	0.2123	0.0864	0.0289 
    0.5569	0.2789	0.1063	0.06233  0.5569	0.2789	0.1063	0.06233
    0.4297	0.3292	0.1476	0.09756  0.4297	0.3292	0.1476	0.09756
    0.3774	0.3457	0.1485	0.1255   0.3774	0.3457	0.1485	0.1255
    0.2149	0.3486	0.1821	0.2526   0.2149	0.3486	0.1821	0.2526
    0.141	0.3254	0.194	0.3401   0.141	0.3254	0.194	0.3401
    0.04921	0.2445	0.1742	0.5277   0.04921	0.2445	0.1742	0.5277
    0.0178	0.1728	0.1732	0.6323   0.0178	0.1728	0.1732	0.6323
    0.006431	0.1091	0.1137	0.7702 0.006431	0.1091	0.1137	0.7702
    0.002595	0.08301	0.08224	0.835 0.002595	0.08301	0.08224	0.835]; 

theta0 = rand(20,1);

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

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = SEIHRDBP(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('Population')
% legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')
% 

function C=SEIHRDBP(theta,t)
    c0=theta(end-3:end);
    [T,Cv]=ode45(@DifEq,t,c0);
   %
    function dC=DifEq(t,c)
        dcdt=zeros(8,1);
        dcdt(1) = theta(1) - (theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) +  theta(20)*c(8))*c(1)-theta(2)*c(1); 
        dcdt(2) =(theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) +  theta(20)*c(8))*c(1) - (theta(2) + theta(3))*c(2);
        dcdt(3) = theta(3)*c(2) - (theta(2) + theta(4) + theta(5))*c(1);
        dcdt(4) = theta(5)*c(3) - (theta(6) + theta(7) + theta(2))*c(4);
        dcdt(5) = theta(7)*c(4) - theta(2)*c(5);
        dcdt(6) = theta(4)*c(3) + theta(6)*c(4) -theta(8)*c(6);
        dcdt(7) = theta(8)*c(6) - theta(9)*c(7);
        dcdt(8) = theta(10) +  theta(14)*c(3) + theta(15)*c(4) + theta(11)*c(6) + theta(12)*c(7) -theta(13)*c(8);
        dC=dcdt;
    end
    C=Cv;
end
        
% Please can some help check it for me.    I have attached the observed caese of all the state variables(SEIHRDBP)  

Answer 1

[Star Strider]

Thank you for quoting my code!

The problem was here:

c0=theta(end-3:end);

use this instead:

c0=theta(end-7:end);

Since my code uses the last elements of the paramter vector for the initial conditions of the differential equations (it fits them as parameters, as well as the rate constants), that needs to be accounted for. Here there are 8 differential equations, so there must be 8 initial conditions.

Try this —

% I want simulate "parameter estimation of SEIHRDBP Model", where
% S=susceptible, E=exposed, I=infected, H=hospitalised, R=recovered,
% D=death, B=buried, P=pathogen. I tried Matlab codes from the Matlab
% community but couldn't get right.  The parameters here are assumed values

        % Parameter value
        % Delta =  params(1);   mu    =  params(2); delta =  params(3); gamma =  params(4);  gamma1 = params(5); 
        % eta  =   params(6);   xi   =   params(7);   d    =   params(8); b    =   params(9); phip =   params(10);
        % d1   =   params(11);  b1   =   params(12); alpha =  params(13);  sigma =  params(14);
        % eta1  =  params(15); beta  =  params(16); beta1  = params(17);  beta2  = params(18);  beta3 =  params(19); 
        % betap =  params(20);
       % params=  [10, 0.5, 0.5, 0.005, 0.12, 0.12, 0.6, 0.8, 0.5, 0, 0.04, 0.04, 0.03, 0.04, 0.04, 0.006, 0.002, 0.012, 0.0012; 0.01];


        % dxdt =   zeros(8,  1);    % initializing the state variables
        % params = zeros(20, 1);    % initializing the parameters
        % 
        % dxdt(1) = params(1) - (params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) +  params(20)*x(8))*x(1)-params(2)*x(1); 
        % dxdt(2) =(params(16)*x(3) + params(17)*x(4) + params(18)*x(6)+ params(19)*x(7) +  params(20)*x(8))*x(1) - (params(2) + params(3))*x(2);
        % dxdt(3) = params(3)*x(2) - (params(2) + params(4) + params(5))*x(1);
        % dxdt(4) = params(5)*x(3) - (params(6) + params(7) + params(2))*x(4);
        % dxdt(5) = params(7)*x(4) - params(2)*x(5);
        % dxdt(6) = params(4)*x(3) + params(6)*x(4) -params(8)*x(6);
        % dxdt(7) = params(8)*x(6) - params(9)*x(7);
        % dxdt(8) = params(10) +  params(14)*x(3) + params(15)*x(4) + params(11)*x(6) + params(12)*x(7) -params(13)*x(8);

T1 = readtable('observed.csv')

T1 = 35×8 table
    susceptible    exposed    infected    Hospitalized    Recovered    deaths    Buried    Pathogens
    ___________    _______    ________    ____________    _________    ______    ______    _________

        3444          323        164           164            0           55        55         0    
        3444          426          1             1            0            1         1         0    
        3444          300          5             5            0            5         5         0    
        3444          335         11            11            0            9         9         0    
        3444          386         12            12            0            6         6         0    
        3444          350         23            23            0           12        12         0    
        3444          404        130           130            0           55        55         0    
        3444          416       3470          3470            0         2287      2287         0    
        3444          449         54            54            0           33        33         0    
        3444          446          8             8            0            4         4         0    
        3444          374         66            66            0           49        49         0    
       34440        38646      28646         28646            0        11323     11323         0    
        3444          377          6             6            0            3         3         0    
        3444          491         36            36            0           13        13         0    
        3444          491         11            11            0            4         4         0    
        3444          437          1             1            0            1         1         0    

VN = T1.Properties.VariableNames;
L = size(T1,1);
t = linspace(0, L-1, L).';                              % Time Vector (NECESSARY) Units = ?
        
%     % This is code I tried. I basically used this help site: https://uk.mathworks.com/matlabcentral/answers/1710260-need-help-with-ode-parameter-estimation
%     % The c1---c8 are the state variables, while the theta1---theta20 are the parameter values.   
%     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.902	0.06997	0.02463	0.00218
%     0.8072	0.1353	0.0482	0.008192 0.8072	0.1353	0.0482	0.008192
%     0.6757	0.2123	0.0864	0.0289   0.6757	0.2123	0.0864	0.0289 
%     0.5569	0.2789	0.1063	0.06233  0.5569	0.2789	0.1063	0.06233
%     0.4297	0.3292	0.1476	0.09756  0.4297	0.3292	0.1476	0.09756
%     0.3774	0.3457	0.1485	0.1255   0.3774	0.3457	0.1485	0.1255
%     0.2149	0.3486	0.1821	0.2526   0.2149	0.3486	0.1821	0.2526
%     0.141	0.3254	0.194	0.3401   0.141	0.3254	0.194	0.3401
%     0.04921	0.2445	0.1742	0.5277   0.04921	0.2445	0.1742	0.5277
%     0.0178	0.1728	0.1732	0.6323   0.0178	0.1728	0.1732	0.6323
%     0.006431	0.1091	0.1137	0.7702 0.006431	0.1091	0.1137	0.7702
%     0.002595	0.08301	0.08224	0.835 0.002595	0.08301	0.08224	0.835]; 

c = table2array(T1);

theta0 = rand(28,1);

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

Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.

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

	Rate Constants:

for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

		Theta(1) =  2.40444
		Theta(2) =  0.03125
		Theta(3) =  1.41359
		Theta(4) =  1.11936
		Theta(5) =  1.11908
		Theta(6) =  0.30996
		Theta(7) =  0.06624
		Theta(8) =  0.60836
		Theta(9) =  0.41146
		Theta(10) =  0.29405
		Theta(11) =  0.29017
		Theta(12) =  0.13046
		Theta(13) =  1.72754
		Theta(14) =  0.69723
		Theta(15) =  0.38384
		Theta(16) =  1.29961
		Theta(17) =  0.99986
		Theta(18) =  0.78319
		Theta(19) =  0.54491
		Theta(20) =  0.44378
		Theta(21) =  0.41763
		Theta(22) =  1.15486
		Theta(23) =  1.41118
		Theta(24) =  0.89700
		Theta(25) =  0.14331
		Theta(26) =  0.62479
		Theta(27) =  0.24113
		Theta(28) =  0.66698

tv = linspace(min(t), max(t));
Cfit = SEIHRDBP(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('Population')
% legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')
legend(hlp, VN, 'Location','best')

% 

function C=SEIHRDBP(theta,t)
    c0=theta(end-7:end);
    [T,Cv]=ode45(@DifEq,t,c0);
   %
    function dC=DifEq(t,c)
        dcdt=zeros(8,1);
        dcdt(1) = theta(1) - (theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) +  theta(20)*c(8))*c(1)-theta(2)*c(1); 
        dcdt(2) =(theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) +  theta(20)*c(8))*c(1) - (theta(2) + theta(3))*c(2);
        dcdt(3) = theta(3)*c(2) - (theta(2) + theta(4) + theta(5))*c(1);
        dcdt(4) = theta(5)*c(3) - (theta(6) + theta(7) + theta(2))*c(4);
        dcdt(5) = theta(7)*c(4) - theta(2)*c(5);
        dcdt(6) = theta(4)*c(3) + theta(6)*c(4) -theta(8)*c(6);
        dcdt(7) = theta(8)*c(6) - theta(9)*c(7);
        dcdt(8) = theta(10) +  theta(14)*c(3) + theta(15)*c(4) + theta(11)*c(6) + theta(12)*c(7) -theta(13)*c(8);
        dC=dcdt;
    end
    C=Cv;
end
        
% Please can some help check it for me.    I have attached the observed caese of all the state variables(SEIHRDBP)  

That works, however it is not an escellent fit. I would keep funning it with random initial parameter vectyors until you get a good fit. Using the ga (genetic algorithm) function might produce better results.

.

[University]

Thank you so very much for you help. Please how does the ga (genetic algorithm) works? And why is it better?

I tried to run it my laptop, it was sucessful but subsequently, it returns this error:

Warning: Failure at t=5.304862e+00. Unable to meet integration tolerances without reducing the step size below the

smallest value allowed (1.421085e-14) at time t.

> In ode45 (line 350)

In sample>SEIHRDBP (line 58)

In lsqcurvefit (line 256)

In sample (line 29)

Error using lsqcurvefit

Function value and YDATA sizes are not equal.

Error using lsqcurvefit

Function value and YDATA sizes are not equal.

Error in sample (line 56)

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

Error using helpview

optim not found.

[Star Strider]

Here is the genetic algorithm approach to this problem. It takes longer than the allotted 55 seconds to run, so I cannot run it here. However if you have the Global Optimization Toolbox, it should give a better result than the gradient descent approach. Also, change the solver to ode15s, since the system is likely ‘stiff’.

Also, if you run it on your own computer on in MATLAB Online, use the ‘opts’ options structure instead of ‘optsAns’, and change that in the ga call as well (the last argument).

The error results in the ODE integrator encountering a singularity (±Inf) that then stops it because it can go no further. Start the code again and keep running the code, since it will eventually choose a set of parameters that will allow it to run completely and produce a useful result.

Try this —

T1 = readtable('observed.csv')

VN = T1.Properties.VariableNames;
L = size(T1,1);
t = linspace(0, L-1, L).';                              % Time Vector (NECESSARY) Units = ?

c = table2array(T1);

ftns = @(theta) norm(c-SEIHRDBP(theta,t));
PopSz = 500;
Parms = 28;

optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-2, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10);                                                                  % Options Structure For 'Answers' Problems
% opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
% optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);   % With 'HybridFcn'

t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end

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','N')
legend(hlp, VN, 'Location','best')

function C=SEIHRDBP(theta,t)
    c0=theta(end-7:end);
    [T,Cv]=ode15s(@DifEq,t,c0);
   %
    function dC=DifEq(t,c)
        dcdt=zeros(8,1);
        dcdt(1) = theta(1) - (theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) +  theta(20)*c(8))*c(1)-theta(2)*c(1); 
        dcdt(2) =(theta(16)*c(3) + theta(17)*c(4) + theta(18)*c(6)+ theta(19)*c(7) +  theta(20)*c(8))*c(1) - (theta(2) + theta(3))*c(2);
        dcdt(3) = theta(3)*c(2) - (theta(2) + theta(4) + theta(5))*c(1);
        dcdt(4) = theta(5)*c(3) - (theta(6) + theta(7) + theta(2))*c(4);
        dcdt(5) = theta(7)*c(4) - theta(2)*c(5);
        dcdt(6) = theta(4)*c(3) + theta(6)*c(4) -theta(8)*c(6);
        dcdt(7) = theta(8)*c(6) - theta(9)*c(7);
        dcdt(8) = theta(10) +  theta(14)*c(3) + theta(15)*c(4) + theta(11)*c(6) + theta(12)*c(7) -theta(13)*c(8);
        dC=dcdt;
    end
    C=Cv;
end

.

[University]

Thank you so very much for your help. Why do you choose theta0 = (1, 28)? Is it because we have 8 state variables and 20 parameters. Since there 20 parameters, why don't theta0 = (1, 20)?

[Star Strider]

That is because I have my code set up to estimate the initial conditions of the differential equations as separate parameters, rather than simply setting them as the initial concentrations in the data (although it is perfectly acceptable to use the first data entries as the initial estimates of the initial conditions).

The parameter vector ‘theta’ is set up so that the first 20 elements (in this instance) are the kinetic parameters, and the last 8 (in this instance) are the initial conditions. Since everything is measured and therefore subject to the same sort of measurement uncertainty, it actually makes sense to estimate the initial conditions as parameters.

The genetic algorithm will probably return the best estimates of the parameters, however it may take it a while. In situations such as yours, when you are estimating 28 parameters, I set ga to execute in a loop for several iterations (usually between 10 and 25, although occasionally as many as 100), saving the parameters and final fitness value (and execution times and the number of generations and other statistics that interest me) in each iteration in a cell array or numerical matrix, and then go off and do something else for a few minutes.

[University]

Oh, I see. Thank you

[Star Strider]

As always, my pleasure!

EDIT — (3 Mar 2024 at 17:18)

The data you are currently using are not appropriate for estimating the parameters of your system, because the data does not demonstrate any change. These data are essentially steady-state. You need a different data set.

.

[University]

Hi Strider,

I am back please. I have gotten the data correct data with date. The whole idea is that I want the parameter estimation for parameters for SEI. I.e., run the parameter estimation these classes (SEI). These are only the active human class. The parameters associated with classes are theta16, theta17, theta18, theta19, theta20, theta3, theta7 and theta8.

I tried to run the new code use the "ga" algorithm but I got the error:

Arrays have incompatible sizes for this operation.

Error in sample>@(theta)norm(c-SEIHRDBP(theta,t)) (line 9)

ftns = @(theta) norm(c-SEIHRDBP(theta,t));

Error in createAnonymousFcn>@(x)fcn(x,FcnArgs{:}) (line 11)

fcn_handle = @(x) fcn(x,FcnArgs{:});

In fact, this is my new corrected codes and new data is attached.

SEIHRDBP_MODEL_PARAMETER

Arrays have incompatible sizes for this operation.

Error in solution>SEIHRDBP_MODEL_PARAMETER (line 59)
opts = optimoptions('ga','PopulationSize',PopSz,'InitialPopulationMatrix',randi(1E+4,PopSz,Parameter).*[ones(1,12)*1E-3 ones(1,7)*10],'MaxGenerations',2E3)%,'PlotFcn','gaplotbestf');

function SEIHRDBP_MODEL_PARAMETER

%7 ODEs to describe the model (S,E,I,H,R,D,B, P)
    function M = SEIHRDBP(p,t)
        %Initial conditions with population 
        M0 = p(21:28);
        [~,m] = ode15s(@DifEq,t,M0);
        function dM = DifEq(~,k)
            dMdt = zeros(8,1);
            dMdt(1) = p(1) - (p(16)*k(3) + p(17)*k(4) + p(18)*k(6)+ p(19)*k(7) +  p(20)*k(8))*k(1)-p(2)*k(1);
            dMdt(2) = (p(16)*k(3) + p(17)*k(4) + p(18)*k(6)+ p(19)*k(7) +  p(20)*k(8))*k(1) - (p(2) + p(3))*k(2);
            dMdt(3) = p(3)*k(2) - (p(2) + p(4) + p(5))*k(1);
            dMdt(4) = p(5)*k(3) - (p(6) + p(7) + p(2))*k(4);
            dMdt(5) = p(7)*k(4) - p(2)*k(5);
            dMdt(6) = p(4)*k(3) + p(6)*k(4) -p(8)*k(6);
            dMdt(7) = p(8)*k(6) - p(9)*k(7);
            dMdt(8) = p(10) +  p(14)*k(3) + p(15)*k(4) + p(11)*k(6) + p(12)*k(7) -p(13)*k(8);
            dM = dMdt;
        end
        M = m;
    end
%From day Aug to nov 20
t = [1
    2
    3
    4
    5
    6
    7
    8
    9
    10];
%Official data for I, D for 10 days
k = [24 23 
    62 35 
    68 41 
    70 42 
    70 42 
    70 43 
    68 49 
    67 49
    67 49
    66 49]; 
% [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

function fv = ftns(p,t)
    SEIHRDBPout = SEIHRDBP(p,t);
    fv = norm(k-SEIHRDBPout(:,[3 6]));
end
% ftns = @(p) norm(k-EIDAGG2(p,t));
PopSz = 50;

% Parameter = 12;
Parameter = 28;

opts = optimoptions('ga','PopulationSize',PopSz,'InitialPopulationMatrix',randi(1E+4,PopSz,Parameter).*[ones(1,12)*1E-3 ones(1,7)*10],'MaxGenerations',2E3)%,'PlotFcn','gaplotbestf');
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n',t0)

[p,fval,exitflag,output] = ga(@(p)ftns(p,t),Parameter,[],[],[],[],zeros(Parameter,1),[ones(20,1); Inf(8,1)],[],[],opts);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n',t1)
GA_Time = etime(t1,t0);
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

fprintf(1,'\n\tRate Constants:\n')
for k1 = 1:12
    fprintf(1,'\t\tp(%2d) = %11.8f\n',k1,p(k1))
end
fprintf(['\n\tInitial Conditions:\n'])
for k1 = 1:8
    fprintf('\t\tdMdt(%d) = %11.8f\n',k1,p(k1))
end

tv = linspace(min(t),max(t));
Mfit = SEIHRDBP(p,tv);

figure(1)
plot(t,k,'p')
hold on
hlp = plot(tv,Mfit);
hold off
grid
xlabel('time')
ylabel('Population')
legend(hlp,'S','E','I','H','R','D','B','P',Location','best')
end

[Star Strider]

I posted a reply to this in your other question Parameter estimation of SEIHRDBP model a few minutes ago.

[University]

Thank you so much. I was about to reply you. I have gotten it.

[Star Strider]

As always, my pleasure!