function CellCulture

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                       %
% This is a companion code for the results shown in Fig.5  "Collective  %
% behavior and spatial organization of cells with different motilities" %
% from amanuscript: Kim M, Reed D, Rejniak KA. "The formation of tight  %
% tumor clusters affects the efficacy of cell cycle inhibitors:a hybrid %
% model study", published in the Journal of Theoretical Biology. 2014;  %
% 352:31-50, doi: 10.1016/j.jtbi.2014.02.027. PubMed PMID: 24607745.    %
% Author: MunJu Kim, Moffitt Cancer Center                              %                                                    
%                                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



random_force_max_factor=0;
while ~((random_force_max_factor==1)||(random_force_max_factor==10)||(random_force_max_factor==100))
    random_force_max_factor=input('Choose the motility. Type in 1, 10, or 100 and hit enter:  ');
end

randseed=1316;

filesaveon=0;
showfigures=0;
figuresave=0;
domain_size=480;

endtime=1680; % simulation time [0.1hour] == 60hrs

onecellstart=0;  %1; % will show around only one cell
initial_cell_num=10;
randomize_position=1;

r_c1i=0;
r_c2i=0;
r_w1i=0;


dir_descrip=['_CDK1_',num2str(r_c1i),'_CDK2_',num2str(r_c2i),'_WEE1_',num2str(r_w1i),'_control_Motility_',num2str(random_force_max_factor),'_random_seed_',num2str(randseed)];


edge_cell_iter=10;
radius_min=5/sqrt(2);
radius_max=5*sqrt(2);
file_save_freq=10;
dt=0.1;
Q0=0.01;

S=1;
synthesis_per_hour=1/8;
ds=dt*synthesis_per_hour;
t=0;

timestamp=clock;

if filesaveon==1
    path=['initcells_',num2str(initial_cell_num),dir_descrip];
    mkdir(path)
end

rand('seed',randseed)            % random number seed OLD STYLE

big_margin=10;
marg =10;

xsize=domain_size+big_margin;                           % xsize
ysize=domain_size+big_margin;                           % ysize

dx=10;                                           % discretization size
dy=10;                                           % discretization size

Nx=xsize/dx;
Ny=ysize/dy;

if ~((mod(Nx,1)==0)&&(mod(Ny,1)==0))
    
    disp('Discretize carefully!!')
    
end


t_iter=10;  % frequency of drawing figure

if onecellstart==0
    close all
end


cell_radius=[];   % These are alll in the .mat file
cells=[];         % These are alll in the .mat file



diam =10;                     % corresponds to ~ 10mu
rad  =diam/2;                    % cell radius
% domain marigin for periodic translocation
dt=0.1;               % in hours ~ 6 minutes
nu=1;                 % damping constant
% these num bers are not used, but are left for program consistenscy



repulse_spring=5;     % parameters for cell mechanics
k=repulse_spring;


a0=5*sqrt(2);
a_default=5;

r_a=0.19;
r_w1=0.48;


r_c1=0.48;


r_c2=0.48;


r_mf=0.48;
r_c1y=2.4;
r_c1yd=1;

r_p=0.44;


a_th=a_default*sqrt(2)*0.95;
MF_th=0.99;
CDK1Y_th=0.1;

r_bd=0.01;
r_br=0.3;
r_n=0.1;
r_h=0.4;
r_cr=1;


r_cp=1.5*r_bd;
Q_th=0.02;

k_s=0.1;



P_th=2-0.01;





% define intial cell population

initial_cell_num_idx=randi(480,initial_cell_num,1);





cells(1:initial_cell_num,3)=[23.7000
    23.9000
    2.8500
    1.9000
    17.0500
    8.4500
    14.7000
    18.2000
    22.4500
    6.1000];


cell_phase(1:initial_cell_num)=[3
    3
    1
    1
    2
    1
    2
    2
    3
    1];

cell_radius(1:initial_cell_num)=[7.0329
    7.0344
    5.0295
    4.6235
    6.9353
    6.3701
    6.8585
    6.9621
    7.0227
    5.9733];

Q(1:initial_cell_num)=[0.0157
    0.0136
    0.0179
    0.0152
    0.0358
    0.0224
    0.0277
    0.0382
    0.0387
    0.0217];

CDK1(1:initial_cell_num)=[1.0000
    1.0000
    0.8778
    0.8061
    0.9999
    0.9920
    0.9996
    0.9999
    1.0000
    0.9748];

CDK1Y(1:initial_cell_num)=[0.2129
    0.1785
    0.0965
    0.0718
    0.5876
    0.3027
    0.5052
    0.6278
    0.6815
    0.2134];

CDK1h(1:initial_cell_num)=[1.0000
    1.0000
    0.8778
    0.8061
    0.9999
    0.9920
    0.9996
    0.9999
    1.0000
    0.9748];


CDK2(1:initial_cell_num)=[1.0000
    1.0000
    0.8778
    0.8061
    0.9999
    0.9920
    0.9996
    0.9999
    1.0000
    0.9748];
WEE1(1:initial_cell_num)=[1.0000
    1.0000
    0.8778
    0.8061
    0.9999
    0.9920
    0.9996
    0.9999
    1.0000
    0.9748];

