Which optimizer to use for SIR Epidemic Model?

Hi,

Some background: The SIR compartmental model is used to model the spread of an epidemic. The population is split into 3 bins: Susceptible, Infectious, Recovered, and the transition from S --> I and from I --> R are governed by a set of differential equations. There are two parameters in this system of ODEs: beta (transmission rate) and gamma (recovery rate). I've written a code which discretizes the ODEs and it predicts the number of people in each bin during each timestep (1 day), it runs for a set number of timesteps/days (e.g. 200) and produces a plot with 3 curves showing the evolution of S,I,R numbers throughout the simulation.

I now have some data for S/I/R and want to fit the model to the data by varying the 2 parameters (beta and gamma) and minimizing a cost function (which is just a simple least squares - sum of (predicted - actual)^2 )

Is there a black box optimizer (no functions and no derivatives) I can use? Basically just need something to iterate through 2 arrays of parameters (for example beta = [0.2:0.01:0.4] and similarly for gamma) to minimize the cost function, I've read the documentation on fmincon, fminunc, patternsearch, lsqnonlin etc but am quite lost, seems like alot of re-coding? Also, I intend to expand the model (SEIRD) to have 5 bins, 5 ODEs and 5 parameters so I need something which works well in higher dimensions too.

Any recommendations would be greatly appreciated!

Thanks.

Answer 1

[Bjorn Gustavsson]

One recommendation that you're not likely to like: The constant transition-probability of recovering/dieing has to be unbiological and should be changed to something more adequate. A constant transition-rate leads to exponential response - which cannot be correct (sure everyone and their aunts use this all over the place but that is not an excuse but rather an indictment on lazy/arrogant mathematicians) which implies that the largest number of infected recover the first day - which seems optimistic, and then the tail is most likely way too long with far too many infected remaining after 2, 4, 8 etc "half-recovery-times". Before you can do a parameter fitting and present the results with a straight back and face you'll have to make the model at least sensibly good.

For fitting parameters to systems of ODEs you can have a look at these links:

nonlinear-data-fitting-for-ode

fitting-parameters-in-ode-for-a-kinetic-reaction

parameter-estimation-for-a-system-of-differential-equations

monod-kinetics-and-curve-fitting

In essence once you've written a function that integrates the ODE-system for a set of initial conditions and parameters and calculate the sum-of-squared-residuals you're good to go with any of the optimization-functions.

HTH

[Star Strider]

@Bjorn Gustavsson — Thank you!

[Bjorn Gustavsson]

@Star Strider, you're welcome (though I'm not exactly sure for what it is good to have done something making another happy enough to say "thank you"...)

[Star Strider]

@Bjorn Gustavsson — You cited to several of my Answers!

[Shaun Tan]

Hi, this is the code I have so far:

So what I'm trying to do is to optimize beta, Ti, delta, and kappa so that the output matches the real data as closely as possible. (I've expanded it to SIRD instead of SIR)

@ Star Strider: I've also tried to implement one of your suggestions on an earlier thread

https://uk.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

but I'm not sure it's the right approach because I don't see where I can input the real values for comparison and also uncertain about the sizes of B, Sv, t etc

The codes are attached below.

------- MAIN CODE -------------

% Model Constants

N = 6.722*10^7; % Total population N = S + I + R + D

I0 = 93693; % Initial number of infected

R0 = 4215757;

D0 = 127767;

T = 100; % Trial Period

dt = 1; % Time interval (of a day)

% Model Variables

beta = 1.37e-9; % infection rate

Ti = 15; % recovery rate

delta = 0; % immunity loss rate

kappa = 0.01; % mortality rate

r0 = N*beta*Ti

% Calculate the model

tt = 0:dt:T-dt;

[S,I,R,D] = sird_d1(beta,Ti,delta,kappa,N,I0,R0,D0,T,dt);

% Real Values

