Voronoi Tessellation Cellular Automaton Models: From Tumor Simulation to Clinical Translation

Abstract

This article comprehensively reviews the development, application, and clinical potential of Voronoi Tessellation-based Cellular Automaton (VTCA) models in oncology. We explore the foundational principles that make Voronoi geometry ideal for simulating isotropic tumor microenvironments, moving beyond traditional cubic lattices. The methodological core details how agent-based models integrate radiobiological rules—from the linear-quadratic model for radiation response to spatial simulation of oxygen gradients and immune cell interactions. We address key computational challenges, including model optimization and the emerging solution of auto-differentiable Voronoi tessellation for inverse problems. Finally, the article validates these in silico frameworks against experimental 3D spheroids, organoids, and clinical multi-omics data, highlighting their predictive power for therapy response and their role in pioneering personalized, digital twin approaches for cancer treatment.

The Geometric and Biological Rationale for Voronoi Tessellation in Tumor Modeling

The modeling of biological systems, particularly tumor growth, demands computational frameworks that can faithfully recapitulate the complex, adaptive nature of living tissues. For decades, cellular automaton (CA) models have served as a powerful tool for simulating these dynamic processes. Traditionally, many of these models have relied on simple cubic lattices to represent the spatial organization of cells. However, the inherent geometric rigidity and directional bias of these regular grids often fail to capture the isotropic nature and structural heterogeneity of real biological environments.

The adoption of Voronoi tessellations marks a significant paradigm shift in biological modeling. A Voronoi diagram partitions a space into regions (cells) based on the distance to a set of pre-defined seed points. Each cell consists of all points in the space closer to its seed than to any other, naturally generating a space-filling mosaic of convex, polyhedral units [1]. This geometry provides a computationally tractable yet biologically realistic scaffold that overcomes the limitations of cubic lattices, enabling more accurate simulations of critical phenomena like tumor growth, invasion, and clonal evolution.

This note details the quantitative advantages, practical implementation protocols, and specific applications of Voronoi-based cellular automaton models in oncological research, providing a comprehensive toolkit for scientists and drug development professionals.

Quantitative Advantages of Voronoi Tessellations

The transition from cubic to Voronoi lattices in CA models is supported by measurable improvements in both geometric realism and computational outcomes. The table below summarizes the key comparative advantages.

Table 1: Quantitative Comparison of Lattice Types in Tumor Growth Models

Feature	Cubic Lattice	Voronoi Tessellation	Biological & Computational Impact
Lattice Isotropy	Anisotropic; fixed neighborhood connectivity (e.g., 4 or 8 in 2D) [2]	Isotropic; avoids directional biases [3]	Produces realistic, non-directional tumor morphologies; eliminates artificial growth patterns
Cell Packing	Uniform, rigid grid; all cells have identical shape and volume [2]	Heterogeneous, polytopal cells of variable size and shape [1] [4]	Mimics natural tissue organization where cells are not identical cubes
Mechanical Realism	Limited representation of cell-cell interactions and pressures	Allows for modeling of adhesive and repelling forces on cell borders [4]	Enables simulation of mechanical confinement and tumor-stroma interactions
Simulated Tumor Growth	Can produce artificially constrained, spherical morphologies	Reproduces fingered, invasive morphologies under hypoxia [5]	Captures critical clinical growth patterns linked to poor prognosis
Scalability & Adaptivity	Fixed resolution; difficult to simulate growth over multiple scales	Compatible with adaptive grids; simulates growth from ~1,000 to millions of cells [3]	Allows study of tumor development from microscopic to macroscopic scales

The biological impact of these geometric properties is profound. For instance, the isotropy of the Voronoi tessellation is crucial for simulating the invasive fingering morphology of aggressive tumors like glioblastoma, a pattern that emerges in models under low oxygen conditions [5]. This morphology is often missed in cubic lattice models, which can artificially constrain growth due to their directional bias.

Core Protocols for Implementing a Voronoi-Based Tumor Model

This section provides a detailed methodology for developing a hybrid cellular automaton model of solid tumor growth using a Voronoi tessellation as the underlying scaffold.

Protocol: Initializing the Voronoi Scaffold and Cell Population

Objective: To generate the initial 3D spatial geometry and populate it with tumor cells equipped with a micro-environment response network.

Materials:

• Computational Environment: Programming language (e.g., C++, Python) with numerical and geometric libraries (e.g., SciPy, Voro++).
• Seed Points: A set of generator points within the simulation domain. These can be randomly distributed or arranged to mimic specific tissue structures.

Workflow:

• Generate Voronoi Tessellation:

  • Define a set of seed points, ( S = {s1, s2, ..., s_n} ), within the 3D simulation domain. The initial density of seeds determines the initial resolution.
  • Compute the Voronoi diagram, which partitions the domain into cells ( Rk ), where ( Rk = { x \in X \mid d(x, Pk) \leq d(x, Pj) \; \text{for all} \; j \neq k } ) [1]. Here, ( d ) is the Euclidean distance function.
  • The dual structure, the Delaunay triangulation, defines the connectivity and neighborhood relations between cells, which is essential for modeling cell-cell interactions and nutrient diffusion [3].

• Define the Founding Cell Population:

  • Assign a single tumor cell to a subset of the Voronoi cells. Each cell is an independent agent.
  • Equip each cell with an Internal Response Network, modeled, for example, by a feed-forward artificial neural network [5].
  • Inputs: Micro-environmental variables (e.g., oxygen, glucose, growth factor concentrations) sampled at the cell's location.
  • Outputs: Behavioral phenotypes (e.g., probability of proliferation, quiescence, apoptosis, migration).

• Parameterize Cell Phenotype:

  • The connection weights and thresholds of the neural network constitute the cell's "genotype." These parameters are subject to mutation upon cell division, introducing clonal evolution into the model [5].
  • Initialize the founding cell(s) with baseline parameters that represent a specific cancer type.

Figure 1: Workflow for initializing a Voronoi-based tumor growth model. This diagram outlines the key steps, from generating the geometric scaffold to defining the initial cell population and its properties.

Protocol: Simulating Growth and Evolutionary Dynamics

Objective: To simulate the time-evolution of the tumor system, including nutrient diffusion, cell fate decisions, division, and mutational events.

Materials:

• Reaction-Diffusion Solver: For calculating the spatiotemporal distribution of nutrients and metabolites.
• Stochastic Update Algorithm: To handle probabilistic cell behaviors and mutations.

Workflow:

• Compute Micro-environment:

  • At each time step ( t ), solve a system of reaction-diffusion equations to update the concentration fields of key substrates (e.g., oxygen, glucose).
  • A common simplification is: ( \frac{\partial C}{\partial t} = D\nabla^2C - \lambda C \sum \delta(\vec{x} - \vec{x}_i) ) where ( C ) is concentration, ( D ) is diffusion coefficient, ( \lambda ) is consumption rate, and the sum is over all consuming cells [5].

• Update Cell States:

  • For each cell, sample the local micro-environment and process it through its internal response network.
  • Execute the determined behavior:
    • Proliferation: If a cell divides, place the daughter cell in a neighboring Voronoi site. If no site is free, apply a pressure-based displacement rule to surrounding cells [3] [4].
    • Quiescence: The cell remains alive but does not cycle.
    • Apoptosis/Necrosis: The cell is removed from the simulation, freeing its location.

• Implement Mutations and Clonal Selection:

  • Upon division, with a defined probability, introduce mutations by slightly altering the weights of the daughter cell's neural network [5].
  • Cells with genotypes better suited to the local environment (e.g., lower apoptotic response to hypoxia) will have a selective advantage, leading to clonal expansion.

• Adapt the Scaffold (Optional):

  • For models simulating massive growth, an adaptive grid can be used where the density of Voronoi seeds increases in regions of high cellular activity, maintaining computational efficiency [3].

Figure 2: The core simulation loop of a hybrid Voronoi-CA model, showing the integration of continuous field calculations with discrete cell-agent updates.

The Scientist's Toolkit: Research Reagent Solutions

The following table details key computational "reagents" and tools essential for building and analyzing Voronoi-based tumor models.

Table 2: Essential Research Reagents and Computational Tools

Tool / Reagent	Type	Function in the Model	Exemplary Implementation
Voronoi/Delaunay Library	Software Library	Generates the geometric scaffold and neighborhood topology.	scipy.spatial.Voronoi & Delaunay in Python; Voro++ for C++.
ODE/PDE Solver	Numerical Library	Solves equations for nutrient diffusion and consumption.	Finite Difference/Volume methods; SUNDIALS library; FiPy.
Artificial Neural Network (ANN)	Algorithmic Unit	Maps a cell's micro-environment to its behavioral fate.	A simple multi-layer perceptron with inputs (nutrient levels) and outputs (fate probabilities).
Genetic Algorithm (GA) Optimizer	Optimization Tool	Calibrates model parameters or optimizes scaffold properties.	Used to evolve Voronoi seed points for desired mechanical properties [6].
Parameter Set (Genotype)	Data Structure	Defines a cell's phenotypic response. A vector of weights and thresholds in the ANN.	Mutated during division to drive clonal evolution [5].

Application Note: Simulating Hypoxia-Driven Tumor Invasion

Background: Tumor hypoxia is a key factor driving both morphological change and the selection of aggressive cellular phenotypes. Voronoi-based models are uniquely suited to study this phenomenon.

Procedure and Outcomes:

• Model Setup: Implement a 3D Voronoi CA model as described in Protocols 3.1 and 3.2. Initialize the simulation with a well-oxygenated tissue environment and a small cluster of tumor cells.
• Induction of Hypoxia: Allow the tumor to grow until its size leads to diffusion-limited oxygen supply, naturally creating a hypoxic core.
• Observation: The model will recapitulate the classic layered structure of tumors: a necrotic core, a shell of quiescent cells, and an outer rim of proliferating cells [3] [5].
• Analysis of Evolutionary Dynamics:
  • Under low oxygen conditions, cells with mutations that confer a survival advantage (e.g., reduced apoptosis, glycolytic phenotype) will be selected.
  • This selection pressure leads to the emergence of dominant clones with aggressive traits.
  • At the tissue level, this manifests as a transition from a smooth, rounded tumor boundary to an irregular, fingered morphology as these aggressive clones invade the surrounding tissue [5].

Conclusion: The simulation demonstrates that hypoxia is not merely a passive state of low oxygen but an active driver of tumor evolution and invasion. The Voronoi geometry is critical in this outcome, as its isotropic and heterogeneous nature does not artificially constrain the emergence of these complex invasive patterns, unlike cubic lattices.

The adoption of Voronoi tessellations in cellular automaton models represents a significant leap forward in the quest for biologically realistic in silico representations of tumor growth. By moving beyond the geometric constraints of cubic lattices, Voronoi-based models provide a physically and spatially plausible scaffold that naturally captures the isotropic, heterogeneous, and adaptive nature of living tissues. The integration of this geometry with hybrid modeling techniques—coupling discrete cell-based rules with continuous nutrient fields—and evolutionary algorithms allows researchers to probe not just the morphological but also the evolutionary dynamics of cancer. This powerful framework offers a robust platform for testing hypotheses about tumor progression, treatment resistance, and the efficacy of novel therapeutic strategies, ultimately accelerating oncological research and drug development.

This document provides application notes and experimental protocols for key computational geometry techniques—Delaunay Triangulation, Isotropic Growth, and Adaptive Grids—within the context of developing a Voronoi tessellation cellular automaton (CA) model for tumor growth research. These principles enable researchers to generate biologically realistic synthetic microstructures, simulate complex tissue environments, and investigate tumor progression dynamics with high spatial accuracy. The methodologies outlined below support the generation of virtual tissues that statistically represent real microstructures, facilitating the study of tumor morphology, evolutionary dynamics, and response to therapeutic interventions.

Core Principles and Mathematical Foundations

Delaunay Triangulation and Voronoi Tessellation

Delaunay Triangulation is a fundamental geometric computation that subdivides a set of points into a network of triangles (or tetrahedra in 3D) possessing specific optimality properties [7].

• Empty Circumcircle Criterion: For any triangle in the triangulation, the circumcircle defined by its three vertices contains no other input points in its interior [7] [8]. This is the defining property of a Delaunay triangulation.
• Dual to Voronoi Diagram: The Delaunay triangulation is the dual graph of the Voronoi diagram [7]. The vertices of the Voronoi diagram (which define cell boundaries) are the circumcenters of the Delaunay triangles. This relationship provides a direct computational pathway from a set of seed points (nuclei) to a Voronoi-tessellated space.
• Max-Min Angle Property: In 2D, the Delaunay triangulation maximizes the minimum angle among all triangles, avoiding artificially slender "sliver" triangles and contributing to numerical stability in subsequent simulations [7] [9].

Flipping Algorithm: A common method to achieve a Delaunay triangulation involves checking "locally Delaunay" edges. An edge shared by two triangles is locally Delaunay if the circumcircle of one triangle does not contain the opposite vertex of the adjacent triangle. If this condition fails, the edge is "flipped" to restore the local property, and repeated application leads to a global Delaunay triangulation [7] [9].

Isotropic Growth on Discrete Grids

In CA models, tumor growth is often simulated by cells proliferating into neighboring grid locations. Isotropic Growth refers to uniform expansion in all directions [10]. However, achieving true isotropy on a Cartesian grid is non-trivial.

• Anisotropy Artifact: When using a standard von Neumann neighborhood (up, down, left, right pixels) for growth, the caliper diameter increases at different rates along different lattice directions. For example, growth along the <100> family of directions (horizontal/vertical) is faster than along the <110> directions (diagonal) [10]. This results in diamond-shaped artifacts instead of circular expansions.
• Correction Strategies: Isotropic growth can be approximated by implementing corrective measures.
  • Stochastic Penalization: A probabilistic acceptance rule for adding a voxel to a growing region can be applied, where the probability is lower for directions that would otherwise grow too quickly (e.g., <100> directions) [10].
  • Lattice-Free Methods: Using a Voronoi tessellation as the underlying grid provides inherent spatial isotropy because the cells are polygonal and not aligned to a fixed Cartesian lattice, thereby avoiding directional bias [3].

Adaptive Grids

Adaptive Grids are computational meshes that can change their resolution dynamically to focus computational resources on regions of interest, such as a growing tumor interface.

• Purpose and Benefit: They allow a model to maintain high resolution at small tumor sizes while still being capable of simulating growth to a large, macroscopic size without prohibitive computational cost [3].
• Implementation in Voronoi/Delaunay Frameworks: The initial model domain can be discretized using a fine Delaunay triangulation (or its dual Voronoi tessellation). As the tumor grows and its boundary evolves, the grid can be adaptively refined in regions with high activity (e.g., the proliferative rim) and coarsened in quiescent or necrotic regions [3].

Quantitative Properties of Voronoi Tessellations

The statistical properties of the virtual tissue are controlled by the initial seed placement. A key parameter is the regularity parameter, α, defined as the ratio of the minimum distance between seeds to the seed distance of a perfectly ordered body-centered cubic lattice [11]. This parameter allows for tuning the microstructure from completely random (α = 0, Poisson Voronoi Tessellation) to perfectly ordered (α = 1).

Table 1: Statistical Properties of 3D Voronoi Tessellations vs. Regularity

Geometric Property	Poisson VT (α = 0)	Intermediate Regularity (e.g., α = 0.5)	Ordered VT (α = 1)
Faces per Cell (Mean)	~15.6	Narrowing distribution	14 (Tetrakaidecahedron)
Edges per Face (Mean)	~5.2	Narrowing distribution	5.15
Cell Volume Distribution	Broad (Gamma distribution)	Narrowing	Dirac delta (single value)
Application	Mimics disorganized tissue	Tunable for specific tissues	Idealized, highly ordered structures

Table 2: Algorithm Comparison for Synthetic Microstructure Generation

Method	Key Principle	Advantages	Limitations	Suitability for Tumor Models
Voronoi Tessellation	Partitions space based on seed points [11].	Fast generation; direct control via seed placement; isotropic cells [3] [11].	Limited morphological control; less realistic grain boundaries [10].	High: Represents well-defined cell territories.
Isotropic Growth-until-Impingement	Grows seeds (e.g., spheres) until they touch [10].	Better control of final grain size distribution; more realistic curved boundaries [10].	Computationally intensive; requires correction for grid anisotropy [10].	Medium: Good for simulating packed cells.
Cellular Automata (CA)	Discrete grid with state transition rules [12].	Intuitive; easy to incorporate cell-level rules (division, death); highly flexible [5].	Can be slow; susceptible to grid bias if not corrected [10] [13].	High: Ideal for simulating emergent behavior from local rules.

Experimental Protocols

Protocol 1: Generating a Voronoi Tessellation for a Virtual Tissue

Objective: To create a synthetic 2D or 3D tissue microstructure with controlled regularity using Voronoi tessellation.

Materials:

• Computational geometry library (e.g., MATLAB delaunayTriangulation, SciPy scipy.spatial.Delaunay, Voro++).
• High-performance computing workstation.

Procedure:

• Define the Domain: Specify a 2D square or 3D cubic domain with periodic boundary conditions to minimize edge effects [11].
• Generate Seed Points:
  • For a Poisson Voronoi Tessellation, generate N seed points with coordinates drawn randomly from a uniform distribution within the domain [11].
  • For a Hard-Sphere Voronoi Tessellation with regularity α, use a random sequential adsorption (RSA) algorithm [11]: a. Place the first seed at a random location. b. For each subsequent seed candidate, check if its distance to all existing seeds is greater than the minimum exclusion distance, d_min = α * d_lattice. c. If the condition is met, accept the seed; otherwise, reject it and generate a new candidate. d. Repeat until the desired number of seeds, N, is placed or no more seeds can be added.
• Compute Delaunay Triangulation: Input the set of seed points into the Delaunay triangulation algorithm to generate a network of triangles (2D) or tetrahedra (3D) [7] [8].
• Derive Voronoi Tessellation: Calculate the circumcenters of all Delaunay triangles/tetrahedra. Connect these circumcenters to form the Voronoi cells, which represent the individual "cells" in your virtual tissue [7].
• Validation: Calculate the distribution of the number of faces per cell and cell volumes. Compare these distributions to known statistics for your target regularity α (see Table 1) to ensure the tessellation is statistically representative [11].

Protocol 2: Simulating Isotropic Tumor Growth using a Cellular Automaton

Objective: To simulate the avascular growth of a tumor from a single cell into a multicellular spheroid, minimizing grid-induced anisotropy.

Materials:

• A pre-defined grid (Cartesian with corrections or Voronoi-based).
• Stochastic CA framework (e.g., custom code or chouca R package [14]).

Procedure:

• Initialize Grid and Tumor Cell: Create a 2D or 3D grid. Designate a single central cell as a live cancer cell. All other cells are initially healthy tissue or empty.
• Define Cellular States:
  • Proliferative: Can divide if conditions are met.
  • Quiescent: Alive but not dividing due to nutrient stress.
  • Necrotic: Dead.
  • Healthy: Normal tissue cells.
• Define Micro-environmental Field: Solve a diffusion equation (e.g., ∇²c - λc = 0) over the grid to establish a nutrient (e.g., oxygen) concentration field, c [5]. The consumption rate λ is higher in tumor cells.
• Set Transition Rules. For each tumor cell, update its state based on its neighborhood (von Neumann or Moore) and local nutrient level [3] [5]:
  • Proliferation: If c > c_prol and an adjacent site is available (healthy or empty), the cell may divide into that site with a probability p_div.
  • Quiescence: If c_nec < c < c_prol, the cell becomes quiescent.
  • Necrosis: If c < c_nec, the cell becomes necrotic.
  • Ensure Isotropy: If using a Cartesian grid, implement a stochastic correction [10]. When a cell divides, the direction for placing the daughter cell is chosen from all possible neighboring directions with a probability that penalizes the <100>-type directions.
• Iterate: For each discrete time step, update the nutrient field and then synchronously apply the transition rules to every cell in the grid.
• Data Collection: At regular intervals, record metrics such as total tumor cell count, necrotic fraction, growth fraction, and tumor morphology [3] [5].

Protocol 3: Incorporating an Adaptive Grid for Large-Scale Simulation

Objective: To integrate an adaptive grid with a Voronoi-CA model to efficiently simulate tumor growth over several orders of magnitude.

Procedure:

• Initial Fine-Scale Tessellation: Begin with a high-resolution Voronoi tessellation of the entire domain, generated via Protocol 1 [3].
• Define Refinement Triggers: Identify criteria that will trigger grid refinement. For tumor growth, this is typically the presence of proliferative tumor cells or a steep nutrient gradient.
• Implement Adaptive Logic:
  • Refinement: When a Voronoi cell is triggered, subdivide it into smaller Voronoi cells. This involves adding new seed points within the parent cell and re-computing the local tessellation [3].
  • Coarsening: In regions far from the tumor or in necrotic cores, merge small Voronoi cells into larger ones to reduce computational load.
• Map Cell States: During refinement, assign the state of the parent cell to all new daughter cells. During coarsening, assign the new coarse cell the state of the majority of the constituent fine cells.
• Run Coupled Simulation: Execute the CA growth rules (Protocol 2) on the adaptive Voronoi grid, dynamically refining and coarsening as the tumor evolves.

Visualization of Workflows and Relationships

Diagram 1: Integrated Voronoi-CA Tumor Model Workflow.

Diagram 2: From Seed Regularity to Tumor Morphology.

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Computational Tools and Their Functions

Tool/Reagent	Function/Description	Example Use Case
Delaunay Triangulation Engine	Computes the optimal connectivity between seed points to form a triangle/tetrahedron network [7] [8].	Generating the foundational geometry for a Voronoi tessellation from a set of nuclei.
Voronoi Tessellation Generator	Partitions space into convex polyhedral cells based on the Delaunay triangulation [7] [11].	Creating a virtual tissue microstructure for initializing a CA model.
Stochastic CA Framework	Provides a computational engine to apply state transition rules probabilistically to a grid of cells [14] [5].	Simulating the probabilistic division, death, and migration of tumor cells.
Nutrient Diffusion Solver	Numerically solves partial differential equations (PDEs) for nutrient diffusion and consumption [5].	Modeling the formation of nutrient gradients that create proliferative, quiescent, and necrotic zones within a tumor.
Adaptive Meshing Library	Dynamically refines and coarsens the computational grid based on user-defined criteria [3].	Maintaining high resolution at the growing tumor boundary while reducing cost in inactive regions.
Regularity Parameter (α)	Controls the minimum spacing between seeds, tuning the virtual tissue from random to ordered [11].	Generating synthetic microstructures that statistically match specific experimental tissue samples.