P(1:initial_cell_num)=[1.9902
    1.9902
    1.0000
    1.0000
    1.8943
    1.0000
    1.6992
    1.9366
    1.9902
    1.0000];

MF(1:initial_cell_num)=[0.9954
    0.9954
    0.6439
    0.4741
    0.9954
    0.9737
    0.9954
    0.9954
    0.9954
    0.9191];


CDK2i(1:initial_cell_num)=0;
WEE1i(1:initial_cell_num)=0;

cells(1:initial_cell_num,1:2)=domain_size*rand(initial_cell_num,2);
















num_cells=numel(cells(:,1));

% randomize cell location if 1
if randomize_position==1
    
    cells(:,1)=domain_size*rand(num_cells,1);
    cells(:,2)=domain_size*rand(num_cells,1);
    
end

cells(:,5)=zeros(num_cells,1);  % initialize cell damage


% don't want to see cell 2:4
if onecellstart==1
    cells(2:4,2)=domain_size;
end

% to see only cell 1
axislimit=[cells(1,1)-domain_size/3, cells(1,1)+domain_size/3, cells(1,2)-domain_size/3, cells(1,2)+domain_size/3];


cellgeneration=zeros(num_cells,1);      % record cell generation

dr=0.0354;
random_force_max=random_force_max_factor*dr*k;  % size of the random move
iteration_number=ceil(random_force_max_factor/2);  % handles one random move in many steps
everynstep=5;                                    % during the calculation of random move re-triangulization everynstep






xmin=0;
xmax=domain_size;
ymin=0;
ymax=domain_size;

need_to_reorganize=0; %test for re-triangulization





%%%% organize neighborhood
for organize_neighborhood=1:1
    neighbor_tri=DelaunayTri(cells(:,1),cells(:,2));
    num_cells=length(neighbor_tri.X(:,1));
    [tri_size_m,tri_size_n]=size(neighbor_tri.Triangulation);
    tetrahedron_tracker=zeros(num_cells,1);
    tetrahedron_number=zeros(num_cells,1);
    neighbor_cells=zeros(num_cells,1);
    neighbor_cells_number=zeros(num_cells,1);
    
    for idx=1:num_cells
        temp=find(neighbor_tri.Triangulation==idx);
        temp=rem(temp,tri_size_m);
        temp(find(temp==0))=tri_size_m;
        tetrahedron_tracker(idx,1:length(find(neighbor_tri.Triangulation==idx)))=temp;
        tetrahedron_number(idx)=length(temp);
        neighbor_cells_temp=neighbor_tri.Triangulation(temp(1),:);
        
        for jdx=2:tetrahedron_number(idx)
            neighbor_cells_temp=union(neighbor_cells_temp,neighbor_tri.Triangulation(temp(jdx),:));
        end
        
        neighbor_cells(idx,1:length(neighbor_cells_temp)-1)=setdiff(neighbor_cells_temp,[idx]);
        neighbor_cells_number(idx)=length(neighbor_cells_temp)-1;
    end
end
% calculate neighborhood_cells
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


ttt_count=0;
ttt_total=endtime/10;
total_grow_cell=zeros(ttt_total+1,1);
total_dying_cell=zeros(ttt_total+1,1);
total_cell=zeros(ttt_total+1,1);


if filesaveon==1
    
    fname=[path,'/zave',num2str(0),'.mat'];
    save(fname,'*')
    
end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
%%%%%  SIMULATION STARTS

