update code

CRA_Updated_Version/CRA_Toolbox_bash/Examples/parameters_EGFR_IGF1R_model.m

+% 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_MIRI.m

@@ -29,21 +29,21 @@
 
 
 
-function compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,AllDerivedParam,folder)
+function compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,folder)
 
 nominal_parameters=model.nominal_parameters;
-%num_parameters=length(nominal_parameters);
-nominal_parameters_name=model.nominal_parameters_name;
-
-if(isfield(model,'derived_parameters'))
-   derived_parameters=model.derived_parameters;
-   derived_parameters_name=model.derived_parameters_name;
-   num_parameters=length(nominal_parameters)+length(derived_parameters);
-   parameters_name=[nominal_parameters_name derived_parameters_name];
-else
-   num_parameters=length(nominal_parameters);
-   parameters_name=[nominal_parameters_name];
-end
+num_parameters=length(nominal_parameters);
+parameters_name=model.parameters_name;
+
+% if(isfield(model,'derived_parameters'))
+%    derived_parameters=model.derived_parameters;
+%    derived_parameters_name=model.derived_parameters_name;
+%    num_parameters=length(nominal_parameters)+length(derived_parameters);
+%    parameters_name=[nominal_parameters_name derived_parameters_name];
+% else
+%    num_parameters=length(nominal_parameters);
+%    parameters_name=[nominal_parameters_name];
+% end
 
 obs_name=model.observables_name;
 
@@ -124,7 +124,7 @@ try
         Results.currentTailMethod=current_tm;
         
         [XiMax,XiMin]=Results.compute_tail(AllPerturbations{k,1},tail_size);
-        [XiMaxD,XiMinD]=Results.compute_tail(AllDerivedParam{k,1},tail_size);
+        %[XiMaxD,XiMinD]=Results.compute_tail(AllDerivedParam{k,1},tail_size);
         %[XiMax,XiMin]=handles.Results.compute_tail(handles1.perturbation,handles.tail_size);
         
         NcloudMax=size(XiMax,1);
@@ -138,17 +138,17 @@ try
         %[all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=obj_MIRI.estimate_parampdfs(handles1.nominal_p, handles1.perturbation, Ncloud, handles.XiMax, handles.XiMin);
         [MIRI_param,CloudCondMin, CloudCondMax]=obj_MIRI.MIRI_computation(all_ks_max,all_ks_min,all_xbin_max,all_xbin_min);
         
-        if(isfield(model,'derived_parameters'))
-           obj_MIRI_derived=MIRI();
-           [all_xbin_max_d,all_xbin_min_d,all_ks_max_d,all_ks_min_d]=obj_MIRI_derived.estimate_parampdfs(derived_parameters, AllDerivedParam{k,1}, XiMaxD, XiMinD);
-           %[all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=obj_MIRI.estimate_parampdfs(handles1.nominal_p, handles1.perturbation, Ncloud, handles.XiMax, handles.XiMin);
-           [MIRI_param_d,CloudCondMin_d, CloudCondMax_d]=obj_MIRI_derived.MIRI_computation(all_ks_max_d,all_ks_min_d,all_xbin_max_d,all_xbin_min_d);
-        
-        
-           MIRI_param=[MIRI_param MIRI_param_d];
-           CloudCondMax=[CloudCondMax CloudCondMax_d];
-           CloudCondMin=[CloudCondMin CloudCondMin_d];
-        end
+%         if(isfield(model,'derived_parameters'))
+%            obj_MIRI_derived=MIRI();
+%            [all_xbin_max_d,all_xbin_min_d,all_ks_max_d,all_ks_min_d]=obj_MIRI_derived.estimate_parampdfs(derived_parameters, AllDerivedParam{k,1}, XiMaxD, XiMinD);
+%            %[all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=obj_MIRI.estimate_parampdfs(handles1.nominal_p, handles1.perturbation, Ncloud, handles.XiMax, handles.XiMin);
+%            [MIRI_param_d,CloudCondMin_d, CloudCondMax_d]=obj_MIRI_derived.MIRI_computation(all_ks_max_d,all_ks_min_d,all_xbin_max_d,all_xbin_min_d);
+%         
+%         
+%            MIRI_param=[MIRI_param MIRI_param_d];
+%            CloudCondMax=[CloudCondMax CloudCondMax_d];
+%            CloudCondMin=[CloudCondMin CloudCondMin_d];
+%         end
         
         MIRIT(k,:)=MIRI_param;
 
@@ -164,14 +164,14 @@ try
             ks_param{2,k}(ip,:)=all_ks_min(ip,:);
         end   
         
-        if(isfield(model,'derived_parameters'))
-            for ip=1:length(derived_parameters)
-               xbin_param_d{1,k}(ip,:)=all_xbin_max_d(ip,:);
-               xbin_param_d{2,k}(ip,:)=all_xbin_min_d(ip,:);
-               ks_param_d{1,k}(ip,:)=all_ks_max_d(ip,:);
-               ks_param_d{2,k}(ip,:)=all_ks_min_d(ip,:);
-            end
-        end
+%         if(isfield(model,'derived_parameters'))
+%             for ip=1:length(derived_parameters)
+%                xbin_param_d{1,k}(ip,:)=all_xbin_max_d(ip,:);
+%                xbin_param_d{2,k}(ip,:)=all_xbin_min_d(ip,:);
+%                ks_param_d{1,k}(ip,:)=all_ks_max_d(ip,:);
+%                ks_param_d{2,k}(ip,:)=all_ks_min_d(ip,:);
+%             end
+%         end
         
    end
    
@@ -210,9 +210,9 @@ try
    disp('saving probability density functions of all parameters')
    save(fullfile(folder,variable_name,'pdf_param.mat'),'xbin_param','ks_param');
    
-   if(isfield(model,'derived_parameters'))
-      save(fullfile(folder,variable_name,'pdf_param_derived.mat'),'xbin_param_d','ks_param_d');
-   end
+%    if(isfield(model,'derived_parameters'))
+%       save(fullfile(folder,variable_name,'pdf_param_derived.mat'),'xbin_param_d','ks_param_d');
+%    end
       
 catch ME
     

CRA_Updated_Version/CRA_Toolbox_bash/main_EGFRIGF1Rmodel.m

@@ -83,15 +83,15 @@ 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;
+Nr=3;
 
 LBpi=0.1;
 UBpi=10;
-Ns=10000;
+Ns=1000;
 
 variable_name='Erk_active';
 current_func=Area(); %current evaluation function
-tail_size=1000;  %number of samples for the lower and upper tail
+tail_size=10;  %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);
 
@@ -108,7 +108,7 @@ 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');

CRA_Updated_Version/CRA_Toolbox_bash/main_PulseGeneratorNetwork.m

@@ -64,21 +64,21 @@ parameters_name={'k1','K1','lambda2','k12','K2','lambda'};
 x0=[0 0];
 u=470;
 num_observables=2;
-observables_name={'R1','Y'};
+observables_name={'R2','Y'};
 
-model=struct('name',model_name,'odesolver',ode_solver,'time',time_axis,'stop',stop_time,'step',step_size,'nominal_parameters',nominal_parameters,'nominal_parameters_name',{parameters_name},'num_observables',num_observables,'observables_name',{observables_name},'initial_conditions',x0,'input',u,'total_proteins',0);
+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',0);
 
-Nr=2;
+Nr=10;
 
 folder='2014sortedprovaparfor_PulseGeneratorNetwork_results';
 
 LBpi=0.1;
 UBpi=10;
-Ns=3000;
+Ns=1000;
 
 variable_name='Y';
-current_func=Area(); %current evaluation function
-tail_size=1000;  %number of samples for the lower and upper tail
+current_func=TimeOfMaximum(); %current evaluation function
+tail_size=100;  %number of samples for the lower and upper tail
 step_size=0.1;  %this parameter needs to be defined only when current_tm=tmp_sum(step_size)
 %current_tm=tmp_sum(0.1);
 current_tm=sorted();
@@ -92,8 +92,7 @@ disp('All done! Model simulation completed!');
 %variable
 disp('Starting computation of the MIRI for each parameter...');
 try
-  AllDerivedParam=cell(Nr,1);
-  compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,AllDerivedParam,folder)
+  compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,folder)
 catch ME
     %break  
     return

