High resolution microfluidic assay and probabilistic modeling reveal cooperation between T cells in tumor killing

Gustave Ronteix, Shreyansh Jain, Christelle Angely, Marine Cazaux, Roxana Khazen, Philippe Bousso, Charles Baroud

<figure><img src="image_1.png" /></figure>

► To cite this version:

Gustave Ronteix, Shreyansh Jain, Christelle Angely, Marine Cazaux, Roxana Khazen, et al.. High resolution microfluidic assay and probabilistic modeling reveal cooperation between T cells in tumor killing. Nature Communications, 2022, 13 (1), pp.3111. 10.1038/s41467-022-30575-2 . pasteur-03775301

HAL Id: pasteur-03775301

https://pasteur.hal.science/pasteur-03775301

Submitted on 12 Sep 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
ARTICLE

https://doi.org/10.1038/s41467-022-30575-2

OPEN

<figure><img src="image_3.png" /></figure>

Check for updates

# High resolution microfluidic assay and probabilistic modeling reveal cooperation between T cells in tumor killing

Gustave Ronteix<fnref n="1" />,²,⁴, Shreyansh Jain<fnref n="1" />,²,⁴, Christelle Angely<fnref n="1" />,², Marine Cazaux³, Roxana Khazen³, Philippe Bousso<fnref n="1" />,³✉ & Charles N. Baroud<fnref n="1" />,²✉

Cytotoxic T cells are important components of natural anti-tumor immunity and are harnessed in tumor immunotherapies. Immune responses to tumors and immune therapy outcomes largely vary among individuals, but very few studies examine the contribution of intrinsic behavior of the T cells to this heterogeneity. Here we show the development of a microfluidic-based in vitro method to track the outcome of antigen-specific T cell activity on many individual cancer spheroids simultaneously at high spatiotemporal resolution, which we call Multiscale Immuno-Oncology on-Chip System (MIOCS). By combining parallel measurements of T cell behaviors and tumor fates with probabilistic modeling, we establish that the first recruited T cells initiate a positive feedback loop to accelerate further recruitment to the spheroid. We also provide evidence that cooperation between T cells on the spheroid during the killing phase facilitates tumor destruction. Thus, we propose that both T cell accumulation and killing function rely on collective behaviors rather than simply reflecting the sum of individual T cell activities, and the possibility to track many replicates of immune cell-tumor interactions with the level of detail our system provides may contribute to our understanding of immune response heterogeneity.

<footnote n="1">¹Physical microfluidics and Bioengineering, Institut Pasteur, Université Paris Cité, 25-28 Rue du Dr. Roux, Paris, France. ²LadHyX, CNRS, École Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France. ³Dynamics of Immune Responses Unit, Equipe Labellisée Ligue Contre le Cancer, Institut Pasteur, Université Paris Cité, INSERM U1223, 75015 Paris, France. ⁴These authors contributed equally: Gustave Ronteix, Shreyansh Jain.
✉email: philippe.bousso@pasteur.fr; charles.baroud@pasteur.fr</footnote>
The capacity of cytotoxic T lymphocytes (CTL) to eliminate tumor cells is the basis for the development of important tumor immunotherapies such as immune checkpoint inhibitors (e.g., anti-CTLA-4, anti-PD1, or anti-PD-L1 mAbs) and the development of cellular therapies such as CAR T cells¹,². However, patient response to these therapies can be highly variable. While many parameters are known to influence patient response to immunotherapies, the number, phenotype, and distribution of CTLs can have a strong predictive value in several types of cancer³.

These observations underscore the need to better understand how a successful T cell attack proceeds and what are the critical parameters associated with CTL behavior and function that favor tumor regression. In this respect, several key questions remain unanswered. For example, how do CTLs encounter tumor cells and accumulate within the tumor microenvironment (TME)? What are the dynamics of CTL killing and what CTL density is needed for tumor eradication? Are individual CTLs acting autonomously in the TME or do they interact together? Understanding the basic principles that dictate whether a tumor mass is regressing or not is in fact essential to design, optimize, and evaluate tumor immunotherapeutic strategies.

Multiple approaches are available to evaluate T cell cytotoxicity against tumors. In vitro assays in cell suspension have been used to measure CTL killing capacity and, when performed at the single-cell level, provide information on the extent of functional heterogeneity within a T cell population⁴. These in vitro assays, however, lack the complexity of the 3D tumor microenvironment, which strongly impacts T cell behavior and function. At the other end of the spectrum, intravital imaging offers direct insights into the dynamics, signaling, and killing behavior of single T cells within a developing tumor⁵⁻⁸. Limitations of these approaches, however, include the fact that they provide a view of the interactions in a limited spatial and temporal window. Indeed continuous observation periods are generally limited to a few hours, precluding a full understanding of T cell histories in the TME.

An interesting emerging platform comes from advanced in vitro models that recapitulate some aspects of the TME while providing access to the system dynamics⁹,¹⁰. These include organoids¹¹,¹², where cells are allowed to organize in three dimensions (3D), or organ-on-a-chip devices¹³, where the microfluidic device represents the organ geometry and the microfluidics enable temporal control of the flows and physical conditions. Recent work has also dealt with combining the advantages of both approaches to produce organoids-on-a-chip¹⁴. Gathering general rules from these systems that would explain the outcome of a CTL attack remains challenging as it requires to link quantitative measurements of T cell behavior, which is inherently stochastic, with tumor cell fate⁹,¹⁴.

Here, we introduce a microfluidic-based approach for the multiplexed analysis of tumor spheroid fate in the presence of defined immune cell populations. This methodology is based on the parallel formation, manipulation, and observation of hundreds of tumor spheroids within stationary microfluidic droplets. When associated with mathematical models, the quantity and quality of spatiotemporally resolved data allow us to pinpoint key behaviors leading to spheroid destruction and to detect and understand heterogeneity of tumor outcomes.

## Results

**Parallelized immune challenge on an integrated microfluidic chip.** The immunogenic rejection of 3D cancer models is studied using the classic model of Ovalbumin-expressing mouse B16 melanoma (B16-OVA), which are challenged by OVA-specific CD8+ cytotoxic T lymphocytes (CTLs) bearing the OT-1

transgenic TCR (see “Methods”)¹⁵⁻¹⁹. Further, the results obtained with B16-OVA are compared with non-OVA expressing B16 Wild-Type cells (B16-WT) as a control. These two cell types were previously reported to exhibit similar proliferation dynamics and do not differ in their in vivo immunogenicity¹⁸. The experiments rely on a microfluidic device that consists of a droplet generating region followed by a droplet trapping region that serves to culture the cells and observe them (Fig. 1a, b)²⁰,²¹. This trapping region is patterned with 234 microfluidic anchors²², that allow the droplets to be held in place even in the presence of an external flow (see “Methods”). The anchors are diamond-shaped (see Fig. 1b) in order to allow for multiple droplet pairings²³.

The experiment begins by producing aqueous droplets (volume = 50 nl) containing Matrigel and a suspension of B16 cells at a concentration of 1.5 × 10⁶ cell/ml. To obtain single spheroids in each droplet, a concentration of 2.0 mg/ml of Matrigel was used, with higher Matrigel concentrations leading to the formation of multiple spheroids per droplet (Supplementary Fig. 1a) (see “Methods”). Once the droplets are anchored, the device is placed in an incubator at 37 °C overnight, which allows a single B16 spheroid to form in each droplet (Fig. 1b). At these cells and Matrigel concentrations we obtain spheroid radii in most cases ranging between 35 to 45 µm (Fig. 1c). A live-dead staining shows that less than 3% of the cells were dead after 48 h in the chip (Fig. 1d).

After overnight incubation, the CTLs are brought to the Matrigel droplets by generating a group of secondary smaller droplets (volume = 10 nL) that contain a broad distribution of CTLs (Supplementary Fig. 1b). These secondary droplets are trapped in the triangular sections of the anchors and then merged with the spheroid-containing Matrigel, thus bringing the two cell populations into the same Matrigel droplet (Fig. 1e)²³. The interactions between the CTLs and spheroids are observed by time-lapse microscopy, typically over 24 h. The whole process from spheroid preparation to CTL addition and imaging is done on a single chip (Fig. 1f) and each experiment yields up to 234 individual replicates, of which we typically obtain 50 time-lapse movies, due to the small image acquisition time-intervals (2min/ frame). For higher time-intervals, more data points can be collected from the same chip. The movies are then analyzed using our custom-made scripts (see Supplementary Software 1).

Three stages in the cell-cell interactions can be identified: the CTL exploration of the gel, their accumulation on the spheroid, and the killing of B16 Ova-expressing cells (see Supplementary Movie 1 for representative cases). These stages are studied in detail below.

## CTL migration in micro-device reproduces in vivo behavior.

