Main

Physiological processes involve cell populations that change over time1. Some cell populations expand and others are removed, as occurs in development and in the immune response to pathogens9. A clinically important example of changing cell populations is the evolution of the cancer microenvironment, in which the growing tumour recruits stromal and immune cells that are crucial for tumour survival10,11. Understanding cell population dynamics and the underlying cell–cell communication circuits is a major goal of tissue biology, and can enable new treatment strategies based on sculpting cell populations in desired ways7,12.

However, measuring the dynamics of cell populations in human tissues is currently very difficult. Biopsies provide a single snapshot, and taking multiple biopsies from a single patient is infeasible and does not provide longitudinal evidence from the same cells.

Current strategies to measure tissue dynamics do not apply to human biopsies. For example, cell lines, mice, organ-on-a-chip13 or organoid models14 are treated as replicas and are analysed at different time points. Ex vivo tissues can be followed over time but lack the native physiological context. Intravital fluorescence microscopy has made it possible to follow living cells within animal models15,16. Recent advances use synthetic biology to engineer cells to record their activity or lineage17,18. These approaches are restricted to animal models or in vitro settings, limiting their ability to capture the complexity of in vivo dynamics in humans.

The emergence of single-cell technologies offers new opportunities for studying tissues at high resolution. Notable approaches use single-cell data for understanding cell–cell communication at a single time point19,20,21,22. Other approaches attempt to infer dynamics of processes such as transcription within individual cells on a timescale of hours. Examples include RNA velocity, ergodic rate analysis and Zman-seq23,24,25. These methods do not address the challenge of understanding how cell populations change on the tissue level, processes that could take days to weeks.

Here we present an approach to estimate cell population dynamics based on a spatial biopsy snapshot (Fig. 1a). This approach, OSDR, is based on using a cell division marker to determine division rates as a function of neighbourhood composition, providing a dynamical model of how cell populations change over time. We applied OSDR to human breast cancer spatial proteomics samples from three large cohorts2,3,4. OSDR reconstructs a fibroblast–macrophage circuit with two steady states of hot and cold fibrosis5,7, in agreement with co-culture experiments that measure dynamics directly8. We then used OSDR to discover a pulse-generating excitable circuit of T and B cells, suggesting that cancer surveillance in the tumour microenvironment operates in temporal flares, as opposed to the steady-state picture implicit in current literature. Finally, we validated OSDR using longitudinal data from patients with triple-negative breast cancer who received either chemotherapy or chemotherapy and immunotherapy. In both treatment regimes, OSDR predicted the collapse of the tumour cell population in responders but not in non-responders, based on biopsies at treatment initiation. The present approach opens the way to infer cell population dynamics from spatial snapshots of patient biopsies.

Fig. 1: Overview of OSDR.

a, OSDR uses spatial proteomics of a biopsy to infer cell population dynamics. The illustration was created by Nigel Orme. b, The division probability is learned based on the number of cells of each type in the neighbourhood. c, Cell division rates are inferred from a cell division marker (Ki67). d, Rate of change in the number of cells (X) is the difference between divisions and removals as described by a stochastic model. e, Simulations can take an initial spatial arrangement of cells and propagate the populations over time. f, OSDR also provides a phase portrait showing the direction of change of cell populations as a function of their neighbourhoods, revealing fixed points.

Tissue dynamics from a spatial snapshot

The aim of this study was to introduce a method for inferring a dynamical model of cell populations based on a single tissue sample (Fig. 1a). We considered biopsy sections with cell-type markers and cell division markers. From this data, we generated a list of cell coordinates and division status for each cell type. We used imaging mass cytometry (IMC) samples2,3,4 and the division marker Ki67 (Methods), but in principle, other spatial proteomic or transcriptomic assays can be used (Methods).

The key idea is to use a spatial omics snapshot to estimate the division and removal rates of cells as a function of the composition of the neighbourhood of the cell. We call this approach OSDR. Its essential advance is the use of division rates measured at the level of single cells at one time point to produce dynamics at the tissue level.

