function modelWhAM(IsResistance,KindOfResistance,PreExist_Location,...
                   delta_death_rate,p_dna_repair_rate)      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                         %% 
%% A Research Collaboration Workshop for Women in Applied Mathematics,     %%
%% Dynamical Systems with Applications to Biology and Medicine.            %%
%% IMA-WhAM!-2013-2015 Research Project:                                   %%
%% Anti-Cancer Drug Resistance: A Pre-Exisitng or Emerging Phenomenon?     %%
%% published as Chapter I of "Applications of Dynamical Systems in Biology %%
%% and Medicine", edited by A. Radunskaya and T. Jackson, IMA Volumes in   %%
%% Mathematics and Its Applications, Springer-Verlag, 2015,                %%
%% ISBN 978-1-4939-2781-4;                                                 %%
%%                                                                         %% 
%% Katarzyna Rejniak - Moffitt Cancer Center & Research Institute, IMO     %%
%% Jana Gevertz   - The College of New Jersey, Dep. of Math & Stats        %%
%% Kerri-Ann Norton - Johns Hopkins Uinversity                             %%
%% Judith Perez-Velazquez - Helmholtz Zentrum, Muenchen, Germany           %%
%% Zahra Aminzare - Rutgers University                                     %%
%% Alexandria Volkening - Brown University                                 %%
%%                                                                         %% 
%%                                                                         %% 
%%   parameters:                                                           %% 
%% IsResistance      : 0=no resistance, 1=resistance to drug level,        %%
%% KindOfResistance  : 0=acquired resistance; 1=pre-existing resistance,   %%
%% PreExist_Location : 1=two cells closest to vessel are resistant,        %% 
%%                     2=two cells at intermediate distance are resistant, %%
%%                     3=two cells furthest from vessel are resistant,     %%
%% delta_death_rate  : rate of increase of death threshold,                %%
%% p_dna_repair_rate : DNA repair rate,                                    %%
%%                                                                         %% 
%% for acquired   : change delta_death_rate,   keep p_dna_repair_rate=0    %%
%% for pre-exising: keep   delta_death_rate=0, change p_dna_repair_rate    %%
%% to reconstruct Figure 1 simulations use:                                %%
%% Figure 1c modelWhAM(0,0,0,0,0.00015)                                    %%
%% Figure 1e modelWhAM(1,1,2,0,0.00015)                                    %%
%% Figure 1f modelWhAM(1,0,0,0.000059,0.00015)                             %%
%% Figure 1g modelWhAM(1,0,0,0.000045,0.00015)                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Table 1 parameters
CellDiam       =  5;         % cell diameter [microns]
F_spr          =  1;         % spring stiffness    
nu             = 15;         % mass/viscosity term 
neighborMax    = 14;         % number of neighbors when cell is overcrowded
maturConst     =360;         % age at which the cell will divide [min] 
xmin=-75;  xmax= 75;         % domain [microns]
ymin=xmin; ymax=xmax; 
hb             =  2;         % grid width [microns]
dt             =  0.5;       % time step [min]


% Table 2 parameters
Source_oxy =1;               % influx of oxygen from the vessels
DiffOxy    =0.5;             % oxygen diffusion coefficient
rho_oxy    =5e-5*Source_oxy; % oxygen uptake rate
Thr_hypo   =0.05*Source_oxy; % hypoxia threshold

Source_drug=1;               % influx of drug from each vessel
DiffDrug   =0.5;             % drug diffusion coefficient  
d_drug     =1e-4;            % decay of drug 
rho_drug   =1e-4*Source_drug;% dose taken up by a single cell 
Thr_death  =0.5;             % threshold for cell death

omega      =0.45;            % boundary outflux rate

drug_exp   =0.01;  % value of drug concentration to increment cell_death
time_exp   =5*dt;  % amount of time of drug exposure to increment cell_death



% Switches for drawing the images and saving the data during a simualtion
% ShowImages        : 0=don't show images; 1=show images,                  
% SaveData          : 0=don't save data; 1=save data.                      
ShowImages = 1;    % the default option is to show the iamges
SaveData   = 1;    % the default option is to do not save the data 


% Define iteration parameters 
Niter      =25000;           % max number of iterations
Nsave      =  100;           % save every Nsave iteration    
Ndraw      =  500;           % draw every Ndraw iteration


% Define domain parameters 
Ngx=1+floor((xmax-xmin)/hb); % number of grid points - x axis
Ngy=1+floor((ymax-ymin)/hb); % number of grid points - y axis
xgg=xmin:hb:xmax;            % data for drawing
ygg=ymin:hb:ymax;            % data for drawing



% create a directory to save all data
if (SaveData==1)
  path=MakeDir; 
end


% Fixed random seed for comparison between runs with differnet parameters
seed=2; rng(seed); fprintf('\tUsing seed = %d\n',seed);