The CTL migration is tracked in the time-lapse movies as the cells perform 3D migration within the gel or on the surface of the spheroid (see Fig. 2a). The recorded velocities display alternating periods of motility and arrest phases, as seen by the high and low velocities in Fig. 2b. This behavior, as well as the value of the velocities, correspond well to previously reported CTL velocities in collagen gels in vitro²⁴ or within tissues in vivo²⁵,²⁶. We can infer the 3D motility properties from the acquired microscopy data (see “Methods”). Note that the value of the measured velocities depends on the image acquisition frequency (see Supplementary Fig. 2), so we maintain the sampling at 2min/frame for all the experiments (see “Methods”).

In order to evaluate the influence of the spheroid presence in the droplet on the motility of the CTLs, the migration statistics without a spheroid present are compared with the statistics in the presence of a spheroid during the first 500 min of an experiment. We restrict the analysis to cells that are not in contact with the
<figure><img src="image_4.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Chamber top view</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_9.png" /></figure>

**Fig. 1** *Microfluidic immuno-oncology chip and protocol.* **a** Microfluidic chip on a standard glass slide. **b** Expanded view of the trapping region of the chip (dashed box) showing an array of 234 trapped droplets. Each droplet contains a single B16 spheroid in Matrigel, as shown in the inset. **c** Distribution of spheroid radii within a single chip (N = 215). **d** Viability measurements using live-dead staining after 24 and 48 h (N = 54). **e** Schematic showing a primary droplet with a tumor spheroid, followed by the addition and fusion of a secondary droplet containing GFP-labeled CTLs, eventually leading to tumor cell killing and spheroid fragmentation. Scale bar is 200 µm. **f** Schematic representation of the complete experimental protocol.
<figure><img src="image_10.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>i</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>g</figcaption></figure>

<caption>h</caption>

**Fig. 2 CTL migration in droplets recapitulates in vivo behavior.** **a** (left) Representative image of CTLs with instantaneous velocity vectors inside Matrigel droplet. (right) CTL tracks in one droplet over 24 h, each color represents an individual cell track. The dashed black circle outlines the spheroid boundary. **b** Representative velocities as a function of time for three different T cells. **c** Probability distribution of a cell to migrate by a given distance ($\Delta r$) during a fixed time step $\Delta t = 1$ min ($n = 67965$ points without spheroid and $n = 34072$ individual points for CTLs in presence of the B16 spheroids). **d** Mean-square displacement (MSD) of CTL migration with ($N = 20$ droplets) and without ($N = 26$ droplets) spheroids. Error bars represent the SEM. **e** Time sequence showing the initial CTL approach and contact with a spheroid. **f** Track of a single CTL as it migrates in the matrigel and on the spheroid surface. Colormap represents the instantaneous velocity of the cell. **g, h** Average velocity and mean square displacement exponent ($\alpha$) of cells migrating in the gel and on the spheroid. Each data point is the average velocity in a given droplet ($N_{gel} = 55$, $N_{spheroid} = 54$, respective $p$-values of $1.3 \times 10^{-10}$ and $1.2 \times 10^{-15}$). **i** Mean-square displacement of cells migrating in the matrigel and on the spheroid. Bold and dashed lines represent the best fits for the MSD of CTLs on the spheroid and in the matrigel, with respective exponents of 1.1 and 1.4 (measurement conducted over $N = 54$ droplets). Error bars represent the SEM. Source data are provided as a Source data file.

spheroid. The displacement distributions (Fig. 2c) and the mean-square displacements (MSD) (Fig. 2d) do not show any significant difference between the two conditions. In both cases, the CTLs undergo super-diffusive random walks (Fig. 2d) with MSD ~ $\tau^\alpha$, where $\tau$ is the time between two observations and $\alpha = 1.6$, in agreement with what was reported in vitro and in vivo$^{24,26}$.

After some time, one CTL comes in contact with the spheroid. This contact generally leads the T cell to adhere to the spheroid and explore its surface over the course of a few hours (Fig. 2e, f and Supplementary Movie 2). We select individual tracks with segments both on and off the spheroid to investigate more precisely the CTL motility change upon reaching the spheroid. We observe that the CTL behavior is strongly modified: They display lower mean velocity (Fig. 2g and Supplementary Movie 2)

with a lower MSD exponent (Fig. 2h, i). The average MSD exponent goes from 1.4 when the cells move in the gel to 1.1 after the same cells have reached the spheroid.

The CTL migration in the gel therefore recapitulates behaviors that have been reported in vivo$^{24-27}$, with the current data highlighting the switch in motility before and after the CTL contacts the spheroid surface.

**A positive feedback loop drives CTL accumulation on the spheroid.** We now investigate the contact time statistics of the CTLs on the spheroids. The distribution of first-passage times is consistent with the distribution of randomly migrating particles in an enclosed environment, indicating that the initial contact is random and that there is no attraction from the spheroid on the
<figure><img src="image_17.png" /><figcaption>**Fig. 3** CTL accumulation is enhanced by a positive feedback-loop after first contact. **a** Experimental distribution of first CTL-spheroid contact times and theoretical distribution for randomly migrating CTLs. **b** Number of CTLs detected on each spheroid and **c** the fraction per droplet as a function of time. Each thin line represents a single tracked spheroid, in bold is the averaged value. In red is the accumulation for B16-OVA spheroids and in blue for B16 WT spheroids [84 individual B16-OVA spheroids and 81 B16 WT spheroids tracked]. **d** Number of CTLs as a function of time on two representative spheroids showing the detection of attachment/detachment events. **e** Schematic of the stochastic accumulation model: CTLs can switch from the gel to the spheroid with different probabilities. $p_{in}(\Delta t)$ (conv. $p_{out}(\Delta t)$) is the probability for a cell to attach to (conv. detach from) the spheroid during a time interval $\Delta t$. Counting attachment and detachment events in the experiments allow us to infer the rates $\lambda_{in}$ and $\lambda_{out}$. **f** Estimates for the attachment rates ($\lambda_{in}$), detachment rates ($\lambda_{out}$), and affinity ratio ($\lambda_{in}/\lambda_{out}$) for B16-WT (blue) and B16-OVA (red) cells. The box plots are obtained by using a bootstrapping method with 50 repetitions as described in the methods (respective p-values of $7 \times 10^{-18}$, $1.2 \times 10^{-16}$, and $7 \times 10^{-18}$). **g** The affinity ratio as a function of the number of CTLs detected on the spheroid for B16 WT and B16 Ova spheroids (respective p-values of $6 \times 10^{-9}$, $3 \times 10^{-15}$, and $1.3 \times 10^{-13}$). **h** Normalized attachment rate $\lambda_{in}$ (white) and detachment rate $\lambda_{out}$ (red) as a function of the number of CTLs attached to the spheroid. $\lambda_{in}$ and $\lambda_{out}$ are normalized by their mean values for 0 and 1 CTL on the spheroid, respectively. Source data are provided as a Source data file.</figcaption></figure>

CTLs (Fig. 3a). After the initial CTL contact with the B16 cells, the arrival of successive CTLs leads to an accumulation of T cells on the spheroid in the case of the B16-OVA spheroids. This accumulation is shown in absolute numbers (Fig. 3b) and also by computing the ratio of CTLs in each droplet that are present on the spheroid as a function of time (Fig. 3c). However, the accumulation is not observed in the case of wild-type B16 cells, which do not express Ova (Supplementary Movie 3).

At this stage, an important question is whether this accumulation results from cells reaching the spheroid randomly, as they explore the droplet volume, or if the accumulation rate is enhanced due to cell-cell signaling. We address this question by analyzing the accumulation of the CTLs at the spheroid level. At this scale, the CTL accumulation is not homogeneous but stochastic, with CTLs both attaching and detaching over time (Fig. 3d, Supplementary Movie 2, Supplementary Movie 4). A

time series of the number of CTLs present on each spheroid is obtained by counting the CTLs that are detected within the spheroid region on the images, with a 2 min time resolution. The statistics of these time series on all of the parallel realizations are then analyzed using a Markov chain method (see “Methods”). This in turn yields the attachment rate $\lambda_{in}$ and detachment rate $\lambda_{out}$ for each individual CTL (Fig. 3e).

The value of $\lambda_{in}$ is found to be significantly higher in the B16-OVA spheroids compared the B16-WT spheroids, while the opposite is true for $\lambda_{out}$ (Fig. 3f). These measurements indicate that the arrival rate of CTLs increases when the cells composing the spheroid express the cognate antigen recognized by the CTLs and that they stay attached for longer periods of time. Therefore the accumulation of CTLs on the spheroids is mediated by two independent phenomena: first the increase in arrival frequency and second by the decreased leaving frequency. The net effect of
the attachment/detachment dynamics can then be summarized by the affinity ratio, $\lambda_{in}/\lambda_{out}$, which accounts for the net effective accumulation of CTLs on target. This ratio is found to be significantly higher with B16-OVA spheroids when compared to the B16-WT spheroids (Fig. 3f).