The change in cell counts within a tissue is a balance between division, death and migration. We separated the dynamics into two parts: one part that results from cell division and death within a tissue, and a second part that results from influx of cells from circulation. Our work focused on inferring cell division and death. We then used this estimated component to quantify the net contribution of migration from external sources of cells (Methods; Supplementary Figs. 3m,n and 4m,n). In our analyses, local proliferation and removal were sufficient to explain the cell population dynamics, but this might not be the case in settings with a massive influx of cells.

To estimate the effect of each cell type on the growth and removal of other cells, we assumed that the neighbourhood surrounding each cell contains growth factors secreted by nearby cells. We thus tabulated, for each cell in a slide, its type, as well as the number and types of cells within a radius (r) (Fig. 1b). We used r = 80 µm based on measurements of in vivo cell–cell interaction ranges26 (Methods). The ability of neighbourhood composition to predict division rates, as well as the inferred dynamics, are not sensitive to this choice of radius (r) (Supplementary Figs. 2e, 3f and 4g).

We defined division events using Ki67 thresholds based on experiments in human cell lines by Uxa et al.27 (Methods; Fig. 1c). By using data on many different neighbourhoods, available from the spatial heterogeneity of the samples, we obtained the probability of division as a function of neighbourhood cell composition (Methods; Fig. 1d).

This approach is precise when cells across the sample obey the same underlying dynamics. In reality, varying proximity to blood vessels or tumour mass can create zones with different levels of metabolites, hypoxia, inflammation and other factors that can affect the dynamics29. To control for these effects, we show below that the estimated dynamics are preserved across such zones by including terms for density of endothelial or tumour cells in our model (Supplementary Figs. 3h and 4h). We also show that the estimated dynamics are preserved across various patient subgroups defined by external factors such as tumour genetics or stage (Supplementary Figs. 3i and 4i,j).

Ideally, one would also like a marker for cell death, but currently available markers are not considered to be sufficient to quantify the diverse forms of cell death28. To make progress, we bypassed the death marker by assuming that the death rate is a constant for each cell type, namely, that the removal rate is not affected by the composition of the neighbourhood of the cell. We found that the results are robust to wide variations in the value of this constant removal rate (Supplementary Figs. 3g and 4e). We therefore approximate the removal rate by the mean division rate for each cell type, an assumption that is exact in the case of (quasi) steady-state tissues in which mean removal and division rates are equal.

We used the estimated models for cell division and death to produce dynamics of cell populations at the scale of a tissue. To produce a trajectory of tissue composition over time, we performed the following steps. We computed the probabilities for division and death for each cell in the tissue. We then used the computed probabilities to sample division and death events for each cell. We added a new cell next to each dividing cell and removed cells that died. This produces the composition of the tissue after one time step ∆T. By repeating this process, we obtained a trajectory of tissue composition over time (Fig. 1e).

In addition to these detailed trajectories, our approach allowed a second perspective to study the circuitry of cell interactions in a manner that does not depend on the spatial organization of a specific tissue. To do so, we analysed the dynamical system of neighbourhood composition dynamics (Fig. 1f). Each neighbourhood composition was mapped to a location in a state-space where each axis corresponds to the number of cells of one cell type (Fig. 2a). Plotting the direction of change for each possible neighbourhood composition produced a phase portrait (Fig. 2b,c), in which arrows mark the direction and magnitude of the rate of change. The phase portrait provides an overview of the dynamics, including its stable and unstable fixed points. Figure 2a–c denotes two cell types for ease of visualization, yet this approach can provide phase portraits with numerous cell types as shown below.

Fig. 2: OSDR recovers known dynamical systems from simulations.