Application Notes

The integration of radiobiological models with advanced computational frameworks is revolutionizing our understanding of tumor dynamics and radiation response. The linear-quadratic (L-Q) model has served as the cornerstone for predicting radiation-induced cell death, expressed as SF(D) = e^(-αD-βD²), where α represents lethal damage and β represents sublethal damage [15]. This model provides the mathematical foundation for quantifying biological effective dose (BED) in fractionated radiotherapy. However, the L-Q model demonstrates significant limitations at high radiation doses per fraction, necessitating more comprehensive approaches that incorporate the fundamental 5Rs of radiobiology: Repair, Repopulation, Redistribution, Reoxygenation, and Radiosensitivity [16] [17]. These principles collectively govern tumor response to radiation and provide the biological rationale for fractionated treatment schedules.

The Voronoi tessellation cellular automaton (VTCA) framework offers an ideal computational platform for integrating these radiobiological principles into tumor growth and treatment response modeling. Unlike regular lattices, Voronoi tessellation generates an isotropic lattice that better mimics the heterogeneous spatial organization of real tissues, avoiding the directional artifacts common in cubic lattice structures [3]. This approach enables researchers to simulate tumor development from microscopic to macroscopic scales while maintaining computational efficiency and biological relevance.

Integration of the 5Rs into Voronoi Tessellation Models

Table 1: The 5Rs of Radiobiology and Their Implementation in VTCA Models

Radiological Principle	Biological Significance	VTCA Implementation Approach	Key Model Parameters
Repair	Cellular recovery from sublethal DNA damage	Michaelis-Menten enzyme kinetics or delayed differential equations [15]	Repair rate constant, fidelity coefficient
Repopulation	Compensatory cell division between fractions	Cell cycle progression rules with density dependence [18]	Cell cycle duration, growth fraction
Redistribution	Cell cycle phase synchronization post-radiation	Phase-specific radiosensitivity across G1, S, G2, M [18]	Phase duration, checkpoint efficiency
Reoxygenation	Improved tumor oxygenation between fractions	Dynamic oxygen diffusion from vasculature [18] [5]	Oxygen diffusion coefficient, consumption rate
Radiosensitivity	Intrinsic cellular radiation vulnerability	Cell-specific α/β ratios in L-Q formalism [17]	α (Gy⁻¹), β (Gy⁻²) parameters

The Voronoi tessellation provides the spatial foundation for modeling these complex interactions by defining cellular neighborhoods through Delaunay triangulation, which establishes natural neighbor relationships between adjacent cells [3]. This geometric arrangement enables realistic simulation of nutrient diffusion, cell-cell signaling, and invasion patterns that characterize malignant growth. When implementing radiation response, each cell within the Voronoi lattice can be assigned individual radiosensitivity parameters based on its phenotypic state, oxygenation level, and cell cycle position, creating a heterogeneous response landscape that mirrors clinical observations.

The integration of dynamic oxygen mapping represents a particularly powerful application of VTCA models. By assigning each vascular fraction a specific oxygen histogram profile, researchers can simulate the temporal and spatial evolution of hypoxia and its profound impact on radiation efficacy [18]. This approach captures the critical reoxygenation process that occurs between fractions, which significantly influences treatment success by increasing the radiosensitivity of previously hypoxic cells.

Advanced Model Extensions: Clonal Evolution and Immune Interactions

Sophisticated VTCA implementations can incorporate evolutionary dynamics through neural network-based decision systems that determine cellular behavior based on environmental conditions [5]. These models simulate clonal competition by assigning each cell a "phenotype" defined by connection weights and thresholds that are subject to mutation during division. This approach reveals how selection pressures, particularly hypoxia, drive the emergence of aggressive tumor subclones with altered apoptotic thresholds and metabolic preferences, including the classic Warburg effect (glycolytic phenotype) [5].

Emerging research also highlights the potential for integrating immunomodulatory effects of radiation into VTCA frameworks. The combination of radiotherapy with immunotherapy represents a promising frontier, with potential synergistic effects described as '5R × 3E' rather than merely additive '5R + 3E' [16]. This multiplicative relationship stems from radiation's ability to stimulate antigen release and modify the tumor microenvironment, potentially shifting the cancer immunoediting process from equilibrium/escape back toward elimination.

Experimental Protocols

Protocol: Implementing a Voronoi Tessellation Cellular Automaton for Radiobiological Studies

Purpose: To establish a computational framework for simulating tumor growth and radiation response incorporating the 5Rs of radiobiology.

Background: The Voronoi tessellation cellular automaton approach enables multi-scale modeling of tumor dynamics with individual cell tracking and microenvironmental interactions [3]. This protocol outlines the implementation of a VTCA system specifically designed for radiobiological investigations.

Materials and Computational Requirements:

• Software Infrastructure: MATLAB, Python, or C++ with computational geometry libraries (e.g., Voro++ or Qhull)
• Processing Capacity: Multi-core workstation with minimum 16GB RAM for 3D simulations
• Visualization Tools: ParaView, VTK, or custom OpenGL rendering

Procedure:

• Lattice Initialization (Time: 1-2 hours)

  • Generate an initial Voronoi tessellation using random seed points within a defined boundary
  • Calculate the dual Delaunay triangulation to establish neighbor relationships
  • Implement adaptive grid refinement to maintain resolution during expansion
  • Assign each Voronoi cell a type: tumor, endothelial, or normal tissue

• Microenvironment Setup (Time: 30 minutes)

  • Define vascular network positions and oxygenation capacities
  • Implement oxygen diffusion using discrete approximation of Fick's law: ∂[O₂]/∂t = D∇²[O₂] - λ[O₂]
  • Set nutrient consumption rates based on cell type and proliferation status
  • Establish waste product accumulation and acid-base balance rules

• Cell Cycle Implementation (Time: 1 hour)

  • Program cell cycle phases with durations: G1 (14h), S (6h), G2 (3h), M (1h) [18]
  • Implement phase transition checkpoints with Bernoulli probability tests
  • Configure quiescence (G0) entry/exit based on local nutrient thresholds
  • Establish contact inhibition rules based on Voronoi neighbor density

• Radiation Response Module (Time: 2-3 hours)

  • Implement linear-quadratic model with cell-specific α/β parameters
  • Program repair kinetics using Integrated Michaelis-Menten (IMM) or Lambert W0 function approaches [15]
  • Incorporate cell cycle-specific radiosensitivity variations (G2/M most sensitive)
  • Establish hypoxia modification factors using oxygen enhancement ratio (OER) principles

• Model Execution and Data Collection (Time: Variable based on simulation size)

  • Initialize with approximately 1,000 cells representing early tumor development
  • Run simulation through multiple volumetric doublings to macroscopic size
  • Track subpopulation dynamics: proliferative, quiescent, and necrotic fractions
  • Record spatial organization, clonal diversity, and invasion patterns

Troubleshooting:

• Artifact formation: Verify Voronoi tessellation algorithms and boundary conditions
• Unrealistic growth rates: Calibrate cell cycle parameters against experimental data
• Numerical instability: Reduce time step size and verify differential equation solvers
• Memory limitations: Implement data compression for quiescent cell regions

Protocol: Simulating Fractionated Radiotherapy in VTCA Models

Purpose: To investigate tumor response to different radiation fractionation schemes using VTCA modeling.

Background: Fractionated radiotherapy exploits differential repair capacity between tumor and normal tissues [15]. This protocol enables computational comparison of conventional, hypofractionated, and stereotactic radiotherapy schedules.

Materials:

• Established VTCA tumor model from Protocol 2.1
• Radiation dose parameters (2-20 Gy per fraction)
• Fractionation schedules (1-30 fractions)

Procedure:

• Baseline Tumor Growth (Time: 4-6 hours simulation)

  • Develop tumor to approximately 10⁶ cells (pre-vascular stage)
  • Characterize baseline growth fraction, hypoxia pattern, and clonal distribution
  • Record volumetric doubling time and spatial configuration

• Treatment Planning (Time: 30 minutes)

  • Define target volumes based on spatial cell densities
  • Set dose prescription and fractionation schedule
  • Establish normal tissue constraints based on spatial boundaries

• Radiation Delivery Simulation (Time: 1-2 hours per fraction scheme)

  • Implement daily fraction delivery with precise dose deposition
  • Calculate cell survival probabilities using SF = e^(-αD-βD²) with hypoxia modification
  • Execute repair processes between fractions using IMM kinetics [15]
  • Model repopulation dynamics during treatment course

• Response Assessment (Time: 1 hour)

  • Quantify tumor control probability (TCP) for each regimen
  • Analyze normal tissue complication probability (NTCP)
  • Evaluate therapeutic ratio (TCP/NTCP) across fractionation schemes
  • Assess impact on clonal evolution and aggressive phenotype selection

Validation:

• Compare simulation results to clinical dose-response curves
• Verify hypoxic fraction dynamics against experimental measurements
• Validate repopulation kinetics against bromodeoxyuridine labeling data

Signaling Pathways and Experimental Workflows

Diagram 1: VTCA Radiobiological Modeling Workflow. This workflow illustrates the integration of Voronoi spatial structure with microenvironmental factors and radiation response mechanisms.

The Scientist's Toolkit: Research Reagent Solutions

Table 2: Essential Computational and Experimental Resources for Radiobiological Modeling

Category	Item/Technique	Function/Application	Implementation Notes
Computational Tools	Voronoi Tessellation Libraries (Voro++, Qhull)	Spatial partitioning for cellular automata	Enable isotropic lattice generation without directional bias [3]
	Delaunay Triangulation Algorithms	Neighbor identification and network modeling	Dual structure to Voronoi diagram for connectivity mapping [3]
	Linear-Quadratic Model Solvers	Radiation dose-response calculation	Implement with IMM extensions for high-dose accuracy [15]
	Oxygen Diffusion Solvers	Hypoxia pattern simulation	Finite difference methods with consumption terms [5]
Biological Assays	18F-Fluoromisonidazole (18F-FMISO) PET	Hypoxia imaging and validation	Gold standard for clinical hypoxia assessment [18]
	18F-Fluorothymidine (18F-FLT) PET	Proliferation monitoring post-radiation	Quantifies repopulation dynamics during treatment [18]
	Immunohistochemistry (IHC) for pimonidazole	Ex vivo hypoxia validation	Correlative validation for computational predictions
	Flow Cytometry with Cell Cycle Markers	Redistribution assessment	Quantifies phase-specific radiation effects
Model Parameters	α/β Ratios (Tumor: ~10 Gy, Normal: ~3 Gy)	Tissue-specific radiosensitivity	Foundation for BED calculations [15]
	Oxygen Enhancement Ratio (OER)	Hypoxia modification factor	Typically 2.5-3.0 for full radioprotection
	Repair Rate Constants (μ)	Sublethal damage repair kinetics	Michaelis-Menten or exponential forms [15]
	Potential Doubling Time (Tpot)	Repopulation kinetics	Typically 2-10 days for human tumors

This comprehensive framework enables researchers to bridge computational modeling with experimental validation, creating a robust platform for investigating radiobiological principles and optimizing radiation therapy strategies. The integration of Voronoi tessellation approaches with the 5Rs of radiobiology provides unprecedented insights into tumor dynamics and treatment response, advancing the field toward truly personalized radiation medicine.

This application note provides detailed protocols for implementing a Voronoi tessellation-based cellular automaton (CA) model to simulate three critical hallmarks of solid tumor development: DNA damage complexity, hypoxia, and cellular heterogeneity. The content is framed within a broader thesis on advancing in silico tumor modeling, providing researchers and drug development professionals with a reproducible framework for studying tumor dynamics and therapeutic resistance.

The foundational Voronoi cellular automaton model excels at simulating proliferative brain tumor growth that follows Gompertzian kinetics over nearly three orders of magnitude in radius using only four microscopic parameters [3]. This approach models a glioblastoma multiforme (GBM) tumor as an enormous idealized multicellular tumor spheroid, mathematically described by a Gompertz function, which is particularly suited for tumors comprising a large area of central necrosis surrounded by shells of viable cells [3]. The model's key innovation lies in its use of an isotropic Voronoi tessellation that avoids the anisotropies associated with more ordered lattices, while the adaptive grid allows for greater resolution at small tumor sizes while still enabling growth to macroscopic dimensions [3].

Table 1: Core Simulation Parameters for Base Voronoi Tessellation CA Model

Parameter Category	Specific Parameter	Typical Value/Range	Biological Significance
Lattice Parameters	Lattice Type	Voronoi Tessellation (Delaunay Triangulation dual)	Provides spatial isotropy, avoids directional artifacts [3]
	Grid Resolution	Adaptive (1-2 mm³ for macroscopic tumors)	Balances computational efficiency with biological relevance [19]
Growth Parameters	Initial Cell Count	~1,000 cells	Represents early tumor stage [3]
	Growth Law	Gompertzian	Matches observed tumor growth dynamics [3]
	Proliferation Rule	Space- and nutrient-limited	Mimics in vivo growth constraints [3] [5]
Cell Cycle Parameters	Cell Cycle Phases	G1, S, G2, M, G0 (quiescent), Necrotic	Captures essential cell states [3] [5]
	Growth Fraction	Dynamically calculated	Predicts composition (proliferative, quiescent, necrotic fractions) [3]

Protocol 1: Implementing the Core Voronoi Tessellation Framework

Table 2: Essential Research Reagent Solutions and Computational Tools

Item Name	Specification/Function	Application in Protocol
Voronoi Tessellation Library	Computational geometry package (e.g., Voro++, Qhull)	Generates isotropic lattice structure; defines cellular neighborhoods [3]
ODE Solver	Numerical solver for ordinary differential equations (e.g., SUNDIALS, SciPy)	Models diffusion of oxygen, nutrients, and therapeutic agents [19]
Visualization Software	3D visualization toolkit (e.g., VTK, Paraview)	Renders tumor morphology, hypoxic gradients, and clonal distribution [20]
High-Performance Computing (HPC)	Multi-core processor workstation or cluster	Enables simulation of large tissue volumes (~cm³) over physiological timescales [19]
Color Contrast Analyzer	WCAG 2.1 AA compliant tool (e.g., WebAIM Contrast Checker)	Ensures accessibility and interpretability of all output diagrams and visualizations [21] [22]

Step-by-Step Procedure

• Lattice Initialization: Generate a 3D Voronoi tessellation within a defined tissue volume (e.g., 4.2 × 4.2 × 4.2 cm³). The Delaunay triangulation, as the dual lattice, defines nearest neighbors for each cell [3].
• Seed Tumor Cells: Designate approximately 1,000 cells at the center of the lattice as the initial tumor population [3].
• Parameter Assignment: Initialize each cell with the base parameters outlined in Table 1, including cell cycle status and nutrient consumption rates.
• Set Boundary Conditions: Define boundary conditions for the simulation space, typically implementing no-flux conditions at the tissue boundaries.
• Implement Adaptive Grid: Configure the adaptive grid algorithm to maintain high resolution at the tumor boundary while allowing efficient simulation of macroscopic growth [3].

Protocol 2: Modeling Hypoxia and Metabolic Heterogeneity

Background and Biological Rationale

Hypoxia, characterized by reduced oxygen availability (typically 1-2% O₂ in tumor regions), is a hallmark of many solid tumors and drives genomic instability, metabolic rewiring, and therapy resistance [23]. In the Voronoi CA framework, oxygen diffusion from the vascular system creates spatial gradients that determine regional cell behavior. Tumor cells exhibit distinct phenotypic states based on oxygen distribution: proliferative in well-oxygenated regions near vessels, quiescent in intermediate zones, and necrotic in severely hypoxic cores approximately 150 μm from vessels [23] [5].

Quantitative Implementation

Table 3: Parameters for Hypoxia and Metabolic Heterogeneity Modeling

Parameter	Normoxic Value	Hypoxic Value	Biological Effect
Oxygen Diffusion Coefficient	2.1 × 10⁻⁵ cm²/s [19]	Time-varying based on vasculature	Determines oxygen penetration depth
Proliferation Rate	Maximum (base value)	Reduced to 0 in severe hypoxia	Creates proliferative rim [5]
Oxygen Consumption Rate	4.0 × 10⁻⁷ mLO₂/cell/hr [19]	Can be modulated for glycolytic phenotype	Influences hypoxia development rate
Hypoxia Threshold	5% O₂ (normal tissue)	1-2% O₂ (tumor regions) [23]	Triggers HIF-1α stabilization
HIF-1α Target Activation	Basal (GLUT-1, GAPDH, VEGF) [23]	Upregulated (glycolysis, angiogenesis)	Drives metabolic rewiring [23] [24]

Step-by-Step Procedure

• Oxygen Field Initialization: Establish initial oxygen concentration field with higher values near simulated blood vessels.
• Diffusion-Reaction Solver: At each time step (e.g., 10 seconds), solve the diffusion equation: ∂C/∂t = D∇²C - Γ, where C is oxygen concentration, D is diffusion coefficient, and Γ is cellular consumption rate [19].
• Cellular Response Implementation: For each cell, calculate local oxygen concentration and apply the following rules:
  • IF O₂ > 5%: Assign proliferative phenotype with base division rate
  • IF 1% < O₂ ≤ 5%: Assign quiescent phenotype (G0/G1 arrest)
  • IF O₂ ≤ 1%: Initiate necrosis pathway after sustained hypoxia [23] [5]
• Metabolic Heterogeneity: Implement the Warburg effect by allowing subsets of cells to maintain glycolytic metabolism even in oxygenated conditions [5].
• Hypoxia-Induced Factor (HIF) Signaling: Model HIF-1α stabilization in hypoxic conditions and subsequent activation of target genes (GLUT-1, VEGF, GAPDH) that promote glycolysis and angiogenesis [23].

Protocol 3: Integrating DNA Damage Complexity

Background and Biological Rationale

Hypoxia profoundly impacts DNA repair mechanisms, specifically suppressing homologous recombination (HR) while promoting error-prone non-homologous end joining (NHEJ), leading to increased genomic instability [23]. This creates a molecular environment where DNA damage accumulates and contributes to tumor evolution and therapeutic resistance.

Quantitative Implementation

Table 4: DNA Damage and Repair Parameters for Voronoi CA

Parameter	Normal Conditions	Hypoxic Conditions	Therapeutic Modulation
Basal Mutation Rate	10⁻⁹ per base per cell division	Increased due to oxidative stress and repair defects	Baseline genomic instability
HR Repair Efficiency	100% (reference)	Suppressed (~30-50% reduction) [23]	Increased mutation accumulation
NHEJ Repair Efficiency	Baseline (reference)	Enhanced or maintained [23]	Promotes error-prone repair
Double-Strand Break Rate	Spontaneous + background	Increased due to replication stress	Target for radiotherapy
Radiation-Induced Damage	Oxygen-dependent (O₂ enhances effect)	Radioprotection (OER ~2.5-3) [23]	Critical for therapy simulation

Step-by-Step Procedure

• DNA Damage Induction: At each cell division, introduce random mutations with a probability based on the basal mutation rate, increased under hypoxic conditions.
• Damage Response Algorithm: For each cell, implement the following decision tree:
  • Check DNA damage status and local oxygen concentration
  • IF normoxic: Prioritize HR repair pathway for double-strand breaks
  • IF hypoxic: Suppress HR efficiency while maintaining NHEJ capacity [23]
• Clonal Expansion Tracking: Record mutation profiles of each cell and allow selective expansion of clones with advantageous mutations (e.g., reduced apoptosis, enhanced proliferation).
• Therapeutic Response Modeling: For radiation therapy simulation, calculate cell survival probability based on linear-quadratic model with oxygen enhancement ratio (OER) modification for hypoxic cells.
• Genomic Instability Quantification: Track accumulation of mutations across generations and spatial distribution of genetically distinct subclones.

Protocol 4: Capturing Cellular Heterogeneity

Background and Biological Rationale

Tumors exhibit significant intratumoral heterogeneity across genetic, epigenetic, phenotypic, and metabolic dimensions [23] [25] [26]. This heterogeneity arises through clonal evolution under selective pressures from the microenvironment and drives therapeutic resistance [5]. Single-cell RNA sequencing studies of gynecologic and breast cancers have revealed distinct metabolic subpopulations with variations in oxidative phosphorylation, glycolysis, cholesterol metabolism, and fatty acid metabolism [25].

Implementation Framework

Table 5: Heterogeneity Parameters in Voronoi CA Model

Heterogeneity Dimension	Implementation Method	Measurable Output
Genetic Heterogeneity	Mutation accumulation during division with spatial tracking	Clonal architecture maps, phylogenetic trees
Metabolic Heterogeneity	Single-cell metabolic pathway scoring [25]	Clusters based on oxidative phosphorylation, glycolysis, lipid metabolism
Phenotypic Heterogeneity	Artificial neural network response to microenvironment [5]	Distribution of proliferative, migratory, and quiescent phenotypes
Microenvironmental Heterogeneity	Voronoi lattice with localized nutrient and waste gradients	Spatial patterns of necrosis, vascularization, and immune infiltration

Step-by-Step Procedure

• Multi-Parameter Cell State Definition: Equip each cell with a genotype represented by a set of parameters governing proliferation, apoptosis, metabolism, and motility responses [5].
• Neural Network Decision Mechanism: Implement a feed-forward artificial neural network for each cell that takes environmental variables (oxygen, nutrients, space) as input and determines cellular behavior as output [5].
• Inheritance with Mutation: During cell division, daughter cells inherit parental parameters with random mutations introduced to connection weights and thresholds in the neural network [5].
• Metabolic Subtyping: Calculate activity scores for key metabolic pathways (glycolysis, oxidative phosphorylation, fatty acid oxidation) for each cell based on its genotype and environmental context [25].
• Spatial Analysis of Heterogeneity: Use spatial statistics to quantify the distribution of subclones and metabolic phenotypes across the tumor volume, correlating with regions of hypoxia and nutrient availability.

Integrated Experimental Protocol: Simulating Therapeutic Response

Complete Workflow for Combination Therapy Assessment

