Analysis

f_analysis

Code

function rst = f_analysis(stg,analysis,mmf,analysis_options)

if contains(analysis,analysis_options(1))
    disp("Starting " + analysis_options(1))
    rst.diag = f_diagnostics(stg,mmf);
    disp(analysis_options(1) + " completed successfully")
end

if contains(analysis,analysis_options(2))
    disp("Starting " + analysis_options(2))
    rst.opt = f_opt(stg,mmf);
    disp(analysis_options(2) + " completed successfully")
end

if contains(analysis,analysis_options(3))
    disp("Starting " + analysis_options(3))
    rst.gsa = f_gsa(stg,mmf);
    disp(analysis_options(3) + " completed successfully")
end
end

Calls the proper analysis functions depending on the analysis that was chosen on the settings file. The supported analysis right now are:

• Model diagnostics functions

• Optimization

• Global sensitivity analysis

• Inputs

  • stg

  • analysis - (string) analysis being run (stg.analysis)

• Outputs - rst

Diagnostics

f_diagnostics

Code

function rst = f_diagnostics(stg,mmf)

% Run the model and obtain scores for fitness Multi Core
if stg.optmc
    disp("Running the model and obtaining Scores (Multicore)")

    pa = stg.pa;
    % Iterate over the parameter arrays to be tested
    parfor n = stg.pat
        [~,rst(n),~] = f_sim_score(pa(n,:),stg,mmf);
    end
    
    % Run the model and obtain scores for fitness single Core
else
    disp("Running the model and obtaining Scores (Singlecore)")
    
    % Iterate over the parameter arrays to be tested
    for n = stg.pat
        [~,rst(n),~] = f_sim_score(stg.pa(n,:),stg,mmf);
    end
end
end

Used to understand the effects of different parameters sets on model behaviour or in comparing different parameters sets.

It loads the user defined configurations, performs all the needed simulations, and calculates scores of the error functions either per experimental output, per experiment, or in total (check results).

• Inputs

  • stg - stg.optmc , stg.pat

• Outputs - rst (diagnostics results)

Optimization

f_opt

Code

function rst = f_opt(stg,mmf)
% Call function to run fmincon optimization algorithm if chosen in settings

if stg.fmincon
    rst(1) = f_opt_fmincon(stg,mmf);
end

% Call function to run simulated annealing optimization algorithm if chosen
% in settings
if stg.sa
    rst(2) = f_opt_sa(stg,mmf);
end

% Call function to run pattern search optimization algorithm if chosen in
% settings
if stg.psearch
    rst(3) = f_opt_psearch(stg,mmf);
end

% Call function to run genetic algorithm optimization if chosen in settings
if stg.ga
    rst(4) = f_opt_ga(stg,mmf);
end

% Call function to run Particle swarm optimization algorithm if chosen in
% settings
if stg.pswarm
    rst(5) = f_opt_pswarm(stg,mmf);
end

% Call function to run Surrogate optimization algorithm if chosen in
% settings
if stg.sopt
    rst(6) = f_opt_sopt(stg,mmf);
end
end

Calls the correct optmizer or optimizers that have been chosen in the settings file.

• Inputs

  • stg - stg.fmincon, stg.sa, stg.psearch, stg.ga, stg.pswarm, stg.sopt

• Outputs - rst (optimization results)

f_opt_start

Code

function [spoint,spop] = f_opt_start(stg)

% Set the randomm seed for reproducibility
rng(stg.rseed);

% Optimization Start method 1
if stg.osm == 1
    
    % Get a random starting point or group of starting points, if using
    % multistart, inside the bounds
    spoint = lhsdesign(stg.msts,stg.parnum).*(stg.ub-stg.lb)+stg.lb;
    
    % Get a group of ramdom starting points inside the bounds
    spop = lhsdesign(stg.popsize,stg.parnum).*(stg.ub-stg.lb)+stg.lb;
    
    % Optimization Start method 2
