%{
This script runs the simulations used for the publication

Parameter references
[1] Sgro et al, From intracellular signaling to population oscillations: bridging size- and time-scales in collective behavior
[2] Hoefer et al, Dictyostelium discoideum: Cellular Self-Organization in an Excitable Biological Medium
[3] Vidal-Henriquez et at, Spontaneous center formation in Dictyostelium discoideum
[HF] From experiments by Hugh Ford

Written 12/2021 by Angelika Manhart
%}

%% Set base parameters

p.Lx=4; % size of domain in x in mm
p.Ly=4; % size of domain in y in mm
p.maxTime=90; % in min
p.nX=round(p.Lx*50); % nr of grid points in x
p.nY=round(p.Ly*50); % nr of grid points in y

% cell movement
cellDens=8000; % per mm^2: 4000-5000 [HF]
p.N=round(p.Lx*p.Ly*cellDens); % nr of cells
p.movNoise=0; % positional noise
p.repRad=5e-3; % cell repulsion range in mm (to do with cell size) in mm [HF]
p.repStrength=0.1; % cell repulsion strength

% chemotaxis parameters
p.chemDesens=2.67*1e-5; % [HF] & numerics
p.chemTimeScale=6.7*1e-3; % [HF] & numerics
p.turnRate=3*1e-6; % [HF] & numerics
p.speedModelP1=1e-8; % [HF] & numerics ("gamma")
p.speedModelP2=0.6; % [HF] & numerics ("beta")

% internal cell dynamics
p.epsi=0.1;   %[1] ratio between activator and inhibitor timescales
p.preA=0.058;  % [1] magnitude of response of pre-processing unit
p.preKd=1e-5;  % [1] threshold for response of pre-processing unit
p.inhDecay=0.5; % [1] inhibitor decay
p.inhBasal=1.2; % [1] inhibitor parameter
p.intTimeScale=1/6; % (estimated from data) 1/n means roughly n spike every 30mins in a constant ext cAMP scenario

% signal dynamics
p.sigDiff=2.4e-2; % signal diffusion constant in mm^2/min: 2.4e-2 [3]
p.sigDecay=5; % signal decay rate in /min: 1-10 [3]
p.sigSecret=1e3; % % signal secretion rate per active cell
p.sigSpreadRad=2*0.05; % the radius over which the secreted signal is initially spread in mm;

p.dt=1/24; % time step in minutes

% simulation setting
p.timeDepParas=0;
p.optBC=1; % 0 .. PBC, 1.. no flux
p.optSigDisc=1; % 0 .. explicit, 1.. implcit
p.saveEvery=2;
p.plotEvery=2;
p.isVideo=1;

addpath('Program_Files\')

%%
%{
1) cell repulsion without time dependent internal dynamic and with cell movement
2) cell repulsion with time dependent internal dynamics and with cell movement
3) cell repulsion without time dependent internal dynamics and without cell movement (I understand that this is the same as without cell repulsion)
4) no cell repulsion without time dependent internal dynamics and with cell movement
%}

%% 1) cell repulsion without time dependent internal dynamic and with cell movement
load('ICs/ICs_naturalSpiral.mat');

p.repStrength=0.1; % cell repulsion strength
p.timeDepParas=0;
p.speedModelP1=1e-8; % [HF] & numerics 

results=main_Dicty(p,ICs);


%% 2) cell repulsion with time dependent internal dynamics and with cell movement
% changes p.intTimeScale dynamically from 1/6 to 1/15 (in code)
load('ICs/ICs_naturalSpiral.mat');

p.repStrength=0.1; % cell repulsion strength 
p.timeDepParas=1; % dynamically change this parameter
p.speedModelP1=1e-8; % [HF] & numerics
p.folderName='forPaper_2';

results=main_Dicty(p,ICs);

save(strcat('Results/forPaper_2/results.mat'), 'results','p', '-v7.3')

%% 3) cell repulsion without time dependent internal dynamics and without cell movement

load('ICs/ICs_naturalSpiral.mat');

p.repStrength=0; % cell repulsion strength -> doesn't matter, but speeds up code to set it to zero
p.timeDepParas=0;
p.speedModelP1=0; % [HF] & numerics 
p.folderName='forPaper_3';

results=main_Dicty(p,ICs);

save(strcat('Results/forPaper_3/results.mat'), 'results','p', '-v7.3')

%% 4) no cell repulsion without time dependent internal dynamics and with cell movement
load('ICs/ICs_naturalSpiral.mat');

p.repStrength=0;
p.timeDepParas=0;
p.speedModelP1=1e-8;

results=main_Dicty(p,ICs);

save(strcat('Results/forPaper_4/results.mat'), 'results','p', '-v7.3')