Today I’ve decided to show how to efficiently code a cancer stem cell driven tumor growth model in less then 70 lines of code in MATLAB.
Quick guide to the model: A cell, either cancer stem cell (CSC) or non-stem cancer cell (CC), occupies a single grid point on a two-dimensional square lattice. CSCs have unlimited proliferation potential and at each division they produce either another CSC (symmetric division) or a CC (asymmetric division). CCs are eroding their proliferation potential at each division and die when its value reaches 0. Moreover, at each proliferation attempt, CCs may undergo spontaneous death and then be removed from the system.
We start with defining initial settings of a domain and simulation timespan.
N = 1000; %square domain dimension nSteps = 6*30*24; %number of simulation steps
Next we define the values of all parameters used in the simulation.
pprol = 1/24; %probability of proliferation pmig = 10/24; %probability of migrating pdeath = 1/100; %probability of death pmax = 10; %initial proliferation capacity of CC ps = 3/10; %probability of symmetric division
Now we can intialize domain and place initial CSC in its center. Our computational domain is represented by N x N Boolean matrix L (a true value indicates the lattice point is occupied). An additional integer vector cells is maintained to store indices for all viable cells present in the system (to avoid costly search for cells on the lattice).
L = false(N,N); L([1:N 1:N:N*N N:N:N*N N*(N-1):N*N]) = true; %boundary L(N*round(N/2)+round(N/2)) = true; cells = int32(N*round(N/2)+round(N/2)); cellsIsStem = true; cellsPmax = uint8(pmax);
To reduce the amount of used memory we utilized int32 and unit8 types as they occupy less memory than double (double is default type in MATLAB). We also set all the boundary values of L to true in order to make sure that we will not address any site out of the lattice. Typically one would use an if statement when addressing the lattice site to make sure that index to a site is within the domain. Here, because we have set the boundaries to true without adding those sites to cells vector, we don’t need to do that – boundary will be treated as an occupied site, but not as a viable cell.
We can now define auxiliary variables that will be utilized further.
aux = int32([-N-1 -N -N+1 -1 1 N-1 N N+1])'; %indices to heighborhood Pms = perms(uint8(1:8))'; %permutations nP = size(Pms,2); %number of permutations
Variable aux simply defines cell’s neighborhood and Pms is a variable in which we store all possible permutations of variable aux.
Now we can start the main loop of the program and randomly shuffle cells at its beginning.
for i = 1:nSteps sh = randperm(length(cells)); cells = cells(sh); cellsIsStem = cellsIsStem(sh); cellsPmax = cellsPmax(sh);
Now few lines in which we first create indices to neighborhoods of all of the cells already in random order (variable S), then we flag by 0 all of the spots that are occupied and finally select only indices to those cells that have at least one free spot in the neighborhood (variable indxF).
S = bsxfun(@plus,cells,aux(Pms(:,randi(nP,1,length(cells))))); S(L(S)) = 0; %setting occupied spots to 0 indxF = find(any(S)); %selecting cells with at least one free spot nC = length(indxF); %number of cells with free spot
Now we can decide about the faith of the cells that can perform action (are not completely surrounded by other cells).
P = rand(1,nC)<pprol; %proliferation Ps = P & rand(1,nC)<ps & cellsIsStem(indxF); %symmetric division De = P & (cellsPmax(indxF) == 0); %proliferation capacity exhaution D = P & (rand(1,nC)<pdeath) & ~cellsIsStem(indxF); %death at proliferation attempt M = ~P & (rand(1,nC)<pmig); %go when no grow
In the additional loop we perform update of the system.
del = D | De; %cells to delete act = find((P | M) & ~del); %indices to the cells that will perform action for ii = act %only for those that will do anything ngh = S(:,indxF(ii)); %cells neighborhood ngh(ngh==0) = ; %erasing places that were previously occupied indO = find(~L(ngh),1,'first'); %selecting free spot if ~isempty(indO) %if there is still a free spot L(ngh(indO)) = true; %updating occupancy if P(ii) %proliferation cells = [cells uint32(ngh(indO))]; %adding new cell if Ps(ii) %symmetric division cellsIsStem = [cellsIsStem true]; cellsPmax = [cellsPmax cellsPmax(indxF(ii))]; else cellsIsStem = [cellsIsStem false]; cellsPmax = [cellsPmax cellsPmax(indxF(ii))-1]; if ~cellsIsStem(indxF(ii)) cellsPmax(indxF(ii)) = cellsPmax(indxF(ii))-1; end end else %migration L(cells(indxF(ii))) = false; %freeing spot cells(indxF(ii)) = uint32(ngh(indO)); end end end
At the end we need to erase cells that died.
if ~isempty(del) %updating death L(cells(indxF(del))) = false; cells(indxF(del)) = ; cellsIsStem(indxF(del)) = ; cellsPmax(indxF(del)) = ; end end
This is it. Whole CA is coded and we are ready to run simulations!
In order to visualize simulations results we can implement short function.
function visualize( N, cells, cellsIsStem, cellsPmax, pmax ) M = ones(N,N,3); %matrix for image color = hot(3*pmax); %colormap %drawing CCs M(cells(~cellsIsStem)) = color(cellsPmax(~cellsIsStem)+1,1); M(cells(~cellsIsStem)+N*N) = color(cellsPmax(~cellsIsStem)+1,2); M(cells(~cellsIsStem)+2*N*N) = color(cellsPmax(~cellsIsStem)+1,3); %drawing CSCs, we want to draw them as 3x3 points aux = int32([-N-1 -N -N+1 -1 1 N-1 N N+1])'; CSCs = cells(cellsIsStem); plusSurr = bsxfun(@plus,CSCs,aux); M(plusSurr) = color(2*pmax,1); M(plusSurr+N*N) = color(2*pmax,2); M(plusSurr+2*N*N) = color(2*pmax,3); figure(1) imshow(M); end
Below is the visualization of the tumor simulated using the above code. It is consisted of 369,504 cells in total (8563 CSCs). The whole simulation took about 12 minutes on my laptop, what taking into account the lattice size and number of cells is a reasonable time.
12 thoughts on “Cancer stem cell driven tumor growth model in less than 70 lines of code”
[…] we will modify the code from the previous post in order to simulate tumor in 3D. There are only a few changes: redefinition of the domain and the […]
[…] cell driven tumor growth model in C++. It is almost the same model as implemented in my previous post (I’ll explain why it is “almost” the same at the end of this […]
[…] the previous posts (CA in MATLAB and C++) I’ve shown how to implement cancer stem sell driven tumor growth model. In both […]
[…] weeks ago I posted a MATLAB code for the very similar model, which was quite fast but a lot slower than C++ implementation. On the […]
[…] shown how to implement agent-based model of cancer stem cell driven tumor growth, both in MATLAB and C++. Having the implementations we can go one step further and perform some analysis of the […]
[…] will use the codes for basic ABM model posted before as a template (MATLAB version here and C++ version here). In both cases we considered rectangular lattice with Moore neighborhood. If […]
[…] Python doesn’t seem to be the first programming language people go to when developing cellular automata models. However, given that Python is an object-orientated language that is easy to read and write, it might actually be ideal for such models, especially if you prefer to think from the perspective of the agent (if you’d rather model using matrices you can do that too by using Python’s NumPy package). As we’ll see, it can also be pretty fast. Below I’ll cover one way to implement the Cancer Stem cell CA described in the Cancer Stem Cell CA in Matlab post. […]
[…] MATLAB (post with the code here) […]
[…] do is to update speed comparison post with Julia. First, I’ve decided to rewrite the partially vectorized MATLAB code in Julia without making any essential changes. After browsing the web I found couple of useful […]
Hi How can I cite your work if I want to use it?
I’m glad that you are considering to use the model. You can either cite the blog post directly (follow standard guidelines) or cite one of my papers with Heiko Enderling (e.g. https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004025 ).
Thanks. I finally used the C code you had which was so fast. Thanks again.