elseif stg.osm == 2
    
    % Get a random starting point or group of starting points, if using
    % multistart, near the best point
    spoint = stg.bestpa - stg.dbpa +...
        (stg.dbpa*2*lhsdesign(stg.msts,stg.parnum));
    
    % Get a group of ramdom starting points near the best point
    spop = stg.bestpa - stg.dbpa +...
        (stg.dbpa*2*lhsdesign(stg.popsize,stg.parnum));
end
end

Creates the starting parameter set or sets of the optimizations, if single or multistart selected in settings file.

It supports two different random distributions for the starting points.

• Inputs

  • stg - stg.rseed, stg.osm, stg.msts, stg.parnum, stg.ub, stg.lb, stg.popsize, stg.bestpa, stg.dbpa

• Outputs

  • spoint - (double) starting parameter set for the optimization

  • spop - (double) Starting parameter sets for multiple start optimizations

f_opt_fmincon/sa/psearch/ga/pswarm/sopt

Code

f_opt_fmincon

function rst = f_opt_fmincon(stg,mmf)

% Set the randomm seed for reproducibility
rng(stg.rseed);

% Set the starting point for the optimization
[startpoint,~] = f_opt_start(stg);

% Get the optimization options from settings
options = stg.fm_options;
options.UseParallel = stg.optmc;

% Display console messages if chosen in settings
if stg.optcsl
    options.Display = 'iter-detailed';
end

% Display plots if chosen in settings
if stg.optplots
    options.PlotFcn = ...
        {@optimplotx,@optimplotfunccount,@optimplotfval,...
        @optimplotstepsize,@optimplotfirstorderopt};
end

% Optimize the model with multiple starting points if chosen in settings
if stg.mst
    parfor n = 1:stg.msts
        disp(string(n))
        [x(n,:),fval(n),exitflag(n),output(n)] =...
            fmincon(@(x)f_sim_score(x,stg,mmf),startpoint(n,:),...
            [],[],[],[],stg.lb,stg.ub,[],options);
    end
    
    % Optimize the model
else
    [x(1,:),fval(1),exitflag(1),output(1)] =...
        fmincon(@(x)f_sim_score(x,stg,mmf),startpoint(1,:),...
        [],[],[],[],stg.lb,stg.ub,[],options);
end

% Save results
rst.name = 'fmincon';
rst.x = x;
rst.fval = fval;
rst.exitflag = exitflag;
rst.output = output;
end

f_opt_sa

function rst = f_opt_sa(stg,mmf)

% Set the randomm seed for reproducibility
rng(stg.rseed);

% Set the starting point for the optimization
[startpoint,~] = f_opt_start(stg);

% Get the optimization options from settings
options = stg.sa_options;
options.MaxTime = stg.optt;
options.InitialTemperature = ones(1,stg.parnum)*2;

% Display console messages if chosen in settings
if stg.optcsl
    options.Display = 'iter';
end

% Display plots if chosen in settings
if stg.optplots
    options.PlotFcn =...
        {@saplotbestf,@saplottemperature,@saplotf,@saplotstopping,@saplotx};
end

% Optimize the model with multiple starting points if chosen in settings
if stg.mst
    parfor n = 1:stg.msts
        disp(string(n))
        [x(n,:),fval(n),exitflag(n),output(n)] =...
            simulannealbnd(@(x)f_sim_score(x,stg,mmf),startpoint(n,:),...
            stg.lb,stg.ub,options);
    end
    
    % Optimize the model
else
    [x(1,:),fval(1),exitflag(1),output(1)] =...
        simulannealbnd(@(x)f_sim_score(x,stg,mmf),startpoint(1,:),...
        stg.lb,stg.ub,options);
end

% Save results
rst.name = 'Simulated annealing';
rst.x = x;
rst.fval = fval;
rst.exitflag = exitflag;
rst.output = output;
end

f_opt_psearch

function rst = f_opt_psearch(stg,mmf)

