matlab code

Area.m

+%Class computing the area of the temporal behavior of a model variable. It extends
+%the abstract class EvaluationFunction. The method compute_evalfunc takes in input
+%an object TimeBehavior and returns in output the value of the area. 
+
+classdef Area < EvaluationFunction    
+    
+    properties
+    end
+   
+    methods
+               
+        function EvalFuncValue=compute_ef(obj,obj_TimeBehavior)
+            EvalFuncValue=sum(obj_TimeBehavior.values,2);            
+        end
+        
+    end
+    
+end
\ No newline at end of file

EvaluationFunction.m

+%Abstract class representing a generic evaluation function to compute
+
+classdef (Abstract) EvaluationFunction
+   
+    methods (Abstract)
+        values=compute_ef(obj,obj_TimeBehavior);        
+    end
+    
+end
\ No newline at end of file

LatinHypercube.m

+%Class representing a latin hypercube generated through the function
+%lhsdesign. The object LatinHypercube has the following properties:
+% -Nsample is the number of samples of the hypercube
+% -LBpi and UBpi are, respectively, the lower and upper boundaries of the
+% interval in which the hypercube is centered
+
+classdef LatinHypercube
+    
+    properties
+       LBpi
+       UBpi
+       Nsample
+    end
+    
+    methods
+       
+        function obj=LatinHypercube(lb,ub,ns)
+            obj.LBpi=lb;
+            obj.UBpi=ub;
+            obj.Nsample=ns;            
+        end
+        
+        function perturbation=generate_sample(obj,len_p)
+            perturbation=obj.LBpi+(obj.UBpi-obj.LBpi).*lhsdesign(obj.Nsample,len_p); 
+        end
+        
+        
+    end
+    
+end
\ No newline at end of file

MIRI.m

