ode_covid19_v3.m 2.44 KB
%ODE model to describe the spread and clinical progression of COVID-19

function dx=ode_covid19_v3(t,x,p,Tlock,s00,s11,s22)%,Tlock2,s02,s12)

be = p(1);
b0=p(2);
b1=p(3);
b2=p(4);
b3=p(5);
        
Norm=10^5;

be=be/(Norm);
b0=b0/(Norm);
b1=b1/(Norm);
b2=b2/(Norm);
b3=b3/(Norm);
        
FracSevere = p(6);
FracCritical = p(7);
FracMild = 1-FracSevere-FracCritical;
FracAsym = p(8);
               
IncubPeriod = p(9);
DurMildInf = p(10);

DurAsym = p(11);
   
DurHosp = p(12);
TimeICUDeath = p(13);
ProbDeath=p(14);

PresymPeriod=p(15)*IncubPeriod;

           
CFR=ProbDeath*FracCritical/100; 
            
a1=1/PresymPeriod; %presymptomatic period of transmission
a0=1/(IncubPeriod-PresymPeriod); %true latent period  

f=FracAsym;
    
g0=1/DurAsym;

g1=(1/DurMildInf)*FracMild;
p1=(1/DurMildInf)-g1;
    
p2=(1/DurHosp)*(FracCritical/(FracSevere+FracCritical));
g2=(1/DurHosp)-p2;
    
u=(1/TimeICUDeath)*(CFR/FracCritical);
    
g3=(1/TimeICUDeath)-u;


s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3)) * (t<Tlock(4)) + s00(3)*(t>=Tlock(4));
s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s11(2)*(t>=Tlock(3))* (t<Tlock(4)) + s11(3)*(t>=Tlock(4));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
%s3 = (t<Tlock(1)) + s33 * (t>=Tlock(1));

%susceptibles individuals (not infected)
S = x(1);
%exposed individuals, infected but no symptoms and no transmission
E0 = x(2);
%exposed individuals, infected but no symptoms and can transmit
E1 = x(3);
%infected individuals but asymptomatic infection
I0 = x(4);
%infected individuals, mild infection
I1 = x(5);
%infected individuals, severe infection
I2 = x(6);
%infected individuals, critical infection
I3 = x(7);
%recovered
R = x(8);
%dead
D = x(9);


%dS = - (be *(1/(1+s0^n))* E1 + b0*(1/(1+s0^n))*I0 + b1*(1/(1+s1^n))*I1 + b2*(1/(1+s2^n))*I2 + b3*(1/(1+s2^n))*I3)*S; % + k_e0s * E0;
%dE0 = (be *(1/(1+s0^n))* E1 + b0*(1/(1+s0^n))*I0 + b1*(1/(1+s1^n))*I1 + b2*(1/(1+s2^n))*I2 + b3*(1/(1+s2^n))*I3)*S -a0*E0; % - k_e0s * E0; 
dS = - ((be * s0)* E1 + (b0*s0)*I0 + (b1*s1)*I1 + (b2*s2)*I2 + (b3*s2)*I3)*S; % + k_e0s * E0;
dE0 = ((be *s0)* E1 + (b0*s0)*I0 + (b1*s1)*I1 + (b2*s2)*I2 + (b3*s2)*I3)*S -a0*E0; % - k_e0s * E0; 
dE1 = a0 * E0 - a1*E1; % -g4*E1;
dI0 = f * a1 * E1 - g0*I0;
dI1 = (1-f) * a1 * E1 - g1*I1 - p1*I1;
dI2 = p1 * I1 - g2*I2 - p2*I2; % -u2*I2;
dI3 = p2 * I2 - g3 *I3 - u * I3;
dR = g0*I0 + g1*I1 + g2*I2 + g3*I3; % +g4 * E1;
dD = u * I3; % + u2 * I2;



dx=[dS;dE0;dE1;dI0;dI1;dI2;dI3;dR;dD];