function Call_Hybrid_Model

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Authors:
% Turner, Leo; Burbanks, Andrew; Cerasuolo, Marianna
% 
% Affiliation:
% University of Portsmouth, UK
%
% Date:
% 13/05/2022
%
% Purpose:
% Call the simulation code for each schedule and treatment in batches.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

%% Input Parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User sets the environment here
% For O,G,H these will be copied into the BC in the main function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

InitialO=0.014;
InitialG=0.05;
InitialH=10^(-7.2);
InitialATP=1;  %--- New line

% List of available treatments (do not edit unless defining new treatments)
treatments = ["none"; "enza"; "caba"; "ever"; "enza_then_caba"; "enza_then_ever"; "enza_alt_caba"; "enza_alt_ever"; "enza_with_caba"; "enza_with_ever"];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User selects the treatments here - a separate run will be done for each
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% prefix to name this particular batch of runs
% batch_name = "params20220513" % for the preprint
batch_name = "params20220801_0000"

start_weeks = 5;
treatment_weeks = 12; %12; %18, %52;
end_weeks = 0;

%%%ADB[
all_batch_names = [
    "20220906_0000",
    "20220906_0001",
    "20220906_0002", 
    "20220906_0003", 
    "20220906_0004",
    "20220906_0005",
    "20220906_0006",
    "20220906_0007",
    "20220906_0008",
    "20220906_0009",
    "20220906_0010",
    "20220906_0011",
    "20220906_0012", 
    "20220906_0013", 
    "20220906_0014",
    "20220906_0015",
    "20220906_0016",
    "20220906_0017",
    "20220906_0018",
    "20220906_0019",
    "20220906_0020",
    "20220906_0021",
    "20220906_0022", 
    "20220906_0023", 
    "20220906_0024",
    "20220906_0025",
    "20220906_0026",
    "20220906_0027",
    "20220906_0028",
    "20220906_0029",
    "20220906_0030",
    "20220906_0031",
    "20220906_0032", 
    "20220906_0033", 
    "20220906_0034",
    "20220906_0035",
    "20220906_0036",
    "20220906_0037",
    "20220906_0038",
    "20220906_0039",
    "20220906_0040",
    "20220906_0041",
    "20220906_0042", 
    "20220906_0043", 
    "20220906_0044",
    "20220906_0045",
    "20220906_0046",
    "20220906_0047",
    "20220906_0048",
    "20220906_0049"
    ]
dummy = [
];
all_treatment_weeks = [52, 18, 12];
for i_batch_name=1:length(all_batch_names)
    batch_name = all_batch_names(i_batch_name)
    for i_treatment_weeks=1:length(all_treatment_weeks)
        treatment_weeks = all_treatment_weeks(i_treatment_weeks)
%%%ADB]

selected_treatments = 1:10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Drug Doses. Treatments are switched on/off by the choice of treatment,
% above, so these should *not* be changed between runs unless the user
% is recalibrating dosages.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

EnzaBC=0.02;            %mM
EverBC=11.2;            %nuM
%CabaBC=1.85; %20220513 paper figure 18, 20220725 stats   %nuM
CabaBC=1.79; %20220509 paper other figures, 20220801 stats   %nuM

do_plot_probabilities = true; % if set to true, probability curves will be plotted once for the batch.
do_independent_sensitivity = true; % generate indep rand() for Sp<->Rp choice within proliferation.
do_dry_run = false; % switch off chemical integration, to see drugs and CA only
do_figures = false; % plot some figures directly after each simulation
do_reproduce_bug = false; % reproduce the issue with transpose errors

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of user-selected parameters - do not edit beyond this point
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ii=selected_treatments
treatment = treatments(ii);

% total simulation time, can extend beyond all treatments
T = 7*(start_weeks + treatment_weeks + end_weeks);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Automatic naming convention used for data files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% simulation name
schedule = start_weeks+"_"+treatment_weeks+"_"+end_weeks;
simulation_name = batch_name+"_"+schedule+"_"+treatment;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute treatment details from the information selected above
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set default values for treatment period, computed from the above
T_start = 7*start_weeks;
T_stop = T_start + 7*treatment_weeks;

% by default, no treatments begin until after the simulation ends,
% in other words, all treatments are switched off by default.
T_enza = T+1; % first enza dose
T_caba = T+1; % first caba dose
T_ever = T+1; % first ever dose

% default for enza; repeat every stop-start steps, on constantly.
enza_period = T_stop-T_start;
enza_length = enza_period;

% caba is given once every 2 weeks
caba_period = 14;
caba_length = 1;

% ever is given once every other day
ever_period = 2;
ever_length = 1;

% gap timings
mid_week = 4;
next_day = 1;