Evidence of a positive feedback loop for the attraction among the CTLs can be obtained by calculating the change of the affinity ratio, defined as $\lambda_{in}/\lambda_{out}$, as a function of the number of T cells present on the spheroid. Indeed, the depth of the experimental data allows us to obtain a value of $\lambda_{in}/\lambda_{out}$ before the first contact and then successively track the change of this ratio after every contact, in each droplet (Fig. 3g). Again, the data for B16-OVA show a significant difference with the WT case. More interestingly, the higher the number of CTLs present on the B16-OVA spheroid, the higher the ratio, and thereby, the faster the accumulation rate: this is a hallmark of a positive feedback. The increase in the accumulation rate is driven by an increase in $\lambda_{in}$ for a constant value of $\lambda_{out}$ (Fig. 3h). It demonstrates that the accelerated accumulation is mediated by the increasing attraction of the CTLs to the spheroid, which starts from the very first CTL attached to the spheroid (Fig. 3h).

The current results confirm recently published observations¹⁵ that show a CCR5-mediated swarming of T cells in vitro, as well as in vivo studies that report the accumulation of T cells on targets²⁸,²⁹. Figure 3 shows that these effects occur even for cell populations consisting of a few individuals and that a single contact can trigger the beginning of the positive feedback. Below we go beyond the CTL accumulation to address the relationship between the accumulation of CTLs and their capacity to kill the B16 spheroids.

**Killing of B16 cells by CTLs is heterogeneous.** After focusing on the behavior of the T cells we now turn to the response of the cancer spheroids upon CTL accumulation. The time-lapse movies allow us to identify individual cell death events in the spheroids at the molecular level by detecting the activation of Caspase 3/7, which provides an early marker of apoptosis³⁰ (Fig. 4a and Supplementary Movie 5). In the brightfield image, we observe instances of rapid shedding of cellular material and debris from the spheroids (Fig. 4a, Supplementary Movies 1, 3, 5), which we refer to as “spheroid fragmentation”. Combining this information with the position of CTLs relative to the spheroid, it is thus possible to record a detailed chronology of the key events taking place in each droplet by timing successive CTL contacts with the spheroid and the apparition of fragmentation events and Caspase 3/7 signals, as shown in Fig. 4b.

We observe that the timing of the first Caspase event post-CTL contact is well correlated with the first fragmentation event (Fig. 4c), indicating that the two observations are closely related. For this reason, we will hereafter use brightfield images to quantify spheroid killing, which simplifies the analysis pipeline. Furthermore, the timing of these killing events is highly variable, ranging from a few hours to beyond 24 h. For some spheroids, no fragmentation or Caspase events are observed over the course of an experiment. In the following analysis, we will label “successful killing” the cases when the first fragmentation event is observed before $t = 14$ h. The statistics of such events are summarized in Fig. 4d, which shows that 44% ($N_{OVA} = 84$ spheroids) of the OVA expressing spheroids are successfully killed by the CTLs. This contrasts with the B16-WT spheroids, where we do not observe any fragmentation events (Fig. 4d, Supplementary Movie 3). An analysis of the statistical impact of each of the problem parameters on a successful killing shows that the dominant parameters are the number of CTLs in the droplet and the number of CTLs that reach the spheroid.

**Tumor spheroid killing involves collective effects.** We now consider the relationship between the spatiotemporal dynamics of the CTLs and the tumor spheroid outcomes (successful or unsuccessful killing). An indication of the relevance of this link is obtained first by observing that the CTL accumulation rates are faster on the spheroids that display fragmentation than in the opposite case, both in absolute numbers (Supplementary Fig. 3a) and as a fraction of the total number of cells per droplet (Supplementary Fig. 3b). This indicates that a faster accumulation is correlated with efficient killing.

Moreover observing the CTLs on the spheroid reveals that fragmentation events are associated with the presence of several CTLs in the vicinity (Fig. 5a, b). This local effect is quantified by counting the number of T cells present within a 30 µm radius of the first cell fragmentation event, giving a mean number of 3.4 cells (median at 3 cells), with fragmentation very rarely occurring with only one cell present at the fragmentation site (Fig. 5b). Indeed, CTLs sometimes appear to besiege a salient B16 cell, causing it to burst after a few minutes of attack. This suggests that the CTLs tend to cluster at particular sites on the spheroid and that their clustering enhances their ability to induce spheroid fragmentation (Fig. 5c).

The above observations can be analyzed in greater depth to demonstrate that the cooperative effect of the CTLs enables the killing of tumor cells. Here, we first infer the fragmentation probability of cancer cells as a function of the number of CTLs present on the spheroid (n). If we consider a time-interval of length $\Delta t$, the probability of a cancer cell fragmentation event happening is $\Gamma_{frag}(n)\Delta t$ and conversely, the probability of the cell not fragmenting is $1 - \Gamma_{frag}(n)\Delta t$ (Fig. 5d) (see “Methods”). Comparing the different droplets we see that the fragmentation events can happen at different time-points as is schematically represented in Fig. 5e. From the experimentally observed fragmentation times $\Gamma_{frag}$ can be inferred: $\Gamma_{frag}(n) = \frac{K_{Death,n}}{K_{Death,n}+K_n}$. Here $K_{Death,n}$ is the number of fragmentation occurrences with exactly n CTLs on the spheroid and $K_n$ the number of instances with n CTLs on the spheroid before the first fragmentation event.

The experimental measurements are compared with a probabilistic model where the CTLs behave independently from each other as is schematically represented in Fig. 5f (see “Methods”). This independent model leads to the parabolic probability of killing (blue line in Fig. 5g) which poorly matches with the experimental data. On the other hand, the experimental evolution of $\Gamma_{frag}$ is well-fitted by an exponential increase (red dashed line), demonstrating that a model of independent interactions between the CTLs and the target cells fails to capture the underlying mechanism of spheroid fragmentation.

**Long and short-range interactions combine to determine probability of CTLs to kill a tumor spheroid.** The probability for successful killing to occur for any particular spheroid can now be explained as a combination of the effect of CTL accumulation on the spheroid and their cooperative killing behavior (Fig 6a). This probability is a function of the total number of CTLs in the droplet, which also correlates with the maximum number observed on the spheroid (see Supplementary Fig. 4 and Table 1). The spheroid fate can be simulated as a two-step branching process. First, we simulate the evolution of the number of CTLs on the spheroid at each time step based on the experimentally derived parameters $\lambda_{in}$ and $\lambda_{out}$, shown in Fig. 3. Second, we account for the possible spheroid fragmentation depending on the number of CTLs on the spheroid and $\Gamma_{frag}$, as described in Fig. 5. We then repeat the process for each time-interval, updating the number of T cells on the spheroid and the spheroid state (fragmented or intact) iteratively until the time of the simulated
<figure><img src="image_18.png" /><figcaption>**Fig. 4** CTLs kill tumor cells within spheroids. **a** Representative sequence showing CTL (OT1) positions on the spheroid, caspase 3/7 fluorescent death marker, and B16 cell fragmentation. The white arrow at 3h indicates the appearance of caspase signal next to a CTL. At 9 h it represents a fragmented dead cell. **b** Representative chronology showing the key events for a given spheroid interaction with CTLs: Contact times of CTLs on spheroids, detection of caspase 3/7 signal, detection of fragmentation events. **c** Time of first caspase signal vs. first observation of cell fragmentation. **d** Percentage of WT (black) and OVA (red) spheroids that show at least one fragmentation event in under 14 h. N equals 54 and 84 spheroids, respectively. The error bar is the SEM. Source data are provided as a Source data file.</figcaption></figure>

experiment has lapsed. The model process is illustrated in Fig. 6b (see “Methods”).

This computation is performed for a number of CTLs per droplet (N) ranging from 0 to 20 cells and repeated 50 times for each value of N. This allows us to first compare the experimental first fragmentation times with the simulated ones (Fig. 6c). In both cases, the first fragmentation time decreases as the number of CTLs increases in the droplet and the simulated values closely fit the experimentally observed fragmentation times. We then compare the simulated fragmentation probability curve as a function of the number of CTLs in the droplet to the experimental data (Fig. 6d). Not only does the simulation confirm the key role of CTL number in causing spheroid fragmentation, but the close match between simulated and experimental measurements indicates that the spheroid fragmentation process is well recapitulated from these two mechanisms: collective accumulation and cooperative killing of CTLs at the spheroid site.

## Discussion

