Skip to content

SysBiOThe / CRA-Matlab

Go to a project
Toggle navigation
Toggle navigation pinning
Commit a14857d7 by Chiara Antonini
Options

code paper

1 parent 2ef245c7
Inline Side-by-side
Showing 20 changed files with 0 additions and 853 deletions
CRA_Updated_Version/CRA_Toolbox_bash/Examples/ITALIA/init_italia.m deleted 100644 → 0
region_name='Italy_pred';
%define initial conditions
N = 60e6; %population Italy
%define time
%time = 0:1:110;
E00 = (1/N)*10^5;
S0 = 10^5-E00;
x0 = [S0 E00 0 0 0 0 0 0 0];
%inizio epidemia: 27 Gennaio
%dopo 24 giorni (20 Febbraio): paziente 1
%dopo 25 giorni (21 Febbraio): zona rossa in Lombardia e Veneto
%dopo 28 giorni (24 Febbraio): chiusura scuole in Emilia-Romagna, Friuli-Venezia Giulia,
%Liguria, Lombardia, parte delle Marche, Piemonte, Veneto
%dopo 25 giorni (21 Febbraio) DPI
%dopo 38 giorni (5 Marzo): chiusura scuole nel resto d'Italia
%dopo 40 giorni (7 Marzo): fuga verso sud
%dopo 41 giorni (8 Marzo): lockdown nord Italia + social distancing +
%eventi vietati
%dopo 43 giorni (10 Marzo): total lockdown
%RILASCIO DEL LOCKDOWN:
%dopo 98 giorni (4 Maggio 2020)
%dopo 112 giorni (18 Maggio 2020)
Tlock=[25 28 38 43 98 112];
%s0 = [0.8 0.75 0.7 0.2]; %reduction of 30%, 70%, 80%
%s1=[0.7 0.65 0.6 0.15]; %reduction of 40%, 75%, 85%
%s2=0.4; %reduction of 60%
\ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/Examples/compute_intervention.m deleted 100644 → 0
function [s0,s1,s2]=compute_intervention(t,Tlock,s00,s11,s22,region_name)
if(strcmp(region_name,'Italy'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3))*(t<Tlock(4)) + s00(4)*(t>=Tlock(4));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Italy_pred'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3))*(t<Tlock(4)) + s00(4)*(t>=Tlock(4))*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Umbria'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Campania'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3)); %*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Lazio'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Lombardia'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(1)) + s11 * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Veneto'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11 * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'EmiliaRomagna'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11 * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Calabria'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Abruzzo'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Basilicata'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Marche'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Piemonte'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Toscana'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Puglia'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Liguria'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Sicilia'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Sardegna'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Friuli'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(2)) + s22 * (t>=Tlock(2));
end
if(strcmp(region_name,'Aosta'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(2)) + s22 * (t>=Tlock(2));
end
if(strcmp(region_name,'Molise'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Trentino'))
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
end
CRA_Updated_Version/CRA_Toolbox_bash/Examples/ode_covid19_v2.m deleted 100644 → 0
%ODE model to describe the spread and clinical progression of COVID-19
function dx=ode_covid19_v2(t,x,p,Tlock1,s01,s11)%,Tlock2,s02,s12)
be = p(1);
b0=p(2);
b1=p(3);
b2=p(4);
b3=p(5);
N=882000;
be=be/(N);
b0=b0/(N);
b1=b1/(N);
b2=b2/(N);
b3=b3/(N);
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);
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;
% be = p(1);
% b0 = p(2);
% b1 = p(3);
% b2 = p(4);
% b3 = p(5);
% a0 = p(6);
% a1 = p(7);
% f = p(8);
% g0 = p(9);
% g1 = p(10);
% p1 = p(11);
% g2 = p(12);
% p2 = p(13);
% g3 = p(14);
% u = p(15);
%seas_amp = p(16);
%seas_phase = p(17);
%s0 = p(18);
%s1=p(19);
%se=p(20);
%k_e0s=p(21);
%g4=p(22);
%u2 = p(23);
%seas=(1 + seas_amp*cos(2*pi*(t-seas_phase)/365));
% if t<Tlock
% s0 = 1;
% s1 = 1;
% elseif t>=Tlock && t<Tlock2
% s0 = int_s0;
% s1 = int_s1;
% else
% s0 = v11
% s1 = v2;
% end
s0 = (t<Tlock1) + s01 * (t>=Tlock1); %* (t<Tlock2) + s01(2)*(t>=Tlock2);
s1 = (t<Tlock1) + s11 * (t>=Tlock1); %* (t<Tlock2) + s11(2)*(t>=Tlock2);
%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 *s0* E1 + b0*s0*I0 + b1*s1*I1 + b2*I2 + b3*I3)*S; % + k_e0s * E0;
dE0 = (be *s0* E1 + b0*s0*I0 + b1*s1*I1 + b2*I2 + b3*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];
CRA_Updated_Version/CRA_Toolbox_bash/Examples/ode_covid19_v3_NO.m deleted 100644 → 0
%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))* (t<Tlock(5)) + s00(4)*(t>=Tlock(5))* (t<Tlock(6)) + s00(5)* (t>=Tlock(6));
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];
CRA_Updated_Version/CRA_Toolbox_bash/Examples/ode_covid19_v3_OK.m deleted 100644 → 0
%ODE model to describe the spread and clinical progression of COVID-19
function dx=ode_covid19_v3_OK(t,x,p,Tlock,region_name) %s00,s11,s22,region_name)%,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;
s00=[p(16) p(17) p(18) p(19) p(20) p(21)];
s11=p(22);
s22=p(23);
[s0,s1,s2]=compute_intervention(t,Tlock,s00,s11,s22,region_name);
%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];
CRA_Updated_Version/CRA_Toolbox_bash/Examples/parameters_EGFR_IGF1R_model.m deleted 100644 → 0
% Model parameters
gamma_EGFR=0.02;
gamma_IGF1R=0.02;
kd_PIK3_active=0.005;
k_p90Rsk_Erk=0.0213697;
KM_p90Rsk_Erk=763523;
k_SOS_EGFR=694.731;
KM_SOS_EGFR=6086070;
k_Ras_SOS=32.344;
KM_Ras_SOS=35954.3;
k_ERK_MEK=9.85367;
KM_ERK_MEK=1007340;
k_DSOS_p90Rsk=161197;
KM_DSOS_p90Rsk=896896;
k_SOS_IGF1R=500;
KM_SOS_IGF1R=100000;
k_PIK3_IGF1R=10.6737;
KM_PIK3_IGF1R=184912;
k_PIK3_EGFR=10.6737;
KM_PIK3_EGFR=184912;
k_Akt_PIK3=0.0566279;
KM_Akt_PIK3=653951;
kd_Akt=0.005;
k_Erk_PP2A=9.85367;
KM_Erk_PP2A=1007340;
k_PIK3_Ras=0.0771067;
KM_PIK3_Ras=272056;
k_Raf_Ras=0.884096;
KM_Raf_Ras=62464.6;
k_Mek_Raf=185.759;
KM_Mek_Raf=4768350;
k_Raf_Akt=15.1212;
KM_Raf_Akt=119355;
k_Ras_RasGab=1509.36;
KM_Ras_RasGab=1432410;
k_MEK_PP2A=2.83243;
KM_MEK_PP2A=518753;
k_Raf_RafPP=0.126329;
KM_Raf_RafPP=1061.71;
kd_P90Rsk=0.005;
% Input values
RafPP=120000;
PP2A=120000;
RasGapActive=120000;
% Total concentrations
DSOS_tot=120000.0;
Ras_tot=120000.0;
Raf_tot=120000.0;
Mek_tot=600000.0;
Erk_tot=600000.0;
p90Rsk_tot=120000.0;
PIK3_tot=120000.0;
Akt_tot=120000.0;
% Initial conditions for state variables
EGFR_active_0=8000.0;
IGF1R_active_0=650.0;
SOS_0=0.0;
Ras_active_0=0.0;
Raf_active_0=0.0;
Mek_active_0=0.0;
Erk_active_0=0.0;
p90Rsk_active_0=0.0;
PIK3_active_0=0.0;
Akt_active_0=0.0;
CRA_Updated_Version/CRA_Toolbox_bash/compute_R0.m deleted 100644 → 0
function R0=compute_R0(param)%,Tlock2,s02,s12)
be = param(1);
b0=param(2);
b1=param(3);
b2=param(4);
b3=param(5);
Norm=10^5; %882000
%Norm=882000;
be=be/(Norm);
b0=b0/(Norm);
b1=b1/(Norm);
b2=b2/(Norm);
b3=b3/(Norm);
FracSevere = param(6);
FracCritical = param(7);
FracMild = 1-FracSevere-FracCritical;
FracAsym = param(8);
IncubPeriod = param(9);
DurMildInf = param(10);
DurAsym = param(11);
DurHosp = param(12);
TimeICUDeath = param(13);
ProbDeath=param(14);
PresymPeriod=param(15)*IncubPeriod;
%PresymPeriod=param(15);
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<Tlock1) + s01 * (t>=Tlock1); %* (t<Tlock2) + s01(2)*(t>=Tlock2);
%s1 = (t<Tlock1) + s11 * (t>=Tlock1); %* (t<Tlock2) + s11(2)*(t>=Tlock2);
R0=Norm*(((be)/a1)+f*((b0)/g0)+(1-f)*(((b1)/(p1+g1))+(p1/(p1+g1))*(b2/(p2+g2)+ (p2/(p2+g2))*(b3/(u+g3)))));
\ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/main_covid19.m deleted 100644 → 0
%SCRIPT FOR RUNNING THE CRA TOOLBOX WITHOUT USING THE GUI
%The model is written as a Matlab function.
%%%% THE FIRST THING TO DO FOR RUNNING THE SCRIPT IS TO ADD THE FOLDER
%%%% CLASSES AND THE FOLDER Examples TO THE PATH %%%%
%Parameters to be set by the user are:
% 1-model_name: name of the model in .xml format. A struct is created to
% store the following characteristics of an ODE model: nominal values and
% names of the parameters, initial conditions and names of variables, input
% values.
% 2-stop_time: final time point for the model simulation
% 3-step_size: time interval for the vector of time points
% 4- ode_solver: type of ode_solver for simulating the model (possible
% choices are: ode45, ode15s, ode23t, sundials, stochastic, explicit tau
% and implicit tau
% 5- Nr: number of independent realizations to perform
% 6- LBpi: lower boundary of the Latin Hypercube Sampling for perturbation
% of the model parameter space
% 7- UBpi: upper boundary of the Latin Hypercube Sampling
% 8- Ns: number of samples of the Latin Hypercube (parameters n of the
% lhsdesign function)
% 9- variable_name: variable of the model to set as reference node to be
% measured
% 10- current_func: is the type of evaluation function. Currently, it is
% possible to choose among
% three evaluation functions: area under the curve, maximum value and time
% of maximum for the time behavior of the selected variables. The user can
% also define his evaluation function in a .m file by extending the
% abstract class EvaluationFunction
% 11- tail_size: number of samples to include in the upper and lower tail
% when computing the probability density function of the evaluation
% function
% 12- current_tm: is the method for computing the tails of the pdf of the
% evaluation function. Right now, it is possible to choose between two
% methods: sorted() which sorts the values of the evaluation function and
% selects the first and last samples according to tail_size; tmp_sum()
% computes the tails by selecting the upper and lower quartile. When using
% tmp_sum(), the parameter step_size needs also to be specified. Step_size
% is the step for computing the lower and upper quartile of the pdf in an
% iterative way, i.e. when the upper and lower tails do not have the number of
% samples specified by the user, the threshold is increased of a quantity
% equal to step_size and the calculation of the tails is repeated.
% The user can also define his own method for the tails computation by
% extending the abstract class TailMethod().
% 13- folder: name of the folder to create where the results are saved
tic;
model_name='ode_covid19_v3_OK';
stop_time=300;
step_size=1;
ode_solver='ode15s';
%time axis for model simulation
time_axis=[0:step_size:stop_time]';
%parameters and initial conditions of the model
init_italia
load('Examples/ITALIA/italia_moda_nr4.mat');
nominal_parameters=[italia_moda_nr4(1,1:19) 1 1 italia_moda_nr4(1,20:21)];
%nominal_parameters=italia_moda_nr4;
nominal_parameters_name={'be', 'b0', 'b1', 'b2','b3', 'FracSevere', 'FracCritical', 'FracAsym', 'IncubPeriod', 'DurMildInf', 'DurAsym', 'DurHosp', 'TimeICUDeath', 'ProbDeath', 'PresymPeriod','s01','s02','s03','s04','s05','s06','s1','s2'};
s00=[nominal_parameters(16) nominal_parameters(17) nominal_parameters(18) nominal_parameters(19) nominal_parameters(20) nominal_parameters(21)];
s11=nominal_parameters(22);
s22=nominal_parameters(23);
derived_parameters=1;
derived_parameters_name={'R0'};
num_observables=9;
observables_name={'S','E0','E1','I0','I1','I2','I3','R','D'};
model=struct('name',model_name,'odesolver',ode_solver,'time',time_axis,'stop',stop_time,'step',step_size,'nominal_parameters',nominal_parameters,'nominal_parameters_name',{nominal_parameters_name},'derived_parameters',derived_parameters,'derived_parameters_name',{derived_parameters_name},'num_observables',num_observables,'observables_name',{observables_name},'initial_conditions',x0,'Tlock',Tlock,'region_name',region_name);
Nr=10;
%define specific boundaries for each parameter
%LBpi=[0.5 0.25 0.25 0.25 0.25 0.125 0.125 0.25 1 0.4311 0.6988 0.1504 0.9213 0.5 0.2995];
%UBpi=[2 4 4 4 4 32 16 10 2.81 2.8741 6.9881 1.5042 9.2128 2 1.1980];
%LBpi=[0.125 0.25 0.25 0.25 0.25 0.125 0.125 0.25 1 0.5 0.6988 0.5556 0.4 0.25 0.25];
%UBpi=[1.02 4 4 4 4 1.1 1.1 1.5 3 2 1.5 2.7778 4 1.5 1.02];
%LBpi=[0.125 0.25 0.25 0.25 0.25 0.125 0.125 0.25 0.5 0.5 0.6988 0.5556 0.4 0.25 0.25];
%UBpi=[1.02 4 4 4 4 1.1 1.1 1.5 2 2 1.5 2.7778 4 1.5 1.02];
load('Examples/ITALIA/intervalli_figura_bande_90perc_CON_PRCTILE_ITALIA_4Nr.mat')
%LBpi=parameters{1,1}(1,:)./nominal_parameters;
%UBpi=parameters{1,1}(2,:)./nominal_parameters;
LBpi=[parameters{1,1}(1,1:19) 0.1 0.1 parameters{1,1}(1,20:21)]./nominal_parameters;
UBpi=[parameters{1,1}(2,1:19) 0.9 0.9 parameters{1,1}(2,20:21)]./nominal_parameters;
%LBpi=[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.6396 0.7018 0.5403 0.7539 0.6150 0.6177 0.5640 0.6903 0.6754 0.5958 0.5023 0.4850 0.1 0.1 0.6887 0.5429];
%UBpi=[2 2 2 2 2 1.5 1.38 1.2792 1.2281 1.6210 1.2565 1.5376 3.0885 1.0152 1.2426 1.3509 1.3903 1.5068 1.4549 0.9 0.9 2.4105 2.1716];
%LBpi=[0.99 0.5 0.5 0.5 0.5 0.5 0.5 0.6396 0.7018 0.5403 0.7539 0.6150 0.6177 0.5640 0.99 0.6754 0.5958 0.5023 0.4850 0.6887 0.5429];
%UBpi=[1.05 2 2 2 2 1.5 1.38 1.2792 1.2281 1.6210 1.2565 1.5376 3.0885 1.0152 1.05 1.3509 1.3903 1.5068 1.4549 2.4105 2.1716];
Ns=10000;
variable_name='I2';
current_func=Area(); %current evaluation function
tail_size=1000; %number of samples for the lower and upper tail
step_size=0.001; %this parameter needs to be defined only when current_tm=tmp_sum(step_size)
current_tm=sorted();
folder='covid19_area_10Nr_1000sample_time300_ITALIA_4Nr_90perc_predizione';
%model simulation for each sample of the Latin Hypercube
disp('Starting model simulation with perturbed parameters');
[AllResults,AllPerturbations,AllDerivedParam]=start_simulation_v2_apr2020(model,Nr,LBpi,UBpi,Ns);
disp('All done! Model simulation completed!');
%computation of the MIRI for the chosen evaluation function and model
%variable
disp('Starting computation of the MIRI for each parameter...');
try
compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,AllDerivedParam,folder)
catch ME
return
end
%plot and save probability density function of the evaluation function
disp('Plot of the probability density function of the chosen evaluation function');
plotpdf_evalfunc(folder,variable_name);
%plot and save conditional probability density functions of the parameters
disp('Plot of the parameter probability density functions');
plotpdf_param(folder,variable_name,model,Nr);
toc;
\ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/start_simulation_v2_apr2020.m deleted 100644 → 0
%Function for simulating the ODE model using the perturbed parameter
%vector. It takes in input:
%1. model which is the struct containing all the characheristics of the
%model
%2. Nr is the number of independent realizations to perform;
%3. LBpi and UBpi are the lower and upper boundaries of the Latin
%Hypercube. They are arrays since each parameter may have its own range of
%variation;
%4. Ns is the number of samples to generate for creating the Latin
%Hypercube;
%It returns in output the array AllResults containing all the output
%variables simulated using the parameter vectors in AllPerturbations
function [AllResults, AllPerturbations,AllDerivedParam]=start_simulation_v2_apr2020(model,Nr,LBpi,UBpi,Ns)
nominal_parameters=model.nominal_parameters;
num_parameters=length(nominal_parameters);
AllResults=cell(Nr,model.num_observables);
AllPerturbations=cell(Nr,1);
AllDerivedParam=cell(Nr,1);
h = waitbar(0,'Please wait...');
parpool(4)
for k=1:Nr
waitbar(k/Nr,h,['Realization number ',num2str(k)]);
%hypercube=LatinHypercube(LBpi,UBpi,Ns);
%perturbation=hypercube.generate_sample(num_parameters);
perturbation=zeros(Ns,num_parameters);
for b=1:num_parameters
hypercube=LatinHypercube(LBpi(b),UBpi(b),Ns);
perturbation_i=hypercube.generate_sample(1);
perturbation(:,b)=perturbation_i;
end
perturbed_p=repmat(nominal_parameters,length(perturbation),1).*perturbation;
size(perturbed_p)
AllPerturbations{k,1}=perturbation;
disp(strcat('Generated a Latin Hypercube of size',{' '},sprintf('%.0f',size(perturbation,1)),'x',sprintf('%.0f',size(perturbation,2))));
%psim_all=cell(Ns,model.num_observables);
psim_all=cell(Ns,1);
R0_all=zeros(Ns,1);
lastwarn('')
parfor l=1:Ns %parfor
pi=perturbed_p(l,:);
%modelsim=@(t,x)feval(model.name,t,x,pi,model.input,model.total_proteins);
modelsim=@(t,x)feval(model.name,t,x,pi,model.Tlock,model.region_name);
%modelsim=@(t,x)feval(model.name,t,x,pi,model.Tlock1,model.s01,model.s11);
[tsimi,psimi]=feval(model.odesolver,modelsim,model.time,model.initial_conditions);
R0i=compute_R0(pi);
%plot(tsimi, psimi(:,6))
%hold on;
% for t=1:model.num_observables
% psim_all{i,t}=psimi(:,t);
%end
msgstr=lastwarn;
if isequal(msgstr,'')
psim_all{l,1}=psimi;
R0_all(l,1)=R0i;
else
psim_all{l,1}=NaN(length(model.time),model.num_observables);
R0_all(l,1)=NaN;
end
lastwarn('')
end
AllDerivedParam{k,1}=R0_all;
for j=1:model.num_observables
for w=1:Ns
AllResults{k,j}{w,1}=psim_all{w,1}(:,j);
end
end
%AllResults{k,1:end}=psim_all;
disp(strcat('Completed one realization. The array of results has size',{' '},sprintf('%.0f',size(psim_all,1)),'x',sprintf('%.0f',size(psim_all,2))));
end
delete(gcp) %close parallel pool
close (h) %close the waitbar
end
\ No newline at end of file
Styling with Markdown is supported
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!