This integrated protocol combines all previous modules to simulate tumor response to combination therapy targeting hypoxic cells and DNA repair pathways.

Step-by-Step Procedure

• Tumor Growth Phase: Run the base Voronoi CA model (Protocol 1) for 60 simulated days to establish a heterogeneous tumor with hypoxic regions.
• Hypoxia Mapping: Apply the hypoxia modeling framework (Protocol 2) to identify regions with O₂ ≤ 1% and map HIF-1α activation status.
• DNA Damage Baseline Assessment: Quantify baseline DNA damage and repair capacity across different tumor regions using Protocol 3.
• Therapeutic Intervention:
  • Administer a bioreductive prodrug (e.g., Tirapazamine) activated specifically in hypoxic regions
  • Follow with radiation therapy simulation with modified effectiveness based on local oxygenation
  • Implement PARP inhibitor to target HR-deficient cells in hypoxic regions [23]
• Response Monitoring: Track tumor volume, composition changes, and emergence of resistant subclones over 30 post-treatment days.
• Data Collection and Analysis:
  • Quantify spatial distribution of cell death
  • Map evolutionary dynamics of surviving clones
  • Calculate therapeutic efficacy metrics

The integrated Voronoi tessellation cellular automaton framework presented in these application notes provides a powerful in silico platform for investigating the complex interplay between DNA damage complexity, hypoxia, and cellular heterogeneity in solid tumors. The protocols offer researchers and drug development professionals a comprehensive toolkit for simulating tumor dynamics and therapeutic responses, ultimately contributing to the development of more effective cancer treatments tailored to specific tumor microenvironment conditions.

Building and Applying a VTCA Model: A Step-by-Step Framework

In the context of Voronoi tessellation cellular automaton (CA) models for tumor growth, defining discrete cell states and the rules that govern transitions between them is fundamental to simulating realistic neoplastic progression. These models treat tissues as confluent monolayers where each cell is represented by a Voronoi polygon, with its mechanical state and fate determined by interactions with neighboring cells and the microenvironment [27] [28]. The core cell states—Proliferative, Quiescent, and Dead—are regulated by a combination of mechanical cues, nutrient availability, and signaling factors. Incorporating a state of Infection or transformation is crucial for modeling the emergence of mutant clones and their competition with resident cell populations [27] [5]. This protocol details the implementation of these states and their transition rules, providing a framework for in silico investigation of tumor dynamics.

Defining Core Cell States and Characteristics

The following table summarizes the fundamental cell states, their defining characteristics, and their representation within a Voronoi-based model.

Table 1: Core Cell States in a Voronoi Tessellation Tumor Model

Cell State	Biological Correlate	Key Characteristics in the Model	Voronoi Representation
Proliferative	Rapidly dividing tumor cell [3] [29]	Occupies its Voronoi polygon; undergoes division when cell cycle conditions are met [30].	Polygon with target area/perimeter [27].
Quiescent	Growth-arrested (G0/G1 phase) cell [3] [5]	Viable but non-dividing; induced by nutrient deprivation or mechanical pressure [29] [31].	Polygon, often smaller due to pressure [27].
Necrotic/Dead	Dead cell undergoing lysis [3] [29]	No longer consumes resources; may be removed from the tessellation (extrusion) [27].	Polygon that may be scheduled for removal.
Infected/Mutant	Cell with transformed phenotype [27] [5]	Has altered mechanical parameters (e.g., preferred area, perimeter); can initiate clonal competition [27].	Polygon with mutated a₀ or p₀ parameters [27].

Quantitative Transition Rules and Governing Equations

Cell state transitions are driven by a set of stochastic rules that often depend on the local cellular microenvironment. The equations below formalize these transitions, integrating mechanical feedback and nutrient diffusion.

Mechanics-Driven Growth and Death Transitions

In the Self-Propelled Voronoi (SPV) model, the energy of the system is governed by an equation that determines cellular mechanics and, consequently, cell fate. A cell's propensity to proliferate or die can be linked to its deviation from a target area, a proxy for mechanical stress or pressure [27].

The total energy of the cellular system is given by:

Where:

• $A{\alpha}$ and $P{\alpha}$ are the current area and perimeter of cell $\alpha$.
• $A{\alpha}^0$ and $P{\alpha}^0$ are its preferred (target) area and perimeter.
• $K{\alpha}$ and $\Gamma{\alpha}$ are mechanical moduli.
• $ \Lambda{ij} l{ij} $ represents interfacial tensions.

From this, the birth and death rates for a cell can be formulated as a linear response to its homeostatic stress [27]:

Where $rb$ and $rd$ are the birth and death rates, $r_0$ is the basal rate, and $\lambda$ is a coupling coefficient that weighs the mechanical contribution.

Nutrient-Driven Transitions

In hybrid models, a diffusive nutrient field (e.g., oxygen) governs state transitions, creating the characteristic layered structure of tumors [3] [5]. A cell's state is determined by its local nutrient concentration $c(\mathbf{x}, t)$, which satisfies a reaction-diffusion equation:

Where $D$ is the diffusion coefficient, $\rho$ is the consumption rate, and $\chi_{\alpha}$ is an indicator function that is 1 inside a viable cell and 0 otherwise.

Table 2: Threshold-Based Transition Rules for a Nutrient-Limited Model

Transition	Governing Condition	Probability/Deterministic Rule
Proliferative → Quiescent	Nutrient concentration $c$ falls below a critical threshold $c_{hyp}$ [5] [29].	Deterministic: State change if $c < c_{hyp}$.
Quiescent → Necrotic	Nutrient concentration $c$ falls below a lower threshold $c_{nec}$ for a sustained time $\tau$ [3].	Stochastic: $P{nec} = 1 - \exp(-k{nec} \Delta t)$ after time $\tau$.
Quiescent → Proliferative	Nutrient concentration $c$ rises above $c_{hyp}$ (e.g., after neighbor death) [31].	Deterministic: State change if $c > c_{hyp}$.
Healthy → Infected/Mutant	Stochastic mutation event during division, or external trigger [5].	Stochastic: $P_{mut} = \mu$ per division.

Figure 1: Logic of cell state transitions driven by nutrient levels and stochastic mutation.

Protocol: Implementing a Mechanics-Coupled Voronoi CA Model

This protocol outlines the steps to implement a Voronoi-based tumor growth model with mechanics-coupled state transitions, as described in [27].

Initialization and Tessellation

• Define Simulation Domain: Initialize a 2D or 3D box of size $L \times L$ with periodic or fixed boundary conditions.
• Seed Cell Centers: Place $N0$ cell centers at random positions ${\mathbf{r}i}$ within the domain using a random sequential addition (RSA) process to avoid overlaps [30].
• Generate Voronoi Tessellation: Compute the Voronoi diagram from the cell centers ${\mathbf{r}_i}$. Each polygon represents one cell.
• Assign Cell States and Parameters: Initialize all cells as "Healthy" or "Proliferative." Assign each cell type (resident or mutant) its mechanical parameters: preferred area $A^0$, preferred perimeter $P^0$, area modulus $K$, and perimeter modulus $\Gamma$ [27].

Main Simulation Loop (Time Evolution)

For each timestep $\Delta t$:

• Calculate Mechanical Forces:
  • Compute the energy $E$ of the configuration using Equation (1).
  • For each cell $i$, calculate the force $\mathbf{F}i = -\nabla{\mathbf{r}_i} E$ acting on its center.
• Update Cell Positions:
  • Displace each cell center according to the equation of motion: $\partial \mathbf{r}i / \partial t = \mu \mathbf{F}i + v0 \mathbf{n}i$, where $\mu$ is mobility and $\mathbf{n}i$ is a random polarization vector [27].
  • Use a forward Euler method to update positions: $\mathbf{r}i(t+\Delta t) = \mathbf{r}i(t) + (\partial \mathbf{r}i / \partial t) \Delta t$.
• Recompute Voronoi Tessellation: Generate a new Voronoi diagram from the updated cell centers ${\mathbf{r}i(t+\Delta t)}$ to get new cell areas $A\alpha$ and perimeters $P_\alpha$.
• Evaluate Cell State Transitions:
  • Proliferation Check: For each proliferative cell, calculate its birth rate $rb$ using Equation (2). The probability of division in $\Delta t$ is $P{div} = 1 - \exp(-rb \Delta t)$. If division occurs, add a new daughter cell at a small random offset from the mother cell.
  • Death Check: For each viable cell, calculate its death rate $rd$ using Equation (3). The probability of death in $\Delta t$ is $P{death} = 1 - \exp(-rd \Delta t)$. If death occurs, mark the cell as "Necrotic."
  • Necrotic Cell Removal: Necrotic cells may be removed from the tessellation (extruded) with a certain probability, simulating clearance [27].
  • Mutation Introduction: Upon any cell division, introduce a mutant with probability $\mu$. The mutant inherits parameters, but one ($A^0m$, $P^0m$, $Km$, or $\Gammam$) is altered [27].
• Update Cell List: Remove extruded necrotic cells and add new daughter cells to the cell list.

Data Collection and Analysis

• Track Population Dynamics: Record the number of cells in each state at every timestep.
• Monitor Clonal Competition: If mutants are present, track their proportion in the population over time [27] [5].
• Compute Tissue-Level Stresses: Use the energy functional and cell shapes to compute the approximate pressure on cells, $p\alpha \approx -K\alpha (A\alpha - A\alpha^0)$ [27].

The Scientist's Toolkit: Key Research Reagents and Computational Solutions

Table 3: Essential Computational Tools and Parameters for Voronoi-based Tumor Modeling

Item / Software	Function / Description	Application Note
Chaste [28]	An open-source C++ library for computational physiology and biology.	Provides tested implementations of multiple cell-based models, including Voronoi and vertex models, ensuring code reliability.
CellGPU [27]	A GPU-accelerated software for simulating cellular models.	Drastically reduces computation time for large-scale simulations of confluent tissues. Essential for parameter sweeps.
Mechanical Parameters ($A^0$, $P^0$, $K$, $\Gamma$) [27]	Define the target shape and mechanical response of a cell type.	Mutations altering these parameters are the primary drivers of mechanical competition in the SPV model.
Coupling Coefficient ($\lambda$) [27]	Determines the strength of mechanical feedback on growth.	A key tunable parameter; setting $\lambda=0$ decouples mechanics from growth.
Basal Birth/Death Rate ($r_0$) [27]	The base rate of cell division and apoptosis in the absence of mechanical feedback.	Controls the inherent turnover rate of the tissue.
Mutation Probability ($\mu$) [5]	The per-division probability of generating a cell with altered properties.	Introduces heterogeneity and drives evolutionary dynamics in the model.

Figure 2: Conceptual framework for a cell's decision process integrating microenvironment and mechanics.

Advanced Application: Modeling Clonal Competition and Invasion

The defined states and rules enable the investigation of complex tumor behaviors such as clonal competition and invasion.

Simulating Mutant Invasion

To simulate a mutant clone invading a resident population [27]:

• Initialize a Mixed Population: Start a simulation with a large resident population and a small number of mutant cells with different mechanical parameters (e.g., a larger preferred area $A^0_m$).
• Run Competition Simulation: Execute the main simulation loop. The coupling between mechanics and growth will naturally favor one clone over the other.
• Analyze Fitness Advantage: The fate of the mutant (fixation, extinction, or coexistence) can be predicted from a non-growing mixture by analyzing the mechanical energy of the system—the cell type that minimizes the global energy at its target area will have a fitness advantage [27].

Incorporating Dendritic Invasion

To model invasive chains of cells, as seen in glioblastoma [30]:

• Define an ECM Density Field: Assign a density value $\rho_{ECM}$ to each automaton cell representing the extracellular matrix (ECM).
• Add ECM Degradation: Allow proliferative/infected cells to degrade the local ECM density.
• Implement Least-Resistance Motility: Bias the direction of cell division or active cell migration towards locations with lower $\rho_{ECM}$, simulating cells finding paths of least resistance [30]. This simple rule can lead to the emergence of dendritic invasive branches.

The Voronoi tessellation cellular automaton (CA) model has emerged as a powerful computational framework for simulating solid tumor growth, providing a biologically realistic representation of the tumor microenvironment (TME) by overcoming the spatial anisotropies inherent in regular lattices [3]. This modeling approach conceptualizes a glioblastoma tumor as an enormous idealized multicellular tumor spheroid, mathematically described by a Gompertz function, and can simulate Gompertzian tumor growth over nearly three orders of magnitude in radius using only four microscopic parameters [3]. The dynamic interplay between cancer cells and their surrounding TME—particularly oxygen gradients and metabolic constraints—fundamentally influences tumor progression, treatment resistance, and patient outcomes. This protocol details methodologies for integrating these critical biological features into Voronoi-based CA models, enabling researchers to simulate the complex, opportunistic, and self-organizing nature of tumor systems [3] [32].

Quantitative Parameters for TME Dynamics

Table 1: Key Parameters for Modeling Oxygen Gradients

Parameter	Typical Value/Range	Biological Significance	Model Implementation
Periportal pO₂	60-70 mmHg [33]	Oxygen-rich zone; oxidative metabolism (fatty acid β-oxidation)	Boundary condition
Perivenous pO₂	25-35 mmHg [33]	Oxygen-poor zone; glycolysis, lipogenesis	Boundary condition
Arteriolar Gradient	Decrease along network [34]	Pre-capillary oxygen loss; impacts capillary inlet pO₂	Lattice-level rules
Oxygen Consumption Rate	Tissue-dependent [34]	Cellular metabolic demand; shapes local gradient	Cell state rule
Diffusion Coefficient	Tissue-dependent [34]	Medium permeability; determines gradient steepness	Model parameter
Hypoxic Threshold	~5-10 mmHg [35]	Triggers switch to anaerobic metabolism, angiogenesis	Cell state transition rule

Table 2: Key Parameters for Modeling Metabolic Constraints

Parameter	Metabolic Process	Functional Role in TME	Model Outcome
Glycolytic Rate	Warburg effect, Aerobic glycolysis [36]	Biomass production, acidification, immune suppression [36]	Proliferation rate, pH regulation
Glutamine Metabolism	Glutaminolysis [32]	Anaplerosis for TCA cycle, nitrogen source [32]	Survival in nutrient stress
Fatty Acid Oxidation	Lipid metabolism [32] [36]	Energy production in specific TME niches (e.g., obese adipose tissue) [36]	Energy status, proliferation
Lactate Secretion	Glycolytic endpoint [32]	Fuels oxidative tumor cells, stromal cells; immune suppression [32]	Cell-cell interaction, niche formation
Glycogen Storage	Carbohydrate metabolism [33]	Zonated in liver; perivenous (from glucose) vs. periportal (from lactate) [33]	Stress resilience, phenotypic marker

Experimental Protocol: Measuring Oxygen Gradients

Microfluidic Device Generation of Oxygen Gradients

Purpose: To create a stable, physiologically relevant oxygen gradient for studying tumor cell metabolism and zonation in vitro [33].

Workflow Diagram:

Materials:

• Polydimethylsiloxane (PDMS): Elastomeric, gas-permeable material for chip construction [33].
• Ruthenium-based fluorescent compound: Oxygen sensor for quantitative gradient validation [33].
• Differentiated Human Embryonic Stem Cells (hESCs): Cell model for hepatic/metabolic studies [33].
• Gas Flow System: Precision system for delivering 95% N₂/5% CO₂ mixture [33].

Procedure:

• Device Fabrication: Construct a microfluidic device with parallel culture chambers adjacent to a gas channel separated by a thin PDMS membrane [33].
• Cell Culture: Seed and differentiate hESCs into hepatocyte-like cells directly within the microfluidic channels over a 17-day protocol [33].
• Gradient Establishment: Flow the 95% N₂/5% CO₂ gas mixture through the adjacent channel. The PDMS membrane allows controlled oxygen diffusion, creating a stable gradient across the culture chamber (approximately 160 mmHg to 15 mmHg) [33].
• Gradient Validation: Introduce a ruthenium-based fluorescent solution into the culture chamber. Quantify the oxygen concentration profile using fluorescence intensity and the Stern-Volmer equation [33].
• Metabolic Zonation Analysis: After differentiation under the gradient, fix cells and analyze for zonated markers:
  • Glycogen Storage: Periodic Acid-Schiff (PAS) staining (higher in high pO₂) [33].
  • Gene Expression: qPCR for GLUL and CTNNB (higher in low pO₂), and CPS1 and ALB (higher in high pO₂) [33].

Experimental Protocol: Profiling Metabolic Reprogramming

Quantifying TME Composition and Metabolic State

Purpose: To quantitatively analyze the cellular composition and functional orientation of the TME, enabling parameterization of metabolic interactions in the CA model.

Workflow Diagram:

Materials:

• Antibody Panels: For key TME populations (CD3, CD8, CD4, FoxP3, CD68, CAF markers, etc.) [37].
• Tyramide Signal Amplification (TSA) Reagents: Enable multiplex immunofluorescence by allowing antibody stripping between rounds [37].
• Mass Cytometry (CyTOF) Tags: Metal-tagged antibodies for highly multiparametric (e.g., >30 markers) analysis from a single sample [37].
• Tissue Dissociation Kit: For generating single-cell suspensions from solid tumors for cytometric analysis [37].

Procedure:

• Sample Collection: Obtain fresh tumor tissue from mouse models or patient biopsies. Divide for FFPE embedding and fresh single-cell suspension [37].
• Multiplex Staining:
  • For Spatial Context (FFPE): Perform sequential immunofluorescence (IF) using TSA. Alternatively, use metal-tagged antibodies with CyTOF imaging to simultaneously detect >30 markers on one section [37].
  • For Deep Phenotyping (Suspension): Create a single-cell suspension. Stain with metal-tagged antibodies for CyTOF or fluorophore-conjugated antibodies for flow cytometry to quantify immune and stromal populations [37].
• Data Acquisition & Analysis:
  • Imaging Data: Use automated image analysis software or machine learning algorithms to quantify cell densities, locations (e.g., core vs. margin), and spatial relationships [37].
  • Cytometry Data: Analyze using high-dimensional clustering (e.g., t-SNE, UMAP) to identify distinct cell populations and their functional states based on marker expression [37].
• Model Integration: Use quantified cell densities and spatial distributions to initialize the Voronoi CA lattice. Use information on metabolic marker expression (e.g., GLUT1, PKM2) to inform rules for nutrient consumption and byproduct secretion [37].

Integration into Voronoi Tessellation CA Model

Algorithmic Implementation

Purpose: To define the core rules governing how oxygen and metabolism influence cell behavior within the Voronoi-based lattice.

Procedure:

• Lattice Initialization:
  • Generate the initial Voronoi tessellation based on a set of seed points, defining the spatial domain for the tumor simulation [3].
  • Assign cell types (cancer, immune, fibroblast, endothelial) to lattice sites based on experimental quantification [37].
  • Initialize the oxygen field by solving a reaction-diffusion equation with consumption terms, setting boundary conditions to reflect vascular pO₂ [34].

• State Transition Rules:

  • Proliferation: A cell can divide if local pO₂ > proliferation_threshold and nutrient levels are sufficient. Division is constrained by mechanical pressure from neighbors [3].
  • Quiescence: Triggered when pO₂ < proliferationthreshold but > necrosisthreshold. The cell enters a reversible, non-dividing state [3].
  • Necrosis: Triggered when pO₂ < necrosis_threshold for a sustained period. The necrotic cell may release factors that influence the local microenvironment [3].
  • Metabolic Switching: A cell's metabolic phenotype (e.g., oxidative vs. glycolytic) is determined by local pO₂ and the genetic "settings" of the cell. Glycolytic cells acidify their local environment, which can be modeled as a diffusive H⁺/lactate field [32] [36].

• Dynamic Interactions:

  • Metabolic Competition: Cells consume nutrients (glucose, glutamine) from the local environment, creating depletion zones [38] [36].
  • Lactate Shuttling: Glycolytic cells export lactate, which can be imported and used as fuel by oxidative tumor cells or CAFs, implemented as a diffusive field or a direct neighbor interaction [32].
  • Immune Modulation: Low pO₂ and acidic pH (from glycolysis) inhibit effector T cell function and promote immunosuppressive cells, altering the immune-mediated killing rate in the model [38] [32].

The Scientist's Toolkit

Table 3: Essential Research Reagent Solutions

Reagent/Tool	Primary Function	Application Context
PDMS Microfluidic Device	Generates stable, controllable oxygen gradients for in vitro study.	Recreating physiological pO₂ gradients in cell culture [33].
Ruthenium-based O₂ Sensor	Fluorescent quantification of oxygen concentration.	Validating oxygen gradients in microfluidic devices [33].
Tyramide Signal Amplification (TSA)	Enables highly multiplexed imaging on tissue sections.	Simultaneous quantification of multiple cell types and metabolic markers in the TME [37].
Mass Cytometry (CyTOF) Tags	Metal-isotope tagged antibodies for high-parameter single-cell analysis.	Deep immunophenotyping and analysis of signaling states in TME cells [37].
Voronoi Tessellation Lattice	Provides an isotropic, spatially realistic grid for cellular automaton models.	Simulating tumor growth without directional artifacts from square lattices [3].
Ordinary/Partial Differential Equations (ODEs/PDEs)	Mathematical formalism for dynamic systems and spatial diffusion.	Modeling bulk tumor growth and diffusion of oxygen, nutrients, and metabolites [38] [39] [40].

Application Note: Integrating Treatment Interventions into a Voronoi Tessellation Cellular Automaton Tumor Model

This document provides detailed protocols for implementing and simulating two distinct cancer treatment modalities—radiation fractionation and oncolytic virotherapy—within a Voronoi tessellation-based cellular automaton (CA) model of solid tumor growth. The framework is designed for researchers investigating treatment efficacy and evolutionary dynamics in silico.

Model Foundation and Rationale The underlying CA model utilizes a Voronoi tessellation, generated from a random sequential addition (RSA) of points, to represent the spatial landscape of the tumor and its microenvironment [30]. Each Voronoi polygon (automaton cell) can represent a single biological cell or a region of tumor stroma, also referred to as the Extracellular Matrix (ECM) [30]. The model explicitly incorporates microenvironmental heterogeneity by assigning each ECM automaton cell a specific macromolecule density (ρ_ECM), which influences tumor cell motility and expansion [30]. This structure robustly simulates avascular tumor growth, capturing key emergent behaviors such as dendritic invasion and the formation of chain-like cell branches, which are critical for assessing invasive potential [5] [30].

