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

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
           newSite = returnEmptyPlace(;
           if (newSite) {//if there is a new spot
               newCell = currCell;
      = newSite;
               if ((double)rand()/(double)RAND_MAX < pDiv) {
                   if (currCell.is_stem) {
                       if ((double)rand()/(double)RAND_MAX > ps) {//asymmetric division
                           newCell.is_stem = false;
                   } else if (currCell.p > 0 && (double)rand()/(double)RAND_MAX > alpha) {
                       lattice[newSite] = true;
                   } else {
                       lattice[] = false;
               } else if ((double)rand()/(double)RAND_MAX < pmig) {
                   lattice[] = false;
                   lattice[newSite] = true;
               } else {//doing nothing
           } else {//no free spot

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

int main() {
    srand(time(NULL)); //initialize random number generator
    initialize(); //initialize CA
    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).

Screen Shot 2015-06-06 at 12.24.53 AM

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++

    • 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.


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

      #include <time.h>

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


      • 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.


      • 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;
        	cellsAllegro = al_create_bitmap(N, N);
        	al_clear_to_color(al_map_rgb(0, 0, 0));
        	unsigned char *ptr = (unsigned char *)lr->data; 
        	for (unsigned int j = 0; j < cells.size(); ++j) {
        		ptr[ * 4] = 255; //red
        		ptr[ * 4+1] = 0;   //green
        		ptr[ * 4+2] = 0;   //blue
        		ptr[ * 4+3] = 255; //alpha
        	al_save_bitmap("simResults.bmp", cellsAllegro);

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


  1. 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

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


    • 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

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s