A 3D Cellular Automata Cancer Stem Cell Model using MATLAB and App Designer

Abstract

This is a flexible and user-friendly 3D cellular automata model that uses MathWorks’ MATLAB and App Designer to simulate cancer tumor growth that takes into account certain types of Cancer Stem Cell (CSC) behavior. A Cellular Automata model, or a grid containing cells with rules governing their current state, is ideal for MATLAB (Matrix Laboratory) in that each “cell” location can be a grid point in a 2D or 3D matrix and subsequent element-wise operations can be applied to simulate tumor growth. The simulation begins by taking in several cellular automata inputs: dimensions of the viewing area (where each unit represents a cell), how long (or how many cycles) the simulation will run, and the neighborhood. Each cycle, the tumor growth depends on the maximum lifespan of a non-CSC cell, the probability that a cell divides (proliferates), migrates to a different location, dies, or symmetrically divides (symmetric is when a CSC divides into two CSCs, asymmetric is when a CSC divides into a CSC and a non-CSC cancer cell). Using these inputs, my MATLAB/App Designer model simulates the growth of a tumor from a single CSC and outputs the number of CSCs and the total number of cells and plots the spatial cell distribution in the 3D array. This is an analysis of how each of the factors impact the cancer stem cell and non-cancer stem cell growth rates and I found that the results from the model seem to be consistent with the expected results. The idea is to make sure that any cancer treatment targets cancer stem cells to ensure that the cancer is eradicated. In the near future, I would also like to enhance the model and incorporate angiogenesis and metastasis somehow to make it more accurate, as the ability to migrate to distant sites is something that needs to be considered for successful treatment.

Introduction

This past summer, I stumbled across the concept of Cellular Automata (CA) and how it may be able to model the growth of cancerous tumor cells.^[1][2] A paper titled “Cancer Stem Cell Modeling Using a Cellular Automaton” sounded intriguing.^[3] After reading the paper, I reached out to the authors and they sent me a link to their follow-up research article entitled “Treatment Analysis in a Cancer Stem Cell Context Using a Tumor Growth Model Based on Cellular Automata.”^[4] The paper mentioned how cancerous tumors and other emergent behaviors could be modeled by Cellular Automata and how it differed from previous approaches using differential equations to describe tumor growth. It also highlighted the implications of the treatment of cancer growth in a cancer stem cell context. I wanted to examine Cellular Automata for myself and one of the basic analytical tools that was mentioned during the Stanford EXPLORE Lecture Series was MathWorks’ MATLAB. After completing an online course in MATLAB and App Designer,^[5] I wanted to see if I could create a simple model in MATLAB so I could examine the various parameters that lead to cancerous tumor growth and what parameters might help with the treatment of cancer. This paper is a result of my studies.

What is Cellular Automata?

A Cellular Automata is a collection of “cells” that is arranged on a grid, where each “cell” can change state over time according to a set of rules. These rules and how each “cell” changes can depend on their own state and the state of their neighbors. The singular form of Cellular Automata is Cellular Automaton.
Since many systems can be broken down into the most basic units resembling this grid, Cellular Automata has found uses in many areas such as mathematics, computer science, physics, and biology. There are many applications ranging from bioinformatics to urban flooding to protein coding to modeling snowflake crystallization, etc.
There are two fundamental concepts in Cellular Automata: spatial structure and local interaction. The spatial structure can consist of a one, two, or multi-dimensional grid in a square, hexagonal, or some other regularly shaped configuration. The cells live on a grid and each cell has a neighborhood. The neighborhood can be defined in any number of ways, but it is typically a list of adjacent cells. Microsoft Visio was used to illustrate several two-dimensional examples. See Figure 1 and Figure 2. On a square grid, two common neighborhoods are the von Neumann neighborhood (diamond-shaped) like in Figure 1A and the Moore neighborhood (square-shaped) like in Figure 1B. The black cell is the cell of interest while the red or blue cells are the neighbor cells where r is an indicator for range.

Figure 1 – Square grid neighborhoods: (A) von Neumann (B) Moore

