Wrapping C++ code of agent-based tumor growth model into MATLAB

Last 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 results, doing more complicated analysis, visualization, or performing data fitting. Of course each of those tasks can be performed by writing another large piece of the code.

Two 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 other hand it was straightforward to visualize the tumor and it’s not much of a hassle to perform more complicated analysis or data fitting using variety of MATLABs built-in functions.

Each of those two approaches has their own advantages: C++ is fast, MATLAB has a lot of built-in in tools.

Today I will show how to take the best out of both approaches, i.e. I’ll show how to easily wrap our C++ code into MATLAB. The whole C++ based simulation will be invoked from MATLAB as a function and results will be visible straight from it.

First we need to add necessary libraries.

#include "mex.h"
#include <matrix.h>

Then we need to modify our original definition of global variables (see previous post) so they can be adjusted based on the input parameters (for example lattice has to be dynamically allocated).

int N;
bool* lattice;
vector<cell> cells;

char pmax;
double pDiv, alpha, ps, pmig; 
int* indcNeigh;

Then the only thing that we need to do is to delete old main() function and use mexFunction() instead. In the mexFunction() we need to read first the parameters that will be passed from MATLAB (input parameters to the function). Then we need to assign the memory for the lattice and set other auxiliary variables. Finally, after performing the simulation we need to write variables to the output, so the MATLAB can see the results.

void mexFunction(
		 int          nlhs, //number of output arguments
		 mxArray      *plhs[], //output arrays
		 int          nrhs, //number of input arguments
		 const mxArray *prhs[] //input arrays
		 )
{
    /*READING INPUT PARAMETERS*/
    double *par=mxGetPr(prhs[0]); //model parameters, first input variable
    double *settings=mxGetPr(prhs[1]); //simulation settings, second input variable

    int Nsteps = settings[0]; //number of simulation steps to perform
    N=(int)settings[1]; //size of the lattice
    pmax = par[0]; //divisions before proliferative capacity exhaustion 
    pDiv = par[1]; //division probability 
    alpha = par[2]; //spontaneous death probability 
    ps = par[3]; //probability of symmetric division 
    pmig = par[4]; //probability of migration
    
    /*CREATING LATTICE AND AUXILARY VARIABLES*/
    //creating lattice with a given size
    lattice = new bool [N*N];
    fill_n(lattice, N*N, false);
    
    indcNeigh = new int [8];//cell neighborhood
    //assigning values
    indcNeigh[0] = -N-1; indcNeigh[1] = -N; indcNeigh[2] = -N+1;
    indcNeigh[3] = -1; indcNeigh[4] = 1; indcNeigh[5] = N-1;
    indcNeigh[6] = N; indcNeigh[7] = N+1;

    /*OLD PIECE OF MAIN() FUNCTION*/	
    srand(time(NULL)); //initialize random number generator
    initialize(); //initialize CA
    simulate(Nsteps); //simulate CA
    
    /*WRITING TO OUTPUT VARIABLES*/
    //1) writing the lattice as a logical matrix
	plhs[0] = mxCreateLogicalMatrix(N,N); //lattice
    mxLogical *ptr = mxGetLogicals(plhs[0]);
    for(int i=0; i<N*N; i++) 
        ptr[i] = lattice[i];

    //2) writing the cells vector as a vector
    plhs[1] = mxCreateDoubleMatrix(1,cells.size()*3,mxREAL);
    double* ptr2 = mxGetPr(plhs[1]);
    for(int i=0; i<cells.size(); i++) {
        ptr2[3*i] = cells.at(i).place;
        ptr2[3*i+1] = cells.at(i).p;
        ptr2[3*i+2] = cells.at(i).is_stem;
    }
}

After making all of the changes we can compile our .cpp file using the MATLAB command “mex CA.cpp”. Remember all other pieces of C++ associated with model itself remain the same – there is no need to change anything.

From now on we can perform our simulation in MATLAB environment by invoking the following command

[lattice, cells] = CA([pmax, pdiv, alpha, ps,pmig],[nSteps, N]);

and further use of all MATLABs benefits, especially data fitting features. At the same time all of the simulations will be performed at C++ speed! Imagine using Parallel Toolbox in MATLAB to harness multiple CPUs with only few lines of the code.

Advertisement