a, Neighbourhood composition maps to a location in state-space. b, The phase portrait shows the change in neighbourhood composition. c, The phase portrait corresponds to a set of ordinary differential equations, which are derived from the stochastic model using the mean direction of change. d, OSDR accurately reconstructs neighbourhood dynamics from simulations of known dynamical circuits based on simulated data of 10,000 cells of each type. ‘Known’ panels (left) display the ground-truth phase portrait of the system used to generate the simulated spatial data. The arrows indicate directions of change, and the coloured lines are nullclines in which only one cell type changes. The black and white dots are stable and unstable fixed points, respectively. ‘Inferred’ panels (right) display the estimated OSDR fixed points from ten simulations, for the four two-cell phase-portrait topologies, as well as the inferred phase portrait from the first simulation of the ten.

We thus obtained two outcomes: stochastic simulation starting from the tissue sample initial condition and phase-portrait analysis, which provides a general view of the dynamics. The phase portrait describes the cell-circuit ‘rules’ that govern local neighbourhoods. The simulations propagate these rules starting from a given initial condition. The two views complement each other. The phase portrait allows an overview of the fixed points of the system and basins of attraction; it does not include the effects of initial spatial conditions, which can influence the dynamics (Methods; Supplementary Fig. 1h–m). By contrast, the simulations provide a detailed account of a given tissue sample initial condition, but if the sample is small or not representative of other parts of the tissue, simulation might miss the larger picture provided by the phase portrait.

To test the feasibility of the OSDR approach, we began with simulated data. We asked how many cells are required to reconstruct a known dynamical system. In each simulation, we specified a certain dynamical system in which cells affect the growth rate of each other. We then simulated experimental spatial data by running the dynamics from various initial conditions of cell density, for a period of time that results in distributions of spatial cell concentrations that resemble experimental data. This created various cell compositions. We then fit OSDR to the simulated data and compared the estimated phase portrait with the phase portrait of the known dynamical system used in each simulation.

We simulated all possible phase-portrait topologies of two cell types, in which each cell type can be stable on its own or stable only in the presence of the other cell type. This resulted in four phase-portrait topologies (Fig. 2d). A few thousand cells of each type were sufficient to reliably reconstruct the fixed points and basins of attraction (Supplementary Fig. 2g–o). We conclude that the OSDR algorithm can reconstruct simple 2D dynamics from simulated cell data and recover the underlying phase portrait given enough cells.

OSDR infers fibroblast–macrophage dynamics

To test ODSR, we analysed an extensive spatial proteomics dataset by Danenberg et al.2 (IMC; Methods). Data included biopsies from 715 patients with breast cancer (Methods). The samples included various breast cancer genotypes and stages. Each sample was an approximately 500 µm × 500 µm tissue section. The samples included a total of 859,710 cells. Cell types were identified using standard markers and include fibroblasts, macrophages, endothelial cells, adaptive immune cells and epithelial cells. We adopted the cell-type definitions from the original publication2 (Supplementary Fig. 2a). The data included the Ki67 division marker.

As OSDR infers division probability (Ki67 > threshold) based on the identity of neighbouring cells, we first asked whether variation in Ki67 is indeed explained by the composition of the neighbourhoods of the cells in the data. This is not obvious, because this variation may stem from unmeasured factors such as gradients of nutrients, hypoxia or inflammation29. We performed logistic regression with cell counts as features (number of neighbouring cells of each type) and Ki67 as target (Ki67 over threshold; Methods). We found that cell counts accurately capture a wide range of division rates (Fig. 3a and Supplementary Fig. 2b). All fits were significant (log-likelihood ratio test P < 10−13; Supplementary Fig. 2c,d). We conclude that the counts of neighbouring cells are strong predictors of division rates.

Fig. 3: OSDR reconstructs breast cancer fibroblast and macrophage dynamics in agreement with in vitro experiments.