% Set the randomm seed for reproducibility
rng(stg.rseed);

% Set the starting point for the optimization
[startpoint,~] = f_opt_start(stg);

% Get the optimization options from settings
options = stg.psearch_options;
options.MaxTime = stg.optt;
options.UseParallel = stg.optmc;

% Display console messages if chosen in settings
if stg.optcsl
    options.Display = 'iter';
end

% Display plots if chosen in settings
if stg.optplots
    options.PlotFcn = ...
        {@psplotbestf,@psplotfuncount,@psplotmeshsize,@psplotbestx};
end

% Optimize the model with multiple starting points if chosen in settings
if stg.mst
    parfor n = 1:stg.msts
        disp(string(n))
        [x(n,:),fval(n),exitflag(n),output(n)] =...
            patternsearch(@(x)f_sim_score(x,stg,mmf),startpoint(n,:),[],[], ...
            [],[],stg.lb,stg.ub,[],options);
    end
    
    % Optimize the model
else
    [x(1,:),fval(1),exitflag(1),output(1)] =...
        patternsearch(@(x)f_sim_score(x,stg,mmf),startpoint(1,:),[],[], ...
        [],[],stg.lb,stg.ub,[],options);
end

% Save results
rst.name = 'Pattern search';
rst.x = x;
rst.fval = fval;
rst.exitflag = exitflag;
rst.output = output;
end

f_opt_ga

function rst = f_opt_ga(stg,mmf)

% Set the randomm seed for reproducibility
rng(stg.rseed);

% Get the optimization options from settings
options = stg.ga_options;
options.MaxTime = stg.optt;
options.UseParallel = stg.optmc;
options.PopulationSize = stg.popsize;

% Display console messages if chosen in settings
if stg.optcsl
    options.Display = 'iter';
end

% Display plots if chosen in settings
if stg.optplots
    options.PlotFcn = {@gaplotbestf,...
        @gaplotexpectation,@gaplotrange,@gaplotscorediversity,...
        @gaplotstopping,@gaplotscores,@gaplotdistance,...
        @gaplotselection};
end

% Set the starting population for the optimization
[~,startpop] = f_opt_start(stg);
options.InitialPopulationMatrix = startpop;

% Optimize the model with multiple starting points if chosen in settings
if stg.mst
    parfor n = 1:stg.msts
        disp(string(n))
        [x(n,:),fval(n)] = ga(@(x)f_sim_score(x,stg,mmf),stg.parnum,...
            [],[],[],[],stg.lb,stg.ub,[],options);
    end
    
    % Optimize the model
else
    [x(1,:),fval(1)] = ga(@(x)f_sim_score(x,stg,mmf),stg.parnum,...
        [],[],[],[],stg.lb,stg.ub,[],options);
end

% Save results
rst.name = 'Genetic algorithm';
rst.x = x;
rst.fval = fval;
rst.exitflag = [];
rst.output = [];
end

f_opt_pswarm

function rst = f_opt_pswarm(stg,mmf)

% Set the randomm seed for reproducibility
rng(stg.rseed);

% Get the optimization options from settings
options = stg.pswarm_options;
options.MaxTime = stg.optt;
options.UseParallel = stg.optmc;
options.SwarmSize = stg.popsize;

% Display console messages if chosen in settings
if stg.optcsl
    options.Display = 'iter';
end

% Display plots if chosen in settings
if stg.optplots
    options.PlotFcn = @pswplotbestf;
end

% Set the starting population for the optimization
[~,startpop] = f_opt_start(stg);
options.InitialSwarmMatrix = startpop;

% Optimize the model with multiple starting points if chosen in settings
if stg.mst
    parfor n = 1:stg.msts
        disp(string(n))
        [x(n,:),fval(n),exitflag(n),output(n)] =...
            particleswarm(@(x)f_sim_score(x,stg,mmf),...
            stg.parnum,stg.lb,stg.ub,options);
    end
    
    % Optimize the model