CRA_Updated_Version/CRA_Toolbox_bash/main_modelPeng2016.m

@@ -66,11 +66,11 @@ observables_name={'A','X1','X2','D_m','C2','R2','I2','D_C','D_R','C1','R1','I1'}
 
 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',0,'total_proteins',0);
 
-Nr=100;
+Nr=10;
 
 LBpi=0.1;
 UBpi=10;
-Ns=10000;
+Ns=100;
 
 variable_name='C1';
 current_func=Area(); %current evaluation function
@@ -91,7 +91,7 @@ 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');

CRA_Updated_Version/CRA_Toolbox_bash/plotpdf_param.m

@@ -7,10 +7,7 @@
 function plotpdf_param(folder,variable_name,model,Nr)
 
 nominal_parameters=model.nominal_parameters;
-nominal_parameters_name=model.nominal_parameters_name;
-
-derived_parameters=model.derived_parameters;
-derived_parameters_name=model.derived_parameters_name;
+nominal_parameters_name=model.parameters_name;
 
     try
         load(fullfile(folder,variable_name,'pdf_param.mat'));
@@ -26,20 +23,6 @@ derived_parameters_name=model.derived_parameters_name;
             legend('Conditional upper density','Conditional lower density');
             savefig(fullfile(folder,variable_name,nominal_parameters_name{1,i}));
         end
-            
-         load(fullfile(folder,variable_name,'pdf_param_derived.mat'));
-         for i=1:length(derived_parameters)
-            figure
-            for k=1:Nr
-                plot(xbin_param_d{1,k}(i,:),ks_param_d{1,k}(i,:),'b',xbin_param_d{2,k}(i,:),ks_param_d{2,k}(i,:),'r');
-                hold on
-            end
-            xlabel(derived_parameters_name{1,i},'Interpreter','none');
-            ylabel('pdf');
-            legend('Conditional upper density','Conditional lower density');
-            savefig(fullfile(folder,variable_name,derived_parameters_name{1,i}));
-        
-        end
     catch ME
         errordlg(ME.message)
     end