a, Statistical inference using IMC data from 715 human breast cancer biopsies2 captures a wide range of division probabilities based on cell counts in the neighbourhood. Each dot corresponds to a subset of 4,000 cells (859,710 cells in total). b, The breast cancer microenvironment includes interactions between cancer-associated fibroblasts and macrophages. The illustration was created by Nigel Orme. c, The in vitro co-culture experiments of Mayer et al.8 followed fibroblast and macrophage cells in breast cancer conditioned medium (red arrows) or standard medium (grey arrows) to establish a phase portrait with several fixed points (coloured circles).The arrows are changes in cell concentration over 4 days. Panel c was reproduced from ref. 8, Springer Nature Ltd, under a Creative Commons licence CC BY 4.0. d, Ki67 marker in fibroblasts as a function of neighbourhood composition (Methods). Each point corresponds to one fibroblast (69,873 cells), and its position is the composition of the neighbourhood of that cell. The blue line separates inferred regions of rising and falling cell numbers (inferred nullcline). The black contours represent density. e, Same as d but for macrophages (3,761 cells). f, OSDR inferred phase portrait for breast cancer-associated fibroblasts and macrophages. Black, black–white and white circles are stable, semi-stable and unstable fixed points, respectively. g, Phase portraits of fibroblasts and macrophages in an OSDR model that includes tumour cells. The panels show sections of the 3D phase portrait at low (left) or high (right) tumour density (tumour cells fixed at zero or at their mean density of 64 cells per neighbourhood, respectively). Circles are as in g. h, Kaplan–Meier curves for patients whose biopsies indicate hot or cold fibrosis (defined by a macrophage density of more than the halfway point between hot and cold fibrosis fixed points). The hot-fibrosis state is associated with poor prognosis (log-rank test P = 0.0046, n = 607 patients with associated clinical data). The shading indicates 95% confidence intervals based on Greenwood’s exponential formula38.

The predictive ability of neighbouring cell counts could be explained by a number of factors (Supplementary Fig. 2f). Growth factor concentrations are influenced by the number of secreting and consuming cells in a neighbourhood26,30. Neighbour counts also reflect the chance of a cell to directly contact another cell, capturing processes such as contact inhibition and signalling via direct contact. Some cell types influence tissue hypoxia (endothelial cells and cancer cells) and inflammation (for example, macrophages and lymphocytes), processes that affect cell proliferation. Counts of such cell types might reflect the extent of hypoxia or inflammation in a region.

In some systems (for example, rapidly proliferating immune or developmental contexts), the assumption might not hold that the proliferation dynamics only depend on the current observed neighbourhood composition, and not on the history of neighbourhood compositions. We therefore recommend testing whether the division rate of each cell is predicted accurately based only on the current number of cells in its neighbourhood (as in Fig. 3a and Supplementary 2b–d). If this is not the case, the current approach is not relevant.

We compared OSDR to recent work that established in vitro dynamics for fibroblasts and macrophages in a breast cancer medium8, supplying a reference phase portrait for the dynamics of these two cell populations. The in vitro phase portrait was obtained by seeding mice mammary fibroblasts and bone marrow-derived macrophages at different initial concentrations in co-cultures and following the changes in the cell populations over several days (Fig. 3b). This is a measurement of the population-level dynamics. In the presence of breast cancer-conditioned medium, fibroblast and macrophage cells showed a phase portrait with several distinct fixed points (Fig. 3c). Fibroblasts and macrophages supported each other in a fixed point called ‘hot fibrosis’, in the sense that both fibroblasts and macrophages coexist. Fibroblasts alone could support themselves in a ‘cold fibrosis’ fixed point. Macrophages were induced by the cancer-conditioned medium to form a macrophage-only fixed point at high macrophage densities.

The OSDR phase portraits from breast cancer samples provide dynamics for fibroblasts and macrophages (Fig. 3f). We included all fibroblast and macrophage subsets defined by Danenberg et al.2 (Supplementary Fig. 2a). Similar results are found when considering neighbourhoods without adaptive immune cells, which were not present in the in vitro experiments, or including all neighbourhoods and adding terms for adaptive immune cells in the model (Supplementary Fig. 4k).