Each cell has a state. The simplest example has two possibilities of either 1 or 0. The cell state can change over time (for modeling and simulation purposes, time would be illustrated in discrete time steps or “generations”) and how each cell state changes depends on its own current state (for example, the black cell in one of the figures above) and the states of the cells in the neighborhood defined by r (for example, the colored cells in the same figure). This is usually dictated by a set of rules. For example, in a two-dimensional r = 1 case, one rule could be “if the cell is ‘1’ and all four cells around it is also a ‘1’, then the state in the next generation will be a ‘0’.” A second rule could be “if the cell is ‘0’ and three of the four cells around it is also a ‘1’, then the state in the next generation will be a ‘1’, and so on. By applying the set of rules on each cell on the grid, the state of each cell is determined for the next generation. Then the same rules are applied for the next generation and for each successive generation. In this way, interesting patterns emerge over time, such as:^[6]

Figure 2 – CA pattern from Wolfram Web Resources

Of course, there are various forms of rules that can be applied, but coupled with spatial structure, CA is capable of producing a wide variety of unexpected behavior.
Stephen Wolfram^[7] stated that empirical studies suggest that four qualitative classes may evolve even when using different initial states and different rules. These are: 1) spatially homogeneous state, 2) sequence of simple stable or periodic structures, 3) chaotic aperiodic behavior, and 4) complicated localized, some propagating.
In the more simple CA models, Wolfram states this can translate to: 1) disappears with time, 2) evolves to a fixed finite size, 3) grows indefinitely at a fixed speed, and 4) grows and contracts irregularly.
The grid, neighborhood definition (von Neumann, Moore, or other), and the range can vary between applications and models. Figure 3 shows an example of a hexagonal grid and neighborhoods. Again, the black cell is the cell of interest while the colored cells are the neighbor cells where r is an indicator for the number of cells that can be considered in a cell’s neighborhood.

Figure 3 – Hexagonal grid neighborhoods

MATLAB (short for Matrix Laboratory) seems to be an ideal tool to use to model CA especially if we assume each cell lives in either a one, two or three-dimensional square grid since the grid can be represented by a simple vector or a two dimensional or three-dimensional matrix.

What is the Cancer Stem Cell Model?

The Cancer Stem Cell (CSC) Model is based on the assumption that within a malignant tumor that is composed of a variety of cancer cells, the cancer is primarily sustained and driven only by a small population of “stem cells”. This is similar to normal stem cells that renew and sustain our bodies, particularly when we are young. In this case, any treatment that does not attack and kill these “stem cells” will ultimately fail.
Figure 3 shows the clonal evolution (sometimes called stochastic) model, which is a non-CSC model. This model assumes that each type of cancer cell (Tumor cell A, Tumor cell B, Tumor cell C, and Tumor cell D) can propagate and produce a new tumor. Over time, mutations can occur and increase the chance of malignancy. The cancer cells can then become invasive, spread to other sites, and therefore become increasingly difficult to treat. In this model, treatment needs to target each type of cancer cell aggressively. We can see examples of CSC Models in Figures 4 and 5.

Figure 4 – General Non-CSC Model

Figure 4 shows the Cancer Stem Cell (CSC) model. This model assumes that only a certain type of tumor cells (CSC) can propagate and drive tumor growth. The non-CSC tumor cells do not necessarily propagate and conventional therapies can kill them. However, the CSCs can renew themselves and can also produce other types of non-CSC cells. Hence if a small portion of the CSCs survives treatment, the tumor can reappear and grow. Therefore, targeting the CSC tumor cells is key to successful eradication of the cancer.

Figure 5 – CSC Model

The CSC model appears to have limits. Researchers at the University of Michigan have reported that not all types of cancer, including melanoma, follow this model.^[8][9] For most of the reported cancer types that follow the classic CSC model, approximately 5% or less of the cancer cells are CSCs, whereas for melanoma it is estimated at 25%. Additional research and studies need to be done to see if the types of cancers that do not follow the model can be modeled using a modified CSC model. In any case, my paper and model will only consider the CSC types of cancer because, as mentioned before, CSC types of cancer are primarily driven by the propagation of cancer stem cells.
Over the past few years, there have been various CA models applied to Cancer Stem Cell growth.^[10][11][12] Some of them are very complex with more than a dozen parameters, ranging from uptake rates, binding rates, secretion rates, etc. I plan on using a much simpler model (albeit three dimensional) together with MATLAB and App Designer, which can still be used to visualize the fundamentals of Cancer Stem Cell growth and simplifies the process.