+%Concrete class for computing MIRI of model parameters. The class has two
+%methods:
+% - estimate_parampdfs that estimates the conditional upper and lower pdf
+% of a parameter. It returns in output the array all_xbin_max and
+% all_xbin_min containing the points of the ksdensity in the upper and
+% lower case while all_ks_max and all_ks_min contains the corresponding
+% probability density estimate;
+% - MIRI_computation takes in input the two estimated pdfs of all
+% parameters and computes MIRI for each of them. MIRI are returned in the
+% array MIRI_param while CloudCondMax and CloudCondMin contain the mode of
+% the parameter vector in the upper and lower case respectively.
+
+classdef MIRI
+   
+    properties
+    end
+    
+    methods
+   
+    function   [all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=estimate_parampdfs(obj,p, Perturbation, XiMax, XiMin)
+   
+      all_xbin_max=[];
+      all_xbin_min=[];
+      all_ks_max=[];
+      all_ks_min=[];
+    
+      pdf_obj=pdfEstimator();      
+      
+      for ip=1:length(p) 
+         
+         %step=p(ip)*(Ncloud^(-1./3)); 
+         step=(max(p(ip).*Perturbation(:,ip)) - min(p(ip).*Perturbation(:,ip)))/size(Perturbation(:,ip),1);
+          
+         BinEdgesMax=min(p(ip).*Perturbation(:,ip)):step:max(p(ip).*Perturbation(:,ip));
+         BinEdgesMin=min(p(ip).*Perturbation(:,ip)):step:max(p(ip).*Perturbation(:,ip)); 
+         
+         [ks_max, xbin_max]=pdf_obj.evaluate_pdf(p(ip).*XiMax(:,ip)',BinEdgesMax);
+         [ks_min, xbin_min]=pdf_obj.evaluate_pdf(p(ip).*XiMin(:,ip)',BinEdgesMin);
+
+         all_xbin_max=[all_xbin_max;xbin_max];
+         all_xbin_min=[all_xbin_min;xbin_min];
+         all_ks_max=[all_ks_max;ks_max];
+         all_ks_min=[all_ks_min;ks_min];
+         
+      end
+    
+        
+    end    
+        
+    function [MIRI_param, CloudCondMin, CloudCondMax]=MIRI_computation(obj, all_ks_max,all_ks_min,all_xbin_max,all_xbin_min)
+
+      len_p=size(all_ks_max,1);
+      MIRI_param=zeros(1,len_p);
+      CloudCondMin=zeros(1,len_p);
+      CloudCondMax=zeros(1,len_p);
+      
+      for ip=1:len_p
+            
+            ks_max=all_ks_max(ip,:);
+            ks_min=all_ks_min(ip,:);
+            
+            xbin_max=all_xbin_max(ip,:);
+            xbin_min=all_xbin_min(ip,:);
+          
+           [ks_max_m, ks_max_I]=max(ks_max);
+           [ks_min_m, ks_min_I]=max(ks_min);
+         
+         
+           CloudCondMin(1,ip)=xbin_min(ks_min_I);
+           CloudCondMax(1,ip)=xbin_max(ks_max_I);
+           s_Xi_ks=sum(abs(ks_max./sum(ks_max)-ks_min./sum(ks_min))); 
+           MIRI_param(1,ip)=s_Xi_ks; 
+
+      end
+
+    end
+   
+           
+ end
+    
+end
\ No newline at end of file

Maximum.m

+%Class computing the maximum value of a model variable curve. It extends
+%the abstract class EvaluationFunction. The method compute_evalfunc takes in input
+%an object TimeBehavior and returns in output the maximum. 
+
+
+classdef Maximum < EvaluationFunction    
+    
+    properties
+    end
+   
+    methods
+               
+        function EvalFuncValue=compute_ef(obj,obj_TimeBehavior)
+            EvalFuncValue=max(obj_TimeBehavior.values,[],2);            
+        end
+        
+    end
+    
+end
\ No newline at end of file

ReadMe.txt

+#################################
+# Files description
+#################################
+
+This is the main code folder. It contains the following files:
+
+1. gui_CRA.m is the function for starting the main GUI. It contains all the functions of the elements in the main GUI.
+   It recalls the following functions:
+   1.1 start_simulation.m: function that executes the CRA algorithm.
+2. gui_CRA.fig is the figure of the graphical interface defined in gui_CRA.m. It allows us to specify the following parameters for
+   performing the CRA:
+   - stop time of the simulation;
+   - type of ode solver and step of the time vector;
+   - number of independent realizations of the procedure to perform;
+   - lower and upper boundaries of the Latin Hypercube to generate;
+   - number of samples of the Latin Hypercube;
+   - variable of the model to set as reference node to measure.
+3. gui_MIRI.m is the function for starting the second GUI for computing MIRI. It contains all the callbacks of the elements in the GUI. 
+   It recalls the following functions:
+   3.1 compute_MIRI.m: function that computes MIRI of the model parameters with respect to the chosen protein
+   3.2 plotpdf_param.m: function for visualizing and saving the conditional probability density function (pdf) of all 
+   the model parameters.
+4. gui_MIRI.fig is the figure of the graphical interface defined in gui_MIRI.m. It allows us to specify the following parameters for 
+   computing MIRI of the output node selected in the previous GUI:
+   - type of evaluation function. It is possible to choose among three evaluation functions: area under the curve, maximum value and time 
+     of maximum of the time behavior of the chosen variable;
+   - number of samples N to include in the upper and lower tail of the evaluation function pdf;
+   - method for computing the tails of the evaluation function pdf. It is possible to choose between two methods: sort and tmp_sum. 
+     The sort method orders all the samples of the evaluation function and selects the N lowest and highest values.
+     The tmp_sum method calculates the lower and upper quartile of the pdf, computing the threshold values in an adaptive manner.
+     The initial lower threshold is set to 0 and it is repeatedly increased of the step size defined by the user until at least N samples 
+     are included in the tail. The same procedure is repeated for the upper tail, starting from a threshold equal to 1. 
+5. start_simulation.m is the function executed when the button START SIMULATION in the main GUI is pushed. It takes in input 
+   the handles of the main GUI in order to retrieve all the variables inserted by the user. For each independent realization, 
+   it generates the Latin Hypercube and simulates the model.
+6. compute_MIRI.m is the function executed when the button COMPUTE MIRI in the second GUI (gui_MIRI) is pushed. It takes in input 
+   the handles of both GUIs and it computes MIRI of the selected protein. It plots boxplot of MIRI values and the pdf of the
+   evaluation function of the chosen protein. It saves MIRI values and the pdfs of both the evaluation function and parameters of the model.
+7. plotpdf_evalfunc.m is the function that plots and saves the estimated pdf of the evaluation function for the chosen protein.
+8. plotpdf_param.m is the function that plots and saves the conditional pdfs of model parameters.
+9. EvaluationFunction.m is the abstract class for a generic evaluation function.
+10. Area.m is a concrete class that extends the abstract class EvaluationFunction. It computes the area under the curve of a 
+   model variable.
+11. Maximum.m is a concrete class that extends the abstract class EvaluationFunction. It computes the maximum value of a variable curve
+12. TimeOfMaximum.m is a concrete class that extends the abstract class EvaluationFunction. It computes the maximum point value of 
+    a variable curve.
+13. TailMethod.m is the abstract class for defining a generic method that computes upper and lower tails of a parameter pdf. 
+    The upper tail is conditioned to high values of the evaluation function while the lower tail is conditioned to low values 
+    of the evaluation function. 
+14. sorted.m is a concrete class that extends the abstract class TailMethod.
+15. tmp_sum.m is a concrete class that extends the abstract class TailMethod.
+16. TimeBehavior.m is the class representing the time behavior of a model variable.
+17. LatinHypercube.m is the class representing the Latin Hypercube generated for perturbing the parameter space.
+18. pdfEstimator.m is the class representing an object for estimating a probability density function through the method ksdensity
+19. MIRI.m is the class for computing the Moment Independent Robustness Indicator (MIRI) of each model parameter   
+
+The folder Models contains some ODE models in SBML format:  
+- PulseGenerator.xml 
+- EGFR_IG1FRmodel.xml
+- model_Peng2016_reduced.xml
\ No newline at end of file

TailMethod.m

+%Abstract class for defining a method to compute the upper and lower tails of the conditional
+%pdfs of parameters
+
+classdef (Abstract) TailMethod
+   
+     methods (Abstract)
+        [XiMax,XiMin]=compute_tailspdf(obj,obj_TimeBehavior,perturbation,tail_size);        
+    end
+    
+end
\ No newline at end of file

TimeBehavior.m

+%Concrete class representing the time behavior of a model variable. It
+%stores the values of a variable curve, the values of the selected
+%evaluation function and the method chosen for tail computation
+
+classdef TimeBehavior < handle
+   
+    properties
+       time   %x-axis of time simulation
+       values  %y-axis of time simulation
+       currentEvalFunc
+       evalFuncValues    
+       currentTailMethod
+    end
+    
+    methods
+       
+        function obj=TimeBehavior(t,v)
+           obj.time=t;
+           obj.values=v;
+        end
+        
+%         function obj=set.currentEvalFunc(obj,func)
+%             obj.currentEvalFunc=func;
+%         end
+        
+        %pattern Strategy for computing values of the evaluation function
+        function EvalValues=computeEvalFunc(obj)
+           EvalValues=obj.currentEvalFunc.compute_ef(obj); 
+           obj.evalFuncValues=EvalValues;
+        end
+        
+        
+        %pattern Strategy 
+        function [XiMax, XiMin]=compute_tail(obj,perturbation,tail_size)
+           [XiMax, XiMin]=obj.currentTailMethod.compute_tailspdf(obj,perturbation,tail_size); 
+        end
+                       
+    end
+    
+end
+
+

TimeOfMaximum.m

+%Class computing the maximum point of a model variable curve. It extends
+%the abstract class EvaluationFunction. The method compute_evalfunc takes in input
+%an object TimeBehavior and returns in output the maximum point. 
+
+classdef TimeOfMaximum < EvaluationFunction    
+    
+    properties
+    end
+   
+    methods
+       
+        function EvalFuncValue=compute_ef(obj,obj_TimeBehavior)
+            [maximum, timeOfMaximum]=max(obj_TimeBehavior.values,[],2);
+            EvalFuncValue=obj_TimeBehavior.time(timeOfMaximum);
+        end
+        
+    end
+    
+end
\ No newline at end of file

compute_MIRI.m

+%Function executed when the button COMPUTE MIRI in the second GUI (gui_MIRI) is
+%pushed. It takes in input the handles of both GUIs and it computes MIRI of
+%the selected protein. It plots boxplot of MIRI values and the pdf of the
+%evaluation function of the chosen protein. It saves the following
+%variables:
+% 1) MIRIT is the array containing MIRI values of all realizations
+% 2) mean_CloudMax is the array containing the mean of values in
+% CloudCondMaxT. CloudCondMaxT is the array of the mode of parameter values
+% that give rise to the upper tail of the evaluation function
+% 3) mean_CloudMin is the array containing the mean of values in
+% CloudCondMinT. CloudCondMinT is the array of the mode of parameter values
+% that give rise to the lower tail of the evaluation function
+% 4) pdf_eval_func contains values of the evaluation function  in all
+% realizations
+% 5) pdf_param contains values of the conditional pdf of parameters in all
+% realizations
+
+
+function compute_MIRI(handles,hObject,handles1)
+
+MIRIT=zeros(handles1.Nr,handles1.num_param);
+CloudCondMaxT=[];
+CloudCondMinT=[];
+xbinT=[];
+ks_y_T=[];
+
+try
+    
+   xbin_param=cell(2,handles1.Nr);
+   ks_param=cell(2,handles1.Nr);
+  
+   h = waitbar(0,'Please wait...');
+   
+   for k=1:handles1.Nr 
+        waitbar(k/handles1.Nr,h,['Realization number ',num2str(k)]);
+        
+        Results_reshape=zeros(handles1.Nsample,length(handles1.time_results));
+        for j=1:handles1.Nsample
+            Results_reshape(j,:)=reshape(handles1.AllResults{k,handles1.num_protein}{j,1},[1 length(handles1.time_results)]);
+        end
+        
+        %creation of an object TimeBehavior where results of the chosen
+        %protein are stored
+        handles.Results=TimeBehavior(handles1.time_results,Results_reshape);
+        
+        %values of the evaluation function are computed and stored in the
+        %TimeBehavior object Results
+        handles.Results.currentEvalFunc=handles.current_func;
+        handles.Results.computeEvalFunc();  
+        guidata(hObject,handles);
+        
+        samples=handles.Results.evalFuncValues;
+        
+        % an evaluation function may have identical values when parameters
+        % do not affect its value. For example, if a protein has a
+        % descending temporal behavior, the maximum and time of maximum
+        % will always be the same, independently of parameter perturbation
+        if (length(unique(samples))==1)
+           ME=MException('Matlab:EvalFuncValuesNull','Array of evaluation function contains identical values. Change evaluation function or protein.');
+           throw(ME)
+           close(h)
+        end
+        
+        %creation of an object pdfEstimator for estimating the pdf of the
+        %evaluation function
+        pdf_obj=pdfEstimator();
+        BinEdges=[min(samples):(max(samples)-min(samples))/length(samples):max(samples)];
+        [ks_y,xbin]=pdf_obj.evaluate_pdf(samples,BinEdges);
+        xbinT=[xbinT;xbin];
+        ks_y_T=[ks_y_T;ks_y];
+        
+        %method for extracting the upper and lower tail of the evaluation
+        %function
+        handles.Results.currentTailMethod=handles.current_tm;
+        
+        [XiMax,XiMin]=handles.Results.compute_tail(handles1.AllPerturbations{k,1},handles.tail_size);
+        %[XiMax,XiMin]=handles.Results.compute_tail(handles1.perturbation,handles.tail_size);
+        handles.XiMax=XiMax;
+        handles.XiMin=XiMin;
+        guidata(hObject,handles);
+  
+        NcloudMax=size(handles.XiMax,1);
+        NcloudMin=size(handles.XiMin,1);
+        Ncloud=min([NcloudMax NcloudMin]);
+ 
+        %creation of an object MIRI for computing MIRI and conditional pdfs
+        %of parameters
+        obj_MIRI=MIRI();
+        [all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=obj_MIRI.estimate_parampdfs(handles1.nominal_p, handles1.AllPerturbations{k,1}, handles.XiMax, handles.XiMin);
+        %[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);
+        
+        MIRIT(k,:)=MIRI_param;
+        CloudCondMaxT=[CloudCondMaxT;CloudCondMax];
+        CloudCondMinT=[CloudCondMinT;CloudCondMin];
+        mean_CloudMax=mean(CloudCondMaxT);
+        mean_CloudMin=mean(CloudCondMinT);
+        
+        for ip=1:length(handles1.nominal_p)
+            xbin_param{1,k}(ip,:)=all_xbin_max(ip,:);
+            xbin_param{2,k}(ip,:)=all_xbin_min(ip,:);
+            ks_param{1,k}(ip,:)=all_ks_max(ip,:);
+            ks_param{2,k}(ip,:)=all_ks_min(ip,:);
+        end   
+   end
+   
+   close (h)
+   
+   %MIRI boxplot are displayed in gui_MIRI
+   boxplot(handles.axes1,MIRIT,'labels',handles1.param_name,'labelorientation','inline');
+   ylabel('MIRI');
+   labelSize = 13; %size of the label
+   set(findobj(handles.axes1,'Type','text'),'FontSize',labelSize); 
+     
+   
+   %saving results
+   mkdir(handles1.folder,handles1.chosen_protein);
+   save(fullfile(handles1.folder,handles1.chosen_protein,'MIRIT.mat'),'MIRIT');
+   
+   
+   save(fullfile(handles1.folder,handles1.chosen_protein,'meanCloudCondMax.mat'),'mean_CloudMax');
+   save(fullfile(handles1.folder,handles1.chosen_protein,'meanCloudCondMin.mat'),'mean_CloudMin');   
+   
+   save(fullfile(handles1.folder,handles1.chosen_protein,'pdf_eval_func.mat'),'xbinT','ks_y_T');
+   save(fullfile(handles1.folder,handles1.chosen_protein,'pdf_param.mat'),'xbin_param','ks_param');
+   
+   %the evaluation function pdf is displayed in gui_MIRI
+   plotpdf_evalfunc(handles,hObject,handles1);
+   
+catch ME
+    
+    if (strcmp(ME.identifier,'MATLAB:badsubscript'))
+        msg='The pdf is a Dirac delta function: change evaluation function or protein';
+        errordlg(msg);
+        close(h)
+    else
+        errordlg(ME.message);
+        close(h)
+    end
+end
+end
\ No newline at end of file

gui_MIRI.m

+function varargout = gui_MIRI(varargin)
+% GUI_MIRI MATLAB code for gui_MIRI.fig
+%      GUI_MIRI, by itself, creates a new GUI_MIRI or raises the existing
+%      singleton*.
+%
+%      H = GUI_MIRI returns the handle to a new GUI_MIRI or the handle to
+%      the existing singleton*.
+%
+%      GUI_MIRI('CALLBACK',hObject,eventData,handles,...) calls the local
+%      function named CALLBACK in GUI_MIRI.M with the given input arguments.
+%
+%      GUI_MIRI('Property','Value',...) creates a new GUI_MIRI or raises the
+%      existing singleton*.  Starting from the left, property value pairs are
+%      applied to the GUI before gui_MIRI_OpeningFcn gets called.  An
+%      unrecognized property name or invalid value makes property application
+%      stop.  All inputs are passed to gui_MIRI_OpeningFcn via varargin.
+%
+%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
+%      instance to run (singleton)".
+%
+% See also: GUIDE, GUIDATA, GUIHANDLES
+
+% Edit the above text to modify the response to help gui_MIRI
+
+% Last Modified by GUIDE v2.5 15-Jun-2018 12:37:24
+
+% Begin initialization code - DO NOT EDIT
+gui_Singleton = 1;
+gui_State = struct('gui_Name',       mfilename, ...
+                   'gui_Singleton',  gui_Singleton, ...
+                   'gui_OpeningFcn', @gui_MIRI_OpeningFcn, ...
+                   'gui_OutputFcn',  @gui_MIRI_OutputFcn, ...
+                   'gui_LayoutFcn',  [] , ...
+                   'gui_Callback',   []);
+if nargin && ischar(varargin{1})
+    gui_State.gui_Callback = str2func(varargin{1});
+end
+
+if nargout
+    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
+else
+    gui_mainfcn(gui_State, varargin{:});
+end
+% End initialization code - DO NOT EDIT
+
+
+% --- Executes just before gui_MIRI is made visible.
+function gui_MIRI_OpeningFcn(hObject, eventdata, handles, varargin)
+% This function has no output args, see OutputFcn.
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+% varargin   command line arguments to gui_MIRI (see VARARGIN)
+
+% Choose default command line output for gui_MIRI
+handles.output = hObject;
+
+% Update handles structure
+guidata(hObject, handles);
+
+% UIWAIT makes gui_MIRI wait for user response (see UIRESUME)
+% uiwait(handles.gui_MIRI);
+
+% --- Outputs from this function are returned to the command line.
+function varargout = gui_MIRI_OutputFcn(hObject, eventdata, handles) 
+% varargout  cell array for returning output args (see VARARGOUT);
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Get default command line output from handles structure
+varargout{1} = handles.output;
+
+
+%function for selecting the desired evaluation function
+function uipanel1_SelectionChangeFcn(hObject, eventdata, handles)
+% hObject    handle to the selected object in uipanel1
+% eventdata  structure with the following fields (see UIBUTTONGROUP)
+%	EventName: string 'SelectionChanged' (read only)
+%	OldValue: handle of the previously selected object or empty if none was selected
+%	NewValue: handle of the currently selected object
+% handles    structure with handles and user data (see GUIDATA)
+switch get(eventdata.NewValue,'tag') 
+    case 'radiobutton1'     
+        handles.current_func=Area();
+    case 'radiobutton2'        
+        handles.current_func=Maximum();
+    case 'radiobutton3'
+        handles.current_func=TimeOfMaximum();
+end
+guidata(hObject,handles);
+
+
+%function for selecting the method for computing the tail
+% --- Executes when selected object is changed in uipanel2
+function uipanel2_SelectionChangeFcn(hObject, eventdata, handles)
+% hObject    handle to the selected object in uipanel13 
+% eventdata  structure with the following fields (see UIBUTTONGROUP)
+%	EventName: string 'SelectionChanged' (read only)
+%	OldValue: handle of the previously selected object or empty if none was selected
+%	NewValue: handle of the currently selected object
+% handles    structure with handles and user data (see GUIDATA)
+
+switch get(eventdata.NewValue,'tag') 
+    case 'radiobutton4'     
+        handles.current_tm=sorted();
+    case 'radiobutton5'        
+        set(handles.text2,'visible','on');
+        set(handles.edit2,'visible','on');
+end
+guidata(hObject,handles);
+
+
+%function for inserting the size of the lower and upper tails of the
+%evaluation function pdf
+function edit1_Callback(hObject, eventdata, handles)
+% hObject    handle to edit1 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of edit1 as text
+%        str2double(get(hObject,'String')) returns contents of edit1 as a double
+handles.tail_size=str2double(get(hObject,'String'));
+guidata(hObject,handles);
+   if isnan(handles.tail_size) || (handles.tail_size<0)
+      errordlg('Input must be a positive number (it should be at least 1000)');
+   end
+
+
+
+% --- Executes during object creation, after setting all properties.
+function edit1_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to edit1 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+
+
+%when the method tmp_sum is selected, the step for finding the threshold of
+%upper and lower tails needs to be specified
+function edit2_Callback(hObject, eventdata, handles)
+% hObject    handle to edit2 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of edit2 as text
+%        str2double(get(hObject,'String')) returns contents of edit2 as a double
+handles.step_thr=str2double(get(hObject,'String'));
+   if isnan(handles.step_thr) || (handles.step_thr<0 && handles.step_thr>1)
+      errordlg('Input must be a positive number less than 1');
+   end
+handles.current_tm=tmp_sum(handles.step_thr);     
+guidata(hObject,handles);
+
+
+% --- Executes during object creation, after setting all properties.
+function edit2_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to edit2 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
+    set(hObject,'BackgroundColor','white');
+end
+
+
+% --- Executes on button press in pushbutton3.
+%pushbutton3 launches the execution of the function for MIRI calculation
+function pushbutton3_Callback(hObject, eventdata, handles)
+% hObject    handle to pushbutton3 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+fig1=handles.fig1;
+handles1=guidata(fig1);
+compute_MIRI(handles,hObject,handles1);
+
+
+% --- Executes on button press in pushbutton2.
+%pushbutton2 launches the execution of the function for visualizing the
+%parameter pdfs
+function pushbutton2_Callback(hObject, eventdata, handles)
+% hObject    handle to pushbutton2 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+fig1=handles.fig1;
+handles1=guidata(fig1);
+plotpdf_param(handles,hObject,handles1);
+
+
+% --- Executes during object creation, after setting all properties.
+function gui_MIRI_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to gui_MIRI (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+

pdfEstimator.m

+%Concrete class representing an object for estimating a probability density function 
+%through the method ksdensity
+
+classdef pdfEstimator
+   
+    properties
+    end
+    
+    methods
+        
+        function [pdf, xbin]=evaluate_pdf(obj,samples,BinEdges)
+           
+           [pdf, xbin]=ksdensity(samples, BinEdges); 
+        end
+        
+    end
+    
+end
\ No newline at end of file

plotpdf_evalfunc.m

+%Function that plots and saves the estimated pdf of the evaluation function
+%for the chosen protein. 
+
+function plotpdf_evalfunc(handles,hObject,handles1)
+    try
+        load(fullfile(handles1.folder,handles1.chosen_protein,'pdf_eval_func.mat'));
+        for i=1:size(xbinT,1)
+            plot(handles.axes2,xbinT(i,:),ks_y_T(i,:));
+            hold on
+        end
+        xlabel(handles1.chosen_protein,'Interpreter','none');
+        ylabel('pdf of evaluation function');
+        
+        fh = figure('Visible','off');
+        h_new=copyobj(handles.axes2, fh);
+        %set(h_new,'Position',get(0,'DefaultAxesPosition'));
+        set(h_new, 'Units', 'Normalized');
+        set(h_new,'OuterPosition',[.1, .1, .85, .85]);
+        set(gcf,'Visible','off','CreateFcn','set(gcf,''Visible'',''on'')')
+        saveas(fh, fullfile(handles1.folder,handles1.chosen_protein,handles1.chosen_protein),'jpeg');
+        close(fh);
+        
+    catch ME
+        errordlg(ME.message) 
+    end
+end
\ No newline at end of file

plotpdf_param.m

+%Function that plots and saves the conditional pdfs of parameters that are
+%used for computing MIRI
+
+
+function plotpdf_param(handles,hObject,handles1)
+
+    try
+        load(fullfile(handles1.folder,handles1.chosen_protein,'pdf_param.mat'));
+        for i=1:length(handles1.nominal_p)
+            figure
+            for k=1:handles1.Nr
+                plot(xbin_param{1,k}(i,:),ks_param{1,k}(i,:),'b',xbin_param{2,k}(i,:),ks_param{2,k}(i,:),'r');
+                hold on
+            end
+            xlabel(handles1.param_name{i,1},'Interpreter','none');
+            ylabel('pdf');
+            legend('Conditional upper density','Conditional lower density');
+            savefig(fullfile(handles1.folder,handles1.chosen_protein,handles1.param_name{i,1}));
+        end
+    catch ME
+        errordlg(ME.message)
+    end
+end
\ No newline at end of file

sorted.m

+%Class that extends the abstract class TailMethod. The method compute_tail
+%computes the upper and lower tail of a pdf by sorting its values and
+%selecting the first and the last N samples. It takes in input
+%an object TimeBehavior containing the values of the evaluation function, 
+%the array perturbation of perturbed parameters and tail_size which specifies 
+%the number of samples in each tail. It returns in output the following variables:
+% - XiMin is the array containing parameter samples in the lower tail
+% - XiMax is the array containing parameter samples in the upper tail
+
+
+classdef sorted < TailMethod
+   
+    properties 
+    end
+    
+    methods
+       
+        function [XiMax,XiMin]=compute_tailspdf(obj,obj_TimeBehavior,perturbation,tail_size)       
+                    
+         samples=obj_TimeBehavior.evalFuncValues;
+         [samples_sorted, indexes_sorted]=sort(samples);                  
+         
+         XiMin=[];
+         XiMax=[];
+         for j=1:tail_size
+             XiMax=[XiMax;perturbation(indexes_sorted(length(samples)+1-j,1),:)];
+             XiMin=[XiMin;perturbation(indexes_sorted(j,1),:)];
+         end             
+        end       
+    end    
+end
\ No newline at end of file

start_simulation.m

+%Function executed when the button START SIMULATION in the main GUI is
+%pushed. It takes in input the handles of the main GUI.
+
+function start_simulation(handles,hObject)
+
+    %cell array where the measured proteins for all independent realizations
+    %are saved
+    AllResults=cell(handles.Nr,length(handles.Model.Species));
+    AllPerturbations=cell(handles.Nr,1);
+    
+    %retrieve parameters and observables of the model    
+    handles.num_param=length(handles.Model.Parameters);
+    handles.param_name=get(handles.Model.Parameters,'Name');
+    handles.nominal_p=cell2mat(get(handles.Model.Parameters,'value'))';
+    obs_name=get(handles.Model.Species,'Name');
+    
+    handles.time_results=[0:handles.step_size:handles.csObj.StopTime]';
+    
+    set(handles.text14,'string','');   
+    h = waitbar(0,'Please wait...');
+    
+    %for each independent realization k
+    for k=1:handles.Nr
+        waitbar(k/handles.Nr,h,['Realization number ',num2str(k)]);
+        
+        %Latin Hypercube generation for perturbing the parameter space
+        hypercube=LatinHypercube(handles.LBpi,handles.UBpi,handles.Nsample);
+        handles.perturbation=hypercube.generate_sample(handles.num_param);
+        perturbed_p=repmat(handles.nominal_p,length(handles.perturbation),1).*handles.perturbation;
+        AllPerturbations{k,1}=handles.perturbation;
+
+        %creation of a SimFunction object for performing model simulations
+        modelsim=createSimFunction(handles.Model, handles.param_name, obs_name,[],'UseParallel',true);
+        [tsim_all,psim_all]=modelsim(perturbed_p,[],[],handles.time_results);
+        
+        %ArrayResults is a cell array for saving measured proteins of each
+        %single realization. Then, these results are all saved in the cell
+        %array AllResults
+        ArrayResults=cell(handles.Nsample, length(handles.Model.Species));
+        
+        for i=1:handles.Nsample
+            for j=1:length(handles.Model.Species)
+                ArrayResults{i,j}=psim_all{i,1}(:,j);
+            end
+        end
+        
+        for j=1:length(handles.Model.Species)
+            AllResults{k,j}=ArrayResults(:,j);
+        end
+        
+    end
+    
+handles.AllResults=AllResults;
+handles.AllPerturbations=AllPerturbations;
+guidata(hObject,handles);
+delete(gcp)
+%save(fullfile(handles.folder,'measuresPeng2016.mat'),'-v7.3');
+
+close (h)  %close the waitbar
+set(handles.text14,'string','All done!');   
+drawnow;
+end

tmp_sum.m

+%Class that extends the abstract class TailMethod. The method compute_tail
+%computes the upper and lower tail of a pdf by selecting the upper and lower quartile. 
+%It takes in input an object TimeBehavior containing the values of the evaluation function, 
+%the array perturbation of perturbed parameters and tail_size which specifies 
+%the number of samples in each tail. It returns in output the following variables:
+% - XiMin is the array containing parameter samples in the lower tail
+% - XiMax is the array containing parameter samples in the upper tail
+
+classdef tmp_sum < TailMethod
+   
+    properties 
+        step_thr
+    end
+    
+    methods
+        
+        function obj=tmp_sum(step)
+           obj.step_thr=step;
+        end
+       
+        function [XiMax,XiMin]=compute_tailspdf(obj,obj_TimeBehavior,perturbation,tail_size)   
+                       
+        % estimation of the pdf of the evaluation function
+         pdf_obj=pdfEstimator();
+         samples=obj_TimeBehavior.evalFuncValues;
+         BinEdges=[min(samples):(max(samples)-min(samples))/length(samples):max(samples)];
+         [ks_y,xbin]=pdf_obj.evaluate_pdf(obj_TimeBehavior.evalFuncValues,BinEdges);
+            
+         [low_index, high_index]=obj.compute_thr(xbin,ks_y,obj_TimeBehavior.evalFuncValues,tail_size);
+      
+         XiMin=[];
+         XiMax=[];
+         for j=1:size(obj_TimeBehavior.evalFuncValues,1)
+           if high_index(j,1)==1
+              XiMax=[XiMax;perturbation(j,:)];
+           end
+           if low_index(j,1)==1
+             XiMin=[XiMin;perturbation(j,:)]; 
+           end
+         end    
+         
+        end
+        
+       %The method compute_thr calculates the upper and lower quartile of a pdf
+       
+       function [low_index,high_index]=compute_thr(obj,BinEdges,ks_y,Results,tail_size)
+           
+         low_thr=0;
+         high_thr=1;
+ 
+         while(low_thr < high_thr)
+          tmp_sum=0;
+          center_pdf=[];
+          for h=1:length(ks_y)
+            tmp_sum=tmp_sum+(BinEdges(2)-BinEdges(1))*ks_y(h); 
+            if(tmp_sum>=low_thr && tmp_sum<=high_thr)
+             center_pdf=[center_pdf, BinEdges(h)];
+           end
+         end
+     
+         low_index=Results<=center_pdf(1);
+         low_tail=Results(low_index);
+         high_index=Results>=center_pdf(end);
+         high_tail=Results(high_index);
+     
+         if(length(low_tail)<tail_size)
+          low_thr=low_thr+obj.step_thr;
+         end     
+         if(length(high_tail)<tail_size)
+          high_thr=high_thr-obj.step_thr;
+         end
+     
+         if(length(low_tail)>=tail_size && length(high_tail)>=tail_size)
+           disp(['THE LOW THRESHOLD IS ',num2str(low_thr)])
+           disp(['THE HIGH THRESHOLD IS ',num2str(high_thr)])
+           break
+         end 
+     
+         if(low_thr>=high_thr)
+           disp('It is not possible to find the low and high thresholds')
+         end     
+        end
+      end 
+    end    
+end
\ No newline at end of file