You must login before you can run this tool.
Introduction: a synthetic multicellular system to deliver a therapeutic "cargo"
This app demonstrates how simple "rules" for chemical- and contact-based communication could be used for active delivery of therapeutic compounds, while overcoming traditional biotransport limits and using hypoxia to target delivery. It is based on a published PhysiCell example in Ghaffarizadeh et al. (2018) . In this model, cancer cells (green) rapidly divide while consuming oxygen, which drives the emergence of hypoxic gradients. Greater availability of oxygen drives faster proliferation, while low oxygen can trigger necrotic death, creating a necrotic core (brown cells) in the center of the tumor. After initial growth, we "inject" a multicellular synthetic therapy.
In this (for now theoretical) synthetic therapy, "worker cells" (red) search for "cargo cells" (blue), form focal adhesions, and haul their cargo towards hypoxic (oxygen-depleted) regions in the tissue. Once they reach a region of sufficiently low oxygen, they release their cargo and resume their initial search for more cargo. Released cargo cells secrete a chemotherapeutic compound that can kill nearby tumor cells. (Cancer cells are shaded according to their level of therapy-induced damage.)
This model and cloud-hosted demo are part of a course on computational multicellular systems biology created and taught by Dr. Paul Macklin in the Department of Intelligent Systems Engineering at Indiana University. It is also part of the education and outreach for the IU Engineered nanoBIO Node and the NCI-funded cancer systems biology grant U01CA232137. The models are built using PhysiCell: a C++ framework for multicellular systems biology .
All chemical signals move by chemical diffusion in the simulated environment, using BioFVM  to solve the reaction-diffusion equations. A key property is the diffusion length scale. Diffusion (with a diffusion coefficient D) helps spread a signal over large distances, while decay (and uptake) (with coefficient λ) eliminates the signal to slow its spread. These effects compete to determine the characteristic distance Lthat a chemical signal travels.
That length scale Lis given by L = sqrt( D/λ ).
Choosing D and λ can help to tune the chemical communication distance in these models.
This is released by the computational boundary (a simplified model of a far-field vasculature), and consumed by tumor cells.
We use a far-field value of 38 mmHg (about 5%), which is a typical physioxic condition. See  and  for more ferences. We use a length scale of 1000 μm, and cancer cell uptake gives a length scale of 100 μm in densely packed regions.
This factor is released by undocked cargo cells to help guide chemotaxis by worker cells.
The chemoattractant is nondimensionalized so that 0 ≤ a ≤ 1, and we initially assign a length scale of 100 μm.
This factor is released by released cargo cells to kill cancer cells.
The therapeutic is nondimensionalized so that 0 ≤ c ≤ 1, and we initially assign a length scale of 80 μm.
We use agent-based models to explore multicellular systems like this one. Each cell is modeled as a software agent with an independent state, and its own rules to change its behavior based on local environmental conditions and communication. We use PhysiCell as our modeling platform .
Key features in the agent-based models:
Birth and death:
Cells can proliferate (divide into two daughter cells of half size), with division rates regulated by cell rules. In our model, cells continuously work to grow towards a "target" mature volume.
Note that birth and death are stochastic events for each cell agent: if such a process occurs at rate r, then between now (t) and soon (t + Δt), the probability of the event occurring for that agent is rΔt.
Secretion and uptake:
Cells can secrete chemical factors, or they can remove them (i.e., consume or uptake). This can contribute to gradients in chemical factors, and it's a key method of communication between cells. Moreover, cells can "sample" the chemical state of their surrounding environment. Secreting is sending a signal. Uptaking and sampling is receiving a signal.
Motile cells perform biased migration. After traveling for Tpersistence time (the mean persistence time), they choose a new migration direction dmigrate. Based on the environmental conditions, the cell chooses a directional bias (intended direction) dbias and a random unit vector r. The motile direction dmigrate is then determined by according to:
dmigrate = b dbias + (1-b) r.
Here, b is a bias parameter (between 0 and 1) that determines how strongly biased migration is towards dbias. If b =1, then migration is deterministic towards dbias. If b = 0, then migration is completely random (Brownian).
Once the migration direction is chosen, it is normalized, and multiplied by the cell's migration speed.
Cell agents can stick to one another within a prescribed interaction distance (some multiple of their radius), and they can exert a pushing force on neighbors. PhysiCell  uses potential functions to implement these simple mechanics. Notably, PhysiCell is an off-lattice model, meaning that cells can have variable sizes, and can move freely without grid artifacts.
They aren't required to move some prescribed number of spatial steps. They also can divide without checking for an open neighbor "site". Instead, they can divide and push their neighbors out of the way.
PhysiCell allows custom functions for any cell agents. We use this to model focal cell-cell adhesions as spring-like forces. If two cells i and j (with positions xi and xj) have a focal adhesion, then cell i is adhered to cell j with a spring-like force that we model as:
Fij = ε (xj - xi),
where ε is the elastic coefficient. Because PhysiCell has an inertialess formatulation (see the method paper in ), this is added directly to the cell's velocity, and so ε has units min-1. To model the finite deformability of cells and these elastic focual adhesions, we detach the cells if their distance (||xj - xi||) exceeds a set maximum distance rmax.
Worker cells (red)
Worker cells do not divide, but they can be assigned an apoptotic death rate to model normal "wear and tear." We model workers as non-adhesive except by spring-like focal adhesions; see focal adhesions. They are made more "repulsive" (by default, fivefold more repulsive) to allow them to more readily move through dense cellular regions. Worker cells consume oxygen, but at a reduced rate (by default, 10% of the rate of cancer cells).
Unattached worker cells use the following as their migration bias direction:
dbias = ∇a,
where a is the chemoattractant released by unattached cargo cells. Note that ∇ is the gradient operator that points towards higher concentrations of a chemical factor.
Unattached worker cells test for near contact with unattached cargo cells. If they detect an unattached cargo cell with active receptor R (R ≈ 1), they form a focal adhesion (see focal adhesions).
Attached worker cells change their migration bias direction to:
dbias = -∇σ,
where σ is the oxygen.
Cargo cells (blue)
Cargo cells do not divide, but they can be assigned an apoptotic death rate to model normal "wear and tear." We model cargo as non-adhesive except by spring-like focal adhesions; see focal adhesions. They are made more "repulsive" (by default, fivefold more repulsive) to allow them to more readily move through dense cellular regions. Cargo cells consume oxygen, but at a reduced rate (by default, 10% of the rate of cancer cells). Note that all cargo cells are non-motile.
Initially, all cargo cells express a receptor R with dimensionless value 1, which is used to control their chemical secretions. If R = 1, they secrete chemoattractant a to help guide worker cells. If R = 0 and they are unattached, then they secrete the therapeutic c.
Cargo cells permanently deactivate their receptor (set R = 0) once they are attached to a worker cell. They detach from worker cells if σ < σrelease.
Cancer cells (green)
Cancer cells divide at a rate proportional to oxygen availability, and they experience necrotic death at a rate that increases as oxygen decreases. We use the default parameter values from  for this sample model; they cannot be adjusted by users in the current demo version. Note that all cancer cells are non-motile in this example.
Beyond oxygen-driven birth and hypoxia-driven death, cancer cells are have a basic damage-repair model. Exposure to the chemotherapeutic compound c causes damage (ω, which can be repaired. We simulate this as
dω ⁄ dt = rdamage c - rrepair ω,
where c is sampled at the cell's position. Cancer cells increase their apoptotic death rate by their therapeutic death rate dtherapy, which we model as
dtherapy = ω ⁄ ωmax dmax ,
where dmax is the maximum death rate that is reach when damage (ω) reaches ωmax.
In this model, we set ωmax by solving the damage model to find the steady state damage at a maximum value of c (which is 1):
ωmax = rdamage ⁄ rrepair .
Modify the parameters in the "config basics" and "user params" tabs. Click the "run" button once ready.
To view the cell plots, click the "cell plots" tab, and slide the bar to advance through simulation frames. Note that as the simulation runs, the "max" field (maximum frame number) will increase, so you can view more simulation frames.
To view the substrate fields, click the "substrate plots" tab, choose a substrate from the drop-down menu, and slide through the saved times. Note that as the simulation runs, the "max" field (maximum frame number) will increase, so you can view more simulation frames.
Note that you can download full simulation data for further exploration in your tools of choice.
This software is powered by PhysiCell [1-2], a powerful simulation tool that combines multi-substrate diffusive transport and off-lattice cell models. PhysiCell is BSD-licensed, and available at:
- GitHub releases: https://github.com/MathCancer/PhysiCell/releases
- SourceForge downloads: https://sourceforge.net/projects/physicell/
It is a C++, cross-platform code with minimal software dependencies. It has been tested and deployed in Linux, BSD, OSX, Windows, and other environments, using the standard g++ compiler.
The Jupyter-based GUI was auto-generated by xml2jupyter , a technique to create graphical user interfaces for command-line scientific applications.
To learn more about our work, please visit MathCancer.org.
Randy Heiland, Research Associate, Intelligent Systems Engineering, Indiana University.
- Lead software architect
Paul Macklin, Ph.D. , Associate Professor, Intelligent Systems Engineering, Indiana University.
- Scientific lead, developed PhysiCell simulation model
The following undergraduate students have contributed to this project:
Jupyter notebook and UI development
- Eric Bower (IU Intelligent Systems Engineering, B.S., 2018-present)
- Daniel Mishler (IU Intelligent Systems Engineering, B.S., 2018-present)
- Tyler Zhang (IU Intelligent Systems Engineering, B.S., 2018)
- NSF EEC-1720625. Network for Computational Nanotechnology - Engineered nanoBIO Node
- Breast Cancer Research Foundation
- Jayne Koskinas Ted Giovanis Foundation for Health and Policy
- NIH U01CA232137
 Ghaffarizadeh A, Heiland R, Friedman SH, Mumenthaler SM, Macklin P (2018) PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems. PLoS Comput Biol 14(2): e1005991. https://doi.org/10.1371/journal.pcbi.1005991
 Ghaffarizadeh A, Friedman SH, Macklin P (2016) BioFVM: an efficient, parallelized diffusive transport solver for 3-D biological simulations. Bioinformatics 32(8):1256-8. https://doi.org/10.1093/bioinformatics/btv730
 Heiland R, Mishler D, Zhang T, Bower E, Macklin P (2019, in preparation) Xml2jupyter: Mapping parameters between XML and Jupyter widgets. J Open Source Software.
Cite this work
Researchers should cite this work as follows:
This tool is still in preparation for submission. For now, please cite:
Ghaffarizadeh A, Heiland R, Friedman SH, Mumenthaler SM, Macklin P (2018) PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems. PLoS Comput Biol 14(2): e1005991. https://doi.org/10.1371/journal.pcbi.1005991