MATLAB and App Designer

As I mentioned earlier, MATLAB seems to be well suited to model CA and therefore the CSC types of cancerous tumor grown if we assume the cells live on either a one, two or three-dimensional square grid. The spatial structure, or the grid the cells exist on, can be visualized as a matrix, with each element representing one cell.
MATLAB is an environment and language for technical computing. The environment is a set of tools and windows that allow you to develop programs and visualize large amounts of data either graphically or through animation. MATLAB’s basic data element is an array. This is an advantage for matrix-oriented computation like CA. It also has specialized toolboxes that have a collection of functions for many different fields of interest. Conveniently, one of these toolboxes is dedicated to Computational Biology, which includes bioinformatics, systems biology, and biostatistics.
App Designer is a development environment within MATLAB that provides a large set of components such as dropdown menus, buttons, various data entry components, tools, and instrumentation controls. This provides a way to customize and create a more interactive dynamic experience. I attempt to use App Designer to create a basic user interface to the CA model I develop in this paper. It will allow a more convenient and user-friendly way to change parameters.

Basic CSC CA Tumor Growth Model

One CSC tumor growth model that seems reasonable to model in MATLAB and App Designer is one that is outlined in several papers.^[13][14] Some of the input parameters they consider are: the probability of symmetric division, proliferation capacity (which includes probability of proliferation and maximum proliferation capacity), migration potential and probability of spontaneous death.
For the rules of the CA model, a tumor population can consist of both CSCs and non-CSCs. Each cell will occupy a single grid point. CSCs are assumed to be immortal and have unlimited proliferation potential, meaning that they can keep reproducing forever. However, the non-CSCs can only survive a limited number of times (p_max) before dying. For every cycle or generation, each CSC has a probability to proliferate (prob_prol) and divide symmetrically to produce two daughter cells. A CSC cell has a chance of either symmetric division producing two new CSCs or asymmetric division (producing one CSC and one non-CSC). Since it is either one or the other, the probability of asymmetric division is (1 – prob_sym), where prob_sym is defined the probability of symmetric division. New non-CSC cells will not reproduce and will die spontaneously after reaching the number of cycles specified by p_max. Since the cells are on a grid, we assume some adjacent cells must be open in order to reproduce and that there is a probability (prob_mig) during each cycle that a cell can migrate into an open neighboring site. If a cell is fully surrounded, we can assume that it will not divide. Because of this migration, I will assume a Moore-type neighborhood (See Figure 1) for my 3D model since it seems unreasonable that cell growth would be restricted only along an axis. I will assume r=1 for simplicity. Lastly, we can assume that in each simulation step, there is a probability that each non-CSC cell can undergo spontaneous death, prob_death, and be removed from the grid. This probability can be set to 0 or some low number that would be realistic, but it can be used to simulate cancer treatments.

Implementation

I was able to successfully put together a three dimensional CA MATLAB model for tumor growth and used App Designer to add a user interface so I could quickly change parameters and re-run simulations. While there is room for improvements and expandability, the basic features are there.
Basically, I started in the MATLAB editor environment to develop the main source code. It started with not only defining all the CSC parameters such as the ones mentioned above, but also the CA parameters like grid size (N x N x N), where N is the number of cell sites along one axis, and the number of steps/cycles of iteration (nSteps). I set all of the cell sites at the edges of the array space to ‘1’s, so the cells cannot propagate out of the space because there are no unoccupied cell sites. One difficulty I had initially was figuring out how to address the nearest neighbors in 3-dimensional space. After searching through the MATLAB documentation and some MATLAB forums, I found a concept called linear indexing. With that issue addressed, I was also able to approximate the center of the 3-dimensional grid array and place a CSC there. The next phase was to define a loop from 1 to N and do the calculations based on the probabilities, keep track of open sites, keep track of cells that may move, remove dead cells, and repeat for N cycles.
After the source code was mostly finished, I used a graphical user interface to prompt the user for the input values: grid size, number of steps, probabilities (as mentioned above), and p_max. When both were nearly complete, I moved the source code into the code that App Designer automatically generates and hooked up all the callback features so that the interface and the source code worked together.
My App Designer interface is shown both in Figure 6 (Design View) and Figure 7 (Code View). For those unfamiliar with App Designer, the Design View is where you can easily drag components for the component library to a canvas, place them where you want, and assign them properties (starting value, min/max/steps, colors, etc.) Many components can be customized and App Designer will automatically generate the code necessary to create the runtime version.