for ttt=1:endtime
    
    
    distance_to_neighbors=zeros(num_cells,max(neighbor_cells_number));
    for idx=1:num_cells
        for jdx=1:neighbor_cells_number(idx)
            distance_to_neighbors(idx,jdx)=norm(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:));
        end
    end
    %Calculate distances to neighbors DONE
    
    
    
    %%%%  CELL GROWTH
    
    
    
    %growing_cells=find(cell_phase==1);
    %growing_cells=1:num_cells;
    cell_radius_new=cell_radius;
    
    %a=a+r_a*(a0-a)*dt;
    cell_radius_new=cell_radius_new+r_a*(radius_max-cell_radius_new)*dt;
    
    %%%% CELL GROWTH DONE
    
    
    %%% CELL SPRING FORCES
    %attract_force=zeros(num_cells,2); %Not In Use Yet
    for calculate_repulforces=1:1
        repulse_force=zeros(num_cells,2);
        for idx=1:num_cells
            for jdx=1:neighbor_cells_number(idx)
                %distance_to_neighbors(idx,jdx)=norm(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:));
                relax_dist=cell_radius_new(idx)+cell_radius_new(neighbor_cells(idx,jdx));
                
                if distance_to_neighbors(idx,jdx)<relax_dist
                    force_dir=(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:))/distance_to_neighbors(idx,jdx);
                    repulse_force(idx,:)=repulse_force(idx,:)+ repulse_spring*(relax_dist - distance_to_neighbors(idx,jdx))*force_dir;
                end
            end
        end
    end
    %%% CELL SPRING FORCES DONE
    
    %%% CELL MOVE DUE TO SPRING FORCE DONE
    neighbor_tri.X=neighbor_tri.X+repulse_force*dt/nu;
    cells(1:num_cells,1:2)=neighbor_tri.X;
    %%% CELL MOVE DUE TO SPRING FORCE DONE
    
    
    %%%% organize neighborhood
    for organize_neighborhood=1:1
        clear neighbor_tri
        neighbor_tri=DelaunayTri(cells(:,1),cells(:,2));
        num_cells=length(neighbor_tri.X(:,1));
        [tri_size_m,tri_size_n]=size(neighbor_tri.Triangulation);
        tetrahedron_tracker=zeros(num_cells,1);
        tetrahedron_number=zeros(num_cells,1);
        neighbor_cells=zeros(num_cells,1);
        neighbor_cells_number=zeros(num_cells,1);
        
        for idx=1:num_cells
            temp=find(neighbor_tri.Triangulation==idx);
            temp=rem(temp,tri_size_m);
            temp(find(temp==0))=tri_size_m;
            tetrahedron_tracker(idx,1:length(find(neighbor_tri.Triangulation==idx)))=temp;
            tetrahedron_number(idx)=length(temp);
            neighbor_cells_temp=neighbor_tri.Triangulation(temp(1),:);
            
            for jdx=2:tetrahedron_number(idx)
                neighbor_cells_temp=union(neighbor_cells_temp,neighbor_tri.Triangulation(temp(jdx),:));
            end
            
            neighbor_cells(idx,1:length(neighbor_cells_temp)-1)=setdiff(neighbor_cells_temp,[idx]);
            neighbor_cells_number(idx)=length(neighbor_cells_temp)-1;
        end
    end
    %done organize neighborhood
    
    
    for overlap_resolving_loop=1:edge_cell_iter
        
        is_overlap=zeros(size(neighbor_cells));
        
        clear idx
        for idx=1:num_cells
            
            for jdx=1:neighbor_cells_number(idx)
                distance_to_neighbors(idx,jdx)=norm(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:));
                if distance_to_neighbors(idx,jdx)<cell_radius_new(idx)+cell_radius_new(neighbor_cells(idx,jdx))
                    is_overlap(idx,jdx)=1;
                end
            end
        end
        clear idx
        
        is_surrounded=ones(num_cells,1);
        
        clear idx
        
        for idx=1:num_cells
            if sum(is_overlap(idx,:))<=2
                is_surrounded(idx)=0;
            end
            if sum(is_overlap(idx,:))>2
                jdx=1;
                while (is_surrounded(idx))&&(jdx<=neighbor_cells_number(idx))
                    if is_overlap(idx,jdx)==1
                        
                        clear oneside_test
                        idx_for_oneside=1;
                        % Line equation of
                        % cells(idx,1:2) and cells( neighbor_cells(idx,jdx) ,1:2)
                        if cells(idx,1)~= cells( neighbor_cells(idx,jdx) ,1)
                            
                            incline=( cells(idx,2)-cells( neighbor_cells(idx,jdx) ,2) )...
                                /(cells(idx,1)-cells( neighbor_cells(idx,jdx) ,1));
                            
                            y_intercept=-cells(idx,1)*incline+cells(idx,2);
                            
                            
                            for kdx=setdiff([1:neighbor_cells_number(idx)],jdx)
                                oneside_test(idx_for_oneside)=incline * cells( neighbor_cells(idx,kdx) ,1)...
                                    +y_intercept - cells( neighbor_cells(idx,kdx) ,2);
                                idx_for_oneside=idx_for_oneside+1;
                            end
                            
                            if max(oneside_test)*min(oneside_test)>=0
                                
                                is_surrounded(idx)=0;
                                
                            end
                            
                        elseif cells(idx,1) == cells( neighbor_cells(idx,jdx) ,1)
                            
                            
                            for kdx=setdiff([1:neighbor_cells_number(idx)],jdx)
                                
                                oneside_test(idx_for_oneside)=cells( neighbor_cells(idx,kdx) ,1) - cells(idx,1);
                                idx_for_oneside=idx_for_oneside+1;
                            end
                            
                            if max(oneside_test)*min(oneside_test)>=0
                                
                                is_surrounded(idx)=0;
                                
                            end
                            
                            
                        end
                        
                    end
                    
                    jdx=jdx+1;
                    
                end
                
            end
            
        end
        
        clear idx
        
        
        
        for calculate_repulforces=1:1
            repulse_force=zeros(num_cells,2);
            
            sum_is_overlap=sum(is_overlap,2);
            which_cell_to_push_more=find( (is_surrounded==0)&(sum_is_overlap>0) );
            
            for idx=which_cell_to_push_more'
                for jdx=1:neighbor_cells_number(idx)
                    %distance_to_neighbors(idx,jdx)=norm(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:));
                    relax_dist=cell_radius_new(idx)+cell_radius_new(neighbor_cells(idx,jdx));
                    
                    if distance_to_neighbors(idx,jdx)<relax_dist
                        force_dir=(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:))/distance_to_neighbors(idx,jdx);
                        repulse_force(idx,:)=repulse_force(idx,:)+ repulse_spring*(relax_dist - distance_to_neighbors(idx,jdx))*force_dir;
                    end
                end
            end
        end
        
        neighbor_tri.X=neighbor_tri.X+repulse_force*dt/nu;
        cells(1:num_cells,1:2)=neighbor_tri.X;
        %%% CELL MOVE DUE TO SPRING FORCE DONE
        
        
        
        %%%% organize neighborhood
        for organize_neighborhood=1:1
            clear neighbor_tri
            neighbor_tri=DelaunayTri(cells(:,1),cells(:,2));
            num_cells=length(neighbor_tri.X(:,1));
            [tri_size_m,tri_size_n]=size(neighbor_tri.Triangulation);
            tetrahedron_tracker=zeros(num_cells,1);
            tetrahedron_number=zeros(num_cells,1);
            neighbor_cells=zeros(num_cells,1);
            neighbor_cells_number=zeros(num_cells,1);
            
            for idx=1:num_cells
                temp=find(neighbor_tri.Triangulation==idx);
                temp=rem(temp,tri_size_m);
                temp(find(temp==0))=tri_size_m;
                tetrahedron_tracker(idx,1:length(find(neighbor_tri.Triangulation==idx)))=temp;
                tetrahedron_number(idx)=length(temp);
                neighbor_cells_temp=neighbor_tri.Triangulation(temp(1),:);
                
                for jdx=2:tetrahedron_number(idx)
                    neighbor_cells_temp=union(neighbor_cells_temp,neighbor_tri.Triangulation(temp(jdx),:));
                end
                
                neighbor_cells(idx,1:length(neighbor_cells_temp)-1)=setdiff(neighbor_cells_temp,[idx]);
                neighbor_cells_number(idx)=length(neighbor_cells_temp)-1;
            end
        end
        
    end
    
    
    
    %Calculate distances to neighbors
    distance_to_neighbors=zeros(num_cells,max(neighbor_cells_number));
    for idx=1:num_cells
        for jdx=1:neighbor_cells_number(idx)
            distance_to_neighbors(idx,jdx)=norm(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:));
        end
    end
    %Calculate distances to neighbors DONE
    
    
    
    
    
    %%%% CELL CYCLE DECISION MAKING
    
    cell_divide=zeros(num_cells,1);
    
    for idx=1:num_cells
        
        %Things happen regardless of the phase
        %cells age
        cells(idx,3)=cells(idx,3)+dt;
        
        %cells grow in size
        clear cell_radius_temp
        cell_radius_temp=zeros(neighbor_cells_number(idx),1);
        
        for iix=1:neighbor_cells_number(idx)
            relax_dist=cell_radius_new(idx)+cell_radius_new(neighbor_cells(idx,iix));
            if distance_to_neighbors(idx,iix)>relax_dist
                cell_radius_temp(iix)=cell_radius_new(idx);
                %cell_radius(1)
            else
                cell_radius_temp(iix)=max(distance_to_neighbors(idx,iix)*cell_radius_new(idx)/relax_dist,cell_radius(idx));
                %cell_radius(1)
            end
        end
        
        cell_radius(idx)=min(cell_radius_temp);
        
        cell_radius(idx)=min(cell_radius(idx),radius_max);
        
        
        %cells grown
        
        %KINASES
        WEE1i(idx)=WEE1i(idx)+r_w1i*(WEE1(idx)-WEE1i(idx))*dt;
        WEE1(idx)=WEE1(idx)+r_w1*(1-WEE1(idx))*dt;%-10*WEE1*dt;
        
        %CDK1h(idx)=CDK1h(idx)+r_c1*(1-CDK1(idx))*dt; %CONTROL
        CDK1h(idx)=CDK1h(idx) + r_c1*(1-CDK1(idx))*dt - r_c1i*CDK1h(idx)*dt; %CDK1 inhibition
        
        
        CDK1(idx) =CDK1(idx)+r_c1*(1-CDK1(idx))*dt;
        
        %CDK2i=CDK2i+r_c2i*(CDK2-CDK2i)*dt;
        CDK2i(idx)=CDK2i(idx)+r_c2i*(CDK2(idx)-CDK2i(idx))*dt; %CDK2 inhibition
        CDK2(idx)=CDK2(idx)+r_c2*(1-CDK2(idx))*dt;
        %CDK2=CDK2+0.48*(1-CDK2)*dt-0.4*CDK2*dt;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % phase 4: M
        % phase 3: G2/G2-checking
        % phase 2: S
        % phase 1: G1
        %              2012 Feb 14 version. Old version is at
        %              the end of codes.
        
        if cell_phase(idx)==4
            Q(idx)=Q(idx)+dt*(r_bd*P(idx)-(r_br+k_s*r_n+r_h*CDK1h(idx))*Q(idx));
            %CDK1Y=CDK1Y+10*Q*WEE1*(CDK1-CDK1Y)*dt;
            if Q(idx)<=Q_th
                cell_divide(idx)=1;
            end
        end
        
        
        if cell_phase(idx)==3
            Q(idx)=Q(idx)+dt*(0*r_bd*P(idx)-(r_br+k_s*r_n+r_h*CDK1h(idx))*Q(idx));
            CDK1Y(idx)=CDK1Y(idx)+r_c1y*Q(idx)*(WEE1(idx)-WEE1i(idx))*(CDK1(idx)-CDK1Y(idx))*dt-r_c1yd*CDK1Y(idx)*dt;
            %CDK1Y=CDK1Y-1*CDK1Y*dt;
            if CDK1Y(idx)<=CDK1Y_th
                
                cell_phase(idx)=4;
            end
            
        end
        
        if cell_phase(idx)==2 %S phase
            
            
            Pc=2*(P(idx)-1);
            Po=P(idx)-Pc;
            
            P(idx)=P(idx)+r_p*Po*dt;
            
            Q(idx)=Q(idx)+dt*(r_bd*Po+r_cp*Pc-(r_br+k_s*r_n+r_h*CDK1h(idx))*Q(idx));
            CDK1Y(idx)=CDK1Y(idx)+r_c1y*Q(idx)*(WEE1(idx)-WEE1i(idx))*(CDK1(idx)-CDK1Y(idx))*dt;
            
            if P(idx)>=P_th
                cell_phase(idx)=3;
            end
        end
        
        
        if cell_phase(idx)==1 %G1 phase
            
            MF(idx)=MF(idx)+r_mf*(CDK2(idx)-CDK2i(idx))*(1-MF(idx))*dt;
            Q(idx)=Q(idx)+dt*( r_bd- (r_br+r_n+r_h*k_s*CDK1h(idx))*Q(idx));
            CDK1Y(idx)=CDK1Y(idx)+r_c1y*Q(idx)*(WEE1(idx)-WEE1i(idx))*(CDK1(idx)-CDK1Y(idx))*dt;
            if (cell_radius(idx)>=a_th)&&(MF(idx)>=MF_th)
                cell_phase(idx)=2;
            end
            
        end
        
    end
    
    %%%% CELL CYCLE DECISION MAKING DONE
    
    
    
    %%% CALCULATE RANDOM MOVE  %%%%%
    
    for random_move_gen=1:1
        random_force=random_force_max*(2*rand(size(neighbor_tri.X))-1);
        left_edge=find(neighbor_tri.X(:,1)<xmin-marg);
        if ~isempty(left_edge)
            for iii=left_edge
                if random_force(iii,1)<0
                    random_force(iii,1)=-random_force(iii,1);
                end
            end
        end
        right_edge=find(neighbor_tri.X(:,1)>xmax+marg);
        
        if ~isempty(right_edge)
            for iii=right_edge
                if random_force(iii,1)>0
                    random_force(iii,1)=-random_force(iii,1);
                end
            end
        end
        up_edge=find(neighbor_tri.X(:,2)>ymax+marg);
        
        if ~isempty(up_edge)
            for iii=up_edge
                if random_force(iii,2)>0
                    random_force(iii,2)=-random_force(iii,2);
                end
            end
        end
        down_edge=find(neighbor_tri.X(:,2)<ymin-marg);
        
        if ~isempty(down_edge)
            for iii=down_edge
                if random_force(iii,2)<0
                    random_force(iii,2)=-random_force(iii,2);
                end
            end
        end
        
        
        for iter=1:iteration_number
            neighbor_tri.X=neighbor_tri.X+random_force*dt/nu/iteration_number;
            cells(1:num_cells,1:2)=neighbor_tri.X;
            
            %%reorganize
            if mod(iter,everynstep)==0
                for reorganize=1:1;
                    clear neighbor_tri
                    neighbor_tri=DelaunayTri(cells(:,1),cells(:,2));
                    num_cells=length(neighbor_tri.X(:,1));
                    [tri_size_m,tri_size_n]=size(neighbor_tri.Triangulation);
                    tetrahedron_tracker=zeros(num_cells,1);
                    tetrahedron_number=zeros(num_cells,1);
                    neighbor_cells=zeros(num_cells,1);
                    neighbor_cells_number=zeros(num_cells,1);
                    
                    for idx=1:num_cells
                        temp=find(neighbor_tri.Triangulation==idx);
                        temp=rem(temp,tri_size_m);
                        temp(find(temp==0))=tri_size_m;
                        tetrahedron_tracker(idx,1:length(find(neighbor_tri.Triangulation==idx)))=temp;
                        tetrahedron_number(idx)=length(temp);
                        neighbor_cells_temp=neighbor_tri.Triangulation(temp(1),:);
                        
                        for jdx=2:tetrahedron_number(idx)
                            neighbor_cells_temp=union(neighbor_cells_temp,neighbor_tri.Triangulation(temp(jdx),:));
                        end
                        
                        neighbor_cells(idx,1:length(neighbor_cells_temp)-1)=setdiff(neighbor_cells_temp,[idx]);
                        neighbor_cells_number(idx)=length(neighbor_cells_temp)-1;
                    end
                end
            end
            
            %done
            
            for idx=1:num_cells
                for jdx=1:neighbor_cells_number(idx)
                    distance_to_neighbors(idx,jdx)=norm(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:));
                end
            end
            
            for calculate_repulforces=1:1
                repulse_force=zeros(num_cells,2);
                for idx=1:num_cells
                    for jdx=1:neighbor_cells_number(idx)
                        %distance_to_neighbors(idx,jdx)=norm(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:));
                        relax_dist=cell_radius(idx)+cell_radius(neighbor_cells(idx,jdx));
                        
                        if distance_to_neighbors(idx,jdx)<relax_dist
                            force_dir=(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:))/distance_to_neighbors(idx,jdx);
                            repulse_force(idx,:)=repulse_force(idx,:)+repulse_spring*(relax_dist - distance_to_neighbors(idx,jdx))*force_dir;
                        end
                    end
                end
            end
            
            neighbor_tri.X=neighbor_tri.X+repulse_force*dt/nu/iteration_number;
            cells(1:num_cells,1:2)=neighbor_tri.X;
            
            %%reorganize
            if mod(iter,everynstep)==0
                for reorganize=1:1;
                    clear neighbor_tri
                    neighbor_tri=DelaunayTri(cells(:,1),cells(:,2));
                    num_cells=length(neighbor_tri.X(:,1));
                    [tri_size_m,tri_size_n]=size(neighbor_tri.Triangulation);
                    tetrahedron_tracker=zeros(num_cells,1);
                    tetrahedron_number=zeros(num_cells,1);
                    neighbor_cells=zeros(num_cells,1);
                    neighbor_cells_number=zeros(num_cells,1);
                    
                    for idx=1:num_cells
                        temp=find(neighbor_tri.Triangulation==idx);
                        temp=rem(temp,tri_size_m);
                        temp(find(temp==0))=tri_size_m;
                        tetrahedron_tracker(idx,1:length(find(neighbor_tri.Triangulation==idx)))=temp;
                        tetrahedron_number(idx)=length(temp);
                        neighbor_cells_temp=neighbor_tri.Triangulation(temp(1),:);
                        
                        for jdx=2:tetrahedron_number(idx)
                            neighbor_cells_temp=union(neighbor_cells_temp,neighbor_tri.Triangulation(temp(jdx),:));
                        end
                        
                        neighbor_cells(idx,1:length(neighbor_cells_temp)-1)=setdiff(neighbor_cells_temp,[idx]);
                        neighbor_cells_number(idx)=length(neighbor_cells_temp)-1;
                    end
                end
            end
            
            %done
            
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    if onecellstart==1
        cell_divide(2:4)=0;
    end
    
    
    
    %%%%%%  CELLS DIVIDE
    
    which_cell_divide=find(cell_divide==1);
    
    for idx=which_cell_divide'
        
        cell_radius(idx)=cell_radius(idx)/2; %size a
        cells(idx,3)=0; %age
        cellgeneration(idx)=cellgeneration(idx)+1;
        cell_phase(idx)=1;
        P(idx)=1;
        Q(idx)=Q(idx)/2;
        CDK1(idx)=CDK1(idx)/2;
        CDK2(idx)=CDK2(idx)/2;
        CDK1Y(idx)=CDK1Y(idx)/2;
        CDK1h(idx)=CDK1h(idx)/2;
        CDK2i(idx)=CDK2i(idx)/2;
        WEE1(idx)=WEE1(idx)/2;
        WEE1i(idx)=WEE1i(idx)/2;
        
        MF(idx)=0;
        cell_divide(idx)=0;
        
        %growth_rate_old=growth_rate(idx);
        %growth_rate(idx)=growth_rate(idx)*(erlot_factor+(1-erlot_factor)*rand(1)^rpower);
        
        new_cell=neighbor_tri.X(idx,1:2);
        
        split=(2*rand(1,2)-1);
        while norm(split)<10^(-5)
            split=(2*rand(1,2)-1);
        end
        
        split=0.4*diam*split/norm(split);
        
        new_cell=new_cell+split;
        neighbor_tri.X(idx,1:2)=neighbor_tri.X(idx,1:2)-split;
        
        num_cells=num_cells+1;
        
        neighbor_tri.X(num_cells,1:2)=new_cell;
        
        
        cells(num_cells,1:2)=new_cell;
        cells(num_cells,3)=0;
        cells(num_cells,4)=cells(idx,4);
        cell_phase(num_cells)=1;
        cell_radius(num_cells)=cell_radius(idx);
        Q(num_cells)=Q(idx);
        Q(idx)=Q(num_cells);
        P(num_cells)=1;
        CDK1(num_cells)=CDK1(idx);
        CDK2(num_cells)=CDK2(idx);
        CDK2i(num_cells)=CDK2i(idx);
        CDK1Y(num_cells)=CDK1Y(idx);
        CDK1h(num_cells)=CDK1h(idx);
        
        WEE1(num_cells)=WEE1(idx);
        WEE1i(num_cells)=WEE1i(idx);
        MF(num_cells)=0;
        cell_divide(num_cells)=0;
        
        
        cellgeneration(num_cells)=cellgeneration(idx);
        
        need_to_reorganize=1;
        %growth_rate(num_cells)=growth_rate_old*(erlot_factor+(1-erlot_factor)*rand(1)^rpower);
    end
    
    
    cells(1:num_cells,1:2)=neighbor_tri.X(1:num_cells,1:2);
    
    %%put them back in
    
    
    for push_iter=1:10;
        
        force_pushin=zeros(num_cells,2);
        
        outsider=find((cells(:,1)<xmin-marg)|(cells(:,1)>xmax+marg)|(cells(:,2)<ymin-marg)|(cells(:,2)>ymax+marg) );
        
        for cell_idx=outsider'
            
            if (cells(cell_idx,1)<xmin-marg)&&(cells(cell_idx,2)<ymin-marg)
                
                force_direction_push=[domain_size/2,domain_size/2]-cells(cell_idx,1:2);
                force_direction_push=force_direction_push/norm(force_direction_push);
                
                force_pushin(cell_idx,1:2)=(nu/2/dt)*max( abs( cells(cell_idx,1)-xmin+marg ), abs(cells(cell_idx,2)-ymin+marg) )*force_direction_push;
                
            elseif ( cells(cell_idx,1)<xmin-marg )&&( cells(cell_idx,2)>=ymin-marg )&&( cells(cell_idx,2)<=ymax+marg )
                
                force_direction_push=[1,0];
                
                force_pushin(cell_idx,1:2)=(nu/2/dt)*abs(cells(cell_idx,1)- xmin+marg )*force_direction_push;
                
            elseif ( cells(cell_idx,1)<xmin-marg )&&( cells(cell_idx,2)>ymax+marg )
                
                force_direction_push=[domain_size/2,domain_size/2]-cells(cell_idx,1:2);
                force_direction_push=force_direction_push/norm(force_direction_push);
                
                force_pushin(cell_idx,1:2)=(nu/2/dt)*max( abs( cells(cell_idx,1)-xmin+marg ), abs(cells(cell_idx,2)-ymax-marg) )*force_direction_push;
                
            elseif ( cells(cell_idx,1)>=xmin-marg )&&( cells(cell_idx,1)<=xmax+marg )&&( cells(cell_idx,2)<ymin-marg )
                
                force_direction_push=[0,1];
                
                force_pushin(cell_idx,1:2)=(nu/2/dt)*abs(cells(cell_idx,2)- ymin + marg )*force_direction_push;
                
            elseif ( cells(cell_idx,1)>=xmin-marg )&&( cells(cell_idx,1)<=xmax+marg )&&( cells(cell_idx,2)>ymax+marg )
                
                force_direction_push=[0,-1];
                
                force_pushin(cell_idx,1:2)=(nu/2/dt)*abs(cells(cell_idx,2)- ymax - marg )*force_direction_push;
                
            elseif ( cells(cell_idx,1)>xmax+marg )&&( cells(cell_idx,2)<ymin-marg )
                
                force_direction_push=[domain_size/2,domain_size/2]-cells(cell_idx,1:2);
                force_direction_push=force_direction_push/norm(force_direction_push);
                
                force_pushin(cell_idx,1:2)=(nu/2/dt)*max( abs( cells(cell_idx,1)-xmax-marg ), abs(cells(cell_idx,2)-ymin+marg) )*force_direction_push;
                
            elseif ( cells(cell_idx,1)>xmax+marg )&&( cells(cell_idx,2)>ymax+marg )
                
                force_direction_push=[domain_size/2,domain_size/2]-cells(cell_idx,1:2);
                force_direction_push=force_direction_push/norm(force_direction_push);
                
                force_pushin(cell_idx,1:2)=(nu/2/dt)*max( abs( cells(cell_idx,1)-xmax-marg ), abs(cells(cell_idx,2)-ymax -marg) )*force_direction_push;
                
            elseif ( cells(cell_idx,1)>xmax+marg )&&( cells(cell_idx,2)<=ymax+marg )&&( cells(cell_idx,2)>=ymin-marg )
                
                force_direction_push=[-1,0];
                
                force_pushin(cell_idx,1:2)=(nu/2/dt)*abs(cells(cell_idx,1)- xmax - marg )*force_direction_push;
                
            end
            
            
        end
        
        cells(1:num_cells,1:2)=cells(1:num_cells,1:2) + force_pushin*dt/nu;
        
        
        for organize_neighborhood=1:1
            clear neighbor_tri
            neighbor_tri=DelaunayTri(cells(:,1),cells(:,2));
            num_cells=length(neighbor_tri.X(:,1));
            [tri_size_m,tri_size_n]=size(neighbor_tri.Triangulation);
            tetrahedron_tracker=zeros(num_cells,1);
            tetrahedron_number=zeros(num_cells,1);
            neighbor_cells=zeros(num_cells,1);
            neighbor_cells_number=zeros(num_cells,1);
            
            for idx=1:num_cells
                temp=find(neighbor_tri.Triangulation==idx);
                temp=rem(temp,tri_size_m);
                temp(find(temp==0))=tri_size_m;
                tetrahedron_tracker(idx,1:length(find(neighbor_tri.Triangulation==idx)))=temp;
                tetrahedron_number(idx)=length(temp);
                neighbor_cells_temp=neighbor_tri.Triangulation(temp(1),:);
                
                for jdx=2:tetrahedron_number(idx)
                    neighbor_cells_temp=union(neighbor_cells_temp,neighbor_tri.Triangulation(temp(jdx),:));
                end
                
                neighbor_cells(idx,1:length(neighbor_cells_temp)-1)=setdiff(neighbor_cells_temp,[idx]);
                neighbor_cells_number(idx)=length(neighbor_cells_temp)-1;
            end
        end
        
        for idx=1:num_cells
            for jdx=1:neighbor_cells_number(idx)
                distance_to_neighbors(idx,jdx)=norm(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:));
            end
        end
        
        for calculate_repulforces=1:1
            repulse_force=zeros(num_cells,2);
            for idx=1:num_cells
                for jdx=1:neighbor_cells_number(idx)
                    %distance_to_neighbors(idx,jdx)=norm(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:));
                    relax_dist=cell_radius(idx)+cell_radius(neighbor_cells(idx,jdx));
                    
                    if distance_to_neighbors(idx,jdx)<relax_dist
                        force_dir=(neighbor_tri.X(idx,:)-neighbor_tri.X(neighbor_cells(idx,jdx),:))/distance_to_neighbors(idx,jdx);
                        repulse_force(idx,:)=repulse_force(idx,:)+repulse_spring*(relax_dist - distance_to_neighbors(idx,jdx))*force_dir;
                    end
                end
            end
        end
        
        %%% CELL SPRING FORCES DONE
        
        %%% CELL MOVE DUE TO SPRING FORCE DONE
        neighbor_tri.X=neighbor_tri.X+repulse_force*dt/nu;
        cells(1:num_cells,1:2)=neighbor_tri.X;
        
    end
    
    
    %%%%%%  CELLS DIVIDE DONE
    
    
    %%%% organize neighborhood
    for organize_neighborhood=1:1
        if need_to_reorganize==1;
            clear neighbor_tri
            neighbor_tri=DelaunayTri(cells(:,1),cells(:,2));
            num_cells=length(neighbor_tri.X(:,1));
            [tri_size_m,tri_size_n]=size(neighbor_tri.Triangulation);
            tetrahedron_tracker=zeros(num_cells,1);
            tetrahedron_number=zeros(num_cells,1);
            neighbor_cells=zeros(num_cells,1);
            neighbor_cells_number=zeros(num_cells,1);
            
            for idx=1:num_cells
                temp=find(neighbor_tri.Triangulation==idx);
                temp=rem(temp,tri_size_m);
                temp(find(temp==0))=tri_size_m;
                tetrahedron_tracker(idx,1:length(find(neighbor_tri.Triangulation==idx)))=temp;
                tetrahedron_number(idx)=length(temp);
                neighbor_cells_temp=neighbor_tri.Triangulation(temp(1),:);
                
                for jdx=2:tetrahedron_number(idx)
                    neighbor_cells_temp=union(neighbor_cells_temp,neighbor_tri.Triangulation(temp(jdx),:));
                end
                
                neighbor_cells(idx,1:length(neighbor_cells_temp)-1)=setdiff(neighbor_cells_temp,[idx]);
                neighbor_cells_number(idx)=length(neighbor_cells_temp)-1;
            end
            
            need_to_reorganize=0;
        end
    end
    %%%% organize neighborhood DONE
    
    
    rand_seed_save=rand('seed');
    
    if filesaveon==1
        if mod(ttt,file_save_freq)==0
            fname=[path,'/zave',num2str(ttt),'.mat'];
            save(fname,'*')
        end
        
        %pause(0.0005)
        
    end
    