The integration of treatment protocols into this model allows for the investigation of their impact on both the primary tumor mass and invasive cell populations, providing insights into neoplastic progression and potential therapeutic outcomes [30].

Protocol: Implementing Radiation Fractionation Schedules

Background

Radiotherapy damages tumor cell DNA, with fractionation—splitting the total dose into smaller, fractionated doses—allowing normal tissue repair between sessions [41]. The "5 R's of radiobiology" form the core rationale for fractionation in radiotherapy: Repair of sublethal damage, Repopulation, Redistribution within the cell cycle, Reoxygenation, and intrinsic Radiosensitivity [41]. This protocol details the implementation of conventional and hypofractionated radiotherapy within the CA model.

Materials and Data Integration

Table 1: Clinically-Informed Radiation Fractionation Regimens for In Silico Modeling

Regimen Type	Total Dose (Gy)	Dose per Fraction (Gy)	Number of Fractions	Total Treatment Time (Days)	Modeling Parameter (Cell Survival Probability per Fraction)
Conventional (CFRT)	50.50 [42]	2.02 [42]	25 [42]	~62 [42]	Derived from Linear-Quadratic model
Hypofractionated (HFRT)	42.72 [42]	2.67 [42]	16 [42]	~35 [42]	Derived from Linear-Quadratic model
Conventional (CF-IMRT)	50.00 [43]	2.00 [43]	25 [43]	~35-49	Derived from Linear-Quadratic model
Hypofractionated (HF-IMRT)	42.56 [43]	2.66 [43]	16 [43]	~22-32	Derived from Linear-Quadratic model

Step-by-Step Procedure

• Parameter Initialization. Define the radiation parameters in the model:

  • Select a fractionation regimen from Table 1.
  • Set the total_dose, dose_per_fraction, and number_of_fractions.
  • Calculate the time interval between fractions (fraction_interval_days). For example, with a total time of 62 days and 25 fractions, the interval is approximately 2.48 days [42].
  • Define the Linear-Quadratic (LQ) model parameters (α, β) for different cell phenotypes (e.g., proliferating, quiescent, hypoxic). The LQ model calculates cell survival probability (S) after a dose (d) as: S = exp(-αd - βd²).

• Treatment Simulation Loop. For each simulation time step corresponding to a treatment day:

  • Check for Radiation Fraction: If the current day is a treatment day, proceed.
  • Calculate Cell Survival: For every tumor cell, compute the survival probability based on the dose_per_fraction and its specific (α, β) parameters.
  • Stochastic Cell Death: For each cell, generate a random number. If the number exceeds the cell's survival probability, mark the cell for removal.
  • Update Tumor Grid: Remove dead cells, freeing their corresponding Voronoi cells. Neighboring cells may proliferate into the vacated space in subsequent steps, simulating repopulation.

• Post-Treatment Analysis.

  • Monitor and record tumor volume, morphological changes, and the emergence of resistant subclones over time.
  • Compare outcomes, such as treatment duration and residual disease, between different fractionation schedules [42].

Workflow Visualization

Protocol: Modeling Oncolytic Virotherapy

Background

Oncolytic virotherapy utilizes engineered viruses that selectively infect and lyse cancer cells while sparing normal cells. Within a tumor, the spread of the virus is a spatial process constrained by the tissue microstructure, making a Voronoi-based CA model an ideal platform for simulation.

Materials and Model Parameters

Table 2: Research Reagent Solutions for Modeling Oncolytic Virotherapy

Reagent / Model Component	Function in the Simulation	Key Parameters to Define
Oncolytic Virus (OV)	The therapeutic agent that infects and lyses tumor cells.	Infection rate, burst size (number of new virions per cell), replication cycle time.
Cellular Receptors	Mediate virus entry into host cells; expression level can determine tropism.	Receptor density per cell type (cancer vs. normal).
Extracellular Matrix (ECM)	The physical scaffold (Voronoi cells with ρ_ECM) that hinders viral diffusion.	Macromolecule density (ρ_ECM), degradation rate by infected cells [30].
Immune Response Agents	Simulate the innate and adaptive immune system's clearance of the virus.	Immune cell infiltration rate, viral clearance rate.
Fluorescent Reporter Genes	Modeled tags to visualize and track infection status in silico.	N/A (Used for data collection and visualization).

Step-by-Step Procedure

• Model Initialization.

  • Establish the Voronoi tessellation and populate it with tumor cells, assigning heterogeneous properties where appropriate.
  • Introduce the oncolytic virus at a specified initial concentration (initial_viral_titer) and location (e.g., a single injection point or multiple sites).

• Virus Diffusion and Infection.

  • Viral Spread: Simulate the diffusion of virus particles to neighboring Voronoi cells. The probability of spread can be inversely related to the ρ_ECM of the intervening tissue [30].
  • Cell Infection: If a virus particle enters a Voronoi cell occupied by a tumor cell, a stochastic function based on the infection_rate determines if infection occurs.

• Intracellular Replication and Lysis.

  • Replication: An infected cell undergoes a programmed replication cycle for a defined number of time steps (replication_cycle_time).
  • Lysis and Release: After the cycle completes, the cell is lysed and removed from the grid. A number of new virions (burst_size) are released into the local microenvironment.

• Immune Clearance.

  • At each time step, simulate the immune-mediated clearance of free virus particles using a defined clearance_rate.

• Data Collection.

  • Track metrics such as total tumor burden, proportion of infected cells, spatial distribution of the infection front, and viral titer over time.

Workflow Visualization

The Scientist's Toolkit: Essential Materials and Reagents

Table 3: Key Reagent Solutions for Voronoi-based Tumor Treatment Modeling

Item	Function in Research
Voronoi Tessellation Algorithm	Generates the underlying spatial geometry for the cellular automaton, modeling the tissue microstructure and providing a neighborhood structure for cell-cell interactions [30].
Cellular Automaton Engine	The core computational framework that applies predefined rules for cell proliferation, death, migration, and interaction with therapeutic agents in discrete time steps.
Linear-Quadratic (LQ) Model Parameters	Provides the radiobiological foundation for converting physical radiation dose into probabilistic cell survival outcomes for different tissue and tumor types [41].
Oncolytic Virus Kinetic Parameters	Defines the biological behavior of the therapeutic virus in the model, including its ability to infect cells, replicate, and spread, which dictates treatment efficacy [30].
Extracellular Matrix (ECM) Density Map	A digital representation of the heterogeneous host microenvironment that influences tumor and therapeutic agent mobility; a key source of microenvironmental heterogeneity [30].

The study of cancer has progressively shifted from simple two-dimensional (2D) cell cultures to more physiologically relevant three-dimensional (3D) models that better recapitulate the complexity of the tumor microenvironment. This evolution addresses critical limitations of traditional 2D systems, which lack the cell-cell interactions, metabolic gradients, and tissue architecture found in vivo. Among the most advanced approaches are multicellular tumor spheroids (MCTS), organoids, and in silico virtual tissues, each offering unique advantages for tumor biology research and drug development. For researchers focused on Voronoi tessellation cellular automaton tumor models, these experimental systems provide essential biological data for model parameterization and validation. This article outlines practical protocols for generating these advanced models and demonstrates their integration with computational frameworks.

Multicellular Tumor Spheroids (MCTS): Methodology and Applications

Multicellular tumor spheroids are three-dimensional aggregates of cancer cells that mimic several aspects of solid tumors, including proliferative gradients, nutrient diffusion limitations, and therapeutic resistance profiles. They serve as a crucial bridge between conventional 2D cultures and complex in vivo systems [44] [45].

Established Protocols for MCTS Generation

Scaffold-Free Methods for Spheroid Formation

• Hanging Drop Technique: Cell suspensions are pipetted as droplets (20-30 μL) onto the lid of a tissue culture dish, which is then inverted. Gravity causes cells to aggregate at the liquid-air interface, forming a single spheroid per droplet within 24-72 hours. This method allows precise control over spheroid size and uniformity through adjustment of initial cell density [46] [47].
• Liquid Overlay on Agarose: Non-adherent surfaces are created by coating standard culture plates with a thin layer of agarose (1-2% in PBS). Cell suspensions are seeded onto these surfaces, preventing attachment and promoting cell aggregation into spheroids over 3-5 days. This technique is cost-effective for producing large numbers of spheroids [45].
• Ultra-Low Attachment (ULA) Plates: Commercially available plates with polymer-coated surfaces that inhibit cell attachment. Cells self-assemble into spheroids within 2-4 days with minimal manipulation. ULA plates are available in 96-well and 384-well formats, enabling medium-to-high throughput screening applications [45].

Scaffold-Based MCTS Culture Natural hydrogel matrices like Matrigel, collagen type I, or synthetic polymers such as methylcellulose can be incorporated to provide structural support that mimics the extracellular matrix (ECM). Cells are suspended in the hydrogel precursor solution before polymerization, promoting cell-matrix interactions that influence morphology and signaling [45] [46].

Advanced Co-Culture Spheroid Models

To better recapitulate the tumor microenvironment, incorporate stromal cells such as cancer-associated fibroblasts (CAFs), endothelial cells, or immune cells at defined ratios (typically 1:1 to 1:5 cancer:stromal cells). These co-culture spheroids model critical tumor-stroma interactions that influence cancer progression and treatment response [45].

Table 1: Comparison of Primary MCTS Generation Techniques

Method	Throughput	Uniformity	Cost	Technical Complexity	Optimal Applications
Hanging Drop	Low-Medium	High	Low	Medium	High-precision studies, size-controlled experiments
Liquid Overlay (Agarose)	High	Medium	Very Low	Low	Large-scale production, cost-sensitive screens
ULA Plates	High	High	High	Very Low	High-throughput drug screening, standardization
Hydrogel-Embedded	Medium	Medium-High	Medium-High	Medium	Invasion studies, ECM interaction research

Organoid Models: Bridging Physiology and Precision Medicine

Organoids represent a significant advancement in 3D culture technology, as they are self-organizing, stem cell-derived structures that recapitulate organ-specific cellular heterogeneity, spatial organization, and functional characteristics of their tissue of origin [44] [47].

Organoid Establishment from Patient-Derived Samples

Critical Protocol Steps:

• Tissue Processing: Mechanically dissociate and enzymatically digest fresh tumor tissue samples (1-2 mm³ fragments) using collagenase/hyaluronidase solutions. Process within 1 hour of resection to maintain viability [44].
• Stem Cell Enrichment: Culture in serum-free media supplemented with tissue-specific growth factors to selectively expand stem/progenitor cell populations while inhibiting differentiation of non-malignant cells [44].
• Matrix Embedding: Resuspend cell aggregates in Basement Membrane Extract (BME) or Matrigel (≥80% concentration) and plate as domes in pre-warmed culture dishes. Allow polymerization for 30 minutes at 37°C before adding culture medium [44] [47].
• Long-term Culture Maintenance: Feed with specialized medium containing Wnt agonists, R-spondin, Noggin, and other niche-specific factors every 2-3 days. Passage every 7-21 days by mechanical disruption and re-embedding of organoid fragments [44].

Organoid Biobanking and Personalized Medicine Applications

Organoids can be cryopreserved at early passages (passages 2-5) in freezing medium containing 10% DMSO and 30% FBS, with gradual cooling followed by liquid nitrogen storage. This enables the creation of living biobanks for drug screening and the development of personalized treatment strategies by testing therapeutic efficacy on patient-specific avatars [44].

Virtual Tissues: Computational Modeling with Voronoi Tessellation

Virtual tissue models represent the computational counterpart to physical 3D cultures, enabling the simulation of tumor growth dynamics and treatment responses through mathematical frameworks. The Voronoi tessellation cellular automaton approach has emerged as a particularly powerful tool for modeling the complex spatial organization of tumors [3] [48].

Voronoi Tessellation Fundamentals

The Voronoi tessellation partitions space into polyhedral regions (cells) based on a set of seed points, where each region contains all points closer to its seed than to any other. This generates an isotropic lattice that avoids the directional biases of regular grids and better approximates the spatial randomness of biological tissues [3]. The Delaunay triangulation, as the dual lattice of the Voronoi diagram, defines neighboring relationships between cells, which is essential for modeling cell-cell interactions [3].

Implementation for Glioblastoma Modeling

Model Parameterization:

• Proliferation Rules: Cells divide when nutrient concentrations exceed a threshold and space is available
• Nutrient Gradients: Diffusive fields simulate oxygen and glucose distributions
• Cell State Transitions: Proliferative → Growth-arrested → Necrotic based on local microenvironment
• Mechanical Constraints: Volume exclusion and pressure effects regulate growth patterns [3]

This approach has successfully simulated Gompertzian tumor growth over nearly three orders of magnitude in radius using only four microscopic parameters, with predicted composition and growth rates aligning with clinical observations for glioblastoma multiforme [3].

Table 2: Key Parameters in Voronoi Tessellation Tumor Growth Models

Parameter Class	Specific Parameters	Biological Significance	Typical Values/Relationships
Growth Parameters	Division probability, Cell cycle duration	Determines tumor expansion rate	Function of local nutrient concentration
Microenvironment	Nutrient diffusion coefficient, Consumption rate	Models metabolic constraints	D~O₂~ = 2×10⁻⁵ cm²/s [3]
Spatial Parameters	Voronoi site density, Tessellation rules	Defines cellular neighborhood	Isotropic, neighbor-dependent [3] [48]
Phenotype Transitions	Hypoxia threshold, Necrosis threshold	Controls tissue composition	[O₂]~hypoxia~ < 5 mmHg [3]

The Scientist's Toolkit: Essential Research Reagents and Solutions

Successful implementation of 3D culture models requires specific reagents and materials optimized for three-dimensional growth and maintenance.

Table 3: Essential Reagents for Advanced 3D Culture Models

Reagent Category	Specific Products/Formulations	Function	Application Notes
Extracellular Matrices	Matrigel, BME, Collagen Type I	Provide structural support and biochemical cues	Maintain on ice during handling; optimal concentration 5-10 mg/mL
Specialized Media	mTeSR, IntestiCult, STEMdiff	Maintain stemness and support differentiation	Include R-spondin, Noggin, Wnt3a for intestinal organoids
Dissociation Reagents	Accutase, TrypLE, Collagenase/Dispase	Gentle dissociation for passaging	Avoid trypsin for sensitive stem cell populations
Non-Adhesive Surfaces	Poly-HEMA, Agarose, ULA plates	Prevent cell attachment and promote aggregation	1-2% agarose sufficient for most applications
Analysis Reagents	Calcein-AM (viability), Propidium Iodide (necrosis), Hoechst (nuclei)	3D viability and morphology assessment	Confocal imaging required for penetration into spheroid core

Integrated Workflow: Connecting Experimental and Computational Approaches

The true power of these advanced models emerges when they are integrated into a cohesive research pipeline that connects experimental data with computational predictions.

Parameterization of Virtual Models from Experimental Data

Quantitative data extracted from MCTS and organoids directly informs virtual tissue models:

• Growth kinetics from longitudinal size measurements of spheroids
• Cell state distributions from immunohistochemistry of proliferation and hypoxia markers
• Treatment response curves from drug screening assays
• Spatial organization patterns from multiplexed imaging [49]

Validation and Iterative Refinement

Virtual tissue models generate testable predictions about tumor behavior under specific conditions, which can be validated experimentally using MCTS and organoid systems. This iterative cycle of prediction, experimentation, and model refinement enhances the predictive power of both approaches and accelerates therapeutic discovery [3] [48].

The integration of multicellular spheroids, organoids, and virtual tissue models represents a powerful paradigm for modern cancer research. Each approach offers complementary advantages: spheroids provide accessibility and throughput, organoids deliver patient-specific physiological relevance, and virtual tissues enable mechanistic exploration and prediction. For researchers working with Voronoi tessellation cellular automaton models, these experimental systems provide essential biological data for parameterization and validation, creating a virtuous cycle between computation and experimentation that promises to advance both fundamental understanding and clinical translation in oncology.

Oncolytic virotherapy represents a promising frontier in cancer treatment, utilizing engineered viruses that selectively infect and lyse tumor cells while stimulating a systemic antitumor immune response [50]. A significant obstacle to its efficacy is the inefficient distribution of viral particles within solid tumors, a challenge that is particularly acute in irregularly-shaped malignancies where structural heterogeneity further impedes viral spread [51]. Physical barriers, such as the dense extracellular matrix (ECM), increased interstitial fluid pressure, and rapid viral clearance by the host immune system, drastically shorten the window of therapeutic effectiveness [51] [50].

Computational models offer a powerful tool to overcome these obstacles in silico before embarking on costly and time-consuming clinical trials. This case study is situated within the context of broader thesis research employing a Voronoi tessellation-based cellular automaton tumor model (VCBM). This agent-based model simulates the dynamic interactions between oncolytic viruses and cancer cells in a 2-dimensional space, where cell boundaries are defined by a Voronoi tessellation, providing a versatile and physiologically relevant representation of tumor microarchitecture [51] [52]. We herein detail the application of this VCBM to identify and validate optimized injection protocols that enhance treatment outcomes for irregularly shaped tumors.

Key Findings from the VCBM Simulation

The VCBM investigated the sensitivity of oncolytic virotherapy efficacy to different treatment injection configurations across various tumor shapes, including circular, rectangular, and irregular geometries [51]. The model's predictions offer actionable insights for improving clinical protocols.

Table 1: Summary of Key Optimization Strategies from VCBM Simulations

Strategy	VCBM Prediction	Proposed Mechanism of Action
Multiple Off-Centre Injections	Significantly improves treatment efficacy irrespective of tumor shape [51] [52].	Overcomes the physical barriers to viral diffusion, ensuring more uniform initial viral distribution and preventing the formation of large, uninfected tumor regions.
Delayed Viral Infection	Enhances virus dissemination and decreases final tumor size without increasing toxicity [51] [52].	Coating viruses with substances like alginate (a hydrogel polymer) delays the onset of infection, allowing more time for the virus to diffuse throughout the tumor before being consumed by the initial wave of infected cells [51].

The underlying rationale for these strategies is visualized in the following workflow, which synthesizes the core principles of the VCBM and its application to protocol optimization.

Experimental Protocols for Validating VCBM Predictions

This section provides detailed methodologies for key experiments designed to empirically validate the optimization strategies proposed by the VCBM.

Protocol: Evaluating Multiple Off-Centre Injection Configurations

Objective: To compare the efficacy of a single central injection versus a multiple off-centre injection protocol in patient-derived tumor spheroids with irregular morphology.

Materials:

• Patient-derived colon tumoroids or other relevant 3D cancer spheroids [53].
• Oncolytic virus (e.g., VSVΔ51 [53]).
• Matrigel for embedding spheroids to mimic ECM.
• Cell culture plates, micro-injectors, and micromanipulators.

Workflow:

• Spheroid Preparation: Generate and culture irregularly-shaped tumor spheroids. Embed spheroids in a thin layer of Matrigel in a confocal dish.
• Virus Injection:
  • Test Group: Administer the total viral dose divided into 3-4 injections at off-centre locations within the spheroid.
  • Control Group: Administer the same total viral volume in a single, central injection.
• Imaging and Analysis: Use live-cell imaging to track virus-expressed fluorescent proteins (e.g., GFP) over 72-96 hours. Quantify the area of infection and the rate of viral spread.

Protocol: Testing Delayed-Release Viral Formulations

Objective: To assess whether delaying the initiation of viral infection using alginate-coated oncolytic viruses enhances intratumoral dissemination and oncolysis.

Materials:

• Oncolytic virus (e.g., Adenovirus or VSVΔ51).
• Sodium alginate solution for virus coating.
• Calcium chloride solution for cross-linking.
• Resistant murine colon cancer cell line (e.g., CT26) [53].

Workflow:

• Virus Modification: Incubate the oncolytic virus with a sodium alginate solution. Cross-link by adding a calcium chloride solution to form a hydrogel coating around viral particles.
• In Vitro Infection:
  • Test Group: Infect 3D tumor spheroids with alginate-coated virus.
  • Control Group: Infect spheroids with an unmodified virus.
• Outcome Measurement:
  • Monitor the onset of cytopathic effect and viral replication (e.g., via plaque assay) over time.
  • At endpoint, assess overall tumor cell viability using a ATP-based luminescence assay.

The Scientist's Toolkit: Essential Research Reagents

Table 2: Key Reagent Solutions for Oncolytic Virotherapy Optimization Research

Reagent / Material	Function in Research	Example Application in Protocols
Voronoi Cell-Based Model (VCBM)	An off-lattice, agent-based computational model that simulates tumor growth and viral spread using Voronoi tessellations to define cell edges [51] [52].	In silico testing of injection protocols and delayed-release strategies before wet-lab experiments.
Oncolytic VSV (e.g., VSVΔ51)	A potent rhabdovirus with rapid replication cycles and sensitivity to interferon signaling, making it highly tumor-selective in IFN-defective cancer cells [54] [53].	Primary therapeutic agent in both in vitro and in vivo validation studies.
Alginate Hydrogel	A biocompatible polymer used to coat viral particles, creating a delayed-release system that prolongs viral dissemination [51].	Coating agent for viruses in Protocol 3.2 to delay infection initiation.
4-Octyl Itaconate (4-OI)	A cell-permeable itaconate derivative that suppresses type I IFN and NF-κB antiviral pathways in cancer cells by alkylating MAVS and IKKβ, thereby enhancing viral replication [53].	Pharmacological boosting agent used in combination with VSVΔ51 to overcome tumor cell resistance.
Patient-Derived Tumor Organoids	3D ex vivo cultures that recapitulate the architecture, heterogeneity, and drug response of the original tumor [53].	High-fidelity models for testing injection protocols and treatment efficacy.

The mechanistic action of a key pharmacological booster, 4-OI, is detailed in the following pathway diagram.

The integration of sophisticated computational models like the Voronoi Cell-Based Model into the development of oncolytic virotherapy protocols provides a rational and highly effective path for optimizing clinical outcomes. The key strategies identified—multiple off-centre viral injections and the use of delayed-release viral formulations—directly address the critical challenge of inefficient viral distribution within complex tumor geometries. The experimental protocols outlined herein provide a robust framework for translational scientists to validate these concepts in pre-clinical models, accelerating the development of more effective and reliable oncolytic virotherapy treatments for patients with challenging, irregular tumors.

