Upload new file

CRA_Updated_Version/CRA_Toolbox_bash/main_EGFRIGF1Rmodel.m

+%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='EGFR_IGF1R_model';
+
+stop_time=60;
+step_size=0.1;
+ode_solver='ode15s';
+
+%time axis for model simulation
+time_axis=[0:step_size:stop_time]';
+
+%parameters and initial conditions of the model
+run parameters_EGFR_IGF1R_model
+nominal_parameters=[gamma_EGFR, gamma_IGF1R,kd_P90Rsk,kd_PIK3_active,kd_Akt,k_SOS_EGFR,KM_SOS_EGFR,k_SOS_IGF1R,...
+KM_SOS_IGF1R,k_DSOS_p90Rsk,KM_DSOS_p90Rsk,k_Ras_SOS,KM_Ras_SOS,k_Ras_RasGab,KM_Ras_RasGab,k_Raf_Ras,KM_Raf_Ras,...
+k_Raf_RafPP,KM_Raf_RafPP,k_Raf_Akt,KM_Raf_Akt,k_Mek_Raf,KM_Mek_Raf,k_MEK_PP2A,KM_MEK_PP2A,k_ERK_MEK,KM_ERK_MEK,...
+k_Erk_PP2A,KM_Erk_PP2A,k_p90Rsk_Erk,KM_p90Rsk_Erk,k_PIK3_Ras,KM_PIK3_Ras,k_PIK3_IGF1R,KM_PIK3_IGF1R,k_PIK3_EGFR,...
+KM_PIK3_EGFR,k_Akt_PIK3,KM_Akt_PIK3];
+
+parameters_name={'gamma_EGFR', 'gamma_IGF1R','kd_P90Rsk','kd_PIK3_active','kd_Akt','k_SOS_EGFR','KM_SOS_EGFR','k_SOS_IGF1R',...
+'KM_SOS_IGF1R','k_DSOS_p90Rsk','KM_DSOS_p90Rsk','k_Ras_SOS','KM_Ras_SOS','k_Ras_RasGab','KM_Ras_RasGab','k_Raf_Ras','KM_Raf_Ras',...
+'k_Raf_RafPP','KM_Raf_RafPP','k_Raf_Akt','KM_Raf_Akt','k_Mek_Raf','KM_Mek_Raf','k_MEK_PP2A','KM_MEK_PP2A','k_ERK_MEK','KM_ERK_MEK',...
+'k_Erk_PP2A','KM_Erk_PP2A','k_p90Rsk_Erk','KM_p90Rsk_Erk','k_PIK3_Ras','KM_PIK3_Ras','k_PIK3_IGF1R','KM_PIK3_IGF1R','k_PIK3_EGFR',...
+'KM_PIK3_EGFR','k_Akt_PIK3','KM_Akt_PIK3'};
+
+x0=[EGFR_active_0,IGF1R_active_0,SOS_0,Ras_active_0,Raf_active_0,Mek_active_0,Erk_active_0,p90Rsk_active_0,...
+PIK3_active_0,Akt_active_0];
+
+num_observables=10;
+observables_name={'EGFR_active','IGF1R_active','SOS','Ras_active','Raf_active','Mek_active','Erk_active',...
+'p90Rsk_active','PIK3_active','Akt_active'};
+
+u=[RafPP PP2A RasGapActive];
+xT=[DSOS_tot Ras_tot Raf_tot Mek_tot Erk_tot p90Rsk_tot PIK3_tot Akt_tot];
+
+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,'input',u,'total_proteins',xT);
+
+Nr=100;
+
+LBpi=0.1;
+UBpi=10;
+Ns=10000;
+
+variable_name='Erk_active';
+current_func=Area(); %current evaluation function
+tail_size=1000;  %number of samples for the lower and upper tail
+step_size=0.01;  %this parameter needs to be defined only when current_tm=tmp_sum(step_size)
+current_tm=tmp_sum(step_size);
+
+folder='modelEGFR_test1';
+
+%model simulation for each sample of the Latin Hypercube
+disp('Starting model simulation with perturbed parameters');
+[AllResults,AllPerturbations]=start_simulation(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    
+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