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)&&(k0)&&(l>0)&&(k0) % 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)=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)(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','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)