Cancer stem cell driven tumor growth model in C++

Because of the feedback that I’ve received after publishing first three posts, I’ve decided to change a tool of interest from MATLAB to C++. Today, I’ll show how to quickly implement a cancer stem 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 post).

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.

First we start with including the necessary headers and setting the namespace. Apart from the standard functions and datatypes, we will use the vector datatype to store the cells and a shuffling function from algorithm library.

```#include <vector>
#include <algorithm>
#include <stdlib.h>
#include <time.h>

using namespace std;
```

Now we can define a cell. It will be defined by three variables: index to the place on the lattice (integer), remaining proliferation capacity (char), and boolean variable defining stemness.

```struct cell {//defining cell
int place;
char p;
bool is_stem;
};
```

Next step is to define the lattice size, the lattice itself, vector that will contain all viable cells, and auxiliary variable defining the cells neighborhood.

```static const int N = 2000; //lattice size
bool lattice[N*N] = {false}; //empty lattice
vector<cell> cells; //vector containing all cells present in the system

static const int indcNeigh[] = {-N-1, -N, -N+1, -1, 1, N-1, N, N+1};//neighborhood
```

The parameters of the model are defined as follows.

```char pmax=10; //proliferation capacity
double pDiv=1./24.; //division probability
double alpha=0.05; //spontaneous death probability
double ps=0.05; //probability of symmetric division
double pmig=10./24.; //probability of migration
```

Having made all of those initial definitions we can finally start coding the main program.
Let us start with writing the function that will initialize the whole simulation, i.e. fill the lattice boundary and put the initial stem cell.

```void initialize() {
for (int i=0; i<N; i++) {lattice[i]=true;};//filling left
for (int i=0; i<N*N; i=i+N) {lattice[i]=true;};//filling top
for (int i=N-1; i<N*N; i=i+N) {lattice[i]=true;};//filling bottom
for (int i=N*(N-1); i<N*N; i++) {lattice[i]=true;};//filling right

lattice[N/2*N+N/2] = true; //initial cell in the middle
cell initialCell = {N/2*N+N/2,pmax,true};//initial cell definition
cells.push_back(initialCell);
}
```

As in the previous post, we set all the boundary values of lattice 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.

The second auxiliary function that we will use returns the index to the randomly selected empty space around a given spot.

```int returnEmptyPlace(int indx) {
int neigh[8], nF = 0;
for(int j=0;j<8;j++) {//searching through neighborhood
if (!lattice[indx+indcNeigh[j]]) {//if free spot
neigh[nF] = indx+indcNeigh[j]; //save the index
nF++; //increase the number of found free spots
}
}
if(nF) {//selecting free spot at random
return neigh[rand() % nF];
} else {//no free spot
return 0;
}
}
```

If there is no free spot the function returns 0.
Finally we can code the part in which all the magic happens – main simulation procedure. It is hard to dissect the whole procedures into parts, so I did my best to explain everything in the comments.

```void simulate(int nSteps) {

vector<cell> cellsTmp;
int newSite;
cell currCell, newCell;

for (int i=0; i<nSteps; i++) {
random_shuffle(cells.begin(), cells.end()); //shuffling cells
while (!cells.empty()) {
currCell=cells.back(); //pick the cell
cells.pop_back();
newSite = returnEmptyPlace(currCell.place);

if (newSite) {//if there is a new spot
newCell = currCell;
newCell.place = newSite;
if ((double)rand()/(double)RAND_MAX < pDiv) {
if (currCell.is_stem) {
lattice[newSite]=true;
cellsTmp.push_back(currCell);
if ((double)rand()/(double)RAND_MAX > ps) {//asymmetric division
newCell.is_stem = false;
}
cellsTmp.push_back(newCell);
} else if (currCell.p > 0 && (double)rand()/(double)RAND_MAX > alpha) {
currCell.p--;
newCell.p--;
lattice[newSite] = true;
cellsTmp.push_back(currCell);
cellsTmp.push_back(newCell);
} else {
lattice[currCell.place] = false;
}
} else if ((double)rand()/(double)RAND_MAX < pmig) {
lattice[currCell.place] = false;
lattice[newSite] = true;
cellsTmp.push_back(newCell);
} else {//doing nothing
cellsTmp.push_back(currCell);
}
} else {//no free spot
cellsTmp.push_back(currCell);
}
}
cells.swap(cellsTmp);
}
}
```

Now we wrap everything in the main function, compile and run.

```int main() {
srand(time(NULL)); //initialize random number generator
initialize(); //initialize CA
simulate(24*30*6);
return 0;
}
```

