Spatial interactions modulate tumor growth and immune infiltration

Abstract

Direct observation of immune cell trafficking patterns and tumor-immune interactions is unlikely in human tumors with currently available technology, but computational simulations based on clinical data can provide insight to test hypotheses. It is hypothesized that patterns of collagen formation evolve as a mechanism of immune escape, but the exact nature of the interaction between immune cells and collagen is poorly understood. Spatial data quantifying the degree of collagen fiber alignment in squamous cell carcinomas indicates that late stage disease is associated with highly aligned fibers. Here, we introduce a computational modeling framework (called Lenia) to discriminate between two hypotheses: immune cell migration that moves 1) parallel or 2) perpendicular to collagen fiber orientation. The modeling recapitulates immune-ECM interactions where collagen patterns provide immune protection, leading to an emergent inverse relationship between disease stage and immune coverage. We also illustrate the capabilities of Lenia to model the evolution of tumor progression and immune predation. Lenia provides a flexible framework for considering a spectrum of local (cell-scale) to global (tumor-scale) dynamics by defining a kernel cell-cell interaction function that governs tumor growth dynamics under immune predation with immune cell migration. Mathematical modeling provides important mechanistic insights into cell interactions. Short-range interaction kernels provide a mechanism for tumor cell survival under conditions with strong Allee effects, while asymmetric tumor-immune interaction kernels lead to poor immune response. Thus, the length scale of tumor-immune interactions drives tumor growth and infiltration.

Figure 1.. Lenia as a cancer model:

(A) Lenia artificial life virtual creatures, reproduced from ref. . (B) List of characteristics possible to produce in Lenia by varying the growth dynamics function (C) or the interaction kernel function (D). (E) adding a fraction of the growth field at each time step to the cell density, forms the Lenia update rule, as seen in eqn. 1. (F) An example snapshot of a simulation shows the density of cells at each lattice location where the interaction kernel (G) specifies the nature of interaction of cells depending on their distance from each other. (H) The density potential, interpreted as a weighted average of interactions at each lattice location, is calculated as the convolution of $\mathbf{A}(\mathbf{x})$ and $\mathbf{K}(\mathbf{x})$. (I) the growth field is calculated by applying a growth map to the density potential.

Figure 2.. Collagen alignment in HNSCC

(A) Pipeline of image analysis to computational model using second harmonic generation to determine collagen density and subsequently the alignment and microenvironmental gradient. (B) Alignment collagen fibers is determined using OrientationJ ImageJ plugin; alignment tends to increase by disease stage.

Figure 3.. Immune cell trafficking model in Particle Lenia

(A) Hypothesized models of the influence of collagen fiber alignment on immune cell trafficking. See corresponding Supplementary Video S1. (B) example simulation showing immune cell infiltration (circles) with track indicating path taken. Cells seed on all four sides, color-coded by initial side (see panel F). (C, D) Simulated immune coverage for perpendicular (C) or parallel (D) immune trafficking model, colored by stage. (E) Representative sample ROI alignment and density of collagen fibers (using OrientationJ ImageJ plugin) ordered by disease stage. (F) Immune coverage, color-coded by immune initial condition of which of the four sides: left, right, top, bottom (see legend). The color of each pixel is determined the most trafficked side. (G) Immune coverage as the percentage of all four sides covered by immune surveillance. For example, if no immune cells reached this pixel it is white, if all four sides reached this pixel, it is dark purple. See corresponding Supplementary Video S2.

Figure 4.. Cancer growth dynamics in Lenia

(A) Average cell density over time for varying kernel sizes, for deterministic and stochastic Lenia. Averages over stochastic runs are shown in dotted lines and individual trajectories are shown in translucent lines. ODE solution (shown in blue) matches simulation with a well-mixed kernel. (B) State of field at half maximal capacity for different kernel sizes for Stochastic Lenia (SL) or Deterministic Lenia (DL). Density potential and Growth distribution shown in inset. (C) State of field at progressive time points for kernel size 4 for deterministic and stochastic Lenia. (D) Kernels used in models in figures (B).

Figure 5.. Short-range interaction kernels are more robust to Allee effects

(A) Average cell density over time for varying kernel sizes, for deterministic and stochastic Lenia. Averages over stochastic runs are shown in dotted lines and individual trajectories are shown in translucent lines. ODE solution shown in blue. (B) Long-term tumor fate for deterministic Lenia where trajectories that eventually result in extinction are shown in red, or growth shown in blue (colorbar indicates change in tumor size at $t^{\mathrm{*}}=3$. (C, D) Spatial maps shown for deterministic Lenia (C) and stochastic Lenia (D) at time $t^{\mathrm{*}}=3$. See corresponding Supplementary Video S3.

Figure 6.. Competition dynamics in Lenia

(A) Schematic representation of the tumor-immune predator-prey model, illustrating various kernel sizes for tumor growth $\left(K_{11}\right)$, predation distance $\left(K_{12}\right)$, immune recruitment $\left(K_{21}\right)$, and density-dependent death distance in the immune system $\left(K_{22}\right)$. (B, C) Impact of different predation distances $\left(K_{12}\right)$ and immune recruitment $\left(K_{21}\right)$ on the integration results of the tumor-immune predator-prey model, with other kernels held well-mixed. (D) Effect of varying predation and recruitment ($K_{12}$ and $K_{21}$) on the time to tumor regression (considered when there are 2% tumor cells in the area) and the peak value in the immune response. (E) Spatial outcomes of the tumor and immune response in the model, showcasing the influence of different kernel sizes. See corresponding Supplementary Video S4. Unless otherwise noted, parameters used are $=5,b=12,g=1.5,d=1,L=0.08$ (see eqns. 6–7).