Need help with ODE parameter estimation

I am trying to develop an SIRD model for modeling covid infections over time. I have developed 4 differential equations that show the change in each population.

1) dS/dt = -beta*S*I (the change in susepctible population)

2) dI/dt= beta*S*I - lambda*I - kd*I (the change in infected population)

3) dR/dt= lambda*I (the change in recovered population)

4) dD/dt= kd*I (the change in dead population)

The only problem is that I do not have the values for the 3 constant parameters (beta, lambda, and kd). beta represents the change from susceptible to infected, lambda represents the change from infected to recovered, and kd represents the transition from infected to dead.

I have data that shows the infected cases reported every day over a series of days. So my question is, how do I estimate the parameters using this data. I have looked at similar links but am still having a hard time putting it together.

Any help is much appreciated, thank you.

Answer 1

[Torsten]

I have data that shows the infected cases reported every day over a series of days.

This will not be enough to separate the different sources for change of I. You must have data for the other groups as well.

[Kevin Burke]

I actually do have data on the recovered and dead as well. My apologies.

[Torsten]

Try to adapt Star Strider's reference code for parameter fitting in ordinary differential equations under

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

according to your needs.

[Kevin Burke]

function S = MonodKinetics1(B, t)
% MONODKINETICS1 codes the system of differential equations describing one
%   version of Monod kinetics:
%       dS/dt = (mmax * X * S)/(Y * (Ks + S));
%       dX/dt = (mmax * X * S)/ (Ks + S) - (b * X);
%   with:
%       Variables:  x(1) = S,  x(2) = X
%       Parameters: mmax = B(1),  Y = B(2),  Ks = B(3),  b = B(4)
x0 = rand(2,1);
[T,Sv] = ode45(@DifEq, t, x0);
      function dS = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = (B(1) .* x(2) .* x(1)) ./ (B(2) .* (B(3) + x(1)));
          xdot(2) = (B(1) .* x(2) .* x(1)) ./ (B(3) + x(1)) - (B(4) .* x(2));
          dS = xdot;
      end
S = Sv(:,1);
end

I am a little confused by this part. The poster had two differentials: (dS/dt and dX/dt) but the imbedded function dS=difEq only accounts for one of them correct? Since i have 4 differntial equations, do I need 4 imbedded functions or just one similar to how Star Strider did it? I am a MATLAB beginner so forgive me if this question doesnt make any sense.

Thanks,

Kevin

[Star Strider]

The Monod Kinetics code only needed to fit one result. To fit all of them (or a subset of them) return all of ‘S’ (or selected columns).

Peruse this search for similar and more general approaches.

[Kevin Burke]

Hello Star Strider, thank you for the response. I have been reading your old comments all day, I feel like a celebrity just responded to me lol. Do you think you could help me to fix my code?

My goal is to solve for the three parameters Beta, lambda and kd. I have data for 31 days in March 2020 (US) that gives the ammount of susepticnle (S), infected (I), recovered (R) and dead (D) people each day. Do you have any insight into how this could be accomplished?

Thanks

Kevin

Call command

clc
clear all

P0=[.01; .01; .01]; %Parameter initial guesses
US_Population = 330000000; %US population
USData = xlsread('USData.xlsx'); %Importing data
IRD=USData([40:70],[2:4]); %Isolating IRD March 2020 data
time=[1:31]'; %Creating a 31 day column vector
Infected = IRD(:,1); %Isolates cases column vector
S_data = US_Population - Infected; %S(t): Population - I(t)
US_March_2020=[S_data IRD];
[P,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Objective,P0,time,S_data,[],[]);

Local objective Function

%Variables: x(1)=S, x(2)=I, x(3)=R, x(4)=D
%Parameters: beta = P(1), lambda = P(2), kd = P(3)

%Adding this duplicate line of code so that X0 is definied in the function
US_Population = 330000000; %US population
USData = xlsread('USData.xlsx'); %Importing data
IRD=USData([40:70],[2:4]); %Isolating IRD March 2020 data
Infected = IRD(:,1); %Isolates cases column vector
S_data = US_Population - Infected; %S(t): Population - I(t)
US_March_2020=[S_data IRD];
x0= US_March_2020(1,:);

opts = odeset('Refine', 10);
[T,Sv] = ode15s(@SIRD, t, x0, opts);

    function dS = SIRD(t,x)
    xdot = zeros(4,1);
    xdot(1) = (-P(1).*x(1).*x(2));
    xdot(2) = (P(1).*x(1).*x(2)-P(2).*x(2)-P(3).*x(2));
    xdot(3) = (P(2).*x(2));
    xdot(4) = (P(3).*x(2));
    dS=xdot;
    end

S = Sv(:,1);
    

end

Edit: Changed initial values. Solver "works" but output parameters are exactly equal to inital values.

Answer 2

[Star Strider]

I appreciate your compliment!

There are a number of SIR models posted that estimate parameters from data.

Modify this kinetic model to work with your model and data (this is an update to previously posted code) —

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 found.

Optimization completed because the size of the gradient is less than
the value of the optimality 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) =  0.76482
		Theta(2) =  0.23492
		Theta(3) =  0.20876
		Theta(4) =  0.49178
		Theta(5) =  0.62211
		Theta(6) =  0.00000
		Theta(7) =  0.90288
		Theta(8) =  0.07146
		Theta(9) =  0.02840
		Theta(10) =  0.00000

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')

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

The initial conditions (the last 4 elements of parameter vector ‘theta’) are estimated as parameters.

It might be necessary to change to ode15s if the parameters vary over a wide range of magnitudes, creating a ‘stiff’ system.

.

[Kevin Burke]

Can you explain why theta0 = rand(10,1)?

[Star Strider]

Yes.

Those are the initial values for the parameter vector.

[Kevin Burke]

Thank you very much for the help. I have just one more question.

I have calcualetd the parameters using the provided data set. Now I want to change one parameter in order to induce a change on the data set. I want to lower the transmission parameter (beta) and I am hoping that will lower the I(t) values.

Thanks.

[Star Strider]

My pleasure!

If you have the solved system, save the parameters to a file (simply to avoid having to re-run the optimsation routine later), change the required parameter in the ‘theta’ vector, and run the ‘kinetics’ function with the new ‘theta’ vector and time vector (that can now be whatever you want it to be, however the time units must be the same).

I did that in my code with:

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

Cfit = kinetics(theta, tv);

to draw the plot of the fitted estimates with respect to the data.

Other options are also the deval and odextend functions.