else
    [x(1,:),fval(1),exitflag(1),output(1)] =...
        particleswarm(@(x)f_sim_score(x,stg,mmf),...
        stg.parnum,stg.lb,stg.ub,options);
end

% Save results
rst.name = 'Particle swarm';
rst.x = x;
rst.fval = fval;
rst.exitflag = exitflag;
rst.output = output;
end

f_opt_sopt

function rst = f_opt_sopt(stg,mmf)

% Set the randomm seed for reproducibility
rng(stg.rseed);

% Get the optimization options from settings
options = stg.sopt_options;
options.MaxTime = stg.optt;
options.UseParallel = stg.optmc;

% Display console messages if chosen in settings
if stg.optcsl
    options.Display = 'iter';
end

% Display plots if chosen in settings
if stg.optplots
    options.PlotFcn = @surrogateoptplot;
end

% Set the starting population for the optimization
[~,startpop] = f_opt_start(stg);
options.InitialPoints = startpop;

% Optimize the model with multiple starting points if chosen in settings
if stg.mst
    parfor n = 1:stg.msts
        disp(string(n))
        [x(n,:),fval(n),exitflag(n),output(n)] =...
            surrogateopt(@(x)f_sim_score(x,stg,mmf),stg.lb,stg.ub,options);
    end
    
    % Optimize the model
else
    [x(1,:),fval(1),exitflag(1),output(1)] =...
        surrogateopt(@(x)f_sim_score(x,stg,mmf),stg.lb,stg.ub,options);
end

% Save results
rst.name = 'Surrogate optimization';
rst.x = x;
rst.fval = fval;
rst.exitflag = exitflag;
rst.output = output;
end

These functions call built in MATLAB® functions that perform parameter optimization . For furher information relating to how these optimizers work please follow the links to the MATLAB® documentation. Optimizers used:

• f_opt_fmincon - fmincon

• f_opt_sa - Simmulated annealing

• f_opt_psearch - Pattern search

• f_opt_ga - Genetic algorihtm

• f_opt_pswarm - Particle swarm

• f_opt_sopt - Surrogate optmization

• Inputs - stg

• Outputs - Optimization results

Global Sensitivity Analysis

f_gsa

Code

function rst = f_gsa(stg,mmf)

rst = f_make_par_samples(stg);

rst = f_make_output_sample(rst,stg,mmf);

rst = f_calc_sensitivities(rst,stg);
end

Calls the global sensitivity analysis functions in the correct order.

f_make_par_samples

Code

function rst= f_make_par_samples(stg)
% Code inspired by Geir Halnes et al. 2009 paper. (Halnes, Geir, et al. J.
% comp. neuroscience 27.3 (2009): 471.)

% MAKE SAMPLE MATRICES
M1 = zeros(stg.sansamples, stg.parnum); % Pre-allocate memory for data
M2 = zeros(stg.sansamples, stg.parnum);
N = zeros(stg.sansamples, stg.parnum, stg.parnum);
rng(stg.rseed)