Of note, the reconstructed phase portrait shares most of the key features of the in vitro portrait (Fig. 3c,f). There is a stable fixed point in which fibroblasts and macrophages support each other, namely, a hot fibrosis fixed point. There is a cold fibrosis fixed point with fibroblasts only. There is a third semi-stable fixed point with high numbers of macrophages without fibroblasts.

We tested the robustness of the OSDR phase portrait. The fixed-point structure of the phase portrait is robust to resampling of cells and samples (Supplementary Fig. 3b–d). It is also consistently found in a wide range of parameter choices (neighbourhood radius and Ki67 threshold; Supplementary Fig. 3e,f) and death rates (Supplementary Fig. 3g), inclusion of endothelial or tumour cells in the statistical model (Supplementary Fig. 3h), or when considering only patient subgroups with specific cancer stage, tumour size, breast cancer genotype, patient survival time or type of treatment (Supplementary Fig. 3i).

We conclude that the reconstructed phase portrait captures the salient features of the in vitro phase portrait. Whereas the in vitro phase portrait used dynamical measurements of cell populations over days, the reconstructed portrait uses a single snapshot with cell density and division information.

We then asked whether cancer cells control the transition between cold and hot fibrosis. We used OSDR to fit a 3D model of cancer cell–fibroblast–macrophage dynamics. We then plotted 2D fibroblast–macrophage phase portraits at increasing cancer cell densities (Fig. 3g). Increasing tumour cell density pushes the stable fixed point towards a ‘hot’ state with abundant macrophages (Fig. 3g).

To ask whether hot fibrosis is clinically important, we partitioned the cohort from Danenberg et al.2 into two groups according to the presence of hot or cold fibrosis (macrophage density > halfway point between hot and cold fibrosis fixed points), and fit a Kaplan–Meier survival curve to each group. The hot fibrosis state was associated with poor prognosis (Fig. 3h; log-rank test P = 0.0046), with a decrease in median survival from 192 to 132 months. This result remained significant when the number of tumour cells was controlled for (Wald test P = 0.01; Supplementary Fig. 3j).

To further test OSDR, we analysed two additional IMC breast cancer datasets, with an additional 1,012 patients and 2.1 million cells3,4. Similar results were found in all three IMC breast cancer datasets (Supplementary Fig. 3k,l).

Excitable dynamics of T and B cells

We next considered two other cell types in the breast cancer microenvironment, T cells and B cells. These cells are components of the adaptive immune system that have a major role in the tumour microenvironment and are crucial for immunotherapy31,32.

We analysed all neighbourhoods in the Danenberg dataset with at least one T and/or B cell. This included a total of 73,961 T cells and 27,642 B cells, with 1,067 and 213 cell division events, respectively. Thus, about 1.4% of the T cells and 0.7% of the B cells showed division events, consistent with previous breast cancer data33.

We estimated the dynamics between T and B cells using OSDR. The inferred phase portrait has a striking feature called excitability. It shows a stable fixed point at zero cells (Fig. 4c). However, if T cells are raised above a threshold, they generate a pronounced pulse in which T cells rise, followed by B cells, and then both populations reduce to zero (Fig. 4d). This pulse is similar to pulses observed in autoimmune diseases such as relapsing–remitting multiple sclerosis34.

Fig. 4: OSDR infers a pulse-generating excitable circuit of T and B cells in the breast cancer microenvironment.