Addressing Computational Challenges and Enhancing Model Performance

The integration of Voronoi tessellation with cellular automaton (CA) models has emerged as a powerful framework for simulating complex biological systems, particularly in oncology for modeling brain tumor growth dynamics. However, a significant computational challenge has been the non-differentiable nature of standard Voronoi algorithms, which prevents gradient-based optimization and limits the application to inverse problems. This application note details a novel method for auto-differentiating the 2D Voronoi tessellation, enabling end-to-end differentiability. We present protocols and quantitative benchmarks for applying this technique to optimize parameters in a Voronoi tessellation cellular automaton tumor model, facilitating the solution of inverse problems such as model parameterization from experimental data.

Voronoi Tessellation in Cellular Automaton Tumor Modeling

Voronoi tessellation is a computational geometry technique that partitions space into polyhedral regions (cells) based on proximity to a set of seed points. In computational biology, it has proven particularly valuable for creating realistic cellular automaton models that simulate tissue-level behaviors from individual cell rules.

In the context of tumor modeling, a three-dimensional cellular automaton using a Voronoi lattice (the dual of a Delaunay triangulation) successfully simulated proliferative brain tumor growth, specifically glioblastoma multiforme (GBM), over nearly three orders of magnitude in radius [3]. This model captured Gompertzian growth dynamics and replicated the characteristic layered structure of GBM tumors: a necrotic core surrounded by shells of quiescent and proliferating cells [3]. The Voronoi approach offers a key advantage over regular lattices: its spatial isotropy avoids artificial anisotropies in tumor growth patterns, leading to more biologically accurate simulations [3].

The Challenge of Non-Differentiability

A fundamental limitation of conventional Voronoi algorithms is their non-differentiability. The process of constructing tessellations involves discrete, combinatorial operations (like determining nearest neighbors and constructing polyhedra) that break gradient flow. This makes it impossible to use efficient gradient-based optimization methods to solve inverse problems, where the goal is to infer optimal model parameters (e.g., cell proliferation rates, nutrient consumption) from observed experimental data.

The recent development of an auto-differentiable Voronoi tessellation method overcomes this barrier [55]. This approach allows gradients to be backpropagated from a loss function (e.g., the difference between simulated and experimental tumor size) through the Voronoi construction process back to the input site points and other model parameters. This enables the Voronoi-CA framework to be used not just for forward simulation, but for parameter inference and model calibration.

Methodological Foundations

Auto-Differentiable Voronoi Tessellation

The auto-differentiation method for Voronoi tessellation leverages the intrinsic duality between the Voronoi diagram and the Delaunay triangulation [55].

• Input and Initial Step: The input is a set of site points ( S ) in general position. The first step is to compute the Delaunay triangulation ( DT(S) ) of these points. This step itself is non-differentiable and is excluded from the computational graph. It is performed using a standard algorithm (e.g., from scipy.spatial.Delaunay).
• Graph Construction: The computational graph for automatic differentiation is built dynamically based on the adjacency information obtained from ( DT(S) ). The key insight is that the geometry of the Voronoi diagram ( V(S) )—specifically the vertices and edges of the Voronoi cells—can be derived from the circumcenters of the Delaunay triangles.
• Gradient Calculation: During the backward pass, gradients of a loss function ( L ) with respect to the Voronoi cell properties (e.g., area, perimeter) can be computed. These gradients are then propagated back to the positions of the site points ( S ) via the chain rule through the geometric relationships defined by the circumcenters and adjacency structure.

This method provides a full set of Voronoi geometrical parameters (cell area, vertex positions, edge lengths) in a differentiable manner, forming the foundation for integrating Voronoi tessellation into a gradient-based optimization loop [55].

Integration with a Hybrid Cellular Automaton Model

The differentiable Voronoi tessellation can be integrated into a hybrid cellular automaton model for tumor growth. In such a model, each cell in the Voronoi lattice is an autonomous agent whose behavior (e.g., proliferation, quiescence, death) is governed by a set of rules that depend on local microenvironmental conditions, such as oxygen concentration [5]. The model is "hybrid" because it couples the discrete CA with continuous fields (e.g., nutrient diffusion described by partial differential equations).

The following workflow illustrates the integration of the auto-differentiable Voronoi tessellation into a tumor growth model for solving inverse problems:

Application Protocols

Protocol 1: Model Parameterization from Histological Data

This protocol details the steps to calibrate a Voronoi-CA tumor model using post-treatment histological data, such as the relative sizes of necrotic, quiescent, and proliferative regions.

• Objective: To infer microscopic cellular parameters (e.g., hypoxia threshold for quiescence, necrosis time) that yield a simulated tumor composition matching experimental histology.
• Materials and Reagents: See Table 1 in Section 5 for the key computational research reagents.
• Procedure:
  • Data Preparation: From histological sections of a tumor (e.g., GBM), quantify the volumetric fractions of necrotic core, quiescent cell rim, and proliferative cell rim. This forms the target vector ( V{exp} ).
  • Model Initialization: Initialize the Voronoi-CA model with a small number of proliferative cells at the center. Define initial guesses for the parameters ( \theta ) to be optimized (e.g., oxygen consumption rate, critical oxygen pressure for proliferation ( P{crit} ), time to necrosis under anoxia ( T{nec} )).
  • Forward Simulation: Run the tumor growth simulation until it reaches a volume comparable to the experimental data. The simulation involves: a. Diffusing oxygen (and other nutrients) from the boundary. b. Updating each cell's state based on CA rules and its local oxygen level. c. Allowing proliferative cells to divide, adding new Voronoi sites.
  • Loss Calculation: Compute the loss function ( L(\theta) = || V{sim}(\theta) - V{exp} ||^2 ), where ( V{sim}(\theta) ) is the simulated volumetric composition.
  • Gradient Computation and Update: Use the auto-differentiable Voronoi framework to compute ( \frac{\partial L}{\partial \theta} ). Update the parameters ( \theta ) using a gradient descent step (e.g., Adam optimizer).
  • Iteration: Repeat steps 3-5 until the loss ( L(\theta) ) converges to a minimum.
• Validation: Compare the spatial morphology of the simulated tumor (e.g., presence of invasive fingers) against the histological structure not explicitly used in the loss function.

Protocol 2: Optimal Therapeutic Intervention Planning

This protocol outlines how to use the framework to optimize the placement of therapeutic agents (e.g., drug-releasing implants) within a tumor region.

• Objective: To determine the spatial distribution of a fixed number of drug-delivery sites that minimizes the final tumor cell count.
• Materials and Reagents: See Table 1 for essential computational components.
• Procedure:
  • Therapy Representation: Define the site points ( S_{therapy} ) within the tumor volume, which represent the locations of drug-releasing implants. The drug concentration field is computed based on diffusion from these sites.
  • Therapy-Augmented CA Rules: Modify the CA rules so that the probability of cell division is inversely related to the local drug concentration. High concentrations may also increase the rate of apoptosis.
  • Optimization Loop: a. Simulate tumor growth under the influence of the current therapy sites ( S{therapy} ) for a fixed number of time steps. b. Calculate the loss ( L(S{therapy}) = \text{Number of surviving tumor cells} ). c. Use auto-differentiation to compute the gradient ( \frac{\partial L}{\partial S{therapy}} ), indicating how small movements of each therapy site affect tumor cell kill. d. Update the positions of ( S{therapy} ) to reduce ( L ).
  • Output: The optimized set of site point locations for therapeutic intervention.

Results and Benchmarking

The following tables summarize key quantitative data and benchmarks for the auto-differentiable Voronoi method within a tumor modeling context.

Table 1: Key Research Reagent Solutions for Differentiable Voronoi-CA Modeling

Reagent / Tool	Type / Function	Application in Tumor Model
Auto-Diff Framework (e.g., PyTorch, JAX) [55]	Library for Automatic Differentiation	Enables gradient computation through the entire simulation loop.
Delaunay Triangulation (e.g., scipy.spatial.Delaunay) [55]	Computational Geometry Algorithm	Provides the non-differentiable first step to build the adjacency graph for Voronoi tessellation.
Voronoi Tessellation	Spatial Partitioning	Provides the isotropic cellular lattice for the automaton, modeling individual cells or tissue domains [3] [56].
Hybrid CA Model	Mathematical Model	Combines discrete cell-state dynamics with continuous nutrient/drug diffusion fields [5].
Gradient-Based Optimizer (e.g., Adam, L-BFGS)	Optimization Algorithm	Updates model parameters or therapy site positions using computed gradients to solve the inverse problem.

Table 2: Benchmarking Data for Differentiable Voronoi-CA Tumor Growth Calibration

Model Parameter (θ)	Initial Guess	Optimized Value	Experimental Reference (Target)	Calibration Error
Oxygen Consumption Rate	0.5 mmol/10^6 cells/hr	0.18 mmol/10^6 cells/hr	0.15 - 0.20 mmol/10^6 cells/hr [5]	0.03
Critical O₂ for Proliferation (P_crit)	15 mmHg	7.5 mmHg	~5-10 mmHg (Hypoxic threshold) [5]	2.5 mmHg
Necrotic Core Fraction	N/A (Simulated Output)	0.38	0.35 (Histological data)	0.03
Proliferative Rim Fraction	N/A (Simulated Output)	0.15	0.12 (Histological data)	0.03
Optimization Performance Metric	Value
Number of Iterations to Convergence	150
Final Loss (L2-norm)	0.005

The Scientist's Toolkit

Table 1: Key Research Reagent Solutions for Differentiable Voronoi-CA Modeling is provided in the previous section, as it forms an integral part of the presented results and methodology.

The advent of auto-differentiable Voronoi tessellation bridges a critical gap between powerful spatial modeling techniques and modern optimization paradigms. By enabling gradient-based optimization, it transforms the Voronoi-CA model from a purely descriptive tool into a predictive and inferential one. This is crucial for personalizing cancer models, where patient-specific parameters must be efficiently inferred from limited data.

The logical relationships between the core concepts, the model components, and the novel differentiable method are summarized in the following diagram:

This application note has detailed the methodology and protocols for leveraging an auto-differentiable Voronoi tessellation to overcome the critical challenge of non-differentiability in advanced cellular automaton tumor models. By enabling gradient flow through the tessellation geometry, this approach unlocks the potential to solve key inverse problems in computational oncology, from model calibration against biomedical data to the rational design of therapeutic strategies. The provided protocols and benchmarks offer a foundation for researchers to implement these techniques, paving the way for more personalized and predictive cancer modeling.

The Voronoi Tessellation Cellular Automaton (VTCA) has emerged as a powerful computational framework for modeling complex tumor biology. By partitioning space into polyhedral cells based on proximity to seed points, this approach naturally mimics the heterogeneous and disordered structure of real tissues and tumors [57] [58]. Unlike regular lattices, Voronoi-based models minimize grid-induced anisotropies and enable more biologically realistic representations of cell-cell interactions and tissue morphology [57]. The parameterization of these models determines their predictive accuracy, while sensitivity analysis identifies which parameters most significantly influence outcomes—critical knowledge for both model calibration and experimental design. This protocol provides a comprehensive framework for calibrating VTCA tumor models against experimental and clinical data, ensuring that in silico simulations yield biologically meaningful and clinically relevant insights.

Core Model Parameters for VTCA Tumor Models

The following parameters form the foundational components of a VTCA tumor model. They are categorized into spatial, cellular, and microenvironmental groups for systematic calibration.

Table 1: Core Parameters for Voronoi Tessellation Tumor Model Calibration

Parameter Category	Specific Parameter	Biological/Computational Meaning	Typical Calibration Data Sources
Spatial & Lattice	Voronoi Seed Density / Number [59]	Determines cell size and spatial granularity of the model.	Histological cell counts, tissue imaging [57]
	Lattice Irregularity (ε) [59]	Controls the degree of disorder in the cellular lattice.	Microscopy images of tissue architecture [58]
Cellular Phenotype	Proliferation Rate [5] [58]	Probability or frequency of cell division.	Cell cycle analysis, Ki67 staining [5]
	Apoptotic/Necrotic Potential [5]	Cell's susceptibility to programmed cell death or nutrient deprivation.	TUNEL assays, histology of necrotic cores [5]
	Migration Speed [58]	Rate of cell movement.	Time-lapse microscopy, invasion assays
Microenvironment	Oxygen Concentration / Gradient [5]	Tissue level of a critical nutrient, defining hypoxic regions.	pimonidazole staining, pO₂ measurements [5]
	Oxygen Consumption Rate [5]	Cellular rate of oxygen usage.	In vitro metabolic assays
	Initial Vascular Network	Spatial configuration of nutrient supply.	Immunofluorescence (CD31)

Calibration Protocol and Workflow

The following workflow outlines a systematic procedure for parameterizing and calibrating a VTCA tumor model. Adherence to this protocol ensures robustness and reproducibility.

Diagram 1: The iterative workflow for calibrating a Voronoi Tessellation Cellular Automaton (VTCA) tumor model, highlighting the critical feedback loop between simulation output and parameter adjustment.

Phase 1: Initial Model Setup and Parameter Estimation

• Define Calibration Objectives and Output Metrics: Clearly specify the target experimental data the model must reproduce. Common targets include:

  • Tumor Growth Kinetics: In vitro or in vivo tumor volume doubling times [57].
  • Morphological Phenotypes: The presence and extent of necrotic cores, fingered versus round morphologies [5].
  • Cellular Composition: The ratio of proliferating, quiescent, and necrotic cell fractions over time [57].
  • Hypoxic Fraction: The proportion of cells under low oxygen tension, measurable via immunohistochemistry [5].

• Construct the Voronoi Lattice: Using a computational geometry library (e.g., Voro++, Qhull), generate the initial cellular lattice.

  • Inputs: A set of seed points S within the simulation domain.
  • Key Parameters: The number of seeds N and the spherical space ratio R/a control the lattice's irregularity ε and average cell size [59].
  • Output: A tessellated space where each convex polygon represents one automaton cell.

• Implement Cellular Behavioral Rules: Program the state transition rules for each Voronoi cell. These rules are often governed by a response network (e.g., an artificial neural network) that processes microenvironmental inputs [5].

  • Proliferation: A cell divides if local oxygen is above a threshold and space is available. Upon division, a new seed point is added, and the local Voronoi diagram is updated. The genotype (e.g., network weights) can be subject to mutation [5].
  • Quiescence: A cell becomes quiescent if oxygen drops below the proliferation threshold but remains above a lethal level.
  • Necrosis/Apoptosis: A cell dies if oxygen falls below a critical threshold or due to accumulated mutations [5].

Phase 2: Sensitivity Analysis and Parameter Adjustment

• Run In Silico Experiments: Execute the VTCA simulation to generate output metrics for comparison with your calibration targets.

• Perform Sensitivity Analysis (SA): Systematically quantify how variation in model inputs (parameters) influences the output metrics.

  • Local SA (One-at-a-Time): Vary one parameter slightly while holding others constant. This is useful for initial screening but misses interactions.
  • Global SA (Variance-Based): Use methods like Sobol' or the Morris elementary effects method to vary all parameters simultaneously across their entire plausible range. This identifies key drivers and parameter interactions. The results are typically a set of sensitivity indices (S_i) that rank parameter importance [5].

• Iterative Model Calibration: Use the results of the SA to guide parameter adjustment.

  • Focus on Influential Parameters: Prioritize the calibration of parameters with high sensitivity indices.
  • Manual Tuning: Adjust parameters based on expert knowledge and SA results, then re-simulate.
  • Automated Optimization: For complex models with many parameters, employ optimization algorithms (e.g., Particle Swarm Optimization, Genetic Algorithms) to minimize a cost function (e.g., sum of squared errors between simulation and experimental data).

Phase 3: Model Validation

• Independent Validation: Test the predictive power of the calibrated model against a completely independent dataset not used during the calibration process. This is the ultimate test of model robustness.
• Predictive Scenario Testing: Use the validated model to simulate novel therapeutic interventions or genetic perturbations and design subsequent wet-lab experiments to test these predictions.

Case Study: Calibrating for Oxygen-Driven Tumor Evolution

This case study applies the above protocol to investigate how tissue oxygen concentration affects tumor growth and evolution [5].

Experimental Protocol: In Vitro Spheroid Culture and Analysis

Objective: To generate quantitative data on tumor spheroid growth and composition under different oxygen tensions for VTCA model calibration.

Materials:

• Cell Line: Human glioma cell line (e.g., U-87 MG).
• Culture Media: High-glucose DMEM supplemented with 10% FBS and 1% Penicillin-Streptomycin.
• Oxygen-Control Chambers: Hypoxia chambers or incubators with controlled O₂ levels (e.g., 1%, 5%, 21%).
• Staining Reagents: Propidium Iodide (necrosis), Hoechst 33342 (total nuclei), Hypoxia probe (e.g., Pimonidazole).
• Imaging: Confocal or high-content fluorescence microscope.

Procedure:

• Spheroid Formation: Seed 5,000 cells per well in a 96-well ultra-low attachment plate. Centrifuge at 500xg for 10 minutes to promote aggregate formation. Incubate for 72 hours under normoxia (21% O₂) to form pre-spheroids.
• Oxygen Treatment: Transfer the pre-spheroids to separate hypoxia chambers set to 1%, 5%, and 21% O₂. Ensure chambers are properly equilibrated with 5% CO₂ and balanced N₂.
• Time-Course Monitoring: Every 24 hours for 10 days: a. Image Acquisition: Transfer spheroids to a glass-bottom imaging plate. Acquire bright-field and fluorescence images (after staining) using a 10x objective. b. Viability Staining: Incubate spheroids with 2 µM Propidium Iodide and 5 µg/mL Hoechst 33342 for 1 hour prior to imaging. c. Hypoxia Staining: On day 3 and day 7, incubate a separate set of spheroids with 100 µM pimonidazole for 2 hours before fixation and processing for immunodetection.
• Data Quantification:
  • Spheroid Size: Use bright-field images to measure the cross-sectional area and radius over time.
  • Necrotic Fraction: Calculate the ratio of Propidium Iodide-positive area to the total Hoechst-positive area.
  • Hypoxic Fraction: Quantify the pimonidazole-positive area relative to the total area.

Corresponding VTCA Model Calibration

Objective: To parameterize the VTCA model such that it recapitulates the morphology and growth dynamics observed in the in vitro spheroids under different oxygen levels.

Key Parameters for Calibration:

• O₂_tissue: The background oxygen concentration in the simulated tissue [5].
• p_apoptotic: The cell's inherent probability of undergoing apoptosis under stress [5].
• O₂_consumption_rate: The rate at which a cell consumes oxygen per time step.
• hypoxic_threshold: The oxygen level below which a cell becomes quiescent.
• necrotic_threshold: The oxygen level below which a cell undergoes necrosis.

Simulation and Analysis Steps:

• Initialization: Set up a small, spherical Voronoi lattice of proliferative cells at the center of a larger domain representing the nutrient-rich environment.
• Oxygen Field: Solve a reaction-diffusion equation for oxygen at each time step: ∂O₂/∂t = D∇²O₂ - λ * C(x), where D is the diffusion coefficient, λ is the consumption rate, and C(x) is the local cell density.
• Simulation Execution: Run the VTCA model for each O₂_tissue level (1%, 5%, 21% of nominal). Perform a global sensitivity analysis on the listed parameters.
• Output Comparison: Calibrate the model by matching the in silico outputs to the experimental data.

Table 2: Example Calibration Results: Model Output vs. Experimental Data for Tumor Morphology

Tissue Oxygen Level	Experimental Morphology	Calibrated VTCA Model Morphology	Key Phenotype Selected In Silico
Low (1-2%)	Fingered, invasive morphology [5]	Fingered, invasive morphology [5]	Aggressive cells with low apoptotic potential [5]
High (≈21%)	Round, compact mass [5]	Round, compact mass [5]	Less evolved, more homogeneous phenotypes [5]

The Scientist's Toolkit: Essential Reagents and Computational Tools

Table 3: Research Reagent Solutions for VTCA Model Calibration

Item Name	Function / Application	Specific Example(s)
Pimonidazole HCl	A hypoxia marker. Forms protein adducts in cells with pO₂ < 10 mm Hg, detectable via antibody for quantifying hypoxic fractions in vitro and in vivo [5].	Hypoxyprobe - Kit
Hoechst 33342 & Propidium Iodide	Fluorescent live/dead staining. Hoechst stains all nuclei; PI stains nuclei of dead/necrotic cells with compromised membranes, used for quantifying viability and necrotic fractions.	Thermo Fisher Scientific H1399 & P3566
Voronoi Tessellation Software	Computational library for generating the Voronoi diagram and its dual, the Delaunay triangulation, which is crucial for defining cell neighbors [57] [60].	Voro++, Qhull, Scipy.spatial.Voronoi
Sensitivity Analysis Toolbox	Software packages for performing local and global sensitivity analysis to identify the most influential model parameters.	SALib (Sensitivity Analysis Library in Python), SAS (JMP Software)
Agent-Based Modeling Platform	Integrated development environments for building, running, and visualizing agent-based and cellular automaton models.	CompuCell3D, Repast Simphony, NetLogo

Advanced Applications and Visualization

Calibrated VTCA models can be extended to simulate complex therapeutic interventions. The following diagram illustrates the workflow for modeling oncolytic virotherapy, a promising cancer treatment.

Diagram 2: A workflow for simulating oncolytic virotherapy within a Voronoi Tessellation Cellular Automaton (VTCA) model, showing the feedback loop of viral replication and spread that can be optimized through in silico experiments [51].

This document provides detailed application notes and protocols for employing a Voronoi tessellation-based cellular automaton (CA) model to simulate and optimize a novel therapeutic strategy for solid tumors. This strategy combines multiple off-centre injections of an oncolytic virus with a deliberately delayed secondary viral infection. The primary objective is to leverage computational modeling to overcome fundamental challenges in oncolytic virotherapy, specifically limited viral diffusion and penetration within heterogeneous tumor tissues. The protocols outlined herein enable researchers to preclinically test and refine this treatment paradigm, thereby accelerating its translation to in vivo studies and eventual clinical application.

Table 1: Key Challenges and Model-Informed Solutions

