Need help with ODE parameter estimation
13 views (last 30 days)
Show older comments
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.
0 Comments
Answers (2)
Torsten
on 2 May 2022
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.
5 Comments
Star Strider
on 3 May 2022
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).
Star Strider
on 3 May 2022
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)));
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 = 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.
.
4 Comments
Star Strider
on 4 May 2022
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.
See Also
Categories
Find more on Biological and Health Sciences in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!