% Define cell properties: 
% page 7: [x,y,age,maturation,oxy,drug,exposure,damage,death,IDcell,IDmother]
[cell_xy,Ncells]=DefineCellPopulation;          % cell coordinates 
cell_age(1:Ncells,1)=rand(Ncells,1)*maturConst; % randomized age between 
for ii = 1:Ncells                               %            [0,maturConst] 
   cell_mature(ii,1)=maturConst*(0.5 + rand(1));% randomized maturation age
end
cell_oxy   =zeros(Ncells,1);         % oxygen sensed by the cells
cell_drug  =zeros(Ncells,1);         % drug accumulated inside the cells
cell_exp   =zeros(Ncells,1);         % time counter for drug exposure 
cell_damage=zeros(Ncells,1);         % damage of the cell (induced by drug)
cell_death=Thr_death*ones(Ncells,1); % threshold of DNA damage for cell death
cell_ID    =[1:1:Ncells;1:1:Ncells]; % cell indices [selfID, motherID]
eps        =0.5*CellDiam;            % distance for placing a new cell [microns] 


% Define vessel properites
[vessel,Nvessel]=DefineVesselPopulation; % vessel locations 
VessRad         =2;                      % vessel radius


% Define oxygen and drug concentration in teh domain
oxyDom=load('oxyinit10.txt');       % read the initial oxygen gradient 
drugDom=zeros(Ngx,Ngy);             % drug distribition inside the domain


% Stability conditions for oxygen and drug
if ((DiffOxy*dt/(hb*hb))>0.5)
  disp('stability condition for oxy is violated'); pause
end
if ((DiffDrug*dt/(hb*hb))>0.5)
  disp('stability condition for drug is violated'); pause
end


% indeces of mother cells, unique cell index and current number of cells
cellsMotherID = 1:1:Ncells;  % motherID for each cell (tracking cell clones)
cellMaxlID=Ncells;           % number of all generated cells (to generate a unique cell ID)
current_num_cells = [];      % for storing number of cancer cells in each iteration



% parameters depending on the chosen kind of resistance
if (KindOfResistance==0)
  Thr_multi    = 1;                  % standard tolerance
  delta_death  = delta_death_rate;   % cell tolerance varied by the user
  p_dna_repair = 0.00015;            % fixed DNA repair rate
elseif (KindOfResistance==1)
  Thr_multi    = 5;                  % increased tolerance for pre-existing
  delta_death  = 0;                  % no changes in cell tolerance 
  p_dna_repair = p_dna_repair_rate;  % DNA repair rate varied by the user
  cell_death   = ChooseResistantCellLocation(cell_xy,vessel,Ncells,...
    PreExist_Location,cell_death,Thr_multi); % chose cells with pre-existing resistance
else  
  disp('wrong Kind of Resistance parameter')
  pause
end



% initiation of imaging
if(ShowImages==1)  
  fg1=figure('position',[ 50,150,650,650]); % image intiation
  fg2=figure('position',[750,150,650,650]);
  pause(0.1)
end

Nclones=Ncells;

   
%%% save all model parameters; separately real and integer
if (SaveData==1)
  paramReal=[xmin,xmax,ymin,ymax,hb,CellDiam,eps,VessRad,dt,F_spr,nu,...
     Source_drug,DiffDrug,rho_drug,d_drug,omega,Thr_death,drug_exp,...
     time_exp,delta_death,Thr_multi,Source_oxy,DiffOxy,rho_oxy,...
     Thr_hypo,p_dna_repair];
  paramInt=[Ngx,Ngy,Nsave,Nvessel,Ncells,Nclones,maturConst,neighborMax,...
     Niter,cellMaxlID,PreExist_Location];
  fname=[path,'/paramReal.txt']; save(fname,'paramReal','-ascii')
  fname=[path,'/paramInt.txt'];  save(fname,'paramInt','-ascii')
end


   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% main loop of the program %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