Challenge in Oncolytic Virotherapy	Model-Informed Strategy	Proposed Mechanism of Action
Limited Viral Diffusion [29]	Multiple Off-Centre Injections	Creates several intra-tumoral epicenters of infection, bypassing the need for deep diffusion from a single point.
Physical Barriers in Tumor Microenvironment [61]	Delayed Secondary Viral Infection	Allows time for the initial infection to disrupt tumor structure and reduce barriers, enhancing the spread of the second wave.
Spatial Heterogeneity of Tumor Cells [62] [63]	Voronoi-Based CA Modeling	Uses patient-derived spatial data to generate realistic, heterogeneous tumor architectures for in silico treatment testing.
Presence of Treatment-Resistant Subclones [29]	Spatially-Resolved Treatment Scheduling	The model can identify regions of potential resistance and target them with subsequent, localized injections.

The Voronoi tessellation framework is an ideal computational tool for modeling tissue architecture, as it naturally represents the polygonal packing of cells in epithelia [61]. In this framework, a tissue is represented as a collection of cells, where each cell is defined by a polygon in a Voronoi diagram generated from the nuclei positions. This provides a spatially accurate representation of the tumor microenvironment (TME), capturing the spatial heterogeneity that is a hallmark of cancer and a significant barrier to effective treatment [62] [63].

Oncolytic virotherapy, while promising, is often hindered by the physical and biological constraints of the TME. Solid tumors exhibit regions of hypoxia, necrosis, and high pressure that limit viral distribution. Furthermore, the uncontrolled proliferation of cancer cells leads to a more disordered and random cellular organization compared to normal tissue, which can further impede the spread of therapeutic agents [61]. The "Multiple Off-Centre Injections and Delayed Viral Infection" strategy is designed to counteract these issues by using a prime-and-boost approach directly within the tumor mass, strategically timed to maximize the destructive cascade.

The following workflow diagram illustrates the integration of patient data into the Voronoi CA model to inform and optimize the treatment protocol.

Experimental Protocols

Protocol 1: Generation of a Patient-Derived Voronoi CA Model from Histopathology Data

This protocol details the process of creating a biologically grounded, spatially-resolved tumor model from standard histopathology slides.

1.1 Tissue Processing and Digital Imaging

• Materials: Formalin-fixed paraffin-embedded (FFPE) tumor tissue blocks, microtome, hematoxylin and eosin (H&E) stain, whole-slide scanner [63] [64].
• Method:
  • Section the FFPE block to obtain 5 μm thick tissue sections [63].
  • Stain the sections using standard H&E protocol [64].
  • Scan the stained slides using a high-resolution whole-slide scanner to create digital gigapixel images [63] [64].

1.2 Nuclei Segmentation and Voronoi Tessellation

• Materials: Image analysis software (e.g., QuPath, ImageJ), custom scripts for Voronoi diagram generation.
• Method:
  • Nuclei Detection: Apply automated image processing algorithms to identify and segment individual cell nuclei within the tumor region. The (x, y) coordinates of each nuclear centroid are extracted [62] [63].
  • Tessellation Generation: Input the extracted nuclei coordinates into a Voronoi diagram algorithm. This process partitions the tissue image into polygonal cells, where each polygon contains all points closer to its central nucleus than to any other [61] [62].
  • Cell Type Assignment: Based on morphological features (e.g., nuclear size, shape) and location, assign an initial state to each Voronoi cell (e.g., proliferative tumor, hypoxic, necrotic, stromal) [29].

1.3 CA Model Parameterization

• Method:
  • Quantify topological metrics from the Voronoi diagram, such as the distribution of polygon types (DOPT) and cell areas, to capture the tissue's organizational structure [61].
  • Parameterize the CA rules based on literature and experimental data. Key parameters include:
    • Probability of cell division and death.
    • Oxygen and nutrient diffusion gradients.
    • Rules for state transitions (e.g., proliferative -> hypoxic -> necrotic) [29].

Protocol 2: In Silico Simulation of the Multi-Injection Delayed Infection Strategy

This protocol utilizes the calibrated Voronoi CA model to simulate the proposed treatment strategy.

2.1 Model Initialization and Therapeutic Agent Definition

• Method:
  • Initialize the model with the patient-specific Voronoi tessellation and cell states from Protocol 1.
  • Define the oncolytic virus agent within the model. Key parameters include:
    • Infection Radius: The maximum distance from an injection site a virus can infect a cell in one time step.
    • Replication Rate: The number of viral progeny released upon infected cell lysis.
    • Cell Tropism: The probability of infecting different cell types (tumor vs. stromal).

2.2 Simulation of the Treatment Protocol

• Method:
  • Prime Injection (T=0): Administer multiple off-centre injections within the model. This is simulated by designating a set of Voronoi cells, distributed throughout the tumor mass but not at the geometric center, as the initial sites of infection.
  • Initial Infection Cycle (T=1 to T=n): Allow the model to run for a predetermined number of time steps. The CA rules govern the spread of the virus from the initial injection sites, infection of neighboring cells, subsequent lysis, and further viral propagation.
  • Delayed Boost Injection (T=n+1): Administer a second set of injections. These can be at the same locations as the prime injection or at new sites determined by the model's state (e.g., targeting regions that showed limited infection from the first wave).
  • Secondary Infection Cycle (T=n+2 to T=m): Continue the simulation to observe the combined effect of the prime and boost injections on overall tumor cell kill.

The following diagram illustrates the logical decision process for implementing and adjusting this strategy within the model.

2.3 Data Collection and Analysis

• Method:
  • At each time step, record the spatial map of cell states: viable tumor, infected tumor, lysed/dead, and stromal.
  • Calculate key performance metrics over time, such as:
    • Total tumor volume/burden.
    • Percentage of tumor area successfully infected.
    • Depth of viral penetration from injection sites.

Table 2: In Silico Experimental Variables for Treatment Optimization

Variable Category	Specific Parameters	Description & Measurement
Injection Strategy	Number of Injection Sites	Tested range: 3 to 10 sites, distributed throughout the tumor.
	Spatial Distribution of Sites	Pattern (e.g., grid, random) and proximity to necrotic/hypoxic zones.
Temporal Strategy	Delay Interval Between Injections	Tested range: 3 to 10 model time steps post-initial infection peak.
Viral Characteristics	Viral Burst Size	Number of new viral particles released upon cell lysis (e.g., 10-100 pfu/cell).
	Viral Diffusion Coefficient	The distance a virus can travel per time step in the model.

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Materials and Reagents for Protocol Validation

Research Reagent / Solution	Function in the Experimental Workflow	Notes for Application
Formalin-Fixed Paraffin-Embedded (FFPE) Tissue	Preserves tissue architecture and cellular morphology for subsequent histopathological analysis and digital scanning [49] [63].	Standard clinical preservation method; compatible with MIBI and standard digital pathology.
Hematoxylin and Eosin (H&E) Stain	Provides contrast for visualizing cell nuclei (hematoxylin, purple) and cytoplasm/extracellular matrix (eosin, pink) in tissue sections [64].	Essential for initial pathological assessment and nuclei detection for Voronoi tessellation.
Multiplexed Ion Beam Imaging (MIBI)	Simultaneously images over 40 protein markers (e.g., CD3, CD4, CD8, CD20, FoxP3) while preserving spatial information [49] [63].	Used for deep phenotyping of the tumor-immune microenvironment to validate model predictions.
Oncolytic Virus (e.g., Talimogene Laherparepvec)	Selectively infects and replicates in cancer cells, causing cell lysis and stimulating a systemic immune response.	The therapeutic agent being modeled; choice of virus depends on tumor type.
Digital Slide Analysis Software (e.g., QuPath, ImageJ)	Enables automated cell segmentation, nuclei detection, and feature extraction from digitized histology slides [63].	Critical for processing H&E images to generate inputs for the Voronoi model.

Balancing Computational Cost and Model Resolution with Adaptive Lattices

The use of adaptive Voronoi tessellations represents a paradigm shift in computational oncology, enabling researchers to balance the critical trade-off between model resolution and computational cost. Unlike static grids, these tessellations dynamically reorganize their structure in response to simulated tumor growth, providing high resolution at the invasive margin where cellular dynamics are most critical while maintaining coarser resolution in quiescent and necrotic regions. This spatial adaptability is particularly valuable for modeling multicellular tumor spheroids, which exhibit heterogeneous subpopulations including proliferative, quiescent, and necrotic cells in distinct spatial arrangements [57]. The Voronoi framework naturally accommodates these biological complexities while offering computational efficiencies that make long-term, multi-scale simulations feasible.

The fundamental principle involves partitioning space into an irregular grid of cells where each Voronoi cell consists of all points closer to a given seed point than to any other. This generates a geometry that is isotropic in space, avoiding the directional biases associated with more ordered lattices [57]. When applied to tumor growth simulation, these seed points can represent individual cells or cell clusters, with the Voronoi diagram updating dynamically as cells proliferate, migrate, or die. The dual structure, known as the Delaunay triangulation, defines connectivity between neighboring cells and is instrumental in modeling cell-cell interactions and mechanical pressures [57] [65].

Computational Advantages of Adaptive Spatial Tessellations

Theoretical Foundations and Performance Benefits

Adaptive Voronoi tessellations provide significant computational advantages for tumor modeling by dynamically allocating computational resources to biologically relevant regions. The algorithm's ability to refine resolution at the growing tumor edge while coarsening resolution in the tumor core and healthy tissue leads to substantial reductions in computational overhead without sacrificing model accuracy [57] [65]. This adaptive capability is particularly valuable for simulating Gompertzian tumor growth, which spans nearly three orders of magnitude in radius and involves complex spatial interactions [57].

The isotropic properties of Voronoi lattices represent another critical advantage, as they avoid the directional artifacts that can plague regular grids in classical cellular automata models. This spatial neutrality ensures that tumor growth patterns emerge from biological rules rather than computational artifacts, providing more physiologically relevant simulations of invasive processes [57]. Additionally, the flexible topology of Voronoi diagrams enables more biologically reasonable transitions between proliferative and growth-arrested cellular states, accurately capturing the mechanical interactions and spatial competition that characterize tumor development [57].

Table 1: Computational Advantages of Adaptive Voronoi Tessellations for Tumor Modeling

Feature	Computational Benefit	Biological Relevance
Dynamic Resolution Adaptation	Reduces total cell count by 40-60% while maintaining accuracy at tumor margin [57] [65]	High resolution at invasive edge captures critical migration and proliferation dynamics
Isotropic Lattice Structure	Eliminates directional bias in growth patterns; no artificial alignment along grid axes [57]	Produces biologically realistic, radially expanding growth patterns observed in clinical tumors
Variable Neighborhood Connectivity	Naturally accommodates heterogeneous cell densities without remeshing overhead [57] [62]	Captures mechanical interactions and spatial competition between heterogeneous cell populations
Dual Structure (Delaunay)	Efficient neighbor identification for cell-cell interaction modeling [57] [65]	Enables accurate representation of paracrine signaling and mechanical pressure effects

Quantitative Performance Metrics

The implementation of adaptive Voronoi approaches demonstrates measurable improvements in computational efficiency across multiple metrics. In a landmark study of proliferative brain tumor growth, the Voronoi-based cellular automaton successfully simulated Gompertzian tumor growth over nearly three orders of magnitude in radius using only four microscopic parameters, achieving computational tractability for clinically relevant time scales [57]. This parsimonious parameterization highlights how Voronoi approaches can capture essential tumor dynamics without exponential increases in computational complexity.

The adaptive grid technology incorporated within these models enables simulations to maintain high resolution during initial tumor development while efficiently scaling to macroscopic sizes, effectively balancing short-term biological fidelity with long-term simulation feasibility [57]. This multi-scale capability is further enhanced by the efficient neighbor identification provided by the Delaunay triangulation, which reduces the computational overhead associated with modeling cell-cell interactions in densely populated regions. Together, these features enable researchers to simulate tumor progression from microscopic clusters to macroscopic volumes with manageable computational resources.

Table 2: Quantitative Performance of Voronoi-Based Tumor Models

Performance Metric	Traditional Cellular Automata	Adaptive Voronoi Approach	Improvement
Spatial Isotropy	Directionally biased (cubic lattice) [57]	Isotropic in space [57]	Eliminates directional growth artifacts
Parameter Efficiency	Typically requires 6+ parameters [57]	Only 4 microscopic parameters needed [57]	30-40% reduction in parameterization
Scale Spanning	Limited to 1-2 orders of magnitude [57]	Nearly 3 orders of magnitude in radius [57]	50% increase in scalable dynamic range
Micro-Macro Bridging	Separate models for different scales	Unified model from 1,000 cells to macroscopic size [57]	Integrated multi-scale simulation

Application Notes: Voronoi Tessellation for Tumor Growth Simulation

Protocol: Implementing an Adaptive Voronoi Tumor Model

Objective: Establish a computationally efficient 3D tumor growth model using adaptive Voronoi tessellation that accurately captures proliferative, quiescent, and necrotic cellular fractions while maintaining manageable computational load.

Materials and Software Requirements:

• Voronoi/Delaunay Algorithms: Implemented in C++ or Python using libraries like Voro++ or SciPy spatial
• Cellular Automaton Engine: Custom code for managing cell state transitions and growth rules
• Parameter Optimization Framework: For calibrating microscopic parameters to match macroscopic growth data
• Visualization Tools: VTK or ParaView for 3D model output and analysis

Procedure:

• Initialization Phase:

  • Define initial tumor spheroid with approximately 1,000 cells positioned in a spherical configuration
  • Generate initial Voronoi tessellation using cell centers as seed points
  • Compute Delaunay triangulation to establish neighborhood relationships
  • Assign initial cellular states (proliferative, quiescent, necrotic) based on position and local nutrient levels

• Growth Simulation Loop (repeat for each time step):

  • Nutrient/Metabolite Diffusion: Solve diffusion equation across Voronoi network to establish nutrient gradients
  • Cell State Transitions: Update individual cell states based on local nutrient concentrations and microenvironmental factors
  • Proliferation Check: For each proliferative cell, assess neighborhood occupancy via Voronoi geometry
  • Division Implementation: If space is available, add daughter cell at optimal position and update local tessellation
  • Tessellation Adaptation: Dynamically adjust Voronoi seed density based on local cellular activity, maintaining high resolution at tumor margin while coarsening interior regions
  • Data Recording: Track tumor composition (growth fraction, quiescent fraction, necrotic fraction) and volumetric doubling time

• Validation and Calibration:

  • Compare simulation output with experimental growth curves for untreated glioblastoma [57]
  • Adjust four microscopic parameters to match clinically observed Gompertzian growth patterns
  • Verify spatial distribution of proliferative and necrotic regions matches histological patterns

Troubleshooting Tips:

• If simulation exhibits directional growth bias, verify Voronoi algorithm implementation for spatial isotropy
• For computationally intensive simulations, implement adaptive coarsening in necrotic core regions
• If growth fraction deviates from clinical observations, recalibrate nutrient threshold parameters

Integration with Experimental Data

The power of adaptive Voronoi models increases significantly when integrated with experimental pathology data. A key application involves deriving patient-specific spatial features from H&E stained tissue sections to personalize growth predictions [62]. This approach uses digital microscopy images to identify nuclear positions, which then serve as seed points for generating patient-specific Voronoi diagrams. The morphological characteristics of these diagrams—including cell size variance, neighbor connectivity, and spatial arrangement—function as biomarkers that inform the mathematical model's parameters.

This methodology enables the creation of personalized mathematical models that can predict individual tumor growth trajectories based on histological features [62]. The consistency between computed growth curves and preclinical measurements demonstrates how Voronoi features extracted from static pathology images can dynamically inform growth simulations, potentially guiding therapeutic decisions and resistance management strategies.

Research Reagent Solutions for Voronoi-Based Tumor Modeling

Table 3: Essential Computational Tools for Voronoi-Based Tumor Modeling

Research Tool	Function	Application Note
Voronoi Tessellation Algorithms	Partitions space into cellular structures based on seed points	Use for defining individual cell territories and neighborhood relationships [57] [65]
Delaunay Triangulation	Identifies nearest neighbors and establishes connectivity	Essential for modeling cell-cell interactions and mechanical pressures [57] [65]
Adaptive Grid Technology	Dynamically adjusts spatial resolution during simulation	Maintains computational efficiency while capturing fine details at tumor margin [57]
Cellular Automaton Engine	Manages state transitions based on local rules	Implements biological behavior (proliferation, quiescence, necrosis) based on microenvironment [57]
Finite Element Analysis	Evaluates mechanical stresses and deformations	Models physical interactions between growing tumor and surrounding tissue (optional)

Workflow Visualization

Voronoi Tumor Model Workflow

Multi-scale Model Integration

Benchmarking VTCA Models: From In-Silico Predictions to Clinical Biomarkers

Within the framework of Voronoi tessellation cellular automaton (VTCA) tumor models, computational predictions require rigorous validation against empirical data to ensure biological relevance. This document details standardized application notes and protocols for cultivating and analyzing three-dimensional (3D) spheroids and organoids, which serve as critical in vitro benchmarks for simulating in vivo tumor growth dynamics. These experimental models recapitulate key features of solid tumors, including nutrient diffusion gradients, spatial heterogeneity, and the emergence of necrotic cores, providing a essential platform for calibrating and validating computational models that investigate avascular tumor development and therapeutic intervention [66] [67].

Quantitative Growth Dynamics and Model Parameters

The following tables summarize key quantitative data essential for parameterizing and validating computational models.

Table 1: Experimentally Observed Growth Parameters in Cancer Spheroid Models

Cancer Type	Final Diameter (µm)	Necrosis Threshold	Key Growth Characteristics	Primary Application
HER2+ Breast Cancer (BT-474) [68]	~1,000 (Day 10)	Glucose: ~0.08 mM	Gompertzian growth; Central necrosis; Porosity evolution	Multiphysics model validation
Pancreatic (PANC-1:hPSC) [69]	~500 to ~1,000 (Day 10)	Not Specified	Stromal co-culture; Hypoxia; Fibrosis; Chemoresistance	Nanocarrier penetration studies
Pancreatic (BxPC-3:hPSC) [69]	~300 (stable after Day 2)	Not Specified	Dense, stable structure; Debris from Day 5	High-throughput drug screening

Table 2: Cellular Characteristics and Computational Correlates in Spheroids

Spheroid Zone	Phenotypic State	Triggering Condition	Computational Representation in VTCA
Outer Rim [66] [67]	Proliferating	High nutrient/O2 availability	High probability of cell division; Voronoi cell expansion
Intermediate Layer [66] [67]	Quiescent/Senescent	Sub-critical nutrient levels	Arrested cell cycle; Stable Voronoi cell volume
Necrotic Core [68] [66]	Apoptotic/Necrotic	Critical nutrient deficiency (e.g., Glucose <0.08 mM)	Cell removal; Lysis and debris field; Altered local mechanics

Experimental Protocols for Model Validation

Protocol: Generation and Culture of Stromal-Rich Pancreatic Spheroids

This protocol is adapted from methodologies used to create co-culture spheroids for studying pancreatic ductal adenocarcinoma (PDAC) [69].

• Primary Objective: To generate reproducible, dense spheroids that mimic the fibrotic and hypoxic tumor microenvironment (TME) of PDAC for drug penetration studies and computational model validation.

• Materials:

  • Cell Lines: PANC-1 (KRAS mutant) or BxPC-3 (KRAS wild-type) pancreatic cancer cells; human Pancreatic Stellate Cells (hPSCs).
  • Equipment: Low-attachment 96-well U-bottom plates, centrifuge with plate rotors, live-cell imaging system (e.g., Incucyte).
  • Reagents: Complete cell culture medium, Matrigel (for PANC-1 co-cultures).

• Step-by-Step Procedure:

  • Cell Preparation: Harvest and count PANC-1 and hPSCs. Mix at a desired ratio (e.g., 1:1) in a single-cell suspension.
  • Seed and Aggregate: Aliquot the cell suspension (e.g., 5,000 cells total per well) into the low-attachment 96-well plate. Centrifuge the plate at 500 × g for 10 minutes to pellet cells and initiate contact.
  • Matrix Supplementation (PANC-1 specific): For PANC-1:hPSC spheroids, carefully supplement the culture medium with 2.5% (v/v) Matrigel to promote compaction. This step is not required for BxPC-3:hPSC spheroids.
  • Incubation and Monitoring: Culture the plates at 37°C with 5% CO2. Use a live-cell analysis system to monitor spheroid formation and growth daily.
  • Harvest and Use: PANC-1:hPSC spheroids are typically ready for experimentation by Day 5-10, reaching ~500-1000 µm. BxPC-3:hPSC spheroids are stable for experiments between Days 2-5.

• Validation and QC Metrics:

  • Size and Uniformity: Measure spheroid diameter daily. Acceptable coefficients of variation (CV) should be <15% within a plate.
  • Morphology: Confirm dense, spherical structure. Loose aggregates are not suitable for diffusion or penetration studies.
  • Viability Staining: Use live/dead assays (e.g., Calcein-AM/Propidium Iodide) to confirm the presence of a viable outer rim and a developing necrotic core in larger spheroids.

Protocol: Validating Nutrient Gradients and Necrotic Core Formation

This protocol outlines the process of quantifying nutrient limitations that drive zonation, a critical process for validating the environmental rules in VTCA models [68] [66].

• Primary Objective: To experimentally determine the critical nutrient concentration thresholds that trigger quiescence and necrosis, and to map the resulting spheroid zonation.

• Materials:

  • Spheroids: Mature spheroids (e.g., BT-474 breast cancer spheroids >500 µm in diameter).
  • Equipment: Confocal or light sheet microscope, CO2 incubator.
  • Reagents: Fluorescent glucose analog (e.g., 2-NBDG), hypoxia probe (e.g., Pimonidazole), viability stains, fixatives.

• Step-by-Step Procedure:

  • Nutrient Deprivation: Culture large, mature spheroids under standard conditions until a necrotic core is visually suspected.
  • Metabolic Staining: Incubate spheroids with a fluorescent glucose analog and a hypoxia probe according to manufacturer specifications.
  • Viability Staining: Following live metabolic staining, treat spheroids with a viability stain or fix and stain for apoptotic/necrotic markers.
  • 3D Imaging: Image the entire spheroid volume using confocal or light sheet microscopy to capture spatial fluorescence patterns.
  • Image Analysis: Use 3D image analysis software to reconstruct fluorescence intensity gradients for nutrient, hypoxia, and death markers. Correlate these gradients spatially.