a, Ki67 marker as a function of neighbourhood composition. Each point corresponds to one T cell (73,961 cells) from the IMC dataset of Danenberg et al.2. b, Same as a but for B cells (27,642 cells). c, OSDR inferred phase portrait indicates an excitable system. The blue highlighted trajectory displays a large pulse of adaptive immune cells. d, Spatial simulation of a 2 mm × 2 mm tissue based on inferred population dynamics shows that a high enough number of T cells triggers a pulse of T cells followed by B cells. The y axis displays the number of cells, normalized to the area of a 80-µm neighbourhood. e, Each pulse is followed by a refractory period, as evidenced in simulations in which additional T cells are introduced at different times (vertical arrows). A new pulse is generated only when T cells are introduced after B cells from the previous pulse have declined. f, 3D phase portrait with CD4 T cells, CD8 T cells and B cells shows an immune flare when the CD4 T cell density crosses threshold. The 2D projections are in grey. g, Distributions of CD4 T cells, CD8 T cells and B cells in the 100 most-proliferative and least-proliferative states of each cell type. The most-proliferative states are characterized by a high density of CD4 T cells and a low density of B cells.

Similar to the fibroblast–macrophage analysis above, the estimated dynamics were robust to resampling at the level of both cells and patients (Supplementary Fig. 4b–d) and were consistent under a wide range of death rates, neighbourhood sizes and Ki67 thresholds (Supplementary Fig. 4e–g). The dynamics are not confounded by unmodelled cell types (Supplementary Fig. 4h).

The pulsatile adaptive immune response in the inferred model is shown in Fig. 4c, in which T cells rise, followed by B cells that inhibit them, and then both cell populations decline. The model further indicates that after a pulse, there is a refractory period, in which B cells are still high and T cells are low. During this period, a new pulse cannot be generated due to the inhibitory effects of the B cells (Fig. 4e). A new pulse is possible only after the recovery period (Fig. 4e).

To study the role of different T cell subsets, we fit a 3D model including CD4 and CD8 T cell subsets as well as B cells. The excitable dynamics apparent in the 2D model were reproduced in the 3D model. The 3D model revealed that a pulse is initiated when the density of CD4 T cells crosses a threshold (Fig. 4f). For all three cell types, proliferation was highest in neighbourhoods with a high density of CD4 T cells and a low density of B cells (Fig. 4g). Proliferation was lowest at a high density of B cells and low density of CD4 cells (Fig. 4g). We conclude that CD4 cells initiate the immune pulse, whereas B cells provide the negative feedback required to terminate the pulse.

We next considered subtypes of breast cancer from the Danenberg cohort. We compared ER+, HER2+, PR+ and triple-negative breast cancer. Fibroblast–macrophage dynamics displayed hot and cold fibrosis fixed points in the phase portrait of all four subtypes (Supplementary Fig. 4j). The T–B immune pulse was similar in the three receptor-positive subtypes: ER+, PR+ and HER2+. However, the triple-negative breast cancer phase portrait showed a difference where the pulse does not end in zero T and B cells, but instead in a fixed point where B cells remain stable (Supplementary Fig. 4j). This result is consistent across breast cancer cohorts, with mixed cohorts displaying a collapse of both cells, and the triple-negative cohort displaying B cell stability. This result is consistent with the high rates of lymphocyte infiltration in triple-negative breast cancer, including tumour-infiltrating B cells at early stages of the disease35,36.

We also asked whether a full OSDR model with all cell types displays the same dynamics as the 2D models. We thus developed a model with the six major cell types in the Danenberg dataset: fibroblasts, macrophages, T cells, B cells, cancer cells and endothelial cells. We found that the 6D model recapitulates the dynamics of the 2D models in which other cell types are at their mean concentrations (Supplementary Fig. 4k). We conclude that the OSDR approach allows reconstruction of the dynamics of the full cell–cell interaction network in the samples.

OSDR predicts response to treatment

Patients with cancer can go through many months of treatment without knowing whether the tumour is growing or shrinking. To adjust treatment when it is ineffective, we need early signs of response. To address this challenge, we applied OSDR to a clinical trial with three longitudinal biopsies: the NeoTRIP trial. A cohort of 279 patients with triple-negative breast cancer were randomly assigned to chemotherapy (n = 141) or chemotherapy + anti-PDL1 immunotherapy (n = 138; Fig. 5a). Treatment consisted of eight 3-week cycles. Biopsies were collected at three time points: before treatment, 3 weeks into treatment (day of the second treatment cycle) and post-treatment (at surgical excision of the tumour following eight 3-week treatment cycles). Pathologists labelled the post-treatment biopsies as pathological complete response or residual disease.