The current study introduces a new paradigm for extracting biological information from in vitro experiments, by treating the parallel realizations as “*Monte-Carlo experiments*” that reach different outcomes in a probabilistic way. This contrasts with existing microfluidic models for cancer-immune interactions, which treat each chip as a single experiment and use traditional biological measurement techniques<sup>10,31</sup>. By comparison the droplet format provides several unique features, including the ability to merge many droplet pairs at a well-defined time, thus providing a common starting time of the parallel experiments<sup>23,32</sup>. Moreover, the encapsulation within droplets allows the conditions in each of the parallel experiments to be well controlled, thus allowing for massively multiplexed experiments on a single device.

Here these technical advantages are associated with probabilistic modeling to infer key biological information about the ability of CTLs to sense and respond to the tumor, by relating the
<figure><img src="image_19.png" /><figcaption>**Fig. 5** CTL number and collective behavior determine probability of killing. **a** Two representative images showing CTL clustering on the spheroid during first fragmentation event. The white circles have a radius of 30 µm around fragmenting cell. **b** Distribution of the CTL numbers within a radius ≤30 µm around the fragmentation areas (N = 31). **c** Sketch summarizing the observed trends: CTLs (green) migrate on the surface of the spheroid and cluster together in particular regions, where fragmentation of B16 cells (red) is observed. **d** During a time-interval Δt, a spheroid with *n* CTLs attached to it has a probability Γ<sub>frag</sub>Δt of fragmenting. **e** Illustration of possible scenarios of the number of CTLs on the spheroids and the apparition of fragmentation. **f** The fragmentation rate Γ<sub>frag</sub> can be modeled as the result of independent CTLs interacting with the spheroid, with an individual fragmentation rate per CTL worth ρ. Conversely, the fragmentation rate can be viewed as the result of a collaborative process. **g** Estimates of Γ<sub>frag</sub> as a function of the number of CTLs on the spheroid *n*. Experimental measurements (N = 84, black dots) are fitted with an exponential (dashed red line) compared with the results of the independent CTL model (dashed blue line). Source data are provided as a Source data file.</figcaption></figure>

spatiotemporal dynamics of the CTLs with the outcome for the tumor spheroid. Specifically, we find that the first CTL-tumor cell contact, which occurs randomly, triggers a positive feedback loop that leads to an accelerated accumulation of CTLs on the spheroid. Later, CTLs form clusters on the spheroid that enhance their ability to kill the target cells, leading to tumor rejection.

Several mechanisms may account for the collaborative CTL accumulation and killing. Chemokines that are both sensed and produced by T cells have the ability to drive their swarming behavior¹⁵. Cooperative killing, in which multiple sublethal cytotoxic hits synergize to induce target cell killing, has been

previously described in the context of viral infection³³ and tumor development¹⁷. Alternatively, initial CTL-tumor cell interactions may facilitate tumor destruction by increasing MHC class I expression through IFN-γ production and diffusion in the tumor microenvironment for example²⁵,³⁴. As illustrated here, our approach helps support and generate new hypotheses that can be subsequently dissected at the molecular level.

Looking ahead, the current platform can now be generalized to include several immune cell types and more realistic tumor models in each droplet. Here again the spatiotemporal resolution and Monte-Carlo approach will be fundamental to understand the
<figure><img src="image_20.png" /><figcaption>**Fig. 6** Combining short and long-range interactions to simulate spheroid fate. **a** The concentration-dependent killing is modeled as the result of two complementary mechanisms: long-distance cooperative attraction of CTLs to the target site and local killing cooperation on the spheroid. **b** The evolution of the spheroid fate in the droplets can be modeled as a branching process: at each time step the spheroid can fragment with a probability $\Gamma_{frag}\Delta t$ or not, and CTLs can either attach to or detach from the spheroid. **c** Simulations of the spheroid fate using the model in (b) and parameters obtained above recover the experimentally derived first-fragmentation times. Shaded area represents the 95% confidence interval of simulated data. The bold line represents the mean. **d** Experimental (black dots) and simulated (blue line) spheroid fragmentation probability as a function of the number of CTLs in the droplet. Shaded area represents the 95% confidence interval of simulated data. The bold line represents the mean. Source data are provided as a Source data file.</figcaption></figure>