% Create a distribution for each parameter acording to settings(Note that
% the sampling is being performed in log space as the parameters and its
% bounds are in log space)
for i=1:stg.parnum
    % Uniform distribution truncated at the parameter bounds
    if stg.sasamplemode == 0
        M1(:,i) = stg.lb(i) +...
            (stg.ub(i)-stg.lb(i)).*rand(1,stg.sansamples);
        M2(:,i) = stg.lb(i) +...
            (stg.ub(i)-stg.lb(i)).*rand(1,stg.sansamples);
        % Normal distribution with mu as the best value for a parameter and
        % sigma as stg.sasamplesigma truncated at the parameter bounds
    elseif stg.sasamplemode == 1
        pd(i) = makedist('Normal','mu',stg.bestpa(i),...
            'sigma',stg.sasamplesigma);
        t(i) = truncate(pd(i),stg.lb(i),stg.ub(i));
        r{i} = random(t(i),stg.sansamples,1);
        r2{i} = random(t(i),stg.sansamples,1);
        M1(:,i) = r{i};
        M2(:,i) = r2{i};
        % Same as 1 without truncation
    elseif stg.sasamplemode == 2
        pd(i) = makedist('Normal','mu',stg.bestpa(i),...
            'sigma',stg.sasamplesigma);
        r{i} = random(pd(i),stg.sansamples,1);
        r2{i} = random(pd(i),stg.sansamples,1);
        M1(:,i) = r{i};
        M2(:,i) = r2{i};
        % Normal distribution centered at the mean of the parameter bounds
        % and sigma as stg.sasamplesigma truncated at the parameter bounds
    elseif stg.sasamplemode == 3
        pd(i) = makedist('Normal','mu',...
            stg.lb(i) + (stg.ub(i)-stg.lb(i))/2,'sigma',stg.sasamplesigma);
        t(i) = truncate(pd(i),stg.lb(i),stg.ub(i));
        r{i} = random(t(i),stg.sansamples,1);
        r2{i} = random(t(i),stg.sansamples,1);
        M1(:,i) = r{i};
        M2(:,i) = r2{i};
        % Same as 3 without truncation.
    elseif stg.sasamplemode == 4
        pd(i) = makedist('Normal','mu',...
            stg.lb(i) + (stg.ub(i)-stg.lb(i))/2,'sigma',stg.sasamplesigma);
        r{i} = random(pd(i),stg.sansamples,1);
        r2{i} = random(pd(i),stg.sansamples,1);
        M1(:,i) = r{i};
        M2(:,i) = r2{i};
    end
end

for i=1:stg.parnum
    % Replace the i:th column in M2 by the i:th column from M1 to obtain Ni
    N(:,:,i) = M2;
    N(:,i,i) = M1(:,i);
end

rst.M1=M1;
rst.M2=M2;
rst.N=N;

Creates parameter sets samples with specific parameter distributions that are used to perform the global sensitivity analysis.

• Inputs

  • stg - stg.sansamples, stg.parnum, stg.sasamplemode, stg.ub, stg.lb

• Outputs - M1, M2, N

Code inspired by Geir Halnes et al. 2009 paper.

f_make_output_sample

Code

function rst = f_make_output_sample(rst,stg,mmf)
% Code inspired by Geir Halnes et al. 2009 paper. (Halnes, Geir, et al. J.
% comp. neuroscience 27.3 (2009): 471.)

nSamples = stg.sansamples;
nPars = stg.parnum;
parameter_array = zeros(nSamples,nPars);
progress = 1;
time_begin = datetime;
D = parallel.pool.DataQueue;

afterEach(D, @progress_track);

for i=1:nSamples
    parameter_array(i,:) = rst.M1(i,:);
end

parfor i=1:nSamples
    [~,~,RM1(i)] = f_sim_score(parameter_array(i,:),stg,mmf);
    send(D, "GSA M1 ");
end
disp("GSA M1 Runtime: " + string(datetime - time_begin) +...
                "  All " + nSamples + " samples executed")


for i=1:nSamples
    fM1.sd(i,:) = [RM1(i).sd{:}];
    fM1.se(i,:) = RM1(i).se(:);
    fM1.st(i,:) = RM1(i).st;
    fM1.xfinal(i,:) = [RM1(i).xfinal{:}];
end

rst.fM1 = fM1;
clear a FM1

for  i=1:nSamples
    parameter_array(i,:)= rst.M2(i,:);
end

progress = 1;
parfor i=1:nSamples
    [~,~,RM2(i)] = f_sim_score(parameter_array(i,:),stg,mmf);
    send(D, "GSA M2 ");
end
disp("GSA M2 Runtime: " + string(datetime - time_begin) +...
                "  All " + nSamples + " samples executed")