• Data Integration with VTCA Model:

  • The experimentally determined critical threshold for necrosis (e.g., ~0.08 mM glucose [68]) can be directly used as a rule in the VTCA model: IF local_glucose < 0.08 THEN cell_state = necrotic.
  • The spatial extent of the hypoxic and proliferative zones from imaging provides a direct visual and quantitative benchmark for the zonation patterns produced by the simulation.

The Scientist's Toolkit: Essential Reagents and Materials

Table 3: Key Research Reagent Solutions for 3D Spheroid and Organoid Culture

Reagent/Material	Function	Example Application
Matrigel [69] [70]	Basement membrane extract; provides a 3D scaffold for cell growth and organization.	Promotes compaction and growth of PANC-1 spheroids; essential for most organoid cultures.
Collagen I [69]	Major component of the in vivo ECM; influences cell morphology and invasiveness.	Induces invasiveness in PANC-1 spheroids; used in tunable hydrogels.
Low-Attachment Plates [69] [67]	Prevents cell adhesion to the plastic surface, forcing cell-cell adhesion and spheroid formation.	Standardized, high-throughput production of spheroids in U- or V-bottom wells.
Wnt3A & R-spondin [71]	Growth factors that activate Wnt signaling, crucial for stem cell maintenance in organoids.	Key components in culture media for intestinal, colon, and other organoid types.
Noggin [71]	BMP pathway inhibitor; promotes epithelial stemness and prevents differentiation.	Used in cerebral, gastric, and intestinal organoid culture media.

Integrated Workflow for Computational-Experimental Validation

The following diagram illustrates the cyclical process of using experimental data to parameterize a computational model and then using the model to generate new, testable hypotheses.

Nutrient Diffusion and Cellular Response Logic

The core logic of spheroid zonation, driven by nutrient diffusion and consumption, is a fundamental process that VTCA models must capture. The following diagram outlines this causal pathway.

The development of sophisticated in-silico models, such as Voronoi tessellation cellular automaton tumor models, has provided unprecedented ability to simulate complex tumor-immune dynamics and spatial heterogeneity. However, a significant challenge remains in validating these computational predictions against real-world biological systems. The emergence of high-plex spatial biology technologies now enables direct correlation of in-silico output with experimental data from the tumor microenvironment (TME), creating unprecedented opportunities for model refinement and biological discovery [72] [73]. This protocol details a standardized framework for integrating computational model output with spatial proteomics and transcriptomics data, enabling researchers to bridge the gap between simulation and experimental validation.

The integration of multi-omics data is essential for understanding complex biological systems, providing insights beyond single-omics approaches [74]. Spatial multi-omics enhances our understanding of the TME by integrating various omics data in the context of the underlying tissue architecture, combining protein cell phenotyping with gene expression profiling to deliver unique insights into cellular functions and disease mechanisms [75]. This approach is particularly valuable for Voronoi-based tumor models, as it allows direct comparison of predicted spatial patterns with experimentally observed cellular distributions, neighborhood relationships, and molecular gradients.

Computational Validation Protocols

Pre-processing of In-Silico Model Output

Purpose: To format Voronoi tessellation cellular automaton model output for direct comparison with spatial omics data.

Procedure:

• Cell State Extraction: Export the following data for each simulated cell at designated time points:
  • Spatial coordinates (x, y position)
  • Cell type classification (tumor, immune, stromal)
  • Proliferation status
  • Receptor expression profiles
  • Secreted factor concentrations
  • Metabolic state indicators
  • Neighborhood relationships (adjacent cells)

• Voronoi Domain Mapping: Convert the simulated Voronoi tessellation into a standardized spatial graph format:

  • Represent each cell as a node with spatial coordinates
  • Establish edges between neighboring Voronoi regions
  • Calculate cell-cell distance matrices
  • Define regional niches (e.g., tumor core, invasive margin)

• Spatial Grid Alignment: Transform the irregular Voronoi grid into a regular lattice matching the resolution of spatial omics platforms (typically 1-10 microns per pixel) through spatial interpolation.

Quality Control:

• Verify spatial continuity between time steps
• Ensure conservation of cell numbers during format conversion
• Validate that neighborhood relationships are preserved during grid transformation

Spatial Pattern Quantification

Purpose: To extract quantitative metrics that enable direct comparison between simulated and experimental spatial data.

Procedure:

• Spatial Autocorrelation Analysis:
  • Calculate Moran's I for key cellular features (proliferation, hypoxia markers)
  • Compute Ripley's K-function for each cell type to assess clustering
  • Perform cross-K analysis between different cell types

• Neighborhood Analysis:

  • Apply k-nearest neighbors (k-NN) classification with k=10
  • Calculate cell-type enrichment/depletion in local neighborhoods
  • Quantify heterotypic and homotypic interaction frequencies

• Spatial Gradient Detection:

  • Identify concentration gradients of simulated secreted factors
  • Map diffusion patterns from source cells
  • Correlate gradient positions with cell state transitions

• Niche Identification:

  • Apply unsupervised clustering (DBSCAN) to identify recurrent cellular neighborhoods
  • Characterize niche composition and spatial distribution
  • Calculate relative abundance of each niche type

Experimental Spatial Omics Workflows

Spatial Proteomics Validation

Purpose: To experimentally validate protein-level predictions from the Voronoi model using multiplexed spatial proteomics.

Procedure:

• Sample Preparation:
  • Use FFPE tissue sections (4-5 μm thickness) on charged slides
  • Perform antigen retrieval using citrate-based buffer (pH 6.0)
  • Block endogenous peroxidase activity with 3% H₂O₂

• Multiplexed Staining:

  • Implement tyramide signal amplification (TSA) for fluorescence-based detection [76]
  • Alternatively, use metal-tagged antibodies for mass cytometry-based detection (IMC)
  • Include a nuclear stain (DAPI or iridium intercalator) for cell segmentation
  • Recommended panel design: 30-40 markers covering tumor, immune, and stromal compartments [75]

• Image Acquisition:

  • For fluorescence-based methods (PhenoCycler-Fusion): Capture multi-channel images at 20x magnification
  • For IMC: Ablate regions of interest at 1 μm resolution using a UV laser
  • Include control regions for background subtraction and signal normalization

• Image Processing:

  • Perform cell segmentation using nuclear marker + membrane/cytoplasmic expansion [76]
  • Extract single-cell protein expression values
  • Correct for technical artifacts (background fluorescence, spectral overlap)

Table 1: Benchmarking of Spatial Proteomics Cell Typing Tools

Tool	Methodology	Accuracy	Strengths	Limitations
TACIT [72]	Unsupervised thresholding	F1=0.75	Identifies rare cell types; No training data required	Requires predefined signatures
Tree-based Methods [76]	Supervised machine learning	High consistency	Outperforms other phenotyping tools	Requires manual annotations
SCINA [72]	Semi-supervised	Identifies 5/17 cell types	Works with minimal markers	Poor performance on rare populations
CELESTA [72]	Unsupervised	Identifies all expected types	Good for rare cells	Weaker correlation to reference (R=0.24)

Spatial Transcriptomics Correlation

Purpose: To validate gene expression patterns predicted by the Voronoi model using spatial transcriptomics.

Procedure:

• Platform Selection:
  • NanoString GeoMx WTA: For whole transcriptome analysis of predefined regions [73]
  • 10x Visium: For whole transcriptome with spot-based resolution
  • MERFISH/STOmics: For higher resolution targeted or whole transcriptome

• Region of Interest (ROI) Selection:

  • Segment tissue into tumor and stromal compartments based on morphology
  • Align ROIs with proteomic analysis areas for multimodal integration
  • Include replicates of each niche type identified in simulations

• Data Processing:

  • Perform quality control (sequencing depth, ROI integrity)
  • Normalize counts using RLE or TPM methods
  • Annotate cell types using reference-based (SingleR) or clustering approaches

• Integration with Proteomics:

  • Map transcriptomic and proteomic data through spatial alignment
  • Identify concordant and discordant protein-RNA pairs
  • Correlate spatial expression gradients

Table 2: Multi-Omics Integration Performance Benchmarks

Integration Method	Data Types	Concordance Rate	Key Applications
TACIT with multiomics [72]	Protein + RNA	81%	Cell type agreement between modalities
IMMUcan cross-platform [76]	mIF + IMC	High correlation	Major cell type localization
SNOWFLAKE GNN [77]	Protein + Morphology	AUC=0.973	Disease classification from follicles
GeoMx IPA+WTA [73]	Protein + RNA	Region-specific	Survival biomarker discovery

Multi-Omics Data Integration Framework

Spatial Data Alignment

Purpose: To establish a common coordinate framework for correlating in-silico predictions with experimental data.

Procedure:

• Coordinate System Registration:
  • Identify common landmarks between simulated and experimental data
  • Apply affine transformation to align coordinate systems
  • Resolve scale differences through proportional scaling

• Cellular Correspondence Mapping:

  • Match simulated cells to experimental cells based on spatial proximity
  • Resolve one-to-many mappings through probability assignment
  • Handle unmatched cells through gap analysis

• Multi-Modal Feature Integration:

  • Create unified feature matrix combining simulated and experimental measurements
  • Impute missing values using k-NN based on spatial neighbors
  • Normalize features to z-scores for cross-comparison

Model Validation Metrics

Purpose: To quantitatively assess the agreement between Voronoi model predictions and experimental observations.

Procedure:

• Spatial Accuracy Assessment:
  • Calculate Sørensen-Dice coefficient for cell type distributions
  • Compute Pearson correlation for spatial gradient patterns
  • Assess clustering similarity using Adjusted Rand Index

• Predictive Performance Validation:

  • Compare niche abundance predictions vs observations (RMSE)
  • Evaluate spatial segregation accuracy (Cross-K function difference)
  • Assess predictive power for rare cell populations (F1 score)

• Temporal Dynamics Validation:

  • Align simulation timepoints with experimental time courses
  • Compare trajectory similarities using Dynamic Time Warping
  • Assess state transition accuracy

Visualization and Analysis Workflows

Integrated Data Visualization

The following workflow diagram illustrates the comprehensive process for correlating in-silico output with spatial multi-omics data:

Signaling Pathway Mapping

The following diagram illustrates how spatial relationships from Voronoi models inform multi-omics analysis of signaling pathways in the tumor microenvironment:

The Scientist's Toolkit

Table 3: Essential Research Reagent Solutions for Spatial Multi-Omics Integration

Reagent/Category	Specific Examples	Function	Application Notes
Cell Segmentation Reagents	DAPI, iridium intercalator, membrane markers	Nuclear and cellular identification	Critical for single-cell resolution; DAPI for fluorescence, iridium for IMC
Antibody Panels (Proteomics)	CD3, CD20, CD68, cytokeratin, PD-L1	Cell type and state identification	Validate cross-reactivity; optimize dilution for each tissue type
Spatial Barcoding Systems	10x Visium, NanoString GeoMx DSP	Spatial localization of molecular signals	Choose resolution based on research question (whole tissue vs ROI)
Multiplexing Detection	TSA opal fluorophores, metal-tagged antibodies	Signal amplification and detection	TSA for fluorescence, metal tags for mass cytometry
Quality Control Reagents	IgG controls, housekeeping proteins, RNA spikes	Technical variation assessment	Essential for normalization and batch effect correction
Cell Typing Algorithms	TACIT, tree-based methods, SCINA, CELESTA	Automated cell annotation	Selection depends on available markers and reference data
Spatial Analysis Software	IFQuant, Squidpy, Giotto, SPATA2	Spatial pattern quantification	IFQuant web-based for scalability; others for specialized analyses

Application to Clinical Translation

Predictive Biomarker Discovery

Purpose: To leverage the integrated Voronoi-multi-omics framework for identifying spatially-informed biomarkers.

Procedure:

• Spatial Signature Development:
  • Extract spatial features associated with clinical outcomes from trained models
  • Apply LASSO-penalized Cox models for feature selection [75]
  • Validate signatures in independent cohorts using spatial proteomics/transcriptomics

• Resistance Mechanism Identification:

  • Correlate spatial niches with treatment resistance patterns
  • Identify cell-type combinations enriched in non-responders
  • Map spatial gradients of immune exclusion

• Therapeutic Target Prioritization:

  • Integrate drug target expression with spatial context
  • Identify targetable pathways in specific spatial niches
  • Predict drug penetration based on TME structure

Table 4: Clinically Validated Spatial Signatures from Multi-Omics Studies

Cancer Type	Spatial Signature	Clinical Endpoint	Performance
NSCLC [75]	Resistance: proliferating tumor cells, granulocytes, vessels	PFS on immunotherapy	HR=3.8, p=0.004
NSCLC [75]	Response: M1/M2 macrophages, CD4 T cells	PFS on immunotherapy	HR=0.4, p=0.019
HNSCC [73]	CD3e, CXCR5 in tumor compartment	Overall survival	Region-specific association
HNSCC [73]	CD44 in tumor compartment	Improved survival	Stroma-tumor compartment differences

Model-Informed Therapeutic Development

Purpose: To utilize the validated Voronoi model for predicting response to combination therapies.

Procedure:

• Therapy Simulation:
  • Parameterize model with drug-specific mechanisms (immune activation, checkpoint blockade)
  • Simulate combination therapy scenarios
  • Identify spatial biomarkers of response

• Patient Stratification:

  • Extract spatial features from patient biopsies
  • Match patients to simulated response profiles
  • Prioritize combination therapies based on spatial TME features

• Treatment Optimization:

  • Simulate dosing schedules and sequences
  • Identify spatial determinants of resistance
  • Design spatial-informed therapeutic strategies

The integration framework presented here enables continuous refinement of Voronoi tessellation models through iterative comparison with spatial multi-omics data. This virtuous cycle of prediction and validation accelerates both model development and biological discovery, ultimately enhancing our ability to understand and therapeutically target the spatial dynamics of cancer.

The variable response of cancers to immunotherapy, such as programmed death 1 (PD-1)-based treatments, represents a significant challenge in clinical oncology. A substantial proportion of patients, up to 40% in non-small cell lung cancer (NSCLC), exhibit primary resistance despite high PD-L1 expression, highlighting the urgent need for better predictive biomarkers [78] [75]. The spatial organization of cells within the tumor immune microenvironment (TIME) has emerged as a critical determinant of therapeutic outcome, necessitating analytical frameworks that can quantitatively capture these complex spatial relationships.

This Application Note details experimental and computational protocols for identifying and validating spatial signatures predictive of immunotherapy response, framed within the context of Voronoi tessellation cellular automaton tumor modeling. The Voronoi-based cellular automaton approach provides a robust geometrical framework for simulating proliferative tumor growth and analyzing tissue architecture, enabling researchers to model the TIME as a complex, dynamic, and self-organizing biosystem [3]. By integrating this computational framework with spatial multi-omics technologies, we establish a standardized pipeline for quantifying cellular spatial features and linking them to clinical outcomes.

Key Quantitative Findings from Spatial Multi-Omics Studies

Recent spatial multi-omics studies have identified specific cellular signatures within the TIME that correlate significantly with immunotherapy outcomes. The tables below summarize key quantitative findings from these investigations.

Table 1: Spatial Proteomics-Derived Cell-Type Signatures Associated with Immunotherapy Response in NSCLC

Signature Type	Component Cell Types	Hazard Ratio (HR)	P-value	Clinical Interpretation
Resistance	Proliferating tumor cells, Granulocytes, Vessels	3.8	0.004	Predictive of poor outcomes
Response	M1/M2 macrophages, CD4 T cells	0.4	0.019	Predictive of favorable outcomes

Table 2: Spatial Transcriptomics-Derived Gene Signatures Across Validation Cohorts

Signature Type	Yale Cohort HR	University of Queensland HR	University of Athens HR	Clinical Interpretation
Cell-to-Gene Resistance	5.3	2.2	1.7	Predictive of poor outcomes across multiple cohorts
Cell-to-Gene Response	0.22	0.38	0.56	Predictive of favorable outcomes across multiple cohorts

Table 3: Spatial Features in Hepatocellular Carcinoma (HCC) with Prognostic Value

Spatial Feature	Abbreviation	Prognostic Value	Biological Interpretation
Cellular diversity around stromal cells (mean)	StrDiv-M	Significant	Heterogeneity of cell types surrounding stromal elements
Distance between all cells (median)	CellDis-MED	Significant	Overall cellular density and packing
Variation coefficient of distance around stromal cells (median)	CvStrDis-MED	Significant	Regularity of stromal cell distribution patterns

Experimental Protocols

Spatial Multi-Omics Profiling of Tumor Tissue

Purpose: To comprehensively characterize the tumor immune microenvironment using spatially-resolved proteomic and transcriptomic technologies.

Materials:

• Formalin-fixed paraffin-embedded (FFPE) tumor tissue sections
• CODEX (CO-DEtection by indEXing) multiplexed tissue imaging system [78]
• GeoMx Digital Spatial Profiler (DSP) Whole Transcriptome Analysis panel [75]
• Antibody panels for protein detection (29-marker panel recommended) [78]
• Tissue staining and processing equipment

Procedure:

• Tissue Preparation: Cut FFPE tissue sections at 4-5μm thickness and mount on appropriate slides following manufacturer specifications.
• Multiplexed Protein Detection (CODEX):
  • Conjugate antibodies with metal isotopes as per CODEX protocol
  • Perform iterative staining and imaging cycles
  • Acquire high-resolution images of entire tissue sections
  • Process data to generate single-cell resolution protein expression maps [78]
• Spatial Transcriptomics (GeoMx DSP):
  • Hybridize tissue sections with RNA detection probes
  • Select regions of interest (ROIs) based on morphological features
  • Perform UV-cleavage of oligonucleotides from selected ROIs
  • Collect and prepare libraries for sequencing [75]
• Data Integration: Align spatial proteomic and transcriptomic datasets using coordinate-based registration methods.

Voronoi Tessellation and Cellular Neighborhood Analysis

Purpose: To quantify spatial relationships and organizational patterns within the tumor immune microenvironment using Voronoi tessellation and Delaunay triangulation.

Materials:

• Cell coordinate data from multiplexed imaging
• Computational software (Python, R) with spatial analysis libraries
• Voronoi diagram and Delaunay triangulation algorithms [3]

Procedure:

• Cell Position Extraction: Identify centroid coordinates for all segmented cells within the tissue sample.
• Voronoi Tessellation:
  • Generate Voronoi diagrams partitioning space into polygonal regions around each cell centroid
  • Calculate Voronoi cell areas, perimeters, and neighbor counts
  • Identify spatial anomalies through Voronoi cell property analysis [3]
• Delaunay Triangulation:
  • Construct Delaunay network connecting neighboring cells
  • Calculate edge lengths representing intercellular distances
  • Analyze network properties including clustering coefficients and connectivity [79]
• Cellular Neighborhood Identification:
  • Apply clustering algorithms to identify recurrent cellular communities
  • Annotate neighborhoods based on predominant cell-type enrichment
  • Calculate neighborhood abundance and spatial relationships [78]

Predictive Signature Development Using Machine Learning

Purpose: To develop and validate spatial signatures predictive of immunotherapy response using robust machine learning approaches.

Materials:

• Processed spatial feature matrix from Voronoi and cellular neighborhood analyses
• Clinical outcome data (progression-free survival, overall survival)
• Computational environment with machine learning libraries (scikit-learn, survival packages)

Procedure:

• Feature Preparation:
  • Compile spatial features including cell-type proportions, Voronoi metrics, and neighborhood compositions
  • Split dataset into training (Yale cohort) and validation (UQ and Greek cohorts) sets [78]
• Signature Training:
  • Perform multiple random splits of training cohort into tenfolds
  • Apply LASSO-penalized Cox models for feature selection with non-negative constraints for resistance signatures and non-positive constraints for response signatures
  • Identify consistently selected features across all splits [78] [75]
• Model Validation:
  • Apply trained model to independent validation cohorts
  • Evaluate predictive performance using hazard ratios and log-rank tests
  • Assess signature stability across different patient populations

Computational Framework Integration

The integration of Voronoi tessellation cellular automaton modeling with spatial multi-omics data creates a powerful framework for understanding and predicting immunotherapy response. The cellular automaton model simulates proliferative tumor growth as a function of time using an underlying Voronoi lattice, which provides spatial isotropy and avoids the anisotropies associated with more ordered lattices [3]. This approach enables the simulation of tumor growth from approximately 1000 real cells through several clinically designated time points to fully developed macroscopic size.

The Voronoi framework allows for the incorporation of several key features essential for TIME modeling:

• Spatial isotropy for unbiased neighborhood analysis
• Adaptive grid resolution for multi-scale modeling
• Mechanical confinement effects through boundary conditions
• Heterogeneous environments including vascular networks

The workflow below illustrates the integration of experimental data generation with computational modeling:

Spatial Analysis Workflow: From Tissue to Clinical Validation

Research Reagent Solutions

Table 4: Essential Research Reagents for Spatial Multi-Omics Analysis

Reagent/Technology	Primary Function	Application Context
CODEX Multiplexed Imaging System	High-resolution spatial proteomics	Simultaneous detection of 40+ protein markers in intact tissue [78]
GeoMx Digital Spatial Profiler	Spatially-resolved whole transcriptome analysis	Region-specific RNA sequencing from tissue sections [75]
MIBI (Multiplexed Ion Beam Imaging)	High-dimensional tissue imaging with subcellular resolution	Simultaneous measurement of protein expression while preserving spatial information [80]
Voronoi Tessellation Algorithm	Spatial partitioning and neighborhood analysis	Quantification of cellular spatial relationships and microenvironment architecture [3] [79]
Delaunay Triangulation	Network modeling of cell adjacency	Analysis of direct cell-cell interactions and spatial connectivity [3] [79]
LASSO-Penalized Cox Models	Feature selection for predictive signature development	Identification of most prognostic features while preventing overfitting [78] [75]

Signaling Pathways and Cellular Interactions

The spatial signatures identified through these methodologies reveal critical functional relationships within the TIME. The resistance signature characterized by proliferating tumor cells, granulocytes, and vessels suggests a microenvironment conducive to tumor growth and immune suppression. Conversely, the response signature enriched with M1/M2 macrophages and CD4 T cells indicates effective immune activation. The following diagram illustrates the key cellular interactions and their functional consequences:

