Using lsqcurvefit with ODE15s (SIR model)

So I am currently working on fitting COVID-19 data to the SIR model. This model consists of a series of differential equations. I am taking a somewhat simplified approach given data constraints and my goals for this project. As such, I have data for the Infected category but unfortunately cannot get accurate data for the Removed population (that incorporates death and recovery). For the Susceptible population, I am simply fixing a country's population and then subtracting off the newly infected. My issue is that I am really struggling using lsqcurvefit to determine the parameters associated with the differential equations. I've read a lot on the function as well as almost this specific application (major shoutout to @Star Strider), yet I still can't get it to work. If anyone could help me out that would be greatly appreciated! I'm trying to avoid using the matrix definition of least squares solution as I'd like to be able to use lsqcurvefit for future projects.

clear;

%Read in complete data set and then label coloumn variables
fullData = readtable('full_data.csv','HeaderLines',1);
fullData.Properties.VariableNames = {'Date' 'Country' 'NewInf' 'NewDea' 'TotInf' 'TotDea'};
usData_nz = fullData(strcmp(fullData.Country, 'United States'),:);
usData = usData_nz(usData_nz.TotInf > 0,:); %Zero the data to when first 100 people effected
usDays = 1:1:size(usData,1);
usDays_T = array2table(usDays', 'VariableNames', {'DaysInf'});
usData = [usData usDays_T];
usN = 329436300; %Population of the United States
usSusc(1) = usN - usData.NewInf(1);
for i=2:size(usData,1)
  usSusc(i) = usSusc(i-1) - usData.NewInf(i);
end
usSusc_T = array2table(usSusc', 'VariableNames', {'Susceptible'});   %Assumes removed are immune
usData = [usData usSusc_T];
usN_vec = ones(size(usData,1),1).*usN;
currInf(1) = usData.TotInf(1);
for i=2:size(usDays,2)
  currInf(i) = currInf(i-1) + usData.NewInf(i) - usData.NewDea(i);
end
usCurrInf_T = array2table(currInf', 'VariableNames', {'CurrentInf'});   %Assumes removed are immune
usData = [usData usCurrInf_T];
usSI = [usData.Susceptible usData.CurrentInf usData.TotDea];
param = knfit(usData.DaysInf,usSI);
B0 = [param; usSI(1,:)'];
usFitSIR = sir_Sys(B0,usData.DaysInf);
usFitSIR = usFitSIR./usN;
%% Plotting the Data

usPercSusc = usData.Susceptible./usN;
usPercInf = usData.TotInf./usN;
figure(1);
plot(usData.DaysInf, usPercInf,'r')%,usData.DaysInf, usPercSusc,'b');
hold on
plot(tv, usFitSIR(:,2), 'r--')%,tv, usFitSIR(:,1), 'b--');
hold off

%% Functions

function S = sir_Sys(B,t)
% Function to calculate derivatives of the SIR model
SI0 = B(3:5);
     [~,Sv] = ode15s(@sirODE, t, SI0);
        function xdot = sirODE(t, x)
            xdot = zeros(3,1);
            xdot(1) = -1.*B(1).* x(1).* x(2) ;
            xdot(2) =  B(1).* x(1).* x(2) - B(2).*x(2);
            xdot(3) = -1.*B(2).*x(2);
        end
%S = Sv(:,1);
S = Sv;
end


function x = knfit(t,f)
% t = time
% f = raw data
objfcn = @(B,t) prefun(B,t);
B0 = [0 0]'; % initial values
x = lsqcurvefit(objfcn,B0,t,f);  
  function S = prefun(B,t)
    x0 = f(1,:)';
    [~,Sv] = ode15s(@DiffEq,t,x0);
    function xdot = DiffEq(t,x)
       xdot = zeros(3,1);
       xdot(1) = -1.*B(1).* x(1).* x(2) ;
       xdot(2) =  B(1).* x(1).* x(2) - B(2).*x(2);
       xdot(3) = -1.*B(2).*x(2);
    end
    S = Sv;
  end
end

[Jincy A]

Dear star strider

The code to estimate the parameters not ending. Will you pl rectify

function abraham

%

function C=model1(B,t)

c0=[ 5078931 2 0];

[T,Cv]=ode45(@DifEq,t,c0);

%

function dS = DifEq(t, x)

xdot = zeros(3,1);

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

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

xdot(3) = B(2).*x(2);

dS = xdot;

end

C=Cv;

end

%%

t=1:66;

ydata=[0.999999606216503,3.93783497439324e-07,0;0.999999212433005,7.87566994878649e-07,0;0.999998818649508,1.18135049231797e-06,1.96891748719662e-07;0.999997637299015,2.36270098463595e-06,3.93783497439324e-07;0.999997046623769,2.95337623079493e-06,3.93783497439324e-07;0.999996849732021,3.15026797951460e-06,3.93783497439324e-07;0.999996455948523,3.54405147695392e-06,1.77202573847696e-06;0.999996259056774,3.74094322567358e-06,1.77202573847696e-06;0.999995668381528,4.33161847183257e-06,1.77202573847696e-06;0.999995274598031,4.72540196927189e-06,1.77202573847696e-06;0.999992321221800,7.67877820006683e-06,2.36270098463595e-06;0.999990549196062,9.45080393854378e-06,2.36270098463595e-06;0.999989761629067,1.02383709334224e-05,2.36270098463595e-06;0.999989170953820,1.08290461795814e-05,3.34715972823426e-06;0.999987005144585,1.29948554154977e-05,3.34715972823426e-06;0.999983461093108,1.65389068924516e-05,3.34715972823426e-06;0.999980507716877,1.94922831232466e-05,3.34715972823426e-06;0.999978538799390,2.14612006104432e-05,4.52851022055223e-06;0.999974207180918,2.57928190822757e-05,4.52851022055223e-06;0.999970072454195,2.99275458053887e-05,4.52851022055223e-06;0.999967119077964,3.28809220361836e-05,4.52851022055223e-06;0.999964756376979,3.52436230208195e-05,5.70986071287020e-06;0.999962196784246,3.78032157541751e-05,6.69431945646851e-06;0.999958652732769,4.13472672311291e-05,6.69431945646851e-06;0.999954518006046,4.54819939542420e-05,1.12228296770207e-05;0.999950383279323,4.96167206773549e-05,1.12228296770207e-05;0.999945460985605,5.45390143953464e-05,1.20103966718994e-05;0.999941326258882,5.86737411184593e-05,1.20103966718994e-05;0.999934828831174,6.51711688262082e-05,1.20103966718994e-05;0.999926953161225,7.30468387749947e-05,1.31917471642174e-05;0.999917502357287,8.24976427135385e-05,1.41762059078157e-05;0.999910020470835,8.99795291648856e-05,2.14612006104432e-05;0.999904704393620,9.52956063803165e-05,2.14612006104432e-05;0.999892497105199,0.000107502894800936,2.14612006104432e-05;0.999882061842517,0.000117938157483078,2.14612006104432e-05;0.999856859698681,0.000143140301319194,2.44145768412381e-05;0.999839927008291,0.000160072991709085,2.55959273335561e-05;0.999820828508665,0.000179171491334893,2.57928190822757e-05;0.999799367308055,0.000200632691945336,3.46529477746605e-05;0.999789522720619,0.000210477279381319,3.46529477746605e-05;0.999767667736511,0.000232332263489201,3.46529477746605e-05;0.999750735046121,0.000249264953879092,4.58757774516813e-05;0.999722382634305,0.000277617365694724,4.68602361952796e-05;0.999703087242931,0.000296912757069251,4.68602361952796e-05;0.999682216717567,0.000317783282433535,4.68602361952796e-05;0.999662133759197,0.000337866240802940,6.04457668569363e-05;0.999647563769792,0.000352436230208195,6.39898183338902e-05;0.999624921218689,0.000375078781310956,6.47773853287689e-05;0.999606610286058,0.000393389713941885,6.55649523236475e-05;0.999596568806873,0.000403431193126588,7.16685965339570e-05;0.999580423683478,0.000419576316521600,7.16685965339570e-05;0.999552268163412,0.000447731836588512,7.16685965339570e-05;0.999537698174006,0.000462301825993767,9.74614156162328e-05;0.999518205890883,0.000481794109117013,9.74614156162328e-05;0.999511117787929,0.000488882212070921,0.000147668811539747;0.999494381989288,0.000505618010712093,0.000147668811539747;0.999480796458626,0.000519203541373749,0.000160663666955244;0.999461501067252,0.000538498932748276,0.000168933120401470;0.999428423253467,0.000571576746533179,0.000174839872863060;0.999417594207287,0.000582405792712761,0.000192953913745269;0.999387272877985,0.000612727122015589,0.000201814042437654;0.999365221002128,0.000634778997872191,0.000210280387632599;0.999330764946102,0.000669235053898132,0.000219928083319863;0.999296505781825,0.000703494218175353,0.000246114685899578;0.999267365803014,0.000732634196985863,0.000246114685899578;0.999208692061896,0.000791307938104322,0.000253793464099645]

beta1=0.005 + (0.01-0.005).*rand(1,1);

gamma1=0.2 + (0.5-0.2).*rand(1,1);

lb=[0,0,0];

ub=[1,1,1];

B0=[beta1; gamma1];

%%

[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jma]=lsqcurvefit(@model1,B0,t(:),ydata);

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

for k1 = 1:length(B)

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

end

tv = linspace(min(t), max(t));

Cfit =model1(B, tv);

figure(1)

plot(t, c, 'p')

hold on

hlp = plot(tv, Cfit);

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')

end

Answer 1

[Star Strider]

When I run your code, It throws:

Warning: Failure at t=7.164307e+01.  Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (2.273737e-13) at time t. 

also citing:

param = knfit(usData.DaysInf,usSI);

and then stops.

The error is likely in the initial conditions or the ODE code itself. What is the original model you¹re using? (Please post a .pdf of it.)

There may be other problems, but the first is to get the ODE to integrate correctly.

[Nick Kirkpatrick]

I've attached the ODE I'm trying to model. in my code: x(1) = S, x(2) = I, and x(3) = R.

[Star Strider]

The ‘xdot(3)’ expression is not negated in the image.

Should it instead be:

xdot(3) = B(2).*x(2);

I wonder?

However, when I make that change, I continue to get a similar Warning to the earlier on, although at a different time.

Checking the initial conditions is likely appropriate at this point.

[Nick Kirkpatrick]

Yes, I forgot to update the code in this post but I did notice that before and made the appropriate change with similar results. Alright I will work on the initial conditions, thank you for the help.

[Star Strider]

My pleasure!

[Nick Kirkpatrick]

Fortunately I was able to get the code to run but no matter how much I change the initial parameter estimate, the model remains at zero. Additionally, for certain parameter values it throws an error that the matrix sizes do not agree. I've checked and the parameter values it finds are highly unrealistic, is this more likely due to poor data or an error in the code? I've attached a plot to show, the data is solid red and the model is the dashed line. It is also telling me a local minum is occurring. I'm assuming th ebest apporach is to continue experimenting with parameters.

[Star Strider]

The matrix sizes not being equal is likely due to the ode15s stopping early, probably because it encounters a singularity, so it doesn’t integrate for the entire time vector. (That’s similar to the original problem that threw the Warning, and it should throw the Warning before it throws the Error.)

I got it to do something with:

B0 = rand(2,1); % initial values

in ‘knfit’ since it is never good practice for initial parameter estimates to be zero, however it clearly does not approximate the curve.

That’s the best I can do.

EDIT — (16 Apr 2020 at 11:59)

If you are using the model: quoted in this post, you need to normalise the first two equations by ‘N’.

[Nick Kirkpatrick]

Hey just wanted to follow back and say thanks for your help! I guess the best way to describe it is that I was able to successfully implement lsqcurvefit, just with unsuccessful results. The model is highly sensitive to the parameter changes and particularly the initial estimates/bounds. In the end, I don't think the data is complete enough to effectively estimate the parameters. I tried other methods (linearized regions, matrix definition of lss, etc.) with similar results. Thanks again for the input!

[Star Strider]

As always, my pleasure!

If you have the Global Optimization Toolbox, see if using the ga (genetic algorithm) function will give better resullts. You may need to run it several times before it converges successfully on the best parameter set. I have attached my prototype code for that to this Comment.