if treatment == "none" %1
    % default values are above

elseif treatment == "enza" %2
    T_enza = T_start;

elseif treatment == "caba" %3
    T_caba = T_start;

elseif treatment == "ever" %4
    T_ever = T_start;

elseif treatment == "enza_then_caba" %5
    T_enza = T_start;
    enza_length = enza_period/2;
    T_caba = T_enza+(enza_length-1)+mid_week;

elseif treatment == "enza_then_ever" %6
    T_enza = T_start;
    enza_length = enza_period/2;
    T_ever = T_enza+(enza_length-1)+next_day;

elseif treatment == "enza_alt_caba" %7
    T_enza = T_start;
    enza_period = 14;
    enza_length = 7;
    T_caba = T_enza+(enza_length-1)+mid_week;

elseif treatment == "enza_alt_ever" %8
    T_enza = T_start;
    enza_period = 2;
    enza_length = 1;
    T_ever = T_enza+(enza_length-1)+next_day;

elseif treatment == "enza_with_caba" %9
    T_enza = T_start;
    T_caba = T_enza+mid_week;

elseif treatment == "enza_with_ever" %10
    T_enza = T_start;
    T_ever = T_enza+next_day;
    
else
    error("Unknown treatment");
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET UP QUANTITIES TO CALL SIMULATION FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define the grid geometry - do not change
h=0.1;
a=0.0;
b=10.0-3*h; %Ensure centre and boundary, divisibility by 10.

% Proportion of cells affected each sub loop
Proportion=1.0/10.0;

% Compute size of the grid; total size dimI x dimI
% We ensure that there is a central cell and boundary cells.
dimI=round((b-a)/h)+3;

% Positions, in normalized coordinates
indexI=1:dimI;
w=(indexI-1)*h+a;
[W,Z]=meshgrid(w,w);

% Initial conditions zeroed
InitialCaba=0;
InitialEnza=0;
InitialEver=0;
InitialPhi=0;

% Defines the various matrices used as they are at T=0.
O=OStart(W,Z);
G=GStart(W,Z);
H=HStart(W,Z);
ATP=ATPStart(W,Z);
Enza=EnzaStart(W,Z);
Ever=EverStart(W,Z);
Caba=CabaStart(W,Z);
Cells=CellStart(W,Z);
Phi=PhiStart(W,Z);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Call The main program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Hybrid_Model(batch_name, schedule, treatment, O,G,H,ATP,Enza,Caba,Ever,Cells,Phi,T,dimI,h,T_start,T_enza,T_caba,T_ever,T_stop,enza_period,enza_length, caba_period,caba_length,ever_period,ever_length, EnzaBC, CabaBC, EverBC,Proportion,do_dry_run,do_figures,do_reproduce_bug,do_independent_sensitivity,do_plot_probabilities)

do_plot_probabilities = false; %We only plot these once per batch.

end %next treatment

%%%ADB[
    end %next treatment_weeks
end %next batch_name
%%%ADB]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sets every cell to Initial O2 condition
    function O=OStart(X,Y)
        O=InitialO+(X+Y).*0;
    end

% Sets every cell to acidity condition
    function H=HStart(X,Y)
        H=InitialH+(X+Y).*0;
    end

% sets every cell to initial glucose conditon
    function G=GStart(X,Y)
        G=InitialG+(X+Y).*0;
    end

% sets every cell to initial ATP conditon %--- New line // Define ATP initial BC
    function ATP=ATPStart(X,Y) 
        ATP=InitialATP+(X+Y).*0;
    end

% Sets the cabazatixal starting concentration for every cell
    function Caba=CabaStart(X,Y)
        Caba=InitialCaba+(X+Y).*0;
    end

% Sets the Enzalutamide starting concentration for every cell
    function Enza=EnzaStart(X,Y)
        Enza=InitialEnza+(X+Y).*0;
    end

% Sets the Ever starting concentration for every cell
    function Ever=EverStart(X,Y)
        Ever=InitialEver+(X+Y).*0;
    end

% Creates Grid for recording prolifs
    function Phi=PhiStart(X,Y)
        Phi=InitialPhi+(X+Y).*0;
    end

% Build initial framework of cells, including boundary region
    function cells=CellStart(W,Z)
        assert(all(size(W)==size(Z)) && all(size(W)==[dimI,dimI]));

        % central 3x3 cluster of Sp cells
        cells = zeros(dimI);
        centrei = floor(dimI/2)+1;
        cells(centrei-1:centrei+1, centrei-1:centrei+1) = 1;
        
        % boundary cells
        cells(:,1)=6;
        cells(:,dimI)=6;
        cells(1,:)=6;
        cells(dimI,:)=6;
    end
end
