# Cancer stem cell driven tumor growth model in less than 70 lines of code

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 (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”

1. […] 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 […]

Like

2. […] 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 […]

Like

3. […] the previous posts (CA in MATLAB and C++) I’ve shown how to implement cancer stem sell driven tumor growth model. In both […]

Like

4. […] 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 […]

Like

5. […] 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 […]

Like

6. […] 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 […]

Like

7. […] 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. […]

Like

8. […] MATLAB (post with the code here) […]

Like

9. […] 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 […]

Like

10. Amir Behjat says:

Hi How can I cite your work if I want to use it?

Like