iter=-1; % Start at -1 so first iteration is numbered 0, not 1
while ((iter<=Niter)&&(Ncells>0))% finish if Niter reached or all cells are dead

   iter=iter+1;                         % increase loop counter
   current_num_cells(iter+1) = Ncells;  % current number of cancer cells

   if(mod(iter,500)==0)          % display current information about the cells
      fprintf('\tAt iteration number %d; cells: %d\n',iter,Ncells);
   end
   
   
   % define oxygen supply from the vessels & oxygen diffusion in teh domain
   oxyDom=Supply(oxyDom,vessel,Nvessel,Source_oxy,VessRad,dt,hb,xmin,ymin);   
   oxyDom=DiffusionBound(oxyDom,Ngx,Ngy,DiffOxy,omega,dt,hb);                

   % define drug supply form the vessel, diffusion in the domain and decay
   drugDom=Supply(drugDom,vessel,Nvessel,Source_drug,VessRad,dt,hb,xmin,ymin);
   drugDom=DiffusionBound(drugDom,Ngx,Ngy,DiffDrug,omega,dt,hb);
   drugDom=Decay(drugDom,Ngx,Ngy,d_drug,dt);                    
    
   
   % update cellular properties
   numNeighbors = Overcrowding(cell_xy,2*CellDiam); % Search for cell neighbors 
   cell_oxy=zeros(Ncells,1);     % cells do not accumulate oxygen
   
   perm_indx = randperm(Ncells); % random permutation of cell indices to avoid bias 
   for ii=1:Ncells          
     iperm = perm_indx(ii);      % use a permuted cell index 
    
     
     % Determine if cell has experienced prolong exposure to drug that 
     % will increase its death threshold.    
     if (IsResistance==1)                 
       cell_exp(iperm)=TimeDrugExposure(drug_exp,iperm,cell_drug,cell_exp,dt);
       if (cell_exp(iperm)>=time_exp)    % Exposure criterion met
         cell_death(iperm,1)=cell_death(iperm,1)+delta_death; % increase death threshold
       end 
     end
     
     
     % Determine cellular uptake of oxygen and drug from grid points inside the cell   
     % the amount of xygen and drug in these points will be dimished by the value
     % sensed by the cell    
     ix=1+floor((cell_xy(iperm,1)-xmin)/hb); % grid point near the cell - x
     iy=1+floor((cell_xy(iperm,2)-ymin)/hb); % grid point near the cell - y
          
     drug_decay_in_cell=d_drug*cell_drug(iperm,1)*dt;  % decay of accumulated drug
     
     drug_uptake  =0;                        % determine total drug uptake by a cell
     oxy_uptakeNum=0;                        % determine total number of grid points the cell is sensing fron 
     
     search = 2;                             % how far the cell can sense the space
     for k=ix-search:1:ix+search
       for l = iy-search:1:iy+search           
         if (k>0)&&(l>0)&&(k<Ngx)&&(l<Ngy)     % is the grid point in inside the domain
           xdist = cell_xy(iperm,1)-(xmin+hb*k); % distanace from (k,l) to
           ydist = cell_xy(iperm,2)-(ymin+hb*l); %    ((cells(iperm,1),cells(iperm,2))
           distance = sqrt(xdist*xdist+ydist*ydist);
           if (distance<0.5*CellDiam)          % is the grid point inside the cell 
             % drug uptake 
             cell_drug(iperm,1)=cell_drug(iperm,1)+dt*min(drugDom(k,l),rho_drug); % drug uptake by the cell
             drug_uptake=drug_uptake+dt*min(drugDom(k,l),rho_drug);  % remember how much is taken up
             drugDom(k,l)=drugDom(k,l)-dt*min(drugDom(k,l),rho_drug);% drug taken from the grid point 

             %oxygen sesing & uptake
             cell_oxy(iperm,1)=cell_oxy(iperm,1)+oxyDom(k,l);% oxygen sensing/uptake by the cell
             oxy_uptakeNum=oxy_uptakeNum+1;    % remeber how many grid point the cell senses  
             oxyDom(k,l)=oxyDom(k,l)-dt*min(oxyDom(k,l),rho_oxy);    % oxygen taken from the grid point       
           end
         end % if (k>0)&&(l>0)&&(k<Ngx)&&(l<Ngy) 
       end % l = iy-search:1:iy+search
     end % k=ix-search:1:ix+search
      
     cell_drug(iperm,1)=cell_drug(iperm,1)-drug_decay_in_cell;% loss of drug in cell due to decay

     drug_uptake_decay = drug_uptake-drug_decay_in_cell;      % calculate damage in a cell which depends on current drug uptake  
     cell_damage(iperm) = cell_damage(iperm)+max(0,drug_uptake_decay)...
                        - p_dna_repair*cell_damage(iperm); % Update damage 

     if (oxy_uptakeNum>0)                      % calculate average amount of oxygen around the cell
       cell_oxy(iperm,1)=cell_oxy(iperm,1)/oxy_uptakeNum;
     end

          
    
     % Cell divides if old enough: age is reset to 0, matur is inherited
     % Cell division also only occurs if cell is not too crowded
     if ((cell_age(iperm,1)>=cell_mature(iperm,1))&&...
         (numNeighbors(1,iperm)<neighborMax))&&(cell_oxy(iperm,1)>=Thr_hypo) 
          
        motherID=cell_ID(1,iperm);  % cell ID become cell mother ID   
         
        
        % first daughter cell
        cell_age(iperm,1)=0;        % age of the first daughter cell set to 0
        cell_drug(iperm,1)=0.5*cell_drug(iperm,1); % drug concentration is set to half
        % note: position, damage, oxy, death and exp the same as mother cell         
         
        cellMaxlID=cellMaxlID+1;    % new unique index for a daughter cell
        cell_ID(1:2,iperm)=[cellMaxlID,motherID]; % remeber ID of the cell and its mother    
        cellsMotherID(cellMaxlID)=motherID;       % remeber ID of the cell and its mother

        
          
        % second daughter cell
        Ncells=Ncells+1;            % add new daughter cell
        cell_age(Ncells,1)=0;       % age of the second daughter cell
       
        cellMaxlID=cellMaxlID+1;    % new unique index for a daughter cell
        cell_ID(1:2,Ncells)=[cellMaxlID,motherID]; % remeber ID of the cell and its mother  
        cellsMotherID(cellMaxlID)=motherID;        % remeber ID of the cell and its mother

        cell_mature(Ncells,1)=cell_mature(iperm,1)+(maturConst/10)*(rand-0.5); 
        % daugher inherits meturation time of the mother with gentle noise term
        % new age in range [19/20*maturConst, 21/20*maturConst]
        cell_drug(Ncells,1)=cell_drug(iperm,1); % drug concentration is set to half
        ang=randi(62);             % random angle to place a new cell
        cell_xy(Ncells,1:2)=[cell_xy(iperm,1)+eps*cos(2*pi/(ang*0.1)),...
        cell_xy(iperm,2)+eps*sin(2*pi/(ang*0.1))];
     
        cell_damage(Ncells,1)=cell_damage(iperm,1); % cell damage is inhetited           
        cell_oxy(Ncells,1)=cell_oxy(iperm,1);       % cell oxy is inherited 
        cell_death(Ncells,1) = cell_death(iperm,1); % dathe threshold is inherited
        cell_exp(Ncells,1) = cell_exp(iperm,1);     % exposure is inherited

     else
        cell_age(iperm,1)=cell_age(iperm,1)+dt;     % cells get older if it didn't just divide   
     end   % if           
   
   end  % for  ii=1:Ncells
    
    
   % algorithm for cell death -- alive cells will be copied to new variables
   Ncells_new=0; cell_xy_new=zeros(1,2); cell_age_new=zeros(1,1); 
   cell_mature_new=zeros(1,1); cell_oxy_new=zeros(1,1); cell_drug_new=zeros(1,1);   
   cell_exp_new=zeros(1,1); cell_damage_new=zeros(1,1); cell_death_new=zeros(1,1); 
   cell_ID_new=zeros(2,1); 
   
   for ii=1:Ncells       
     % cell lives if damage level below its death threshold & cell inside the domain 
     if (cell_damage(ii,1)<cell_death(ii,1))&&...   
        (cell_xy(ii,1)>(xmin+eps))&&(cell_xy(ii,1)<(xmax-eps))&&...
        (cell_xy(ii,2)>(ymin+eps))&&(cell_xy(ii,2)<(ymax-eps))
                
       Ncells_new=Ncells_new+1;   % new alive cell added
       cell_xy_new    (Ncells_new,1:2)=cell_xy(ii,1:2);
       cell_age_new   (Ncells_new,1)  =cell_age(ii,1);
       cell_mature_new(Ncells_new,1)  =cell_mature(ii,1);
       cell_oxy_new   (Ncells_new,1)  =cell_oxy(ii,1);
       cell_drug_new  (Ncells_new,1)  =cell_drug(ii,1);
       cell_exp_new   (Ncells_new,1)  =cell_exp(ii,1); 
       cell_damage_new(Ncells_new,1)  =cell_damage(ii,1);
       cell_death_new (Ncells_new,1)  =cell_death(ii,1);
       cell_ID_new    (1:2,Ncells_new)=cell_ID(1:2,ii);
     end % if                     % if not, omit the cell, cell will be killed
   end % for ii=1:Ncells
   
   % update the old variables
   Ncells=Ncells_new; cell_xy=cell_xy_new; cell_age=cell_age_new;
   cell_mature=cell_mature_new; cell_oxy=cell_oxy_new; cell_drug=cell_drug_new; 
   cell_exp=cell_exp_new; cell_damage=cell_damage_new; cell_death=cell_death_new;
   cell_ID= cell_ID_new;   
   clear cell_xy_new cell_age_new cell_mature_new cell_oxy_new cell_drug_new 
   clear cell_exp_new cell_damage_new cell_death_new cell_ID_new
   
   
   
   % determine forces between neighboring cells to push them apart if necessary
   force=DefineForce(cell_xy,CellDiam,F_spr);
   for ii=1:Ncells     % inpect all cells
     cell_xy(ii,1)=cell_xy(ii,1)+(force(ii,1)*dt/nu);  % new cell positions
     cell_xy(ii,2)=cell_xy(ii,2)+(force(ii,2)*dt/nu);  % due to repulsive fores
   end % for
    


   
   % draw data if multiplicity of Ndraw is reached & ShowImages==1
   if (mod(iter,Ndraw)==0)&&(ShowImages==1)
     DrawTissue(fg1,fg2,cell_xy,cell_oxy,cell_drug,cell_damage,cell_death,...
     cell_ID,Ncells,Nclones,vessel,cellsMotherID,oxyDom,drugDom,xgg,ygg,...
     Source_oxy,Source_drug,xmin,xmax,ymin,ymax,Thr_hypo,dt,iter)
     pause(0.1)
     
     if (SaveData==1)
       figure(fg1)
       fname=[path,'/Fig1_',num2str(iter)];
       print('-djpeg100',fname)
       pause(0.1)

       figure(fg2)
       fname=[path,'/Fig2_',num2str(iter)];
       print('-djpeg100',fname)
       pause(0.1)
     end
   end % if (ShowImages==1)
  
   
   
   % save data if multiplicity of Nsave is reached or all cells are dead
   if ((mod(iter,Nsave)==0)||(Ncells==0))
     if (SaveData==1) 
       SaveAllData(path,cell_xy,cell_age,cell_mature,cell_oxy,cell_drug,...
       cell_exp,cell_damage,cell_death,cell_ID,cellsMotherID,vessel,...
       oxyDom,drugDom,iter)
     end
   end
   
   
   % save data once simulation is done (reaches Niter or no more cells)
   if ((iter == Niter) || (Ncells == 0))
     if (SaveData==1)  
       SaveAllData(path,cell_xy,cell_age,cell_mature,cell_oxy,cell_drug,...
       cell_exp,cell_damage,cell_death,cell_ID,cellsMotherID,vessel,oxyDom,...
       drugDom,iter)
       %save number of cells as a function of time
       fname=[path,'/number_cancer_cells.txt']; save(fname,'current_num_cells','-ascii')
     end
   end 

   
end % main loop while



% display information about tumor death and at which iteration
if(Ncells==0)
    fprintf('Tumor died out after %d iterations\n',iter);
end


% plot tumor size versus time here
if(ShowImages==1)  
  figure; 
  plot(current_num_cells); 
  xlabel('Time','FontSize',12); 
  ylabel('Number of cancer cells','FontSize',12); 
  title('Analyzing tumor size','FontSize',14,'Interpreter','none');

  if(SaveData==1)  
    fname=[path,'/Cancer_size'];
    print('-djpeg100',fname)
  end
  pause(0.1)     
end


end % main program   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function defining coordinates of initial cell population            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[cell_xy,Ncells]=DefineCellPopulation
         
   cell_xy=[0,0;5,5;-5,2;-1,7;-5,6;-10,5;-11,-2;-7,-6;-3,-7;3,-7;7,-2;...
      11,1;9,6;4,8;2,11;3,-1;0,-4;-4,-2;-9,1;-13,2;-13,7;-8,9;-3,12;...
      1,15;5,13;7,11;12,9;13,4;15,0;11,-4;6,-8;2,-10;-4,-10;-10,-10;...
      -14,-3;-11.9366,10.9859;-12.1479,12.4648;-18.9085,4.4366;-18.2746,...
      13.3099;-22.7113,10.1408;-25.4577,2.3239;-17.0070,21.3380;-22.5000,...
      17.5352;-30.1056,13.9437;-31.5845,21.1268;-36.0211,12.8873;-41.5141,...
      10.5634;-37.7113,5.0704;-33.6972,-1.6901;-27.5704,-3.3803;-47.4296,...
      16.4789;-48.4859,24.9296;-41.5141,26.4085;-32.8521,18.1690;-28.2042,...
      42.4648;-39.4014,39.0845;-55.4577,33.8028;-55.8803,26.8310;-55.0352,...
      16.4789;-47.8521,13.3099;-50.1761,30.6338;-42.3592,28.7324;-30.7394,...
      5.0704;-27.5704,7.8169;-40.4577,0.4225];
   Ncells=size(cell_xy,1);
   
%% This part of the code is not used in the simulations, but is left here
%% in order to allow reproduction of the results from our paper. Without 
%% this part using a random number generator the order of ramdomly generated 
%% numbers will not agree with simulations from teh paper
   randStateOrder = randperm(Ncells);
   clear randStateOrder
end % function DefineCellPopulation


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function defining locations of vessels                              % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  [vessel,Nvessel]=DefineVesselPopulation
  vessel=[-20,-40;-40,20;20,-20;60,60];
  Nvessel=4;
end % function DefineVesselPopulation


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function defining the forces between all neighboring cells using   % 
% the Delaunay trangulation algorithm from Matlab                    % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function force=DefineForce(cell_xy,dist,F_spr) 
  Csize=size(cell_xy,1);
  force=zeros(Csize,2);
  if (Csize>=3) 
    neigh=delaunay(cell_xy(:,1),cell_xy(:,2)); % find neighbors
    Nneigh=size(neigh,1);
    
    for ii=1:Nneigh % inspect all triangles
      num1=neigh(ii,1); num2=neigh(ii,2); num3=neigh(ii,3);  

      dx=cell_xy(num1,1)-cell_xy(num2,1);
      dy=cell_xy(num1,2)-cell_xy(num2,2);
      dxy=sqrt(dx*dx+dy*dy);
      force(num1,1:2)=force(num1,1:2)+F_spr*max((dist-dxy),0)*[dx,dy];  
      force(num2,1:2)=force(num2,1:2)-F_spr*max((dist-dxy),0)*[dx,dy];  
       
      dx=cell_xy(num1,1)-cell_xy(num3,1);
      dy=cell_xy(num1,2)-cell_xy(num3,2);
      dxy=sqrt(dx*dx+dy*dy);
      force(num1,1:2)=force(num1,1:2)+F_spr*max((dist-dxy),0)*[dx,dy];  
      force(num3,1:2)=force(num3,1:2)-F_spr*max((dist-dxy),0)*[dx,dy];  

      dx=cell_xy(num3,1)-cell_xy(num2,1);
      dy=cell_xy(num3,2)-cell_xy(num2,2);
      dxy=sqrt(dx*dx+dy*dy);
      force(num3,1:2)=force(num3,1:2)+F_spr*max((dist-dxy),0)*[dx,dy];  
      force(num2,1:2)=force(num2,1:2)-F_spr*max((dist-dxy),0)*[dx,dy];      
    end % for
  end % if
end % function DefineForce


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function determining drug diffusion using the finite difference    %
% method with different boundary conditions                          %  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function conc=DiffusionBound(conc,Ngx,Ngy,DiffConc,omega,dt,hb)
  conc2=zeros(Ngx+2,Ngy+2);
  conc2(2:Ngx+1,2:Ngy+1)=conc;
  
  scale=1-hb*omega;   % conc(x-1)=(1-hb*omega)&conc(x)
  
  % drug boundary conditions
    conc2(1,1:Ngy+2)    =scale*conc2(2,1:Ngy+2);  % von Neumann boundary cond
    conc2(Ngx+2,1:Ngy+2)=scale*conc2(Ngx+1,1:Ngy+2);
    conc2(1:Ngx+2,1)    =scale*conc2(1:Ngx+2,2);
    conc2(1:Ngx+2,Ngy+2)=scale*conc2(1:Ngx+2,Ngy+1);
  
  % drug diffusion
  for i1=2:Ngx+1
    for i2=2:Ngy+1
      conc(i1-1,i2-1)=conc2(i1,i2)+(DiffConc*dt/(hb*hb))*(conc2(i1-1,i2)+...
              conc2(i1,i2-1)+conc2(i1,i2+1)+conc2(i1+1,i2)-4*conc2(i1,i2));  
    end % for
  end % for
  
  clear conc2
end % function DiffusionBound


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function determining drug decay                                    %  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function drugDom=Decay(drugDom,Ngx,Ngy,d_drug,dt)
  % drug decay
  for i1=1:Ngx
    for i2=1:Ngy
      drugDom(i1,i2)=drugDom(i1,i2)*(1-d_drug*dt);  %%% was: *(1-d)  
    end % for
  end % for
end % function Decay


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function determining drug supply from the vessels                  %  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function drugDom=Supply(drugDom,vessel,Nvessel,Source_drug,VessRad,...
  dt,hb,xmin,ymin)
  for ii=1:Nvessel
    ix=1+floor((vessel(ii,1)-xmin)/hb);  % grid point near the vessel - x
    iy=1+floor((vessel(ii,2)-ymin)/hb);  % grid point near the vessel - y
   
    for kk=ix-2:1:ix+2
      for jj=iy-1:1:iy+2
        dx=xmin+kk*hb-vessel(ii,1);
        dy=ymin+jj*hb-vessel(ii,2);
        dxy=sqrt(dx*dx+dy*dy);
        if (dxy<=VessRad)
          drugDom(kk+1,jj+1)=Source_drug*dt; % (kk+1,jj+1)        
        end
      end
    end
  end % for
   
  clear ix iy
end % function Supply


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function determining prolonged exposure to drug                    %  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function new_exp_time=TimeDrugExposure(drug_exp,ii,cell_drug,cell_exp,dt)
  if(cell_drug(ii)>drug_exp) % want to increment exposure time for the cell
    cell_exp(ii,1)=cell_exp(ii,1)+dt;
  else
    cell_exp(ii,1)=0;        % counting time of exposure from 0
  end
  new_exp_time=cell_exp(ii,1);
end % function TimeDrugExposure


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function determining the number of neighboring cells within the   %
%% radius neighDist for each cell in cell_xy                         %   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function neighCount=Overcrowding(cell_xy,neighDist)
  neighbors=pdist2(cell_xy,cell_xy,'euclidean'); % distance between each pair of cells
  neighborsLogic=neighbors<neighDist; % 1 if the cells are nearby 
  neighCount=sum(neighborsLogic,1);   % number of nearby cells
end %function Overcrowding


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function that picks exactly 2 cells to have pre-existing resistance,  %
% based on distance to closest vessel                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cell_death=ChooseResistantCellLocation(cell_xy,vessel,Ncells,...
                PreExist_Location,cell_death,Thr_multi)
    
  minDist_vessels = [];
  for i=1:Ncells
    x = cell_xy(i,:);
    dist_vessels_all = [];
    
    % For cell i, find distance to each vessel
    for j = 1:size(vessel,1)
      dist_vessels_all(j) = pdist([x;vessel(j,:)],'euclidean');
    end
    
    %dist_vessels_all %Test
    minDist_vessels(i)=min(dist_vessels_all); % min distance for cell i
  end

  % Now, sort min distance from smallest distance to largest
  [sort_minDist_vessels,index]=sort(minDist_vessels);

  % Finally, use PreExist_Location parameter to determine which cells
  % have pre-existing resistance to drug

  if(PreExist_Location == 1) % Take 2 cells closest to vessel
    i_start = 1;
  elseif(PreExist_Location == 2) % Take 2 cells at intermediate distance
    i_start = floor(Ncells/2);
  elseif(PreExist_Location == 3) % Take 2 cells at furthest distance
    i_start = Ncells - 1;
  end

  cell_death(index(i_start))  =Thr_multi*cell_death(index(i_start));
  cell_death(index(i_start+1))=Thr_multi*cell_death(index(i_start+1));
end %function ChooseResistantCellLocation


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function create a directory in which all data will be saved         % 
% directory neme structure: WhAM_month_day_hour_minue                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function path=MakeDir
  cl=clock; clN=0;
  for ii=2:5 clN=floor(100*clN+cl(ii)); end % Go to 4 if we want hours
  disp(['your data will be saved in a directory WhAM',num2str(clN)])
  
  path=['WhAM',num2str(clN)];
  mkdir(path)
end % function MakeDir


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function drawing all cells and oxy and drug concentrations          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DrawTissue(fg1,fg2,cell_xy,cell_oxy,cell_drug,cell_damage,...
  cell_death,cell_ID,Ncells,Nclones,vessel,cellsMotherID,oxyDom,drugDom,...
  xgg,ygg,Source_oxy,Source_drug,xmin,xmax,ymin,ymax,Thr_hypo,dt,iter)

  cellSize=8;     % size of the cell marker 
  edge_param = 1.5; % thickness of line around cells to indicate age

  Ncolors=6;
  colors=[0,0,1;1,1,0;1,0,1;0,1,1;1,0,0;0,1,0];  
  Nsymbols=12;
  symbols=['o','>','v','p','d','*','<','h','^','s','x','+'];
 
  
  
  figure(fg1)
   clf(fg1)   % clear the figure
   axis([xmin,xmax,ymin,ymax])
   axis equal
   axis([xmin,xmax,ymin,ymax])
   hold on
  
   contourf(xgg,ygg,oxyDom',[0:0.02:Source_oxy],'edgecolor','none') 
   colormap(copper)
   caxis([0,0.45*Source_oxy])
   colorbar
   
     for ii=1:Ncells
        % inside: shades of green for "how far from death" the cell is 
        sc_damage=min(1,0.75*cell_damage(ii,1)/cell_death(ii,1)); 
        % outside: shades of green for drug levels
        sc_drug=min(1,cell_drug(ii,1)/(0.25*Source_drug)); 
       
        plot(cell_xy(ii,1),cell_xy(ii,2),'o','MarkerEdgecolor',[1,1-sc_drug,1],...
            'MarkerFaceColor',[0,1-sc_damage,0],'MarkerSize',cellSize,...
            'Linewidth',edge_param)
      end % for
      plot(vessel(:,1),vessel(:,2),'ro','MarkerFaceColor','r','MarkerSize',...
           2*cellSize,'LineWidth',2)
        
 
      % plot color legend
      plot(40,-65,'ko','MarkerFaceColor',[0,1.0,0],'MarkerSize',cellSize)
      plot(45,-65,'ko','MarkerFaceColor',[0,0.8,0],'MarkerSize',cellSize)
      plot(50,-65,'ko','MarkerFaceColor',[0,0.6,0],'MarkerSize',cellSize)
      plot(55,-65,'ko','MarkerFaceColor',[0,0.4,0],'MarkerSize',cellSize)
      plot(60,-65,'ko','MarkerFaceColor',[0,0.2,0],'MarkerSize',cellSize)
      plot(65,-65,'ko','MarkerFaceColor',[0,0.0,0],'MarkerSize',cellSize)
      plot(40,-70,'o','MarkerFaceColor','g','MarkerEdgeColor',[1,1.0,1],...
          'MarkerSize',cellSize,'Linewidth',edge_param)
      plot(45,-70,'o','MarkerFaceColor','g','MarkerEdgeColor',[1,0.8,1],...
          'MarkerSize',cellSize,'Linewidth',edge_param)
      plot(50,-70,'o','MarkerFaceColor','g','MarkerEdgeColor',[1,0.6,1],...
          'MarkerSize',cellSize,'Linewidth',edge_param)
      plot(55,-70,'o','MarkerFaceColor','g','MarkerEdgeColor',[1,0.4,1],...
          'MarkerSize',cellSize,'Linewidth',edge_param)
      plot(60,-70,'o','MarkerFaceColor','g','MarkerEdgeColor',[1,0.2,1],...
          'MarkerSize',cellSize,'Linewidth',edge_param)
      plot(65,-70,'o','MarkerFaceColor','g','MarkerEdgeColor',[1,0.0,1],...
          'MarkerSize',cellSize,'Linewidth',edge_param)
    
    title(['Oxygen & Cell State; cell #:',num2str(Ncells),'  Hours:',...
           num2str(iter*dt/60,4)],'FontSize',12)

       
    figure(fg2)
    clf(fg2)   % clear the figure
    axis([xmin,xmax,ymin,ymax])
    axis equal
    axis([xmin,xmax,ymin,ymax]) 
    hold on
    
          
    contourf(xgg,ygg,drugDom',[0:0.05:Source_drug],'edgecolor','none');
    colormap(bone)
    caxis([0,0.35*Source_drug])
    colorbar
   
    cloneID=zeros(size(cell_xy,1),1);
    for ii=1:size(cell_xy,1)
      jj=cell_ID(2,ii);       % motherID
      while (jj>Nclones)
        jj=cellsMotherID(jj); % grand-motherID  
      end
      cloneID(ii)=jj;    
    end
  
    for ii=1:Nclones
      ind=find(cloneID==ii);
   
      ccol=1+mod(ii,Ncolors);
      ssym=1+(ii-mod(ii,Ncolors))/Ncolors;
      plot(cell_xy(ind,1),cell_xy(ind,2),symbols(ssym),'MarkerFaceColor',...
           colors(ccol,1:3),'MarkerEdgeColor',colors(ccol,1:3),...
          'MarkerSize',cellSize)
    end
 
    Nhypo=0;
    for ii=1:Ncells
      if cell_oxy(ii,1)<Thr_hypo
        plot(cell_xy(ii,1),cell_xy(ii,2),'wo','MarkerSize',15)
        Nhypo=Nhypo+1;
      end
    end % for

    plot(vessel(:,1),vessel(:,2),'ro','MarkerFaceColor','r','MarkerSize',...
         2*cellSize,'LineWidth',2)
      
    title(['Drug & Clones;  cell #:',num2str(Ncells),' hypo:',...
            num2str(Nhypo),'  iteration:',num2str(iter)],'FontSize',12)
  
end % DrawTissue  
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function saving all data into a specified directory *path*         % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function SaveAllData(path,cell_xy,cell_age,cell_mature,cell_oxy,cell_drug,...
  cell_exp,cell_damage,cell_death,cell_ID,cellsMotherID,vessel,oxyDom,...
  drugDom,iter)

  fname=[path,'/cell_xy_',num2str(iter),'.txt']; save(fname,'cell_xy','-ascii')
  fname=[path,'/cell_age_',num2str(iter),'.txt'];  save(fname,'cell_age','-ascii')
  fname=[path,'/cell_mature_',num2str(iter),'.txt']; save(fname,'cell_mature','-ascii')
  fname=[path,'/cell_oxy_',num2str(iter),'.txt']; save(fname,'cell_oxy','-ascii') 
  fname=[path,'/cell_drug_',num2str(iter),'.txt'];  save(fname,'cell_drug','-ascii')
  fname=[path,'/cell_exp_',num2str(iter),'.txt']; save(fname,'cell_exp','-ascii')
  fname=[path,'/cell_damage_',num2str(iter),'.txt']; save(fname,'cell_damage','-ascii')
  fname=[path,'/cell_death_',num2str(iter),'.txt']; save(fname,'cell_death','-ascii')     
  fname=[path,'/cell_ID_',num2str(iter),'.txt']; save(fname,'cell_ID','-ascii')
              
  fname=[path,'/cellsMotherID_',num2str(iter),'.txt']; save(fname,'cellsMotherID','-ascii')
  fname=[path,'/vessel_',num2str(iter),'.txt']; save(fname,'vessel','-ascii')
  fname=[path,'/oxyDom_',num2str(iter),'.txt']; save(fname,'oxyDom','-ascii')
  fname=[path,'/drugDom_',num2str(iter),'.txt']; save(fname,'drugDom','-ascii')
end % SaveAllData 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%