<table><caption>**Table 1** Generalized linear model (GLM)⁴² results (number of samples is 96) from the statsmodel API⁴⁵.</caption><thead><tr><th>Variable</th><th>Coeff. β</th><th>Std err</th><th>z</th><th>p</th><th>[0.025</th><th>0.975</th></tr></thead><tbody><tr><td>n</td><td>1.84</td><td>0.786</td><td>2.342</td><td>0.019</td><td>0.300</td><td>3.381</td></tr><tr><td>N</td><td>1.5567</td><td>0.707</td><td>2.203</td><td>0.028</td><td>0.172</td><td>2.942</td></tr><tr><td>Area</td><td>-0.2175</td><td>0.320</td><td>-0.681</td><td>0.496</td><td>-0.844</td><td>0.409</td></tr><tr><td>t<sub>1</sub></td><td>-0.5405</td><td>0.999</td><td>-0.541</td><td>0.588</td><td>-2.498</td><td>1.417</td></tr><tr><td>t<sub>2</sub></td><td>1.1746</td><td>1.732</td><td>0.678</td><td>0.498</td><td>-2.221</td><td>4.570</td></tr><tr><td>t<sub>3</sub></td><td>-2.3092</td><td>2.635</td><td>-0.877</td><td>0.381</td><td>-7.473</td><td>2.854</td></tr><tr><td>t<sub>4</sub></td><td>1.4984</td><td>2.191</td><td>0.684</td><td>0.494</td><td>-2.795</td><td>5.792</td></tr></tbody></table>

n is the maximum recorded number of CTLs on the spheroid and N is the number of T cells detected in the droplet, Area is the spheroid surface and t₁ to t₄ are the first times at which 1 to 4 CTLs are detected on the spheroid surface.

causality of the interactions and the effect of 3D geometry. Finally, working with patient-derived organoids³⁵,³⁶ will have important implications for personalized medicine.

## Methods

### Experiments and analysis

*Statistics and reproducibility.* The box plots represent the quartiles of the distribution in all figures where this data representation method is used (i.e., Figs. 1d,

2g, h, 3f-h, and 4b) as detailed in the Seaborn library documentation³⁷. The whiskers describe the distribution minima and maxima.

The statistical tests used to compare distributions are Mann-Whitney-Wilcoxon tests with two-sided Bonferroni corrections. For all figures concerned (i.e., Figs. 2g, h and 3f, g) we have the following mapping for *p*-value annotations: ns: $5.00 \times 10^{-2} < p \le 1$; *: $1.00 \times 10^{-2} < p \le 5.00 \times 10^{-2}$; **: $1.00 \times 10^{-3} < p \le 1.00 \times 10^{-2}$; ***: $1.00 \times 10^{-4} < p \le 1.00 \times 10^{-3}$; ****: $p \le 1.00 \times 10^{-4}$.

One representative image from a data set of total 84 droplets in the case of B16-Ova and 54 in the case of B16-WT has been included in Figs. 1e and 2a, e, f, respectively. A representative image sequence from a total of 61 droplets has been used for Fig. 4a, b, respectively. One representative image from a total of 200 droplets each for aqueous, 2.025 and 4.05 mg/ml condition is being used for Supplementary Fig. 1a.

The data used to generate Figs. 2d, i, 3g, h, and 6c, d are pooled from three different microfluidic chips. Each droplet in these experiments corresponds to a unique experimental condition (cell number, spheroid size, spatial distribution of cells, initial condition of the T cells, etc.). The stochastic models then provide values of the parameters $\lambda$ and $\Gamma_{frag}$ by tracking the co-evolution of different cellular events within each droplet. As such each droplet in each microfluidic device should be treated as an independent replicate of the stochastic evolution. The bootstrapping analysis of the statistics confirms that the use of different subsets of these data does not modify the results.

*Tumor cells.* B16-WT melanoma and B16-OVA albumin peptide (residues 257–264) expressing cell-lines (B16-OVA) were maintained in RPMI 1640 (Fischer Scientific - 12027599) media containing 10% FBS and 1% penicillin-streptomycin antibiotics and were maintained at 5% CO₂ and 37 °C. B16.F0 (B16) and B16.F0-Ova (B16-OVA) melanoma cells were kindly provided by Claude Leclerc.
### Generation of OVA-specific cytotoxic T lymphocytes (CTLs).
Ubi-GFP Rag1-/- OT-1 TCR mice were bred in our animal facility under specific pathogen-free conditions. The mice were bred and managed by Institut Pasteur's animal facility with a central air conditioning equipment which maintains constant temperature of 22 ± 2 °C. Air is renewed at least 20 times per hour in animal rooms. Fluorescent light is provided with a 14:10-h light:dark cycle. Humidity is monitored but not controlled and in the range of 25-65%. Splenocytes were isolated from Ubi-GFP OT-1 TCR transgenic mice and red blood cells were removed by ammonium-chloride-potassium lysis. One-third of the cells was then pulsed with 50 µM of Ova257-264 peptide (SIINFELK) for 2 h at 37 °C in 1 mL total volume of RPMI medium 1640-GlutaMAX. The rest of the cells were incubated at 37 °C in 15 mL of complete medium (RPMI medium 1640-GlutaMAX supplemented with 10% heat-inactivated fetal bovine serum, 50 µg/mL penicillin, 50 µg/mL streptomycin, 1 mM sodium pyruvate, 10 mM HEPES, and 50 µM β-mercaptoethanol). After 2 h, the two populations were mixed and cultured for 3 days. Cells were then subjected to Ficoll gradient centrifugation to remove dead cells and thus select live OVA-specific CTLs, and cultured in complete RPMI medium, supplemented with human IL-2 (10 ng/mL; R&D) for 2 additional days. All animal studies were approved by the Institut Pasteur Safety Committee in accordance with French and European guidelines (CETEA 2017-0038).

### Microfabrication, microfluidic setup, and droplet formation.
The PDMS-based microfluidic device on which the experiments were conducted is precisely described in refs. 20,23,38. Preceding droplet production, the chip in filled with fluorinated FC40 (3M) oil mixed with 2% (v/v) FluoroSurfactant (Ran Biotechnologies) and cooled at -20 °C for 2 h to prevent Matrigel gelification during the loading. The primary droplets are produced as in ref. 20. The aqueous phase is composed of RPMI media, Corning Matrigel (Dutscher Dominique - 354234) and B16 melanoma cells at a concentration of 1.5 × 10⁶ cells/mL. The primary droplets have a volume of ca 50 nL. The chip is then placed in the incubator at 37 °C leading to Matrigel gelification.

The secondary droplets (volume ca. 10 nL) were produced with an aqueous phase containing RPMI media and OT-1-CD8+GFP cells at a concentration of 1.5 × 10⁶ cells/mL. The secondary droplets were further trapped in the triangular regions adjacent to the individual hexagonal wells already present with primary droplets encapsulated with B16 spheroids (see Fig. 1). To fuse the primary with the secondary droplets, 20% (v/v) of 1H,1H,2H,2H-perfluoro-1-octanol (PFO) (Sigma-Aldrich) was dissolved in NoveCTM-7500 Engineered Fluid (3M) and was perfused in the microfluidic chip. This caused the fusion of adjacent droplets. After the fusion of droplets, a fresh solution of FC40 and Fluorosurfactant was flushed-in to remove the PFO in the microfluidic chamber.

### Spheroid formation in different matrigel concentrations.
We primarily used the Matrigel concentration of 2.0 mg/ml in order to encapsulate the droplets with B16 cells. Twenty-four hours after droplet loading, the B16 cells self-assemble into a single 3D spheroid of B16 tumor cells (Fig. 1b). Using concentrations higher than 2.0 mg/ml (4.05 mg/ml) often leads to multiple spheroids (Supplementary Fig. 1a) located at different droplet heights. Conversely, droplets with pure aqueous media (without Matrigel) resulted in single spheroid formation but lacked the Extra-cellular Matrix (ECM) necessary for the migration of T cells (Supplementary Fig. 1a).

### Viability assay (Fig. 1).
B16-OVA spheroid-containing droplets were made according to the protocol described above with the addition of Propidium Iodide (PI) (Sigma - P4864) at a concentration of 3 µM. The chip was then imaged at 24 and 48 h after seeding. Only spheroids positioned in the center of the droplets were imaged, in order to avoid artifacts due to the microscopy. Using a custom-made Imagej macro, the area of red fluorescent signal by PI was measured and the percentage fluorescent area was calculated when compared to the complete spheroid area.

### Apoptosis assay (Fig. 4).
Caspase-3/7 Red Apoptosis Assay Reagent (Essen Bioscience - 4704) was used at 2 µM concentration (added during primary droplet formation) in order to visualize the apoptotic cells in the spheroids.

### Microscope imaging.
Images were captured using a Nikon Ti2 motorized epifluorescence microscope with a ×20 objective lens. The illumination was produced by a Lumencor LED light source and the images were captured by a Hamamatsu C13440-20CU SCMOS camera. Raw data collection was done using imaging software Nikon Elements (version 5.11.01, Build 1367).

### Image analysis.
A specific image analysis pipeline was developed in order to extract physical and biological variables from the time-lapse movies. The routines were written in Python (version 3) programming language and use several open-source libraries³⁹⁻⁴¹. The code is available on the GitHub repository associated to this publication. Basic data visualization and quantification of spheroid size distribution is done using custom-made macro in Imagej (version 1.53 f51).

### Determining the positions of CTLs with regard to the spheroid.
Each droplet image is multi-channel consisting in brightfield and FITC (510 nm) channels. The former enables well identification and spheroid segmentation, whereas the latter enables us to record CTL positions.

To segment the spheroid in a given droplet we rely upon a border-detection routine based on Laplace filtering. The spheroids properties such as size and position are recorded and stored. This procedure is repeated at every time step for more robustness. The CTLs are detected using the fluorescence channel, and their positions are stored. This information is then crossed with the spheroid positions extracted above to determine the position of the CTLs relative to the spheroid in the droplet (on/off the spheroid). For each time-step an image with the raw image, the mask covering the detected spheroid and the relative positions of the CTLs is generated for manual verification a-posteriori of the algorithm efficacy. Faulty wells (with for example mis-detected spheroids) are then thrown away before data analysis. The full image-analysis algorithm is available on the GitHub repository associated with this publication.

### λ<sub>in</sub> and λ<sub>out</sub> estimation and simulation.
The attachment and detachment rates are estimated using a bootstrapping procedure. From a subset of the experimental accumulation plots (Fig. 3d) we count the attachment and detachment events and estimate the values of λ<sub>in</sub> and λ<sub>out</sub>. We repeat the procedure 50 times selecting for each iteration 70% of the total number of spheroids in each condition (84 B16-OVA and 81 B16-WT spheroids). This gives us a distribution of the estimated parameters λ<sub>in</sub> and λ<sub>out</sub> represented in Fig. 3. The code necessary for this estimation is available on the GitHub repository associated with this publication.

## Mathematical models

### Interpreting the statistics of cell displacements from 2D images.
After the introduction of the T cells into the Matrigel droplets, they explore the space of the droplet in all three dimensions (3D). For practical considerations, however, the imaging was limited to a single plane within the droplet, leading to two-dimensional (2D) slices of the drops. It is therefore important to consider how to interpret the statistics of 2D measurements of cells moving in 3D. Under the hypothesis of a perfectly isotropic medium, we can study the influence of 2D projection of a 3D motion on the mean square displacement (MSD). For a displacement during a time-interval τ, we can write the following relationship for the displacement vector **r**:

$$ \langle \|\mathbf{r}\| \rangle_{\tau} = \left\langle \sqrt{r_{x}^{2} + r_{y}^{2} + r_{z}^{2}} \right\rangle_{\tau} = \sqrt{3} \langle |\mathbf{r}_{k}| \rangle_{\tau}, \quad (1) $$

where $r_i$ is the displacement in the direction i and $\langle \cdot \rangle$ denotes an ensemble average. Therefore, we should have:

$$ \langle \|\mathbf{r}_p\| \rangle_{\tau} = \left\langle \sqrt{r_x^2 + r_y^2} \right\rangle_{\tau} = \sqrt{\frac{2}{3}} \langle \|\mathbf{r}\| \rangle_{\tau}, \quad (2) $$

Correcting the velocity values from experimental 2D measurements yield estimated 3D velocities on the order of 6.3 µm/min, which are in line with observations in vivo²⁵.

Also, the scaling of the mean square displacement as a function of τ is not impacted by the projection from 2D to 3D as demonstrated by the calculation below:

$$ \langle \|\mathbf{r}_p\|^2 \rangle_{\tau} = \langle \|\mathbf{r}_x\|^2 + \|\mathbf{r}_y\|^2 \rangle_{\tau} = \frac{2}{3} \langle \|\mathbf{r}\|^2 \rangle_{\tau} \quad (3) $$

The calculation also shows that the effective diffusion coefficient $D_p$ of a three-dimensional track estimated from a 2D measurement is worth 2/3 of the “true” diffusion coefficient $D$.

### CTL velocity measurements depend on the imaging frame rate.
Measuring the motility properties of particles—or cells—undergoing random motion is heavily dependent on the sampling rate. Indeed, since the cells can move back and forth, we expect the net distance traveled Δr(τ) over a time τ to be sub-linear. By measuring the average distance traveled for different lag times we extract the experimental dependence in Supplementary Fig. 2c, d. We find that the distance traveled scales as a power law with an exponent smaller than 1 (which corresponds to ballistic motion): Δr(τ) ~ τ<sup>0.6</sup>. This translates to an experimental average velocity which scales as: v(τ) ~ τ<sup>-1.4</sup>.

### Testing the statistical power of the different experimental observables on the killing.
In the same chip we can record heterogeneous spheroid outcomes; some spheroids fragment very fast (complete destruction at 8 h), whereas others are left unscathed at 14 h (see Fig. 4d).

Since the secondary droplet contains variable CTL numbers dependent on the initial cell concentration, this leads to a range of CTL numbers in the main droplet after droplet fusion (see Supplementary Fig. 1b). In addition to the number of CTLs in the droplet N, we record several other features: the first to the fourth contact times ($t_1$ to $t_4$, indicating the moment at which the number of cells on the spheroid goes above 1 to 4 cells), the spheroid projected area and the maximum number of CTLs on the spheroid within the 14 h of the experiment duration n.

We conduct the test with a generalized linear model⁴². The total observation sample size is of 96 events. This test enables us to study the influence and statistical power of each variable on visible spheroid fragmentation at 14 h. We see that only two variables have p-value below 0.05: the total number of CTLs in the well ($p_n = 0.019$, $\beta_n = 1.84$) and the maximum number of CTLs on the spheroid during
the time-lapse ($p_N = 0.028$, $\beta_N = 1.56$). These two measures are very correlated (Supplementary Fig. 4); a high number of CTLs in the droplet increases the chances of having a high number of CTLs on the spheroid. The regression coefficients $\beta_n$ and $\beta_N$ were positive in both cases, confirming the positive correlation between these variables and the fragmentation probability. Interestingly, the other variables do not significantly predict the final spheroid state despite varying levels of correlation with spheroid death.

Inferring $\Gamma_{frag}$ from experimental data. At each time frame, the CTLs on the spheroid can either leave the spheroid, cause a fragmentation event or do nothing. Conversely, CTLs in the gel can either stay in the gel or attach to the spheroid. We consider the killing to be independent of the attachment/detachment process and only to depend on the number of CTLs on the spheroid at any given time. From the point of view of the spheroid, we can therefore model the killing process as a Bernoulli process where the probability of seeing a fragmentation event only depends on the number $n$ of CTLs on it at that time, and of a fragmentation rate $\Gamma_{frag}$. $\Gamma_{frag}$ is the probability that the $n$ CTL(s) on the spheroid cause a fragmentation during a single time-interval $\Delta t$. An implicit assumption of this model is that each time-step is independent: the probability of fragmenting during a time-interval only depends on the number of CTLs on the spheroid at the time and not on the previous CTL actions.

From the experimental fragmentation times and the observed evolution of the number of CTLs on the spheroid, we can infer the evolution of $\Gamma_{frag}$ as a function of $n$. The probability that the first fragmentation event happens at time $T_{frag}$ after exactly $k$ time-steps is a function of $\Gamma_{frag}$ and of the number $n$ of CTLs on the spheroid as a function of time. $n$ in bold font is an uni-dimensional vector where the $i$-th component is the number of CTLs on the spheroid at time $i: n_i$. Since we suppose each time-step to be independent of the previous ones, the probability becomes:

$$P(T_{frag} = k | \mathbf{n}, \Gamma_{frag}) = \Gamma_{frag}(n_k) \prod_{t}^{k-1} (1 - \Gamma_{frag}(n_t)). \quad (4)$$

Since each droplet is independent, we can also estimate the probability of observing the set of fragmentation times $\{T_{frag}\}_i$:

$$P(\{T_{frag}\}_i | \{\mathbf{n}\}_i, \Gamma) = \prod_i P(T_{frag,i} | \mathbf{n}_i, \Gamma_{frag}). \quad (5)$$

Now we can use maximum-likelihood estimation (ML) to extract the value of $\Gamma_{frag}$ that maximizes the probability to observe the experimental fragmentation times given the experimentally observed time-series $\mathbf{n}$. In short, we calculate $\Gamma_{frag}$:

$$\Gamma_{frag} = \arg \max_{T_{frag}} P(\{T_{frag}\}_i | \{\mathbf{n}\}_i, \Gamma_{frag}). \quad (6)$$

After some algebra this yields the following value for $\Gamma_{frag}(n)$:

$$\Gamma_{frag}(n) = \frac{K_{Death,n}}{K_{Death,n} + K_n}, \quad (7)$$

where $K_{Death,n}$ is the number of times we observe a fragmentation event happening with exactly $n$ CTLs on the spheroid, and $K_n$ the total number of times $n$ CTLs are on the spheroid before any fragmentation event. We can verify that if there never is any death observed with $n$ CTLs on the spheroid, then $K_{Death,n} = 0$ and $\Gamma_{frag}(n) = 0$ too. Reversely, if there is a fragmentation event as soon as $n$ cells are on the spheroid, then $K_n = 0$ and $K_{Death,n} > 0$, which gives $\Gamma_{frag}(n) = 1$.

Probabilistic modeling of cooperative vs. independent killing of B16 cells on the spheroid by CTLs. In Fig. 5g, we show that the probability of a spheroid to undergo a fragmentation event increases with the number of CTLs on the spheroid surface. Furthermore, the maximum T cell number on the spheroid is a key variable predicting the final spheroid state. However, we do not know if this increase in $\Gamma_{frag}$ is the result of the accumulation of independent random events or if it is the signature of cooperation between T cells.

We consider that the probability for a given T cell to cause fragmentation is independent of the presence of the other T cells. In quantitative terms, the probability of a fragmentation event occurring within the independent cell hypothesis is given by:

$$\Gamma_{frag}(n) = 1 - (1 - \rho)^n, \quad (8)$$

where the total number of CTLs on the spheroid is given by $n$ and $\rho$ is the probability of a single CTL to kill during a time-interval $\Delta t$. Fitting Eq. (8) to the experimental measurements we see that the model does not accurately reflect experimental results (see Fig. 5g).

A first possible improvement for the independent CTL hypothesis is proposed by ref. 33 and consists of accounting for the heterogeneous nature of the CTL population. Indeed, it is well-known that CTLs differ widely in efficacy$^{43}$. In order to test this hypothesis, we now model the CTL population as composed of independent cells, but the probability of being a fragmentation-causing CTL $\rho$ is now itself drawn from a probability distribution $\mathcal{K}$ of probability density $f$. The average value of $\rho$ is given by $\langle \rho \rangle = \int \rho f(\rho) dp$. We choose to not specify $\mathcal{K}$ to

preserve generality. Then Eq. (8) becomes:

$$\begin{aligned} \Gamma_{frag}(n) &= 1 - E \left[ \prod_{i=1}^{n} (1 - \rho_i) \right] \\ &= 1 - \left( \int (1 - \rho) f(\rho) d\rho \right)^n \\ &= 1 - (1 - \langle \rho \rangle)^n. \end{aligned} \quad (9)$$

Provided that the statistics are sufficient, the result is similar to that in Eq. (8). Thus, the expected fragmentation probability is going to display a parabolic profile, whatever the heterogeneity profile of the CTL population is. Therefore the initial CTL population heterogeneity cannot explain the fragmentation profile recovered experimentally. Thus the CTLs cooperate to cause the fragmentation of the spheroid.

Estimating the attachment and detachment rates ($\lambda_{in}$ and $\lambda_{out}$). A dedicated imaging and analysis pipeline was developed for this project in order to segment the spheroids and detect the CTLs in the focal plane within the droplets. The spheroid detection outputs a mask that corresponds to the spheroid position at every frame. Then, by detecting the CTLs that are in contact with this mask, it is possible to label each detected CTL as being on the spheroid or not. This information allows us to track individual arrival and leaving events on each spheroid. The high temporal resolution of the imaging (one image every 2 min) ensures that at most one event happens between two consecutive images (see Supplementary Movie 2).

The statistics of attachment and detachment are then used within a Markov chain model to estimate the attachment and detachment rates ($\lambda_{in}$ and $\lambda_{out}$) on the spheroid as a function of the number of CTLs on the spheroid. The evolution of these rates is then used to detect whether the attraction of the cells takes place through an active process or whether it is due to random chance.

Three hypotheses are required for the purpose of calculating $\lambda_{in}$ and $\lambda_{out}$: First, we suppose that the CTLs have had enough time to distribute themselves in the droplet before they begin to attach/detach. This is justified by noting that the average first contact time is large compared with the time to cross the droplet. Second, each of the attachment/detachment events occurs independently of other attachment/detachment events. Third, We hypothesize that the waiting times for each individual CTL arriving/detaching on the spheroid is exponentially distributed (i.e., the CTL is memory-less), and the waiting time parameter only depends on the number $n$ of CTLs already on the spheroid. Finally, we only consider attachment and detachment events and not the movement of the CTLs in the droplet. Using these main hypotheses it is now possible to derive the model itself.

Consider a droplet with $N$ CTLs and a single spheroid in it. The spheroid can be in one of $N+1$ states, each state representing the number of CTLs on the spheroid (i.e., if the spheroid is in state $n$, then it has $n$ CTLs on it). The experimental time series of discrete observations is denoted $\{n_t\}_t$ where $n_t$ is the number of CTLs on the spheroid at time $t$. The set of transition probabilities is written $\mathbf{p} = \{\mathbf{p}_{k,l}\}_{k,l'}$ where $\mathbf{p}_{k,l}$ is the probability to go from $k$ to $l$ CTLs on the spheroid.

The system and its evolution can be described as a Markov chain. We call $\mathbf{P}(\{n_t|n_{t-1}\})$ the probability of having $n_t$ CTLs on the spheroid at time $t$ knowing that there were $n_{t-1}$ CTLs on it at $t-1$. $\mathbf{P}(\{n_t\}_t)$ is the probability of observing the time-series $\{n_t\}_t$. By defining $a_{k,l}$ as the number of transitions from $k$ to $l$ CTLs on the spheroid, it is possible to link $\mathbf{P}(\{n_t\}_t)$ to the attachment and detachment statistics. Indeed, by virtue of the chain rule we have:

$$\mathbf{P}(\{n_t\}_t | \mathbf{p}) = \prod_t \mathbf{P}(n_t | n_{t-1}) = \prod_{k,l} \mathbf{p}_{k,l}. \quad (10)$$

Equation (10) then allows us to infer the transition probabilities $\mathbf{p}$. We conduct a maximum-likelihood estimate of $\mathbf{p}$, under the constraint that the total out-going probability must add up to one (i.e., $\sum_k \mathbf{p}_{k,l} = 1$). This yields:

$$\mathbf{p}_{k,l} = \arg \max_{\mathbf{p}_{k,l}} \log \left( \mathbf{P}(\{n_t\}_t | \mathbf{p}) \right) = \frac{a_{k,l}}{\sum_j a_{k,j}}. \quad (11)$$

In practice, $\mathbf{p}_{k,l}$ is calculated from Eq. (11) by dividing the number of observed transitions from $k$ to $l$ cells by the total number of transitions starting from $k$ cells (i.e., the number of times there are $k$ cells on the spheroid).

Until now, we have studied the problem from the viewpoint of the spheroid: the probabilities $\mathbf{p}_{k,l}$ do not take into account the combinatorial aspect of the transitions (i.e., there are often more than one CTL in the gel that may attach to the spheroid). It is therefore important to look at the problem from the viewpoint of the individual CTLs and obtain the individual T cell attachment probabilities. Supposing that they attach independently of the other CTLs in the gel, we can calculate the individual probability for each CTL to attach to the spheroid, $p_{in}(n)$, as:

$$p_{in}(n) = 1 - (1 - \mathbf{p}_{n,n+1})^{\frac{1}{n}}. \quad (12)$$

Given the exponential distribution of waiting times for attachment of the CTLs, the probability that a CTL arrives on a spheroid with $n$ CTLs on it during a time-interval $\Delta t$
is $p_{in}(n) = 1 - \exp(-\lambda_{in}(n)\Delta t)$. Then the transition rates can be written as

$$ \lambda_{in}(n) = - \frac{\log(1 - p_{in}(n))}{\Delta t}. \qquad (13) $$

Finally combining Eqs. (11), (12), and (13), provides a formula for $\lambda_{in}(n)$:

$$ \lambda_{in}(n) = - \frac{\log(1 - p_{n,n+1})}{(N - n)\Delta t}. \qquad (14) $$

This formula allows us to directly estimate the attachment rate from the experimentally measured attachment events. To estimate the uncertainty of our estimations we conduct a bootstrapping scheme, where we randomly select droplets, pool them together and calculate the hitting parameters for this experiment subset. We then repeat the procedure over multiple subsets (typically 50) to get a distribution of arrival rates $\lambda_{in}(n)$, which enables us to calculate the mean and the variance of $\lambda_{in}$ for each $n$. The code to conduct the bootstrapping procedure is freely available on the GitHub repository associated with this publication.

The process is repeated while counting leaving events to estimate the leaving rate $\lambda_{out}(n)$ as a function of the number of CTLs on the spheroid, also using Eq. (14). The evolution of $\lambda_{in}$ and $\lambda_{out}$ is shown in Fig. 3.

It is worth noting that the evolution of the number of CTLs on the spheroid can be described using a master equation:

$$ \frac{\partial \mathbb{P}(n, t)}{\partial t} = \mathbf{q}_{n-1,n} \mathbb{P}(n-1, t) + \mathbf{q}_{n+1,n} \mathbb{P}(n+1, t) - (\mathbf{q}_{n,n+1} + \mathbf{q}_{n,n+1}) \mathbb{P}(n, t), \qquad (15) $$

where $\mathbf{q}_{k,l} = \mathbf{p}_{n,n+1}/\Delta t$ is the transition rate from $k$ to $l$ CTLs on the spheroid. Knowing that for small values of $\mathbf{p}_{n,n+1}$, we have $\mathbf{p}_{n,n+1} \approx (N-n)\Delta t\lambda_{in}(n)$, Eq. (15) can be rewritten as

$$ \frac{\partial \mathbb{P}(n, t)}{\partial t} = (N - n + 1)\lambda_{in}(n - 1)\mathbb{P}(n - 1, t) + (n + 1)\lambda_{out}(n + 1)\mathbb{P}(n + 1, t) - ((N - n)\lambda_{in}(n) + n\lambda_{out}(n))\mathbb{P}(n, t). \qquad (16) $$

The master equation would allow an analytical solution to obtain $\lambda_{in}$ and $\lambda_{out}$ if these parameters were constant for all values of $n$. Since this is not the case, obtaining them from Eq. (16) would require a numerical approach. Nevertheless, it is known that the master-equation and the Markov chain represent the same Markov process$^{44}$, so both methods should lead to the same estimated parameters.

*Relating the spheroid fate to the number of CTLs in the droplet.* We know that due to a positive feedback loop, CTLs exhibit different accumulation rates depending on the number of CTLs on the target spheroid. We have estimated the experimental distribution of the hitting probabilities above. We also know that a certain number of CTLs on the spheroid are necessary for a high killing probability (Fig. 5). Combining these two components we have written an algorithm which reproduces the accumulation and fragmentation processes in silico.

First, we simulate the accumulation process. We define the number of CTLs on the spheroid $n$ and the number of CTLs in the gel $n_{gel}$. Thanks to the estimated parameters $\lambda_{in}$ and $\lambda_{out}$ we can calculate the attachment and detachment probabilities $p_{in}$ and $p_{out}$. From Eq. (13) we get the equation for $p_{in}$: $p_{in}(\lambda_{in}(n), \Delta t) \approx 1 - \exp(n_{gel}\lambda_{in}(n)\Delta t)$. We then do a Bernoulli draw with probability $p_{in}$ (resp. $p_{out}$) to simulate the attachment (resp. detachment) of T cells on the spheroid. Then, to simulate the fragmentation of the spheroid, we conduct a Bernoulli draw of parameter $\Gamma_{frag}(n)\Delta t$ to determine whether or not the spheroid fragments at that time-step.

We repeat these two steps iteratively until the experiment is elapsed. This procedure gives us the simulated fate of the spheroid over the course of the experiment. We have access to the evolution of the number of T cells on the spheroid with time, the status of the spheroid (fragmented or intact) and the first fragmentation time. By re-running this procedure 50 times for a given number of CTLs in the droplet $N$ we get a simulated fragmentation probability that is presented in Fig. 6d.

**Reporting summary.** Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

All raw data to reproduce the graphs in the figures are available with the article. Any further data required will be provided upon request. Source data are provided with this paper.

## Code availability

The code to extract spheroid and CTL positions from the image data is available at the following address: https://github.com/BaroudLab/CTL_tracking_in_droplets. Data processing and data analysis codes are available at this address: https://github.com/BaroudLab/CTL_analysis.

Received: 1 April 2021; Accepted: 5 May 2022;

Published online: 03 June 2022

## References

1. Sharma, P. & Allison, J. P. The future of immune checkpoint therapy. *Science* **348**, 56–61 (2015).

2. June, C. H. & Sadelain, M. Chimeric antigen receptor therapy. *N. Engl. J. Med.* **379**, 64–73 (2018).

3. Fridman, W. H., Pagès, F., Sautes-Fridman, C. & Galon, J. The immune contexture in human tumours: Impact on clinical outcome. *Nat. Rev. Cancer* **12**, 298–306 (2012).

4. Vasconcelos, Z. et al. Individual human cytotoxic T lymphocytes exhibit intraclonal heterogeneity during sustained killing. *Cell Rep.* **11**, 1474–1485 (2015).

5. Boulch, M., Grandjean, C. L., Cazaux, M. & Bousso, P. Tumor immunosurveillance and immunotherapies: a fresh look from intravital imaging. *Trends Immunol.* **40**, 1022–1034 (2019).

6. Breart, B., Lemaitre, F., Celli, S. & Bousso, P. Two-photon imaging of intratumoral CD8. T cell cytotoxic activity during adoptive T cell therapy in mice. *J. Clin. Invest.* **118**, 1390–1397 (2008).

7. Cazaux, M. et al. Single-cell imaging of CAR T cell activity in vivo reveals extensive functional and anatomical heterogeneity. *J. Exp. Med.* **216**, 1038–1049 (2019).

8. Marangoni, F. et al. The transcription factor NFAT exhibits signal memory during serial T cell interactions with antigen-presenting cells. *Immunity* **38**, 237–249 (2013).

9. Kumar, V. & Varghese, S. Ex vivo tumor-on-a-chip platforms to study intercellular interactions within the tumor microenvironment. *Adv. Healthc. Mater* **8**, 1–12 (2019).

10. Shelton, S. E., Nguyen, H. T., Barbie, D. A. & Kamm, R. D. iScience II Engineering approaches for studying interactions and immunotherapy. *Science* **24**, 101985 (2021).

11. Bar-Ephraim, Y. E., Kretzschmar, K. & Clevers, H. Organoids in immunological research. *Nat. Rev. Immunol.* **20**, 279–293 (2020).

12. Schnalzger, T. E. et al. 3D model for CAR-mediated cytotoxicity using patient-derived colorectal cancer organoids. *EMBO J.* **38**, 1–15 (2019).

13. Boussommier-Calleja, A., Li, R., Chen, M. B., Wong, S. C. & Kamm, R. D. Microfluidics: a new tool for modeling cancer-immune interactions. *Trend. Cancer* **2**, 6–19 (2016).

14. Park, S. E., Georgescu, A. & Huh, D. Organoids-on-a-chip. *Science* **364**, 960–965 (2019).

15. GaleanoNiño, J. L. et al. Cytotoxic T cells swarm by homotypic chemokine signalling. *eLife* **9**, 10 (2020).

16. Budhu, S. et al. CD8+ T cell concentration determines their efficiency in killing cognate antigen-expressing syngeneic mammalian cells in vitro and in mouse tissues. *J. Exp. Med.* **207**, 223–235 (2010).

17. Weigelin, B. et al. Cytotoxic T cells are able to efficiently eliminate cancer cells by additive cytotoxicity. *Nat. Commun.* **12**, 1–12 (2021).

18. Manaster, Y. et al. Reduced CTL motility and activity in avascular tumor areas. *Cancer Immunol. Immunother.* **68**, 1287–1301 (2019).

19. Nicolas-boluda, A. & Donnadieu, E. Comparative immunology, microbiology and infectious diseases obstacles to T cell migration in the tumor microenvironment. *Comp. Immunol. Microbiol. Infect. Dis.* **63**, 22–30 (2019).

20. Sart, S., Tomasi, R. F.-X., Amselem, G. & Baroud, C. N. Multiscale cytometry and regulation of 3D cell cultures on a chip. *Nat. Commun.* **8**, 469 (2017).

21. Sart, S. et al. Mapping the structure and biological functions within mesenchymal bodies using microfluidics. *Sci. Adv.* **6**, eaaw7853 (2020).

22. Abbyad, P., Dangla, R., Alexandrou, A. & Baroud, C. N. Rails and anchors: guiding and trapping droplet microreactors in two dimensions. *Lab Chip* **11**, 813–821 (2011).

23. Tomasi, R. F. X., Sart, S., Champetier, T. & Baroud, C. N. Individual control and quantification of 3D spheroids in a high-density microfluidic droplet array. *Cell Rep.* **31**, 107670 (2020).

24. Imle, A. et al. Experimental and computational analyses reveal that environmental restrictions shape HIV-1 spread in 3D cultures. *Nat. Commun.* **10**, 2144 (2019).

25. Thibaut, R. et al. Bystander IFN-γ activity promotes widespread and sustained cytokine signaling altering the tumor microenvironment. *Nature Cancer* **1**, 302–314 (2020).

26. Harris, T. H. et al. Generalized Lévy walks and the role of chemokines in migration of effector CD8+ T cells. *Nature* **486**, 545–548 (2012).

27. Fricke, G. M., Letendre, K. A., Moses, M. E. & Cannon, J. L. Persistence and adaptation in immunity: T cells balance the extent and thoroughness of search. *PLoS Comput. Biol.* **12**, 1–23 (2016).
28. Ariotti, S. et al. Subtle CXCR3-dependent chemotaxis of CTLs within infected tissue allows efficient target localization. *J. Immunol.* **195**, 5285–5295 (2019).

29. Cockburn, I. A. et al. In vivo imaging of CD8+ T cell-mediated elimination of malaria liver stages. *Proc. Natl Acad. Sci. USA* **110**, 9090–9095 (2013).

30. Scabini, M. et al. In vivo imaging of early stage apoptosis by measuring real-time caspase-3/7 activation. *Apoptosis* **16**, 198–207 (2011).

31. Miller, C. P., Shin, W., Ahn, E. H., Kim, H. J. & Kim, D. H. Engineering microphysiological immune system responses on chips. *Trends Biotechnol.* **38**, 857–872 (2020).

32. Fradet, E. et al. Combining rails and anchors with laser forcing for selective manipulation within 2D droplet arrays. *Lab Chip* **11**, 4228 (2011).

33. Halle, S. et al. In vivo killing capacity of cytotoxic T cells is limited and involves dynamic interactions and T cell cooperativity. *Immunity* **44**, 233–245 (2016).

34. Hoekstra, M. E. et al. Long-distance modulation of bystander tumor cells by CD8+ T-cell-secreted IFN-γ. *Nat. Cancer* **1**, 291–301 (2020).

35. Jenkins, R. W. et al. Ex vivo profiling of PD-1 blockade using organotypic tumor spheroids. *Cancer Discov.* **8**, 196–215 (2018).

36. Tuveson, D. & Clevers, H. Cancer modeling meets human organoid technology. *Science* **364**, 952–955 (2019).

37. Waskom, M. & the seaborn development team. *mwaskom/seaborn* (2020).

38. Amselem, G., Sart, S. & Baroud, C. N. Universal anchored-droplet device for cellular bioassays. *Methods Cell Biol.* **148**, 177–199 (2018).

39. Van Der Walt, S. et al. Scikit-image: image processing in python. *PeerJ* **2014**, 1–18 (2014).

40. McKinney, W. Data structures for statistical computing in python. in *Proceedings of the 9th Python in Science Conference*, 1697900(Scipy), 51–56 (2010).

41. Harris, C. R. et al. Array programming with NumPy. *Nature* **585**, 357–362 (2020).

42. McCullagh, P. & Nelder, J. A. *Generalized Linear Models*. (1983).

43. Khazen, R. et al. Functional heterogeneity of cytotoxic T cells and tumor resistance to cytotoxic hits limit anti-tumor activity in vivo. *EMBO J.* **40**, e106658 (2021).

44. Van Kampen, N. G. *Stochastic Processes in Physics and Chemistry* Vol. 1 (Elsevier, 1992).

45. Seabold, S. & Perktold, J. Statsmodels: econometric and statistical modeling with python. in *Proc. of the 9th Python in Science Conf.* 57 (2010).

## Acknowledgements

The authors acknowledge enlightening discussions with Clément Roussel, Andrey Aristov, Raphael Tomasi, Gabriel Amselem and all members of the Baroud team. Raphael Voituriez is gratefully acknowledged for his comments on the probabilistic descriptions. P.B. is funded by Institut Pasteur, Inserm and an advanced grant from the European Research Council (ENLIGHTEN). G.R. is funded by the Direction Générale de

l'Armement. The authors acknowledge the support of the Biomaterials and Microfluidics platform and the FabLab at Institut Pasteur.

## Author contributions

G.R., S.J., M.C., P.B., and C.N.B. designed experiments. G.R. and S.J. performed experiments and analysis. C.A., M.C., and R.K. performed experiments. G.R. developed mathematical models, image and data analysis tools. P.B. and C.N.B. supervised the research. G.R., S.J., P.B., and C.N.B. wrote the manuscript. All authors discussed the results and commented on the manuscript.

## Competing interests

C.N.B. is named inventor on several patents related to the technology. C.N.B. is also co-founder of the spinoff company Okomera. The remaining authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-30575-2.

**Correspondence** and requests for materials should be addressed to Philippe Bousso or Charles N. Baroud.

**Peer review information** *Nature Communications* thanks Christine Trinkle and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_21.png" /></figure>

**Open Access** This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022