Integrating digital pathology and mathematical modelling to predict spatial biomarker dynamics in cancer immunotherapy

The study addresses the lack of predictive quantitative methods to guide early-phase oncology trial decisions, particularly the optimal timing of on-treatment biopsies to assess drug mode of action, as well as identifying combination therapies and patient selection. The authors propose an integrated approach that combines spatial biomarkers derived from IHC-based digital pathology with a mechanistic agent-based model to predict the spatiotemporal dynamics of cytotoxic CD8+ T cells in tumor tissue. By leveraging paired baseline and on-treatment biopsies, the goal is to move beyond expert best-guess scheduling to a data-driven, quantitatively justified biopsy timepoint and to enable hypothesis-driven exploration of treatment mechanisms.

Prior work shows spatial biomarkers from IHC images correlate with prognosis, recurrence, and treatment response via statistical and machine-learning methods (e.g., Barua 2018; Corredor 2019; Brown 2014; Saltz 2018; Schwen 2018). Agent-based models (ABMs) of tumor–immune interactions exist (reviewed by Norton 2019; Kather et al. 2017, 2018; Gong 2017), but validation has largely been qualitative or at single timepoints, with limited quantitative spatial comparisons (e.g., use of radial distribution functions by Alfonso 2016). The authors highlight the absence of a widely accepted quantitative measure to compare spatial cell distributions across sets of clinical images, motivating their derivation of the spatial agreement measure (SAM) and an augmented variability metric (VarSAM). They also cite evidence suggesting core-needle biopsies can reflect intra-tumor heterogeneity (Mani 2016).

Data and preprocessing: The dataset comprised 71 patients from two Simlukafusp trials (NCT03063762, NCT02627274) with paired baseline and on-treatment biopsies; 44 pairs contained sufficient tissue for modeling (≥4 tiles per timepoint). Patients were split into training (n=37) and test (n=7). IHC sections stained for CD8 and Ki67 were digitized and analyzed with an in-house ML algorithm (~90% classification accuracy) that outputs cell types (proliferating/non-proliferating CD8+ T cells, proliferating tumor cells) and spatial coordinates. Biopsies were tiled into minimally overlapping 100×100 cell-width fields with ≥90% tissue overlap and ≤10% tile overlap. Because tumor detection focused on Ki67+ cells, empty annotated tumor space was filled to reach ~70% packing density assuming ~5 μm cell width. Each tile’s baseline spatial map of CD8 and tumor cells, plus assigned cell properties (e.g., proliferation capacity, stemness probability, engagement status; see Table 3 in paper), formed the model input. Agent-based model (ABM): Implemented in Matlab, based on Kather et al. The grid-lattice ABM includes tumor and immune (CD8) cells that can proliferate, migrate, die, and interact via probabilistic rules at each time step, with killing requiring proximity. Extensions included periodic boundary conditions (to reflect continuity of tissue) and increased tumor-cell resilience (multiple immune hits required for death). Key immune parameters explored were proliferation probability (IMpprol), killing probability (IMpkill), natural death probability (IMpdeath), randomness of biased random walk (IMrwalk), and influx terms; tumor parameters included proliferation, migration, death, and capacities. Simulations used baseline tiles as initial conditions and generated simulated biopsy snapshots at 24-hour intervals up to each patient’s actual on-treatment timepoint. Spatial agreement measure (SAM) and VarSAM: Spatial distributions were summarized by normalized radial distribution functions (RDFs) computed on the grid (square taxicab paired correlation). For a given patient, SAM quantifies overlap between distributions of RDF values from simulated tiles and observed on-treatment tiles across inter-cell distances, using acceptance bands defined as observed RDF ranges ±20% and a threshold (threshSAM) on the proportion of simulated RDFs within band per distance; SAM is the fraction of distances meeting the threshold. VarSAM captures agreement in variability within the first 15 distance units as the ratio of observed-to-simulated RDF ranges (min of the two ratios), emphasizing short-range clustering features. Cases with <10 CD8 cells per tile were excluded from SAM and instead compared via two-sample t-tests on CD8 counts. Sensitivity analysis: Parameters (13 total; see Table 4) were varied independently over physiologically plausible ranges (linear or logarithmic sampling). Six representative baseline tiles spanning 7–~800 CD8 cells were simulated with five stochastic repeats. Sensitivities of CD8 and tumor counts at day 20 were computed as normalized changes over parameter ranges, revealing which parameters most influenced outcomes. Optimization: 1000 parameter sets for four immune-related parameters (IMpprol, IMpdeath, IMrwalk, IMpkill) were sampled over defined ranges (Table 4) and evaluated on all training tiles with five stochastic repeats, ending simulations at each patient’s actual on-treatment time. SAM served as the objective; parameter sets with mean SAM >0.7 across patients and VarSAM >0.3 were accepted, yielding a population of accepted parameter sets to capture inter-individual variability. Null benchmark (LOCF): As a baseline without modeling, the last-observation-carried-forward approach used baseline RDFs to represent on-treatment RDFs. SAM/VarSAM were computed using the same thresholds; 41% (18/44) of patients were accepted, establishing a reference acceptance rate the ABM must exceed. Validation: On the holdout test set (n=7), simulations using the accepted population parameters (12 sets; five repeats each) were compared to observed on-treatment tiles via SAM/VarSAM. Five independent control groups of randomly selected parameter sets were also tested. Model accuracy was defined as the percentage of simulations accepted (SAM >0.7 and VarSAM >0.3). Applications: Timecourse predictions for CD8 dynamics were generated to inform biopsy scheduling (e.g., identifying time of peak or plateau). Two hypothetical combination scenarios were simulated: scenario 2 doubled CD8 proliferation probability; scenario 3 increased CD8 influx rate eightfold. Predicted spatial metrics and simulated biopsy images can be extracted at arbitrary times to guide decisions.

