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

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

LikeLike

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

LikeLike

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

LikeLiked by 1 person

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

LikeLike

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

LikeLike

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

LikeLike

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

LikeLike

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

LikeLike

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.

LikeLike

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

LikeLike

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.

LikeLike

Thank you. I made it work using #include

LikeLike