Heterogeneous, delayed-onset killing by multiple-hitting T cells: Stochastic simulations to assess methods for analysis of imaging data

Although quantitative insights into the killing behaviour of Cytotoxic T Lymphocytes (CTLs) are necessary for the rational design of immune-based therapies, CTL killing function remains insufficiently characterised. One established model of CTL killing treats CTL cytotoxicity as a Poisson process, based on the assumption that CTLs serially kill antigen-presenting target cells via delivery of lethal hits, each lethal hit corresponding to a single injection of cytotoxic proteins into the target cell cytoplasm. Contradicting this model, a recent in vitro study of individual CTLs killing targets over a 12-hour period found significantly greater heterogeneity in CTL killing performance than predicted by Poisson-based killing. The observed killing process was dynamic and varied between CTLs, with the best performing CTLs exhibiting a marked increase in killing during the final hours of the experiments, along with a “burst killing” kinetic. Despite a search for potential differences between CTLs, no mechanistic explanation for the heterogeneous killing kinetics was found. Here we have used stochastic simulations to assess whether target cells might require multiple hits from CTLs before undergoing apoptosis, in order to verify whether multiple-hitting could explain the late onset, burst killing dynamics observed in vitro. We found that multiple-hitting from CTLs was entirely consistent with the observed killing kinetics. Moreover, the number of available targets and the spatiotemporal kinetics of CTL:target interactions influenced the realised CTL killing rate. We subsequently used realistic, spatial simulations to assess methods for estimating the hitting rate and the number of hits required for target death, to be applied to microscopy data of individual CTLs killing targets. We found that measuring the cumulative duration of individual contacts that targets have with CTLs would substantially improve accuracy when estimating the killing kinetics of CTLs.

