PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems

Abstract

Many multicellular systems problems can only be understood by studying how cells move, grow, divide, interact, and die. Tissue-scale dynamics emerge from systems of many interacting cells as they respond to and influence their microenvironment. The ideal "virtual laboratory" for such multicellular systems simulates both the biochemical microenvironment (the "stage") and many mechanically and biochemically interacting cells (the "players" upon the stage). PhysiCell-physics-based multicellular simulator-is an open source agent-based simulator that provides both the stage and the players for studying many interacting cells in dynamic tissue microenvironments. It builds upon a multi-substrate biotransport solver to link cell phenotype to multiple diffusing substrates and signaling factors. It includes biologically-driven sub-models for cell cycling, apoptosis, necrosis, solid and fluid volume changes, mechanics, and motility "out of the box." The C++ code has minimal dependencies, making it simple to maintain and deploy across platforms. PhysiCell has been parallelized with OpenMP, and its performance scales linearly with the number of cells. Simulations up to 105-106 cells are feasible on quad-core desktop workstations; larger simulations are attainable on single HPC compute nodes. We demonstrate PhysiCell by simulating the impact of necrotic core biomechanics, 3-D geometry, and stochasticity on the dynamics of hanging drop tumor spheroids and ductal carcinoma in situ (DCIS) of the breast. We demonstrate stochastic motility, chemical and contact-based interaction of multiple cell types, and the extensibility of PhysiCell with examples in synthetic multicellular systems (a "cellular cargo delivery" system, with application to anti-cancer treatments), cancer heterogeneity, and cancer immunology. PhysiCell is a powerful multicellular systems simulator that will be continually improved with new capabilities and performance improvements. It also represents a significant independent code base for replicating results from other simulation platforms. The PhysiCell source code, examples, documentation, and support are available under the BSD license at http://PhysiCell.MathCancer.org and http://PhysiCell.sf.net.

Fig 1. PhysiCell and multiple time scales.

PhysiCell uses BioFVM to update the microenvironment at the short green tick marks, corresponding to Δt_diff. It updates cell mechanics (including cell position) less frequently at the medium black tick marks (Δt_mech), and it runs the cell volume and cycle/death models least frequently at the long red tick marks (Δt_cell). Note that the time steps shown are for illustrative purpose; the default step sizes are given in Time steps.

Fig 2. Hanging drop spheroid (HDS) simulations with deterministic necrosis (left) and stochastic necrosis (right), plotted at 4, 8, and 16 days.

Videos are available at S1 and S2 Videos. Legend: Ki67+ cells are green before mitosis (K₁) and magenta afterwards (K₂). Pale blue cells are Ki67- (Q), dead cells are red (apoptotic) and brown (necrotic), and nuclei are dark blue. Bottom: Hanging drop spheroid experiment (HCC827 non-small cell lung carcinoma) showing a similar necrotic core microstructure. PhysiCell is the first simulation to predict this structure arising from cell-scale mechanical interactions. Image courtesy Mumenthaler lab, Lawrence J. Ellison Center for Transformative Medicine, University of Southern California.

Fig 3. HDS growth.

Left: The deterministic and stochastic necrosis models both give approximately linear growth (left), but the HDS with deterministic necrosis model grows faster (∼ 5% difference in diameter at day 18). Right: The HDS with stochastic necrosis has fewer cells than the deterministic model (∼ 26% difference in cell count at day 18), due to its delay in necrosis. The difference in cell count is larger than the difference in tumor diameter because most of the difference lies in the number of necrotic cells, and necrotic cells are smaller than viable cells.

Fig 4. HDS computational cost scaling.

Left: Wall-time vs. cell count for the stochastic (red) and deterministic (blue) necrosis necrosis models on a single HPC compute node. Both models show approximately linear cost scaling with the number of cell agents. Right: Wall time vs. cell count for stochastic necrosis model on the desktop workstation (orange) and the single HPC node (green).

Fig 5. Ductal carcinoma in situ (DCIS) simulations with deterministic necrosis (left) and stochastic necrosis (right), plotted at 10 and 30 days (multiple views).

Videos are available at S3 and S4 Videos. The legend is the same as Fig 2.

Fig 6. DCIS growth.

The deterministic and stochastic necrosis models both result in linear DCIS growth at approximately 1 cm/year (left), even while their cell counts differ by 21% by the end of the simulations (right).

Fig 7. “Biorobots” example.

Director cells (green) release a chemoattractant c₁ to guide worker cells (red). Cargo cells (blue) release a separate chemoattractant c₂. Unadhered worker cells chemotax towards ∇c₂, test for contact with cargo cells, form adhesive bonds, and then pull them towards the directors by following ∇c₁. If c₁ exceeds a threshold, the worker cells release the cargo and return to seek more cargo cells, repeating the cycle. Bar: 200 μm. A video is available at S5 Video.

Fig 8. Anti-cancer “biorobots” example.

By modifying the worker cells in the previous example (Fig 7) to move up hypoxic gradients (along −∇pO₂) and drop their cargo in hypoxic zones, we can deliver cargo to a growing tumor. In this example, the cargo cells secrete a therapeutic that induces apoptosis in nearby tumor cells, leading to partial tumor regression. Bar: 200 μm. A video is available at S6 Video.

Fig 9. Cancer heterogeneity example.

Each cell has an independent expression of a mutant “oncoprotein” p (dimensionless, bounded in [0, 2]), which scales the oxygen-based rate of cell cycle entry. Blue cells have least p, and yellow cells have most. Initially, the population has normally distributed p with mean 1, standard deviation 0.3, and a “salt and pepper” mixed spatial arrangement. The proliferative advantage for cells with higher p leads to selection and enrichment of the most yellow cells. Stochastic effects lead to emergence of fast-growing foci and a loss of tumor symmetry. Bar: 200 μm. A video is available at S7 Video.

Fig 10. Cancer immunology example.

In this 3-D example, each tumor cell secretes an immunostimulatory factor, and its immunogenicity is modeled as proportional to its mutant oncoprotein expression. (See the previous example in Fig 9.) After 14 days, red immune cells perform a biased random walk towards the immunostimulatory factor, test for contact with cells, form adhesions, and attempt to induce apoptosis for cells with greater immunogenicity. The immune cells successfully attack the tumor initially, leading to partial regression; apoptotic cells are cyan. But strong homing towards gradients of the immunostimulatory factor causes immune cells to “pass” some cells at the outer edge, leading to tumor regrowth. Eventually, immune cells leave the necrotic regions and press their attack on the tumor. This highlights the importance of stochasticity in immune cell movement in mixing with the tumor cells for a more successful immune response. A video is available at S8 Video.