%          close all
%         figure
%         set(gcf,'position',[1 1 1000 1000]) 
    if mod(ttt,1)==0
        
%         close all
%         figure
         set(gcf,'position',[1 1 1000 1000])
        
        
        for iii=1:num_cells
            [xx,yy,zz]=cylinder([0,cell_radius(iii)]);
            
            if cell_phase(iii)==1
                facecolorcode=[1,0,0];
            elseif cell_phase(iii)==2
                facecolorcode=[0,1,0];
            else
                facecolorcode=[0,0,1];%[191/255,	62/255,	255/255];
            end
            
            if (cell_phase(iii)==1)&&(cells(iii,3)>24)
                facecolorcode=[191/255,	62/255,	255/255];
            end
            
            
            surf(xx+neighbor_tri.X(iii,1),yy+neighbor_tri.X(iii,2),zz*0-0.01,'facecolor',facecolorcode,'edgecolor','none')
            hold on
 
            
            
        end
        view(2)
        set(gca,'fontsize',24)
        set(gcf,'color','w')
        
        axis equal
        axis([-20  500  -20  500])
        set(gca,'xtick',-200:1000:1000)
        set(gca,'ytick',-200:1000:1000)
        grid off
        xlabel('http://dx.doi.org/10.1016/j.jtbi.2014.02.027')
        %axis tight
        
        set(gca,'fontsize',24,'fontname','times new roman')
        set(gcf,'color','w')
        titlephrase=['time: ',num2str(ttt*dt),' hrs; number of cells: ',num2str(num_cells),';  Red:G1, Green:S, Blue:G2/M'];
        title(titlephrase)
        box on
        
        drawnow
        
        
        hold off
        
        pause(0.01)
        
        
        
        
        
    end
    
    
    
end  %time for loop