i th CTL. We refer to the set of all samples of from a single simulation as , so that . τ X x ij ∈ X For our analysis we consider the subset of target cells which survived until the end of the simulation (containing samples ) separately to those samples where target cells were S ⊆ X s ij killed, (containing samples ). To make inferences about the CTL killing dynamics, we K ⊆ X k ij would ideally like to use both the samples of surviving ( ) and killed targets ( ), thus taking S K into account all the information available. We define a maximum likelihood estimate for the model parameters that maximises the joint probability of all of the samples : where is the probability of observing any killed sample individually, and (k ) P K ij k ij ∈ K (s ) P S ij is the probability of observing any surviving sample individually. The remaining task is to s ij ∈ S find functions expressing and in terms of the parameters for the CTL hitting process, i.e. P K P S and . If CTLs are given sufficient time to kill all targets, the waiting times until death for all η λ the samples forming the set of killed targets are distributed according to the gamma K probability density function , with a shape parameter and a rate parameter (S3A (τ |η, ) f k λ η λ  [1,2] . However for larger numbers of targets and more complex scenarios, analytical solutions become intractable. Instead our aim here was to develop a more general approach that is applicable to a broad range of experimental settings.
Consider the location one would expect a surviving target (with total attention time ) to fall in s ij the distribution of killed targets, had the simulations indeed continued until it was killed.
According to our probability function , this hypothetical killing time should be "spread" (τ |η, ) f k λ over all values greater than . This spreading should be done with a weight following τ s ij w ij , defined as: can be written as follows: Because is constant and for a specified , we may assume there exists a weighting β β > s ij , η λ function , with: which importantly allows the samples to be broken into two components: one expressed in terms of the samples, , and the second, , a correcting function for a specific set of (s ) w ij (τ ) φ η,λ and . With a functional expression of these weightings in hand we can proceed to write the η λ hypothesized probability density function, , in terms of the modified distribution of CTL (τ |η, ) f k λ attention times for killed cells, , plus our gamma distributed forecast based on the (τ |η, ) P K λ surviving cells : Eq S10 Here, is a function containing information about other factors influencing the distribution (τ ) φ η,λ of adjusted cumulative contact times for killed or surviving cells, for example CTL motility, if (τ ) such knowledge is present and one wishes to include it in the model. For our analysis we are seeking a general solution for recovering the hitting parameters and , so we treat as λ η (τ ) φ η,λ a nuisance variable that we aim to remove from the analysis without damaging our fitting results. Note that if we hold constant, we have a function containing all the samples (τ ) φ η,λ depending only on our parameters. Although by discarding any ancillary information contained in we do not estimate the absolute probability of a set of observations, by holding (τ ) φ η,λ (τ ) φ η,λ constant we may still accurately determine the parameters and . That is because, having λ η isolated the components of our model that depend on the unknown parameters, we have obtained a sufficient parametric division of our problem (see also: [3] , appendix 3) . In practise this means that we can hold the residual likelihood as a constant during fitting, thus allowing φ the probability of our samples to be expressed: , Eq S11 with indicating that due to the discarded information we have not computed the ψ (τ ) φ η,λ absolute probability , but only found the parameters and which maximise . Note P K⋂S η λ P K⋂S that this works because during the estimation process the samples and are fixed whilst k i s i η and are varied, therefore the contribution of the residual error remains the same. However, λ the success of the parameter estimation overall does depend on the extent to which the summary statistic depends on the information discarded. τ Fitting procedure for subpopulations of single-hitting CTLs It is desirable to have a procedure to compare the likelihood of multiple-hitting to other hypotheses. An alternative to the multiple-hitting hypothesis, which may also describe the emergence of high-rate killing over time, is the subpopulation model originally suggested by Vasconcelos et. al. [4] . In the subpopulation model high-rate killing is due to a CTL subpopulation that has an intrinsic capability to kill at high rates. In order to compare the subpopulation with the multiple-hitting hypothesis, we followed the strategy of Vasconcelos et.
al. [4] , who sampled the final number of killed targets per-CTL after the 12 hour experiments, then fit a Poisson mixture model to the resultant distribution. However, different to their fitting approach, we utilised a maximum likelihood approach which also integrated the initial number of targets. This takes explicit account of possible censorship of the data (e.g. due to CTLs exhausting the supply of targets). Furthermore, the initial number of targets together with the number of kills per CTL contains information about the number of surviving targets per well, making comparisons with the multiple-hitting model more straightforward.
Our first step was to recover the hitting rate parameter for uniform populations of single-hitting λ

CTLs (
). For such a uniform CTL population, can be inferred by maximising the η = 1 λ likelihood of each number of targets killed per CTL according to an n-truncated Poisson model.
If a CTL had available n targets, then the probability of observing x kills is: , Eq S12 if . When , the probability is given by: which is easily determined since the integral is the Poisson cumulative distribution function.
Note that Eq S13 (for ) is a natural consequence of our basic requirement that the n = x probability distribution should integrate to unity. To extend the Poisson fitting method to a situation with a subpopulation of high-rate killers amongst low-rate killers, we maximise the likelihood of the samples according to a 2-component, weighted, n-truncated Poisson mixture.
The probability of a single sample can be expressed for the i th CTL: where is the probability density function described in Eq S12-13, parameterised by the P relevant rate parameter for either the high-or low-rate killing subpopulation. In the absence of mechanistic insight into the origins of the hypothesised subpopulations, we use a simple binomial model such that each CTL has a fixed prior probability of belonging to the high-rate killing population. This probability is referred to as the mixture parameter m , a constant to be .
The estimation proceeds with randomly selected parameter estimates, which are then iteratively updated using an expectation maximisation procedure. Specifically, estimates for LR and HR were used to produce a posterior estimate for the weight of the i th sample , which is w i the ratio: ,

HR LR
Eq S16 The mean of the weights derived from LR and HR is an estimate m est of the mixture parameter : with m est treated as a realisation of a random variable drawn from the population whose true fraction (of CTLs with = HR ) is the parameter . This has the binomial probability which m (m) B can be used to constrain the parameter in the fit. Our maximum likelihood estimates for the m mixture parameters maximised the log likelihood function: Comparing Multiple-hitting vs subpopulation hypothesis To compare hypotheses we used the Monte-Carlo simulations where all targets shared risk equally throughout the simulations (as in main Text Fig. 1). We retained the previously tested parameters , =10 (hr -1 ) to represent the multiple-hitting hypothesis. In addition, we 0 η = 1 simulated two groups of single-hitting CTLs to represent the subpopulation hypothesis. η ) ( = 1 One group of single-hitting CTLs had a relatively low killing rate LR = 0.2(hr -1 ), the second group a higher killing rate HR =0.7(hr -1 ). For each group we performed 2 rounds of N w =1000 simulations, each round starting with a different mean for the Poisson-distributed initial n number of targets per-CTL. One round of simulations started with mean (S4A Fig, left  n = 8 column), and the other with (S4A Fig, right column).
The condition resulted in heavy censorship of the data for multiple-hitting CTLs and for n = 8 single-hitting CTLs with HR =0.7(hr -1 ), because in both cases the majority of CTLs killed all of their targets (S4B Fig, left column). We simulated the condition to ensure that our n = 8 procedure would handle such censored samples well and these simulations were not n = 8 used for testing the mixed population. Subsequently, we combined the groups with in η = 1 order to generate a mixed dataset representing the subpopulation hypothesis (S4C Fig). As before, the mixing parameter m represents the fraction of the mixed population with LR = 0.2(hr -1 ), such that the remainder (1-m) belongs to the HR =0.7(hr -1 ) group. The mixing parameter m was selected for the conditions such that the mean killing from the resultant 6 n = 1 mixed dataset would be the same as the mean killing for the multiple-hitting CTLs, which was  S4B Fig, bottom right).
We first established that our Poisson estimated rate parameters (denoted ) would be accurate λ ︿ for the individual rate parameters LR and HR when fitted separately (S4D Fig). As expected, for large sample numbers (N w =1000) the Poisson estimators were reliable (S4D Fig, dashed lines).
We also tested a smaller sample number (N w =10; S4D Fig, solid lines), for which the bounds on the estimated values were much less constrained, particularly for the HR =0.7 group λ ︿ compared to the LR =0.2 group (S4D Fig, right versus left column) and for the group with the lower compared to the group with the higher (S4D Fig, top versus bottom row). For n = 8 6 n = 1 N w =10, the differences between simulations with different and were due to the CTLs with n rate HR often eliminating all available targets, especially in those simulations with the low initial only fits with are annotated). However, when the unconstrained Gamma estimator was η > 1 applied to the multiple-hitting-generated testing data, substantial improvements in likelihood resulted (S5B Fig, dark bars). Moreover, satisfactory estimates for the parameters were achieved (S5C Fig, annotated points). To better understand this result, we examined the cumulative interactions, , for all the samples in each generating dataset (S5D Fig). The τ multiple-hitting-generated data was spread over a much smaller range of than was the τ subpopulation-generated data. This reflects the fact that under multiple-hitting, very short interactions are unlikely to lead to enough hits for killing. This explains the reduced sample density of killed targets in the multiple-hitting-generated data in comparison to subpopulation-generated data, for the lowest values of (S5D Fig, "killed"). Very lengthy τ interactions are also unlikely for multiple-hitting CTLs, since targets that have been contacting a CTL for a long time become increasingly likely to have received hits. This explains the reduced sample density of both killed and surviving targets in our datasets, comparing again multiple-hitting-generated to subpopulation-generated data, for the highest values of (S5D τ Fig, both columns). The subpopulation-generated data is in fact spread further over than τ would be the case for a uniform population, as the mixture of exponentials results in a hyperexponential distribution. Thus, the Gamma estimator prefers strongly if there is a η = 1 subpopulation of single-hitting CTLs, since has an opposing effect i.e. it concentrates the η > 1 distribution around the mean.
Considering now the unconstrained Poisson estimator; when applied to the subpopulation-derived data this also recovered satisfactory estimates for the different killing rates, and (S5C Fig, bottom row "subpopulation"). When the Poisson estimator was λ LR λ HR applied to the multiple-hitting-generated data, this produced similar estimates for and λ ︿ LR λ ︿ HR as for the subpopulation-generated data (S5B Fig, "subpopulation"), which was due to the similar distributions (S4B Fig). Thus, in this example there was a 100% false-positive rate.