Figure 6 – App Designer User Interface (Design View)

Figure 7 – App Designer User Interface (Code View)

Results

Figure 8 shows a sample output of the simulation. The first simulation shown in Figure 8(A) is with N=100, nSteps=200, prob_prol=4%, prob_mig=40%, prob_death=1%, prob_sym=30% and p_max = 10. I chose the initial conditions roughly based on what I have seen in a few papers. I added the final number of CSCs, cells and “density” on the plot to give a more quantitative analysis. I will show some plots based on these numbers later. The “density” actually represents an artificial density calculated by dividing the total number of cells by the entire number of sites in the grid space (on NxNxN sites). I would have liked to calculate the volume the simulated tumor encompassed and then divided by the number of cells, but I was unable to figure out an easy way to do this. In any case, this explains why the “density” values are so low as it is dependent on N. Figure 8(B) shows the exact same conditions except that the number of steps (nSteps), or generations, is doubled. The simulated tumor is shown to have grown significantly.

Figure 8 – (A) Simulation with N=100, nSteps =200, prob_prol=4%, prob_mig=40%, prob_death=1%, prob_sym=30% and p_max = 10; (B) Same simulation with nSteps=400

By changing the probabilities, steps, and p_max, we can generate plots to interpret. One thing to note is that even if you run the same exact conditions over and over again, the results will be different (more CSCs, less CSCs, different occupied sites, etc.), since this is a function of probabilities and the final value and position of cells will evolve over time (or steps). For example, Table 1 shows 10 simulation results using the same probabilities and parameters. The variation can be seen in the number of CSCs and number of total cells.

Table 1 – Simulation results for 10 runs with same input values

In order to compare populations as we vary parameters, it may be necessary to carry out multiple runs for each set of parameters and average the results. Table 2 averages the results for 10 runs per nStep value and the results are plotted in Figure 9. The results of the simulation do seem consistent with the expectation that both the number of CSCs and total cells increase with increasing step count. Note the exponential nature of cell/tumor growth as I formatted the y-axis in Figure 9 as logarithmic. Something else that can be noted is the percentage of cells that are CSCs (See Table 2.) Aside from the initial growth phase (nSteps = 100), the percentage is between 1 to 2%, which seems to agree with University of Michigan’s statement^[15] that most of the reported cancer types that follow the classic CSC model will have approximately 5% or less of the cancer cells are CSCs, whereas for melanoma it is estimated at 25% (and is believed not to follow the CSC model).

Table 2 – Average of 10 runs per nStep value

Figure 9 – # of CSCs and # of total cells as a function of nSteps (or time)

Using the same methodology of running multiple simulations for each data point to average out the randomness and yet still show the trend, I created similar plots for the various probability parameters and p_max. This time, though, I will do 100 simulations per point and then average the results. It is fairly easy to do with a loop statement in MATLAB. The only drawback is computing time as the number of cells and CSCs increase since there are already a few loops in the source code. This will put a constraint on the ranges I will simulate.
Figures 10, 11, 12 and 13 are simulations with the baseline condition of N=100, nSteps =200, prob_prol=4%, prob_mig=40%, prob_death=1%, prob_sym=30% and p_max = 10, but will vary only a single parameter as noted in the plot.
Figure 10 shows the effect of increasing the probability of proliferation on the cell count and hence tumor growth. The impact seems exponential for the proliferation of CSCs and the total number of cells but the slope starts to flatten out for the total number of cells for higher probabilities. This can be due to the size limitation of the grid array or the effect of one of the other parameters. We will look into this at another time, but the trend is as expected: higher proliferation probabilities will cause an exponential-like increase in the number of CSCs and total cells.

Figure 10 – Impact of Probability of Proliferation on cell growth

