PhaseODE_IFAC2020.m

%% PhaseODE_IFAC2020.m
%%% MARCH 10, 2020

clear all

close all


%% Set number of cells N
N = 50; % Minimum N = 2.


%% Set network topology
net_top = 'All-to-All'; % net_top = 'All-to-All' | 'Ring' | 'Star-like'


%% 
[x0, L] = PhaseODE_initSIM(N, net_top);


D = 1; % Diffusion rate

kappa = .3; %

h = +Inf; % Hill's Coefficient.


%% Set simulation time and tSPAN
t_i = 0; % initial simulation time (min)

t_f = 4000; %final simulation time (min)


tSPAN = [t_i t_f];


%% Set variables and options for odesolver
options = odeset('RelTol', 1e-03, 'AbsTol', 1e-04, 'MaxStep', 1e-01);


%% Simulate the population of N(t) oscillators
tic
[t_out, x_out] = ode23t(@PhaseODE_NonCycling, tSPAN, x0, options, N, L, D, kappa, h);
toc


%% Retrieve the cells' phase Theta and the cells' Volume
Theta = mod(x_out(:,1:4:end), 2*pi);


%% Compute the mean phase coherence R
R = nan(length(t_out),1);
Psi = nan(length(t_out),1);

for z = 1:length(t_out)
    
    tempTheta = Theta(z,~isnan(Theta(z,:))); % Consider only those cells that exist at time t_out(z).
    
    tempN = numel(tempTheta); % Number of cells at time t_out(z).
    
    KOP = 1 ./ tempN .* sum(exp(1i .* tempTheta));
    
    R(z) = abs(KOP);
    
    Psi(z) = angle(KOP);
    
end


%% Compute the performance index J
J = 1 / (t_f - t_i) .* trapz(t_out, R);


%% Plot cells' phase
makeFIG(t_out, x_out, net_top, N, h)