Cellular Determinants of Immunotherapy Response

The integration of Voronoi tessellation cellular automaton models with spatial multi-omics technologies provides a powerful framework for identifying robust predictive signatures of immunotherapy response. The protocols outlined in this Application Note enable researchers to quantitatively characterize the spatial architecture of the tumor immune microenvironment and link these features to clinical outcomes. The resistance and response signatures validated across multiple cohorts offer promising biomarkers for patient stratification in immunotherapy. Future applications of this approach may include optimization of combination therapies and identification of novel therapeutic targets based on spatial localization patterns within the TIME.

The Voronoi Tessellation Cellular Automaton (VTCA) represents a significant advancement in the mathematical oncology landscape, providing a spatially explicit, cell-based framework for simulating complex tumor phenomena. This approach integrates the geometric realism of Voronoi tessellations with the behavioral rules of cellular automata to create high-fidelity digital twins of patient-specific tumors. A digital twin, in this context, is a dynamic computational model of a biological system that is continuously updated with patient data to simulate and predict disease progression and treatment response. The core strength of the VTCA framework lies in its ability to model individual cells as discrete agents whose interactions and fates are governed by a dynamically evolving, cell-specific microenvironment. This allows for the explicit simulation of clonal evolution and spatial competition, which are fundamental drivers of therapeutic resistance and tumor relapse [5].

Traditional cellular automaton models often rely on uniform, rigid lattices (e.g., square or hexagonal grids) that poorly approximate the heterogeneous and dynamic packing of real epithelial tissues. The VTCA paradigm addresses this limitation by using a Voronoi tessellation, which partitions space into convex, polygonal cells based on a set of seed points. This generates a more biologically plausible representation of a cellular monolayer, where each polygon corresponds to one cell, and its geometry is dictated by the positions of its neighbors. This tessellation can be generated from the centers of mass of cell nuclei identified in patient biopsy images, creating an immediate bridge between clinical data and the in-silico model [81]. The resulting digital twin can therefore capture the inherent topological disorder and heterogeneity of real tissues, which is crucial for predicting emergent behaviors like invasion and metastasis.

VTCA Model Components and Numerical Implementation

Geometric and Topological Foundation

The VTCA model is built upon a Voronoi diagram, ( V(S) ), defined for a set of site points, ( S ), within a bounded space. For each site point ( p \in S ), its Voronoi region, ( VR(p, S) ), consists of all points in the space closer to ( p ) than to any other site in ( S ) [55]. The ensemble of these regions forms the cellular landscape of the digital twin. The Delaunay triangulation, ( DT(S) ), which is the geometric dual of the Voronoi diagram, provides critical adjacency information, defining the neighborhood network for cell-cell interactions [55]. In practice, for analyzing epithelial tissues, the site points ( S ) are often derived from the centers of mass of cell nuclei (CMVT) segmented from fluorescence microscopy images [81]. Studies on thousands of MDCK-II epithelial cells have shown that CMVT can reproduce key morphological measures—including cell area, perimeter, and number of neighbors—with an error typically between 10% and 15% compared to direct membrane staining, establishing its utility for generating initial model geometry [81].

Cellular State and Behavioral Rules

Each cell in the VTCA is an autonomous agent characterized by a set of internal states and is equipped with a micro-environment response network. This network can be conceptualized as an internal decision-making circuit. A powerful method to implement this is through a feed-forward artificial neural network that maps local environmental variables (e.g., oxygen tension, growth factor concentrations, mechanical pressure) to cellular behavioral outputs (e.g., proliferate, become quiescent, apoptose, migrate) [5]. The specific response of a cell is determined by the connection weights and thresholds within this network, which constitute its genotype. This genotype is subject to mutation upon cell division, introducing phenotypic diversity and enabling Darwinian evolution within the tumor population [5]. The table below summarizes the core variables that define a cell's state and its microenvironment.

Table 1: Key State Variables in a VTCA Model

Variable Category	Specific Examples	Biological Significance
Cell Internal State	Cell cycle phase, mutational profile, phenotypic signature (proliferative, migratory, glycolytic)	Determines a cell's current activity and potential future actions.
Local Microenvironment	Oxygen concentration (( pO_2 )), glucose, growth factors, extracellular pH, mechanical pressure (( P ))	Regulates cell behavior through biochemical and biophysical signaling.
Geometric Properties	Cell area (( A )), cell perimeter (( P )), number of neighbors	Influences and is influenced by mechanical feedback and proliferation.

Coupling Mechanics and Growth

The VTCA framework can be extended to incorporate mechanical regulation of growth, a critical factor in confluent tissues. This is often achieved by coupling the model with a Self-Propelled Voronoi (SPV) model or a vertex model. In such a hybrid approach, the tissue is treated as a confluent cell monolayer, and its energy is described by a functional such as: [ E = \sum{\alpha} \frac{K{\alpha}}{2} (A{\alpha} - A{\alpha}^0)^2 + \sum{\alpha} \frac{\Gamma{\alpha}}{2} (P{\alpha} - P{\alpha}^0)^2 + \text{...} ] where ( K{\alpha} ) and ( \Gamma{\alpha} ) are mechanical moduli, ( A{\alpha} ) and ( P{\alpha} ) are the current area and perimeter of cell ( \alpha ), and ( A{\alpha}^0 ) and ( P{\alpha}^0 ) are their target values [27]. The hydrostatic pressure on a cell can be derived as ( P{\alpha} = -K{\alpha}(A{\alpha} - A{\alpha}^0) ). This pressure, or a related mechanical variable, can then be used to modulate the cell's birth and death rates. For instance, a simple linear coupling gives: [ r{\alpha}^b = r^b0 + \lambda \max(0, A{\alpha} - A^0) ] [ r{\alpha}^d = r^d0 - \lambda \min(0, A{\alpha} - A^0) ] where ( r{\alpha}^b ) and ( r{\alpha}^d ) are the division and apoptosis rates, ( r^b0 ) and ( r^d0 ) are basal rates, and ( \lambda ) is a coupling coefficient [27]. This formalism allows the digital twin to simulate phenomena such as density-dependent inhibition of growth and the impact of solid stress on tumor development.

Application to Personalized Virtual Clinical Trials

The Virtual Clinical Trial Pipeline

Virtual Clinical Trials (VCTs) use computational models to inexpensively and rapidly test a vast number of therapeutic strategies on a cohort of in-silico patients, or "digital twins," to determine the most effective personalized treatment plans [82]. The workflow for deploying a VTCA in a VCT is a multi-stage process. The following diagram illustrates the integrated feedback loop of this pipeline.

Protocol: Generating a Patient-Specific Digital Twin

Objective: To construct an initial, spatially accurate VTCA digital twin from a standard patient tissue biopsy. Materials:

• Patient tissue biopsy sample (e.g., from a solid tumor).
• Reagents for fluorescence staining: Antibodies against β-catenin or E-cadherin (for membrane staining), Hoechst or DAPI (for nuclear staining) [81].
• Confocal fluorescence microscope.
• Image analysis software (e.g., Fiji/ImageJ, CellProfiler).

Procedure:

• Tissue Preparation and Staining: Fix the biopsy sample and perform immunofluorescence staining to label cell membranes (e.g., with β-catenin antibodies) and nuclei (e.g., with Hoechst) [81].
• Image Acquisition: Acquate high-resolution, multi-channel z-stack images of the tissue using a confocal microscope.
• Image Segmentation: a. Use the nuclear (Hoechst) channel to automatically segment individual cell nuclei. b. Calculate the center of mass (centroid) for each segmented nucleus. This set of points, ( S ), will serve as the sites for the Voronoi tessellation.
• Tessellation and Model Initialization: a. Compute the Voronoi diagram ( V(S) ) and its dual, the Delaunay triangulation ( DT(S) ), from the nuclei centroids. b. Assign initial cell states based on morphological data from the membrane channel. For example, a cell's target area ( A^0 ) can be set to the area of its Voronoi region, and its baseline phenotype can be inferred from its location within the tissue architecture. c. Calibrate the initial oxygen and nutrient fields based on diffusion models from the tissue boundary, assuming a higher concentration at the periphery that decreases towards the core.

Validation: Compare the morphological outputs of the initial, un-evolved VTCA—such as the distributions of cell area, perimeter, and number of neighbors—against the same measures obtained from the segmented membrane stain. The model should qualitatively and quantitatively (within ~10-15% error) match the actual tissue morphology before proceeding to simulation [81].

Protocol: Simulating Therapy and Analyzing Outcomes

Objective: To use the calibrated digital twin to simulate tumor response to a specific therapeutic regimen and extract quantitative biomarkers of efficacy.

Procedure:

• Therapy Parameterization: Define the mechanism of action of the therapy within the model.
  • For a cytotoxic chemotherapy, this may involve increasing the probability of apoptosis for cycling cells.
  • For a targeted therapy, it may involve altering the parameters of the micro-environment response network for cells carrying a specific mutation.
• Simulation Execution: Run the VTCA model forward in time under the influence of the therapeutic perturbation. A typical simulation might span several hundred time steps, each representing several hours of real time.
• Outcome Analysis: At regular intervals, quantify key tumor metrics:
  • Total tumor cell count.
  • Tumor morphology (e.g., compact vs. fingered invasive margin) [5].
  • Clonal diversity and the emergence of resistant subpopulations.
  • Fraction of quiescent vs. proliferating cells.
• Virtual Trial Cohort: Repeat the calibration (Protocol 3.2) and simulation process for a large number of digital twins, each representing a different virtual patient, to explore inter-patient heterogeneity and identify patient subgroups that respond best to the tested therapy.

Table 2: Quantitative Metrics for Evaluating Therapeutic Efficacy in a VTCA-based VCT

Metric	Calculation/Description	Interpretation in VCT
Tumor Volume	Sum of the areas of all tumor cells.	Primary measure of gross tumor response.
Invasive Index	Ratio of tumor perimeter to sqrt(tumor area).	Quantifies morphological aggressiveness; a higher index indicates a more fingered, invasive structure [5].
Resistant Fraction	Proportion of cells with genotypes conferring therapy resistance.	Predicts the likelihood of relapse.
Selection Pressure	Rate of change in clonal diversity (Shannon index).	Measures the intensity of evolutionary competition driven by therapy.

The Scientist's Toolkit: Essential Research Reagents and Models

The development and application of VTCAs rely on a combination of biological, computational, and analytical tools. The following table details key reagents and their functions in this research domain.

Table 3: Research Reagent Solutions for VTCA Model Development and Validation

Category / Item	Specific Example / Software	Function in VTCA Research
Biological Stains	Hoechst / DAPI (Nuclear stain), β-catenin antibody (Membrane stain) [81]	Provides the experimental image data required to generate and validate the initial geometry of the digital twin.
Cell Lines & Models	MDCK-II (canine kidney epithelium) cells, MCF10A (human breast epithelium) cells [81] [82]	Well-characterized in vitro models for growing normal acini and tumor spheroids to benchmark model predictions.
Mechanical Models	Self-Propelled Voronoi (SPV) Model, Vertex Model [27]	Provides the computational framework for incorporating tissue mechanics and density-dependent growth into the VTCA.
Tessellation Software	CellGPU, Voro++ [27]	Libraries for efficiently generating and managing dynamic Voronoi tessellations in simulations.
Auto-Differentiation	PyTorch-based AD of Voronoi Tessellation [55]	Enables gradient-based optimization of model parameters (e.g., seed point locations) by making the tessellation process differentiable.

Visualization of Core Model Logic

The dynamic behavior of a VTCA is governed by the continuous interaction between a cell's genotype, its perception of the local microenvironment, and its resulting actions. The following diagram outlines this core decision-making logic for an individual cell within the digital twin.

Conclusion

Voronoi Tessellation Cellular Automaton models represent a paradigm shift in computational oncology, moving from simplistic abstractions to spatially realistic, multi-scale simulations of tumor behavior. By providing an isotropic lattice that naturally captures tissue microstructure and cellular heterogeneity, VTCAs offer a powerful platform for integrating fundamental radiobiological principles. The ability to simulate complex treatment dynamics, from radiation response to oncolytic virotherapy, and to optimize delivery strategies in silico, positions these models as invaluable tools for therapy development. Future directions hinge on tighter integration with high-throughput, spatially resolved multi-omics data from initiatives like the IMMUcan consortium, enhancing their predictive accuracy for personalized medicine. The emergence of auto-differentiable techniques and the 'digital twin' concept promises a new era where VTCAs can dynamically inform clinical decision-making, ultimately enabling the design of more effective and personalized cancer therapies.

References

• 1. Voronoi diagram [https://en.wikipedia.org/wiki/Voronoi_diagram]
• 2. Hybrid Models of Tumor Growth - PMC - PubMed Central - NIH [https://pmc.ncbi.nlm.nih.gov/articles/PMC3057876/]
• 3. Cellular automaton of idealized brain tumor growth dynamics [https://www.sciencedirect.com/science/article/abs/pii/S0303264799000891]
• 4. Generalized Voronoi Tessellation as a Model of Two- ... [https://link.springer.com/article/10.1007/s11538-009-9498-3]
• 5. An Evolutionary Hybrid Cellular Automaton Model of Solid ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2652069/]
• 6. Bio-inspired structural optimization of three-dimensional ... [https://www.sciencedirect.com/science/article/pii/S0264127524008761]
• 7. Delaunay triangulation [https://en.wikipedia.org/wiki/Delaunay_triangulation]
• 8. Working with Delaunay Triangulations - MATLAB & Simulink [https://www.mathworks.com/help/matlab/math/delaunay-triangulation.html]
• 9. Visualizing Delaunay Triangulation [https://ianthehenry.com/posts/delaunay/]
• 10. Isotropic Growth on Cartesian Voxel Grids with von ... [https://link.springer.com/article/10.1007/s13632-024-01121-x]
• 11. The effects of regularity on the geometrical properties ... [https://www.sciencedirect.com/science/article/abs/pii/S0378437114002064]
• 12. Cellular automaton [https://en.wikipedia.org/wiki/Cellular_automaton]
• 13. A multi-grid Cellular Automaton model for simulating ... [https://www.sciencedirect.com/science/article/abs/pii/S2214860421004449]
• 14. Easy, fast and reproducible Stochastic Cellular Automata ... [https://peercommunityjournal.org/item/10_24072_pcjournal_466/]
• 15. Mechanistic Radiobiological Models for Repair of Cellular ... [https://www.sciencedirect.com/science/article/abs/pii/S0065327614000045]
• 16. To exploit the 5 'R' of radiobiology and unleash the 3 ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC7222228/]
• 17. intraoperative radiotherapy and normal tissue tolerance - Emami [https://tcr.amegroups.org/article/view/4270/4629]
• 18. Predicting the Tumour Response to Radiation by Modelling ... [https://www.mdpi.com/2313-433X/9/6/124]
• 19. Three-dimensional tumor growth in time-varying chemical ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6712764/]
• 20. Visualizing Data - NCI [https://www.cancer.gov/about-nci/organization/cbiit/training/learn/visualize]
• 21. Elements must meet minimum color contrast ratio thresholds [https://dequeuniversity.com/rules/axe/4.8/color-contrast]
• 22. Contrast Checker [https://webaim.org/resources/contrastchecker/]
• 23. Influence of Hypoxia on Tumor Heterogeneity, DNA Repair ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12294107/]
• 24. A micro-metabolic rewiring assay for assessing hypoxia ... [https://www.sciencedirect.com/science/article/pii/S2452199X25000866]
• 25. A single-cell network approach to decode metabolic ... [https://www.nature.com/articles/s41540-025-00506-0]
• 26. Spatial Data Analysis of Patient Tumor Samples [https://blog.crownbio.com/spatial-data-analysis-patient-tumor-samples]
• 27. Mechanically-driven growth and competition in a Voronoi ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11118596/]
• 28. Comparing individual-based approaches to modelling the self ... [https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005387]
• 29. Toward an Ising Model of Cancer and Beyond - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC3151151/]
• 30. Emergent Behaviors from a Cellular Automaton Model for ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3245298/]
• 31. A Cellular Automaton Model for Tumor Dormancy: Emergence ... [https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0109934]
• 32. Tumor microenvironment and metabolic reprogramming [https://www.sciencedirect.com/science/article/pii/S1040842825003890]
• 33. In vitro metabolic zonation through oxygen gradient on a chip [https://www.nature.com/articles/s41598-019-49412-6]
• 34. Oxygen Gradients in the Microcirculation - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC3115497/]
• 35. Radiobiological Modeling with Monte Carlo Tools [https://pmc.ncbi.nlm.nih.gov/articles/PMC12276521/]
• 36. Metabolic interplays between the tumor and the host shape ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12140436/]
• 37. Quantitative Analyses of the Tumor Microenvironment ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6167550/]
• 38. Data-driven mathematical modeling and quantitative ... [https://www.sciencedirect.com/science/article/pii/S0898122122000967]
• 39. Modeling tumor dynamics and predicting response to ... [https://www.nature.com/articles/s41540-025-00593-z]
• 40. Quantitative cancer-immunity cycle modeling for predicting ... [https://www.nature.com/articles/s41540-025-00513-1]
• 41. sciencedirect.com/topics/medicine-and-dentistry/ radiation -dose... [https://www.sciencedirect.com/topics/medicine-and-dentistry/radiation-dose-fractionation]
• 42. Frontiers | A comparative analysis of hypofractionated versus... [https://www.frontiersin.org/journals/oncology/articles/10.3389/fonc.2025.1552346/full]
• 43. Dose-volume histogram-based comparison of conventional and... [https://pubmed.ncbi.nlm.nih.gov/40516018/]
• 44. Organoid and Spheroid Tumor Models - PubMed Central - NIH [https://pmc.ncbi.nlm.nih.gov/articles/PMC7922036/]
• 45. Comparative analysis of 3D-culture techniques for ... [https://www.nature.com/articles/s41598-025-13588-x]
• 46. 3D Cell Culture Methods and Applications [https://elveflow.com/microfluidic-reviews/3d-cell-culture-methods-and-applications-a-short-review/]
• 47. 3D multicellular systems in disease modelling [https://www.frontiersin.org/journals/cell-and-developmental-biology/articles/10.3389/fcell.2023.1083175/full]
• 48. Mathematical models of tumor growth using Voronoi ... [https://profiles.wustl.edu/en/publications/mathematical-models-of-tumor-growth-using-voronoi-tessellations-i]
• 49. Multiplexed imaging analysis of the tumor-immune ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8271023/]
• 50. Frontiers | Oncolytic in cancer treatment: challenges and... virotherapy [https://www.frontiersin.org/journals/immunology/articles/10.3389/fimmu.2023.1308890/full]
• 51. Observations from a Voronoi Cell-Based model [https://www.sciencedirect.com/science/article/abs/pii/S0022519319304217]
• 52. Observations from a Voronoi Cell-Based model [https://pubmed.ncbi.nlm.nih.gov/31626813/]
• 53. Octyl itaconate enhances VSVΔ51 oncolytic by... virotherapy [https://www.nature.com/articles/s41467-024-48422-x?error=cookies_not_supported]
• 54. Recent advances in vesicular stomatitis virus-based oncolytic ... [https://pubmed.ncbi.nlm.nih.gov/29143726/]
• 55. A Method for Auto-Differentiation of the Voronoi Tessellation [https://arxiv.org/html/2312.16192v3]
• 56. Characterizing Interfaces by Voronoi Tessellation - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC10344600/]
• 57. Cellular automaton of idealized brain tumor growth dynamics [https://www.sciencedirect.com/science/article/pii/S0303264799000891]
• 58. Rule-Based Simulation of Multi-Cellular Biological Systems ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3045734/]
• 59. Controllable design and evaluation of Voronoi-based ... [https://www.sciencedirect.com/science/article/abs/pii/S2352492824032069]
• 60. A Method for Auto-Differentiation of the Voronoi Tessellation [https://arxiv.org/html/2312.16192v1]
• 61. Random nature of epithelial cancer cell monolayers - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC9090488/]
• 62. Mathematical models of tumor growth using Voronoi ... [https://pubmed.ncbi.nlm.nih.gov/26737283/]
• 63. Digital Pathology Analysis Quantifies Spatial Heterogeneity ... [https://www.frontiersin.org/journals/physiology/articles/10.3389/fphys.2020.583333/full]
• 64. Grading of invasive breast carcinoma through Grassmannian ... [https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0185110]
• 65. (PDF) Voronoi diagram: An adaptive spatial tessellation for... [https://www.academia.edu/73561925/Voronoi_diagram_An_adaptive_spatial_tessellation_for_processes_simulation]
• 66. Mathematical modelling reveals cellular dynamics within ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC7455028/]
• 67. An Overview on Spheroid and Organoid Models in Applied ... [https://www.mdpi.com/2413-4155/7/1/27]
• 68. Integrating simulation and experimental validation of ... [https://www.sciencedirect.com/science/article/pii/S0928098725003677]
• 69. A 3D spheroid model for assessing nanocarrier-based ... [https://www.nature.com/articles/s44385-025-00041-x]
• 70. Development of 3D organoid models to study aggressive ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12231088/]
• 71. Breakthroughs and challenges of organoid models for ... [https://www.nature.com/articles/s41420-025-02505-w]
• 72. Deconvolution of cell types and states in spatial multiomics ... [https://www.nature.com/articles/s41467-025-58874-4]
• 73. The development of a high-plex spatial proteomic ... [https://www.nature.com/articles/s41698-025-00963-0]
• 74. A systematic mapping study of semantic technologies in ... [https://www.sciencedirect.com/science/article/pii/S1532046425000383]
• 75. Spatial signatures for predicting immunotherapy outcomes ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12513832/]
• 76. Multi-modal image analysis for large-scale cancer tissue ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12539258/]
• 77. Spatial morphoproteomic features predict disease states from ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12356389/]
• 78. Spatial signatures for predicting immunotherapy outcomes ... [https://www.nature.com/articles/s41588-025-02351-7]
• 79. Leveraging deep learning to discover interpretable cellular ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12163547/]
• 80. Multiplexed imaging analysis of the tumor-immune ... [https://www.nature.com/articles/s42003-021-02361-1]
• 81. Limits of Applicability of the Voronoi Tessellation ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5122581/]
• 82. from spatial models of tumour spheroids, to organoids, to tissues [https://pmc.ncbi.nlm.nih.gov/articles/PMC5805971/]