
%% Setting parameters
clear 
param.bcl = 1000;
param.model = @model_Torord_dynCl;

%%STEP 1: Parameter factors generation
%                  1     2      3     4        5       6       7      
parameter_names={'INa','INaL','Ito','ICaL', 'ICaNa', 'ICaK', 'IKr', ...
                'IKs','IK1','INCX','INaK','INab', 'IKb', 'IpCa', 'ICab', 'IClCa', 'IClb'};
%                 8     9     10     11     12     13      14     15       16       17              

n_parameters = length(parameter_names);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
trials = 1000; % Select the number of initial models
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu=1;
sigma = 0.2; % standard deviation
rng('shuffle');
p=normrnd(mu, sigma, [ trials, n_parameters]); 
%% Simulation in control conditions
% Here, we make an array of parameter structures
params(1:size(p,1)) = param; % These are initially all the default parametrisation
for iParam = 1:size(p,1)
    params(iParam).ICaL_Multiplier = p(iParam,4); 
    params(iParam).INa_Multiplier = p(iParam,1);
    params(iParam).Ito_Multiplier = p(iParam,3);
    params(iParam).INaL_Multiplier = p(iParam,2);
    params(iParam).IKr_Multiplier = p(iParam,7);
    params(iParam).IKs_Multiplier = p(iParam,8);
    params(iParam).IK1_Multiplier = p(iParam,9);
    params(iParam).IKb_Multiplier = p(iParam,13);
    params(iParam).INaCa_Multiplier = p(iParam,10);
    params(iParam).INaK_Multiplier = p(iParam,11);
    params(iParam).INab_Multiplier = p(iParam,12);
    params(iParam).ICab_Multiplier = p(iParam,15);
    params(iParam).IpCa_Multiplier = p(iParam,14);
    params(iParam).ICaNa_Multiplier = p(iParam,5);
    params(iParam).ICaK_Multiplier = p(iParam,6);
    params(iParam).IClCa_Multiplier = p(iParam,16);
    params(iParam).IClb_Multiplier = p(iParam,17); 
end

options = [];
beats = 500;
ignoreFirst = beats - 3;

myCluster=parcluster('local');
poolobj = gcp('nocreate');
if isempty(poolobj)
    poolobj = parpool(myCluster.NumWorkers);
end
poolobj.IdleTimeout= 120;

parfor i = 1:size(params,2) 
    X0 = getStartingState('Torord_endo_dynCl');
    [time{i}, X{i}] = modelRunner(X0, options, params(i), beats, ignoreFirst);
end
delete(poolobj);
s.name='PopulationTorORd_raw';
s.params=params;
s.result=[(time'),(X')];
save(['PopulationTorORd_raw.mat'], 's')
clearvars -except params p options beats ignoreFirst
%% Simulation of the population with a 90% block of IKs
myCluster=parcluster('local');
poolobj = gcp('nocreate');
if isempty(poolobj)
    poolobj = parpool(myCluster.NumWorkers);
end
poolobj.IdleTimeout= 120;

params_IKs90= params;
for iParam = 1:size(p,1)
    params_IKs90(iParam).IKs_Multiplier = params_IKs90(iParam).IKs_Multiplier*0.1; 
end

% Simulation and output extraction
parfor i = 1:size(params_IKs90,2) 
    X0 = getStartingState('Torord_endo_dynCl');
    [time{i}, X{i}] = modelRunner(X0, options, params_IKs90(i), beats, ignoreFirst);
end
delete(poolobj);

s.name='PopulationTorORd_raw_IKs90';
s.params=params_IKs90;
s.result=[(time'),(X')];
save(['PopulationTorORd_raw_IKs90.mat'], 's')
clearvars -except params p options beats ignoreFirst
%% Simulation of the population with a 70% block of IKr
myCluster=parcluster('local');
poolobj = gcp('nocreate'); % If no pool, do not create new one.
if isempty(poolobj)
    poolobj = parpool(myCluster.NumWorkers);
end
poolobj.IdleTimeout= 120;
params_IKr70= params;
for iParam = 1:size(p,1)
    params_IKr70(iParam).IKr_Multiplier = params_IKr70(iParam).IKr_Multiplier*0.3; 
end
% Simulation and output extraction
parfor i = 1:size(params_IKr70,2) 
    X0 = getStartingState('Torord_endo_dynCl');
    [time{i}, X{i}] = modelRunner(X0, options, params_IKr70(i), beats, ignoreFirst);
end
delete(poolobj);

s.name='PopulationTorORd_raw_IKr70';
s.params=params_IKr70;
s.result=[(time'),(X')];
save(['PopulationTorORd_raw_IKr70.mat'], 's')
clearvars -except params p options beats ignoreFirst
%% Simulation of the population with a 50% block of IK1
myCluster=parcluster('local');
poolobj = gcp('nocreate'); % If no pool, do not create new one.
if isempty(poolobj)
    poolobj = parpool(myCluster.NumWorkers);
end
poolobj.IdleTimeout= 120;
params_IK1_50= params;
for iParam = 1:size(p,1)
    params_IK1_50(iParam).IK1_Multiplier = params_IK1_50(iParam).IK1_Multiplier*0.5; % because of this line, the default parametrisation needs to have IKr_Multiplier defined (otherwise Matlab complains about different structure type). 
end
% Simulation and output extraction
parfor i = 1:size(params_IK1_50,2) 
    X0 = getStartingState('Torord_endo_dynCl');
    [time{i}, X{i}] = modelRunner(X0, options, params_IK1_50(i), beats, ignoreFirst);
end
s.name='PopulationTorORd_raw_IK1_50';
s.params=params_IK1_50;
s.result=[(time'),(X')];
save(['PopulationTorORd_raw_IK1_50.mat'], 's')
clearvars -except params p options beats ignoreFirst