for i=1:nSamples
    fM2.sd(i,:) = [RM2(i).sd{:}];
    fM2.se(i,:) = RM2(i).se(:);
    fM2.st(i,:) = RM2(i).st;
    fM2.xfinal(i,:) = [RM2(i).xfinal{:}];
end

rst.fM2 = fM2;
clear b FM2

for i=1:nSamples
    for j=1:nPars
        parameter_array(i,:,j)= rst.N(i,:,j);
    end
end

progress = 1;
parfor i=1:nSamples
    for j=1:nPars
        [~,~,RN{i,j}] = f_sim_score(parameter_array(i,:,j),stg,mmf);
    end
    send(D, "GSA N  ");
end
disp("GSA N  Runtime: " + string(datetime - time_begin) +...
                "  All " + nSamples + " samples executed")

for i=1:nSamples
    for j=1:nPars
        fN.sd(i,:,j) = [RN{i,j}.sd{:}];
        fN.se(i,:,j) = RN{i,j}.se(:);
        fN.st(i,:,j) = RN{i,j}.st;
        fN.xfinal(i,:,j) = [RN{i,j}.xfinal{:}];
    end
end

rst.fN = fN;
clear c FN

    function progress_track(name)
        progress = progress + 1;
        if mod(progress,ceil(nSamples/10)) == 0 && progress ~= nSamples
            disp(name + "Runtime: " + string(datetime - time_begin) +...
                "  Samples:" + progress + "/" + nSamples)
        end
    end
end

For each parameter set given in the matrices M1, M2, and N it runs the function f_sim_score generating new matrices fM1, fM2, and fN respectively.

• Inputs - M1, M2, N, stg.sansamples, stg.parnum,

• Outputs - fM1, fM2, fN

Code inspired by Geir Halnes et al. 2009 paper.

f_calc_sensitivities

Code

function rst = f_calc_sensitivities(rst,stg)

rst = remove_sim_error(rst,stg);

[rst.SiQ,rst.SiTQ,rst.Si,rst.SiT] = bootstrap(rst,stg);
end

function [SiQ,SiTQ,Si,SiT]=bootstrap(rst,stg)
% calculates confidence intervals.
fM1 = rst.fM1;
fM2 = rst.fM2;
fN = rst.fN;

scores_names_list = ["sd","se","st","xfinal"];

[Si,SiT] = bootstrap_h(fM1,fM2,fN,stg,scores_names_list);

if (isempty(stg.gsabootstrapsize))
    stg.gsabootstrapsize=ceil(sqrt(size(fM1.sd)));
end%if

fM1q = cell(stg.gsabootstrapsize,1);
fM2q = cell(stg.gsabootstrapsize,1);
fNq = cell(stg.gsabootstrapsize,1);

parfor j=1:stg.gsabootstrapsize
    [SiQh{j},SiTQh{j}] = bootstrap_hq(fM1,fM2,fN,stg,scores_names_list,j);
end%parfor

for n = 1:size(scores_names_list,2)
    for j=1:stg.gsabootstrapsize
        eval("SiQ." + scores_names_list(n) + "(j,:,:) = SiQh{j}." + scores_names_list(n) + ";")
        eval("SiTQ." + scores_names_list(n) + "(j,:,:) = SiTQh{j}." + scores_names_list(n) + ";")
    end
end
end%function

function [Si,SiT] = bootstrap_h(fM1,fM2,fN,stg,scores_names)
for n = 1:size(scores_names,2)
    eval("fM1h=fM1." + scores_names(n) + ";");
    eval("fM2h=fM2." + scores_names(n) + ";");
    eval("fNh=fN." + scores_names(n) + ";");
    
    [Sih,SiTh] = calcSobolSaltelli(fM1h,fM2h,fNh,stg);
    
    eval("Si." + scores_names(n) + "=Sih;");
    eval("SiT." + scores_names(n) + "=SiTh;");
end
end