real=[62782783 62780680 62778567 62776499 62774415 62772260 62770037 62767766 62765433 62763050 62760611 62758105 62755513 62752766 62749845 62746754 62743519 62740163 62736729 62733019 62728988 62724605 62719906 62714965 62709693 62703856 62697671 62691207 62684495 62677535 62670362 62662968 62655235 62647098 62638565 62629612 62620286 62610655 62600606 62589722 62577909 62565076 62551169 62536158 62520045 62502753 62483776 62463288 62441312 62417933 62393411 62367805 62341264 62313802 62285228 62255299 62224620 62193023 62160808 62127829 62093500 62057547 62019307 61977153 61931853 61885129 61837623 61789653 61742775 61698597 61658705 61622665 61588821 61556939 61527559 61499611 61472416 61445627 61419266 61393236 61367386 61341552 61315608

93693 91691 90571 89860 89792 89101 88687 87947 87379 86975 87296 88419 87599 87397 87381 87263 87750 89613 91476 91977 92750 93792 94898 96003 99261 103822 107555 110378 114087 117796 123095 128844 135481 141774 147550 153380 158638 165359 174065 182926 191786 198913 207529 216639 230621 245353 258457 230345 240576 256081 271585 291453 313670 334742 353207 369523 388957 409075 433332 457892 483753 509613 531837 554746 580321 611643 650286 690959 733534 765151 788013 816169 844325 869078 882449 893745 896999 894132 890588 893388 900282 900338 900393

4215757 4219848 4223070 4225844 4227990 4230829 4233458 4236462 4239357 4242139 4244252 4245628 4249036 4251978 4254906 4258111 4260848 4262333 4263896 4267098 4270348 4273682 4277272 4281097 4283099 4284366 4286813 4290448 4293446 4296689 4298556 4300190 4301280 4303111 4305858 4308972 4313021 4315920 4317257 4319268 4322197 4327880 4333155 4339043 4341167 4343706 4349559 4398142 4409867 4417729 4426723 4432440 4436743 4443105 4453186 4466771 4477988 4489435 4497364 4505745 4514167 4524217 4540200 4559406 4579082 4594442 4603253 4610500 4614751 4627246 4644193 4652019 4657653 4664706 4680634 4697214 4721084 4750666 4780499 4803642 4822519 4848193 4873992

127767 127781 127792 127797 127803 127810 127818 127825 127831 127836 127841 127848 127852 127859 127868 127872 127883 127891 127899 127906 127914 127921 127925 127935 127947 127956 127961 127967 127972 127980 127987 127998 128004 128017 128027 128036 128055 128066 128072 128085 128108 128131 128147 128160 128167 128188 128208 128225 128245 128257 128281 128302 128323 128351 128379 128407 128435 128467 128496 128534 128581 128623 128656 128695 128744 128786 128838 128888 128940 129006 129089 129147 129201 129277 129358 129430 129501 129575 129647 129734 129813 129918 130007

];

% Least Squares Cost Function

vals=[S;I;R;D];

[r,c]=size(real);

for i=1:r

for j=1:c

squares(i,j)=real(i,j)-vals(i,j);

end

end

cost=sum(squares.^2,'all')

% Plot SIRD Visualization

tr=1:c;

figure (1)

plot(tt,S,'-b',tt,I,'-r',tt,R,'-g',tt,D,'-k','LineWidth',1);

hold on

plot(tr,real(1,:),'.b',tr,real(2,:),'.r',tr,real(3,:),'.g',tr,real(4,:),'.k')

grid on; grid minor

xlabel('Days'); ylabel('Number of individuals');

legend('S','I','R','D');

% ylim([0 6e6])

toc

---------- function "sird_d1" which propagates the ODE -----------

function [S,I,R,D] = sird_d1(beta,Ti,delta,kappa,N,I0,R0,D0,T,dt)

% if delta = 0 we assume a model without immunity loss

S = zeros(1,T/dt);

S(1) = N-I0-R0-D0;

I = zeros(1,T/dt);

I(1) = I0;

R = zeros(1,T/dt);

R(1) = R0;

D = zeros(1,T/dt);

D(1) = D0;

for tt = 1:(T/dt)-1

% Equations of the model