- Dataset: 44 usable paired biopsies (71 enrolled), split into 37 training and 7 test patients; tiles per patient ranged from 4 to 65; tiles sized 100×100 cell widths; ML-based cell detection achieved ~90% accuracy. - Sensitivity analysis: Parameters most influencing outcomes were immune-related: CD8 proliferation probability (IMpprol), killing probability (IMpkill), natural death (IMpdeath), randomness of migration (IMrwalk), and immune influx rate. Tumor cell proliferation, migration, and death parameters had weaker effects on CD8 and tumor counts in the analyzed region of parameter space. - Optimization: 12 parameter sets met acceptance thresholds (SAM >0.7; VarSAM >0.3) across patients, representing 1.2% of sampled parameter space; accepted sets tended to have low IMpdeath. - Null benchmark (LOCF): Using baseline spatial distributions to represent on-treatment yielded acceptance in 41% (18/44) of patients. - Validation: On the holdout test set, the ABM with population parameters reproduced spatial features of on-treatment biopsies with mean accuracy 77% and median 100%, versus 35% mean accuracy for randomly selected control parameter sets and the 41% LOCF benchmark. - Timing insight: Model predicts that, at a population level, biopsies taken around day ~30 capture information near maximum CD8 infiltration; however, patient-specific dynamics vary considerably. - Combination simulations: Increasing CD8 proliferation (2×) led to an initial rise followed by a later decrease in CD8 numbers; increasing CD8 influx (8×) produced a more sustained CD8 increase. These emergent differences illustrate the utility of mechanistic modeling for hypothesizing MoA and combination effects. - Personalization: Spatial features and optimal biopsy timing depend strongly on baseline spatial distributions, supporting individualized biopsy scheduling based on patient-specific baseline tiles. - The model demonstrates that a single baseline feature (CD8 spatial distribution) can predict on-treatment spatial distributions with high accuracy (mean 77%).

The integrated digital pathology–mechanistic modeling approach quantitatively addresses a key clinical trial challenge: selecting informative on-treatment biopsy times to evaluate drug mode of action. By validating against paired baseline and on-treatment biopsies and introducing SAM/VarSAM to quantitatively compare spatial cell distributions, the model shows substantial predictive value beyond a LOCF baseline and random parameter controls. The findings suggest that around 30 days may be optimal, on average, to assess maximal CD8 infiltration, but individual variability supports a personalized scheduling strategy. Sensitivity analyses highlight that immune dynamics and migration stochasticity are critical drivers of tissue-level outcomes, aligning with and extending biological intuition. The combination simulations demonstrate how distinct mechanisms (proliferation vs influx) can produce divergent temporal trajectories and spatial features, enabling hypothesis generation for combination therapy design. Overall, the approach can guide trial protocol design, inform combination selection, and complement machine learning by generating interpretable, full timecourses of spatial biomarkers.

This work contributes a quantitatively validated, mechanistic framework that integrates digital pathology with an agent-based tumor–immune model to predict spatiotemporal CD8+ T-cell dynamics from baseline biopsies. The novel SAM/VarSAM metrics enable robust comparison of spatial statistics between simulated and observed biopsies. Trained on real clinical data, the model achieved 77% mean accuracy in reproducing on-treatment spatial features and provided actionable insights for biopsy timing and hypothetical combination effects. Future directions include incorporating pharmacokinetics/pharmacodynamics and intratumoral drug concentration effects on model parameters, expanding datasets to diverse indications, improving computational efficiency to explore broader parameter spaces, and integrating model-generated timecourses with machine learning for hybrid predictive pipelines. Such tools could become routine decision support in drug development and personalized medicine.

- Clinical/data limitations: Small sample sizes typical of early-phase trials; small, fragmented needle-biopsy tissue with some unusable datasets. Despite this, prior evidence suggests small sections can represent intra-tumor heterogeneity. - Modeling/statistical limitations: SAM is less reliable with very low CD8 counts, biasing optimization toward parameter sets that increase CD8 numbers. Computational constraints limited exploration to 1000 parameter sets with five stochastic repeats, potentially missing parameter co-dependencies. Parameter sensitivity was evaluated independently in a localized region of high-dimensional space. One test patient had few on-treatment tiles, reducing evaluability. Drug concentration effects (PK/PD linking to parameter modulation) were not included due to lack of data but planned for future work.