Commit 06a93842 by Chiara Antonini

update: different intervals of variation for the parameters

1 parent c4ede9d3
Showing 30 changed files with 316 additions and 1 deletions
%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);
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];
%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_v2';
stop_time=80;
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
N = 882000;
InitInf=1;
E00 = InitInf;
S0 = N-E00;
x0 = [S0 E00 0 0 0 0 0 0 0];
load('Examples/moda_parametri_inizio14step_LOG_100000_UMBRIA_I2I3DE.mat');
nominal_input_parameters=moda_parametri_inizio14step_LOG;
be = nominal_input_parameters(1);
b0=nominal_input_parameters(2);
b1=nominal_input_parameters(3);
b2=nominal_input_parameters(4);
b3=nominal_input_parameters(5);
be=be/(N);
b0=b0/(N);
b1=b1/(N);
b2=b2/(N);
b3=b3/(N);
FracSevere = nominal_input_parameters(6);
FracCritical = nominal_input_parameters(7);
FracMild = 1-FracSevere-FracCritical;
FracAsym = nominal_input_parameters(8);
IncubPeriod = nominal_input_parameters(9);
DurMildInf = nominal_input_parameters(10);
DurAsym = nominal_input_parameters(11);
DurHosp = nominal_input_parameters(12);
TimeICUDeath = nominal_input_parameters(13);
ProbDeath=nominal_input_parameters(14);
PresymPeriod=nominal_input_parameters(15);
s01=nominal_input_parameters(16);
s11=nominal_input_parameters(17);
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;
nominal_parameters=[be b0 b1 b2 b3 a0 a1 f g0 g1 p1 g2 p2 g3 u];
nominal_parameters=ones(1,15);
parameters_name={'be', 'b0', 'b1', 'b2','b3', 'a0', 'a1', 'f', 'g0', 'g1', 'p1', 'g2', 'p2', 'g3','u'};
Tlock1=23;
num_observables=6;
observables_name={'E1','I0','I1','I2','I3','D'};
model=struct('name',model_name,'odesolver',ode_solver,'time',time_axis,'stop',stop_time,'step',step_size,'nominal_parameters',nominal_parameters,'parameters_name',{parameters_name},'num_observables',num_observables,'observables_name',{observables_name},'initial_conditions',x0,'Tlock1',Tlock1,'s01',s01,'s11',s11);
Nr=3;
%define specific boundaries for each parameter
LBpi=[0.001 0.001 0.001 0.001 0.001 1e-5 1e-5 0.01 2 1 1 1 1 1 1];
UBpi=[3 3 3 3 3 0.5 0.4 0.7 8 20 20 10 30 60 4];
Ns=1000;
variable_name='I0';
current_func=Maximum(); %current evaluation function
tail_size=10; %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_I0_results';
%model simulation for each sample of the Latin Hypercube
disp('Starting model simulation with perturbed parameters');
[AllResults,AllPerturbations]=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,folder)
catch ME
%break
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 \ No newline at end of file
...@@ -34,7 +34,8 @@ for k=1:Nr ...@@ -34,7 +34,8 @@ for k=1:Nr
psim_all=cell(Ns,1); psim_all=cell(Ns,1);
parfor i=1:Ns parfor i=1:Ns
pi=perturbed_p(i,:); pi=perturbed_p(i,:);
modelsim=@(t,x)feval(model.name,t,x,pi,model.input,model.total_proteins); %modelsim=@(t,x)feval(model.name,t,x,pi,model.input,model.total_proteins);
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); [tsimi,psimi]=feval(model.odesolver,modelsim,model.time,model.initial_conditions);
% for t=1:model.num_observables % for t=1:model.num_observables
% psim_all{i,t}=psimi(:,t); % psim_all{i,t}=psimi(:,t);
......
%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]=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);
h = waitbar(0,'Please wait...');
for k=1:Nr
waitbar(k/Nr,h,['Realization number ',num2str(k)]);
perturbation=zeros(Ns,num_parameters);
for i=1:length(LBpi)
hypercube=LatinHypercube(LBpi(i),UBpi(i),Ns);
perturbation_i=hypercube.generate_sample(1);
perturbation(:,i)=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);
parfor i=1:Ns
pi=perturbed_p(i,:);
%modelsim=@(t,x)feval(model.name,t,x,pi,model.input,model.total_proteins);
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);
% for t=1:model.num_observables
% psim_all{i,t}=psimi(:,t);
%end
msgstr=lastwarn;
if isequal(msgstr,'')
psim_all{i,1}=psimi;
else
psim_all{i,1}=NaN(length(model.time),model.num_observables);
end
lastwarn('')
end
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 \ 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!