Figure 11 shows the effect of increasing the probability of migration on the cell count. The impact seems small for the parameter space we all looking at if we look at purely cell count. See Figure 11(A). However, if we look at the cell distribution over the grid, we see a big difference as shown in Figure 11(B). The cells in the blue color use the condition of prob_mig=0% while the cells in the green color use the condition of prob_mig=100%. Therefore, although the cell counts may be the same, an increase in the migration probability leads to an increase in spatial distribution.

Figure 11 – Impact of Probability of Migration on cell growth (A) Cell count (B) Cell distribution

Figure 12 shows the effect of increasing the probability of symmetric division (or reducing asymmetric division). As explained earlier, a CSC cell has a chance of either symmetric division producing two new CSCs or asymmetric division producing one CSC and one non-CSC. Therefore, the higher probability of symmetric division would increase the ratio of CSCs to total cells. This is highlighted by the figure, especially at 100% probability where effectively all the cells are CSCs.

Figure 12 – Impact of Probability of Symmetric Division on cell growth

Figure 13 shows the effect of increasing the probability of spontaneous death for non-CSC cells. Non-CSC cells have two mechanisms for cell death: one is attributed to the proliferation capacity, p_max, which basically specifies the lifetime of the non-CSC cell, and the other mechanism is spontaneous death with the probability defined by prob_death. In this case, the cell is removed immediately. The number of total cells decreases with increasing death rate (prob_death), yet the number of CSCs remains fairly constant. The non-CSC death rate should not have any noticeable impact on CSCs, although one could speculate that CSCs would more easily propagate in a less dense environment. This may be a minor effect and might not be noticeable for these particular sets of parameters.

Figure 13 – Impact of Probability of Death on cell growth

Lastly, we examine the effect of p_max, the proliferation capacity of the non-CSCs. Table 3 shows the impact of p_max on the ratio of CSCs to total cell count.

Table 3 – Impact of p_max on CSC% (averaging over 10 simulation results)

We can see from the table that the CSC ratio to total cell count is dependent on the non-CSC proliferation capacity. This is expected since p_max = 1 means that non-CSCs only live one cycle before dying rather than the 10 cycles that p_max = 10 indicates. Note, I only averaged 10 simulation results for the number of CSCs and total cells. This is due to the longer simulation time needed for N=600 (steps) and p_max=10.

Discussion

I was successfully able to create a three-dimensional Cellular Automata Cancer Stem Cell model using MATLAB and App Designer. The simulation results seem to be consistent with current Cancer Stem Cell theory. Through the analysis of the simulation results I gained a much better understanding of the mechanisms involved in cancer cell proliferation and tumor growth.
From the papers I have read regarding Cancer Stem Cells, targeting these CSCs is significant if we hope to destroy tumors at their ‘roots’. This is because CSCs have the ability of long-term survival and self-renewal which allow them to possibly regenerate whole tumors if not eliminated. By creating a three-dimensional Cellular Automata Cancer Stem Cell model, I have made a major step in my understanding of some of the parameters involved. I hope to do further study and enhance my MATLAB model in the future.
As previously mentioned, a few things that I would have liked to have fixed or enhanced in my MATLAB code such as calculating the true tumor density or using spinners for the entering probabilities. This would have helped me visualize the effect of migration probability more readily. Initially, I thought something was wrong with my code, as I could not understand why everything looked the same when I stepped through the probabilities. It was only when I re-plotted everything during debugging when I realized what was happening. I would also like to create a version that will plot the grid space step by step so that I could visualize the evolution of the tumor over time. I actually tried to do this and create a separate “animation” option to do that, but it may be beyond my present MATLAB skills. I do plan on doing that in the near future.
On the cancer side, I would like to use Design of Experiment (DOE), where I could construct a set of experiments where I can see the interactions between parameters. In this current study, I basically developed the source code and integrated a user interface. I used it to model CSC and cell growth by varying one parameter at a time and running simulations. With DOE, in a future experiment, I can vary several parameters at a time and maybe I can see how the different parameters interact and under which conditions the growth of CSCs can be minimized. I would like to do this in a follow-up paper after I enhance my CSC model. Some ways I can do this is to incorporate several other factors such as the probability of a CSC migrating to another distant site (since this a major factor in being able to “target and destroy” the CSCs),^[16] the type of parameter to model cancer treatment like chemotherapy or radiation, etc. Another way to enhance my model is to correlate the model to actual tumor growth data. The parameter values I used based on the papers I found roughly assume that nStep of 1 is equivalent to a single day. To create a table of parameters for various different types of cancer cells would be ideal since it appears that certain types of cancer, like melanoma, may fall outside the CSC model. Lastly, I would also like to incorporate angiogenesis and metastasis into the model to make it more accurate. To do this, however, I would initially have to take in the location of the tumor, which may affect other parameters. I would have to assign a time (1 hour, 1 day, etc.) to the time it takes for one simulation step. Then I would have to make angiogenesis happen at a random time between selected intervals, and calculate the probability of metastasis after that. The probability of metastasis would definitely be affected by the location of the tumor. The location of the next tumor would be difficult to determine without knowing the previous location and how tumor cells circulate.
Overall, this 3D Cellular Automata Cancer Stem Cell model we made using an accurate service for 3D scanning is successful. It accurately simulates cancer tumor growth that takes into account certain types of Cancer Stem Cell (CSC) behavior. The results from altering each variable align with the expected results.