Fig. 5: OSDR predicts the collapse of the tumour cell population in responders based on early-treatment biopsies.

a, Longitudinal triple-negative breast cancer dataset by Wang et al.3 includes a total of 141 patients in the chemotherapy arm and 138 patients in the chemotherapy + immunotherapy arm. The illustration was created by Nigel Orme. b, Longitudinal predictions for six cell types. OSDR models were fit separately to early-treatment biopsies from responders (47 and 49 patients in the chemotherapy and chemotherapy + immunotherapy arms, respectively) and non-responders (54 and 57 patients in the chemotherapy and chemotherapy + immunotherapy arms, respectively) in each treatment arm (a total of 140,118 and 157,347 cells in responders and non-responders, respectively, in the chemotherapy arm; 135,559 and 144,261 cells in the chemotherapy + immunotherapy arm, respectively). The tumour cell population collapses in responders (dark grey line) but is stable in non-responders (red line). To isolate the effect of the tumour microenvironment dynamics, as opposed to initial tissue composition, trajectories were computed starting from the composition of all biopsies taken at the beginning of treatment. Plots display the average cell count over all trajectories. The shading denotes 95% confidence intervals of the mean (bootstrap)39. c, Fraction of dividing tumour cells in all neighbourhoods (n = 41,117 and 17,932 cancer cells from non-responders and responders, respectively) and in neighbourhoods with the top 10% of T cells. The error bars display 95% confidence intervals of the mean (bootstrap)39. We compared the fraction of Ki67+ cancer cells in responders versus non-responders using chi-squared tests of independence in all neighbourhoods (left) and in high T cell neighbourhoods (right). Both cases had P < 10−7. d, Zoomed in view of the first 90 days of predicted chemotherapy + immunotherapy treatment response dynamics. The T cell density increases in both responders and non-responders, and the tumour cell population collapses only in responders. The shading denotes 95% confidence intervals of the mean (bootstrap)39.

We hypothesized that a response or lack of response to treatment should be evident in the early-treatment dynamics. We applied OSDR to week 3 biopsies from responders and non-responders in each treatment arm (total of four groups). We modelled dynamics of the six major cell types: fibroblasts, macrophages, tumour cells, endothelial cells, T cells and B cells. We then predicted the composition of the tissue over time. To detect differences resulting from dynamics, as opposed to differences resulting from initial tissue composition or tumour burden, we applied both responder and non-responder models to the same initial states, namely, we applied both models to every early-treatment biopsy. Figure 5b shows the predicted cell density over time, averaged over all starting states. In both treatment arms, OSDR predicts the collapse of the tumour cell population in responders but not in non-responders. These results are robust to the choice of patients used to fit each model (Supplementary Fig. 5a; Mann–Whitney U-test P < 10−5 for both treatment arms).

Our analysis showed that the collapse of tumour cells in responders is not a result of the initial tissue state, as we applied the models to all initial states. Moreover, the average proliferation rate of tumour cells at week 3 does not separate between responders and non-responders (Supplementary Fig. 5b). In fact, the mean division rate of tumour cells in responders was higher than non-responders in the immunotherapy treatment arm (Fig. 5c; chi-squared P < 10−7). Thus, the tumour cell population collapses as a result of interactions between cells over time and not initial tissue composition or average proliferation rates.

The predicted trajectories in immunotherapy show a sharp rise in T cell number (Fig. 5b). This rise appears in both responders and non-responders receiving immunotherapy, but not in patients in the chemotherapy arm. Because T cell density increases in both responders and non-responders under immunotherapy, we hypothesized that T cells have a different effect on tumour cells in these two groups. When considering all neighbourhoods, the fraction of dividing cancer ce