dS = (-beta*I(tt)*S(tt) + delta*R(tt)) * dt;

dI = (beta*I(tt)*S(tt) - 1/Ti*I(tt)) * dt;

dR = ((1-kappa)/Ti*I(tt) - delta*R(tt)) * dt;

dD = (kappa/Ti*I(tt)) * dt;

S(tt+1) = S(tt) + dS;

I(tt+1) = I(tt) + dI;

R(tt+1) = R(tt) + dR;

D(tt+1) = D(tt) + dD;

end

end

----------- Adapting Star Strider's Monod Kinetics Code -----------------

function S = MonodKinetics1(B, t)

% MONODKINETICS1 codes the system of differential equations describing one

% version of Monod kinetics:

% dS/dt = -beta*I*S + delta*R;

% dI/dt = beta*I*S - 1/Ti*I;

% dR/dt = (1-kappa)/Ti*I - delta*R;

% dD/dt = kappa/Ti*I

% with:

% Variables: S = x(1), I = x(2), R = x(3), D = x(4)

% Parameters: beta = B(1), Ti = B(2), delta = B(3), kappa = B(4)

x0 = rand(2,1);

[T,Sv] = ode45(@DifEq, t, x0);

function dS = DifEq(t, x)

xdot = zeros(4,1);

xdot(1) = -B(1) .* x(2) .* x(1) + B(3) .* x(3);

xdot(2) = B(1) .* x(2) .* x(1) - 1/B(2).*x(2);

xdot(3) = (1-B(4))./B(2).*x(2) - B(3) .* x(3);

xdot(4) = B(4)./B(2).*x(2);

dS = xdot;

end

S = Sv(:,1);

end

Thanks for your help!

Shaun

[Bjorn Gustavsson]

1, No attention to my comment on the oversimplification of this type of constant-rate SIR-type modelling of biological systems?

2, Write an error-function where you integrate the SIR-model using one or the other of the odeNN-functions for the and calculates the residuals relative to the observed data:

function err = err_of_sir_too_simplified(pars,t_obs,SIR_obs)

SIR0 = pars(1:3); % couldn't be bothered to write this for a SIRD-model
alphabeta = pars(4:5); % if you need that the modification should be obvious.
% Next we model the SIR-evolution from the initial state for the given
% transition-rates:
[t_mod,SIR_mod] = ode45(@(t,SIR) ode_sir_too_simplified(t,SIR,alphabeta),t_obs,SIR0);
% And calculates the sum of the squared residuals...
err = sum(sum((SIR_obs-SIR_mod).^2));
end

Where your SIR-modeling-function ought to look something like this:

function dSIRdt = ode_sir_too_simplified(t,SIR,alphabeta)
S = SIR(1);
I = SIR(2);
R = SIR(3);
alpha = alphabeta(1);
beta  = alphabeta(2);
dSIRdt = [-alpha*S*I;
          alpha*S*I-beta*I;
          beta*I];
end

Then after writing these 2 equations you can run the (pointless due to unbiologicalness - you REALLY should ask your teacher about that. If you're paid for this then just: NO.) parameter fitting with any of the general minimization-functions, for example fminsearch:

SIR0 = [12356,3432,12]; % Adjust this to what fits your population estimates
alphabeta = [0.03,0.1]; % These are your transition-rate-guesses
load T_obs.mat
load SIR_obs.mat % You'll have to do something...

par0 = [SIR0,alphabeta]; % Initial parameters for the fitting:
par_est = fminsearch(@(par) err_of_sir_too_simplified(par,T_obs,SIR_obs),par0);
% That should give you a fit for the initial conditions and the
% transfer-rates that best fits the simplified SIR-model to the
% observations.
SIR0_est = par_est(1:3);
alphabeta_est = par_est(4:5);
[t_mod,SIR_mod] = ode45(@(t,SIR) ode_sir_too_simplified(t,SIR,alphabeta_est),T_obs,SIR0_est);
plot(t_obs,SIR_mod)
hold on
plot(t_obs,SIR_obs,'*')

HTH