function [Si,SiT] = bootstrap_hq(fM1,fM2,fN,stg,scores_names,j)
for n = 1:size(scores_names,2)
    eval("fM1h=fM1." + scores_names(n) + ";");
    eval("fM2h=fM2." + scores_names(n) + ";");
    eval("fNh=fN." + scores_names(n) + ";");
    
    rng(j*stg.rseed)
    I=ceil(rand(size(fM1h,1),1)*size(fM1h,1));
    fM1q = fM1h(I,:);
    fM2q = fM2h(I,:);
    fNq = fNh(I,:,:);
    
    [Sih,SiTh] = calcSobolSaltelli(fM1q,fM2q,fNq,stg);
    eval("Si." + scores_names(n) + "=Sih;");
    eval("SiT." + scores_names(n) + "=SiTh;");
end
end

function [Si,SiT] = calcSobolSaltelli(fM1,fM2,fN,stg)
%Code inspired by Geir Halnes et al. 2009 paper. (Halnes, Geir, et al. J. comp. neuroscience 27.3 (2009): 471.)

[Nsamples,Nvars,Npars]=size(fN);

if(stg.sasubmean) % Makes the model more stable
    fM1 = fM1 - mean(fM1,1);
    fM2 = fM2 - mean(fM2,1);
    for i=1:Npars
        fN(:,:,i) =  fN(:,:,i) - mean(fN(:,:,i),1);
    end
end

EY2 = mean(fM1.*fM2); % Valid definition (see Halnes et. al. Appendix)
VY = sum(fM1.^2)/(Nsamples-1) - EY2;
VYT = sum(fM2.^2)/(Nsamples-1) - EY2;

Si = zeros(Nvars,Npars);
SiT= zeros(Nvars,Npars);

for i=1:Npars
    Si(:,i) = (sum(fM1.*fN(:,:,i))/(Nsamples-1) - EY2)./VY;
    SiT(:,i) = 1 - (sum(fM2.*fN(:,:,i))/(Nsamples-1) - EY2)./VYT;
end
end

function rst = remove_sim_error(rst,stg)
error=[];
error_helper=[];

for n = 1:size(rst.fM1.sd,1)
    if max(rst.fM1.sd(n,:)) == stg.errorscore
        error = [error,n];
    end
    if max(rst.fM2.sd(n,:)) == stg.errorscore
        error = [error,n];
    end
    for m = 1:stg.parnum
        if max(rst.fN.sd(n,:,m)) == stg.errorscore
            error_helper = 1;
        end
    end
    if error_helper == 1
        error = [error,n];
    end
    error_helper = 0;
end

error = unique(error);

if ~isempty(error)
    for n = size(error,2):-1:1
        rst.fM1.sd(error(n),:) = [];
        rst.fM2.sd(error(n),:) = [];
        rst.fN.sd(error(n),:,:) = [];
        
        rst.fM1.se(error(n),:) = [];
        rst.fM2.se(error(n),:) = [];
        rst.fN.se(error(n),:,:) = [];
        
        rst.fM1.st(error(n),:) = [];
        rst.fM2.st(error(n),:) = [];
        rst.fN.st(error(n),:,:) = [];
        
        rst.fM1.xfinal(error(n),:) = [];
        rst.fM2.xfinal(error(n),:) = [];
        rst.fN.xfinal(error(n),:,:) = [];
    end
end
end

Takes the matrices fM1, fM2, and fN and calculates sensitivity indexes. It calculates indexes based on the following oputputs of the f_sim_score function:

• The scores of each experimental output

• The scores of each experiment

• The total score

• The value of each experimental outputs at the end of the simulation

• Inputs - fM1, fM2, fN, stg.sasubmean

• Outputs - SI, STI

Code modified from the Geir Halnes et al. 2009 paper.

References

Halnes, G., Ulfhielm, E., Ljunggren, E.E., Kotaleski, J.H. and Rospars, J.P., 2009. Modelling and sensitivity analysis of the reactions involving receptor, G-protein and effector in vertebrate olfactory receptor neurons. Journal of Computational Neuroscience, 27(3), p.471.