References

1. Poleszczuk, Jan, Philip Hahnfeldt, and Heiko Enderling. \”Evolution and Phenotypic Selection of Cancer Stem Cells.\” PLoS Computational Biology 5, no. 1, 2015.
2. Boondirek, Ankana, Wannapong Triampo, and Narin Nuttavut. \”A Review of Cellular Automata Models of Tumor Growth.\” International Mathematical Forum, 2010.
3. Monteagudo, Angel, and José Santos. \”Cancer Stem Cell Modeling Using a Cellular Automaton.\” International Work-Conference on the Interplay Between Natural and Artificial Computation. 2013.
4. Monteagudo, Angel, and José Santos. Treatment Analysis in a Cancer Stem Cell Context Using a Tumor Growth Model Based on Cellular Automata. July 15, 2015.
5. Azam, Nouman. MATLAB Programming and Problem Solving: Go from Beginner to Pro., https://stackskills.com/. 2018.
6. Wolfram, Stephen. A New Kind of Science (free online!). Wolfram Media, Inc., 2002.
7. Wolfram, Stephen. \”Cellular Automata as models of complexity.\” Nature (Macmillan Journals Ltd.) 311, no. 5985, 1984.
8. University of Michigan. \”More evidence that melanoma does not conform to the cancer stem cell model.\” ScienceDaily, November 2010.
9. Baker, Monya. \”Cancer stem cells, becoming common.\” Nature, December 2008.
10. Monteagudo, Angel, and José Santos. \”Cancer Stem Cell Modeling Using a Cellular Automaton.\” International Work-Conference on the Interplay Between Natural and Artificial Computation. 2013.
11. Poleszczuk, Jan, and Heiko Enderling. \”A High-Performance Cellular Automaton Model of Tumor Growth with Dynamically Growing Domains.\” Applied Mathematics 5, 2014.
12. Li, Fuhai, et al. \”A 3D multiscale model of cancer stem cell in tumor development.\” BMC Systems Biology, December 2013.
13. Poleszczuk, Jan, Philip Hahnfeldt, and Heiko Enderling. \”Evolution and Phenotypic Selection of Cancer Stem Cells.\” PLoS Computational Biology 5, no. 1, 2015.
14. Enderling, Heiko, Alexander R. A. Anderson, Mark A. J. Chaplain, Afshin Beheshti, Lynn Hlatky, and Philip Hahnfeldt. \”Paradoxical Dependencies of Tumor Dormancy and Progression on Basic Cell Kinetics.\” Cancer Research 69, no. 22, 2009.
15. University of Michigan. \”More evidence that melanoma does not conform to the cancer stem cell model.\” ScienceDaily, November 2010.
16. Yusuke Shiozawa, 1 Biao Nie, Kenneth J. Pienta, Todd M. Morgan, and Russell S. Taichman. \”Cancer Stem Cells and their Role in Metastasis.\” Pharmacology & Therapeutics, 2013.

About the Author

Katelyn Chu is a STEM-focused junior at St. Francis High School in Mountain View, California. She aspires to enter the field of cancer research in the future, but currently enjoys a wide variety of subjects, including computer science, mathematics, biology, chemistry, and physics.