Everything in about 100 lines of the code.
What about the speed of the code? It took about 40 seconds to simulate the tumor presented below, which is consisted of about 460,000 cells (on 2000×2000 lattice).

Why is it “almost” the same model as the one presented in previous post? The answer is in details. In the above C++ implementation lattice is updated after every single cell movement/death, i.e. because of the random shuffling a new free spot can be created for a cell that at the beginning of the iteration was quiescent (without a free spot), so it can successfully proliferate in the same iteration. In the MATLAB implementation that kind of behavior was avoided.

19 thoughts on “Cancer stem cell driven tumor growth model in C++”

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

Like

2. […] week I posted a C++ implementation of basic cancer stem cell driven tumor growth model. About 100 lines of a code allowed to perform simulation, but there was nothing about exporting the […]

Like

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

Liked by 1 person

4. […] there is a free spot – nothing else. We can adapt the very basic C++ implementation from the previous post as a baseline: we have a boolean matrix (in which value true indicates that the spot is occupied […]

Like

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

Like

6. […] C++ (post with the code here) […]

Like

7. […] few tricks (in C++) how to efficiently store output of 2D agent-based model presented in one of the previous posts. As a result we will reduce the size of generated simulation output from about 1Gb to reasonable […]

Like

8. How did you obtain the picture at the end of this post?

Like

• I used MATLAB to plot the tumor (see my “Cancer stem cell driven tumor growth model in less than 70 lines of code” post). You can save the data from the C++ sim to a text file and then load it to any other program to plot the results.

Like

9. Hi,Could you please refer to the model behind the code? I think time in line 2 (of last fragment) is not previously defined.

Like

• Hi, that line has nothing to do with the model itself. It is a (pseudo)random number generator initialization (srand) which requires initial seed, which in my case is taken to be current time.

Like

• Thank you. I made it work using #include

Like

10. alexa says:

Hi, i tried to use this code with code block and it does not go … can you help me?

Like

• Hi, alexa. Do you have compilation error or is it something else? If this is a compilation error please provide the error message.

```#include <time.h>
```

at the beginning of the code in order to avoid ‘time’: identifier not found error.

Like

• alexa says:

Hi jpoleszczuk, thank you for answering! I have already add #include but, after build I have this message:
” Build: Debug in … (compiler: GNU GCC Compiler)”, the problem is that I’m new to C ++ and I’m not be able to understand because I only see black screen.

Like

• Hi alexa, you see only black screen, because no output function is implemented in the code.
You can save an image with the tumor visualization using ALLEGRO, by adding

```#include <allegro5/allegro.h>
#include <allegro5/allegro_image.h>
```

at the top of the source file (you need to use Image add-on) and the following lines in the main function (after simulate(); line)

```	ALLEGRO_BITMAP *cellsAllegro = NULL;
al_init();
cellsAllegro = al_create_bitmap(N, N);

al_set_target_bitmap(cellsAllegro);
al_clear_to_color(al_map_rgb(0, 0, 0));

ALLEGRO_LOCKED_REGION *lr = al_lock_bitmap(cellsAllegro, ALLEGRO_PIXEL_FORMAT_ABGR_8888, ALLEGRO_LOCK_WRITEONLY);
unsigned char *ptr = (unsigned char *)lr->data;
for (unsigned int j = 0; j < cells.size(); ++j) {
ptr[cells.at(j).place * 4] = 255; //red
ptr[cells.at(j).place * 4+1] = 0;   //green
ptr[cells.at(j).place * 4+2] = 0;   //blue
ptr[cells.at(j).place * 4+3] = 255; //alpha
}
al_unlock_bitmap(cellsAllegro);

al_save_bitmap("simResults.bmp", cellsAllegro);

al_destroy_bitmap(cellsAllegro);
```

This will write the tumor visualization to simResults.bmp image file.

Like

11. alexa says:

How did you obtain a plot of this post using ALLEGRO?

Like

12. Margriet Palm says:

Thanks for this code! Maybe I’m misreading the code, but it seems to me that a cell that can’t find a spot for proliferation or migration also can’t die:

} else {//no free spot
cellsTmp.push_back(currCell);
}

also, stem cells seem to be immortal. Is there a reason behind this?

Like

• I’m glad that you found the code useful. And you are right – the cell without any empty spot in the surrounding can’t die and cancer stem cells are immortal. Both of the above are modeling assumptions that you can easily modify. The rationale behind the former is that the cell death occurs because of unsuccessful division.

Liked by 1 person