Modelling the dynamics of senescence spread

Abstract Cellular senescence is a cell surveillance mechanism that arrests the cell cycle in damaged cells. The senescent phenotype can spread from cell to cell through paracrine and juxtacrine signalling, but the dynamics of this process are not well understood. Although senescent cells are important in ageing, wound healing and cancer, it is unclear how the spread of senescence is contained in senescent lesions. In the absence of the immune system, senescence could theoretically spread infinitely from one cell to another, but this contradicts experimental evidence. To investigate this issue, we developed both a minimal mathematical model and a stochastic simulation of senescence spread. Our results suggest that differences in the number of signalling molecules secreted between subtypes of senescent cells can limit the spread of senescence. We found that dynamic, time‐dependent paracrine signalling prevents the uncontrolled spread of senescence, and we demonstrate how model parameters can be determined using Bayesian inference in a proposed experiment.


| INTRODUC TI ON
Senescence is a cellular stress response that can be activated by telomere shortening, DNA damage and the up-regulation of oncogenes (d'Adda di Fagagna et al., 2003;Hayflick & Moorhead, 1961;Serrano et al., 1997). Senescent cells display cell cycle arrest and several changes to the cell, including large-scale chromatin rearrangement and the production of secretory proteins, which can induce senescence in nearby cells (the SASP, senescence associate secretory phenotype) (Acosta et al., 2008;Basisty et al., 2020;Copp'e et al., 2008;Kuilman et al., 2008;Narita et al., 2003). Senescence is necessary for many processes in tissues including development, wound healing and defence against tumour formation (Demaria et al., 2014;Kirschner et al., 2020;Muñoz-Espín & Serrano, 2014;Passos et al., 2009). When senescent cells accumulate, the molecules they release can create a pro-tumourigenic and inflammatory microenvironment, which has been linked to both the ageing process and the recurrence of cancerous tumours following treatment (Wang et al., 2020).
SASP molecules can spread the senescence phenotype through paracrine signalling, with molecules diffusing from a secreting senescent cell and binding to non-senescent cells. Several experiments have demonstrated paracrine transmission of senescence, for example, senescent cells contained within a transwell induce senescence in a population of growing cells sharing the same media (Acosta et al., 2013). SASP signalling also recruits immune cells which can clear senescent cells (Prata et al., 2018).
The terms "primary" and "secondary" senescence are used to distinguish between cells in which senescence has been induced cell-intrinsically (e.g., through DNA damage) and extrinsically (through signals from primary senescent cells), respectively (Figure 1b). Recent evidence shows that these subtypes of senescence are transcriptomically distinct Teo et al., 2019).
In addition to paracrine signalling (SASP), senescence can also spread via juxtacrine communication between cells that are in direct contact (NOTCH-mediated) (Hoare et al., 2016;Teo et al., 2019).
Secondary senescent cells created through juxtacrine signalling may differ from secondary senescent cells created through paracrine signalling, with juxtacrine secondary cells producing fewer SASP proteins, alongside other transcriptional differences (Hoare et al., 2016;Rattanavirotkul et al., 2020). Whether these two types of secondary senescent cells play different roles in the spread of senescence is as yet unclear .
Here, we define the controlled spread of senescence as a local (finite) senescence spread over an entire human lifespan. In vivo, a limited and controlled senescence spread is thought to play an important role in at least two scenarios. First, any secondary senescent cells created around a primary senescent lesion may help to amplify the signal to the immune system. Second, as a failsafe mechanism, damaged cells that failed to undergo primary senescence will become senescent in response to the primary cells ( Figure 1a).
However, multiple lines of evidence suggest that an accumulation of senescent cells and uncontrolled spread will ultimately lead to tissue dysfunction (Baker et al., 2011(Baker et al., , 2016Burd et al., 2013;Katzir et al., 2021).
Currently, there are no hypothesised mechanisms to contain the spread of senescence in the absence of the immune system. However, it has been shown that the spread of senescence is contained when immune cells are not present (Acosta et al., 2013).
In vitro experiments have shown that from a seeded circle of senescent cells, the spread of senescence is finite (senescence was found at a distance of approximately 1 mm from seeded cells) (Acosta et al., 2013). Although the medium taken from primary senescent fibroblasts is capable of inducing secondary senescence in growing fibroblasts, it is uncertain if tertiary senescence is induced when the process is repeated (Acosta et al., 2013).
Controlled senescence spread is also observed in vivo, where some senescent lesions persist over time (melanocytic naevi) and are not cleared by the immune system, but the senescence phenotype does not continue to spread (Bartkova et al., 2006). This is not due to SASP production halting, as cells continue to produce a SASP even months after the onset of senescence and enter a deep senescence state, which can lead to chronic senescence in vivo (van Deursen, 2014).
In this paper, we use both a minimal mathematical model and stochastic simulation to investigate the mechanisms that determine the dynamics of senescence spread. We identify that juxtacrine-induced senescent cells can limit senescence spread by acting as a fire break in the wildfire of paracrine signalling and that this effect is more pronounced when SASP-signalling is delayed compared with juxtacrine signalling. We propose an in vitro experiment and demonstrate that the parameters of senescence spread can be inferred using Bayesian statistics. With the emergence of large-scale descriptive efforts such as the SenNet consortium (US) (Lee et al., 2022) and the advancement of spatial transcriptomics (Kiss et al., 2022), our model can serve as a framework for further exploration as new data become available.

| Senescence spread model
We considered senescence spread in a 2D layer of cells covered by a 3D layer of media so that we could compare our results with senescence spread in vitro, as senescence is best characterised in adherent human fibroblasts (Cristofalo et al., 2003;Hayflick & Moorhead, 1961). Despite the 2D nature of our model, we can use it to infer parameters that are unchanged between 2D and 3D models of senescence spread, such as the number of ligands needed to induce senescence in a cell. To simplify the system, we described the cells as circular and non-motile, evenly and randomly distributed on the 2D plane with a specified density. Following work by Batsilas et al. (2003), we considered the cells to be covered by a layer of media of height h (as shown in Figure 1c). The density of the cells in the 2D plane is given by , and the radius of each cell is given by r cell . The number of receptors per cell is R tot . Paracrine signalling is mediated by ligands, and in our model, the ligands are SASP molecules. Ligands are emitted by senescent cells and diffuse through the surrounding media. D L is the diffusion coefficient for the ligand, and a larger diffusion coefficient results in ligands travelling further in a given time frame. Three rates are used to describe a ligand binding to a receptor; on gives the rate at which ligands bind to a receptor, e the rate of endocytosis of bound ligand-receptor complexes and off the rate at which ligands dissociate from receptors. Here, eff is the effective binding rate of ligands to cells (defined in Table 2). As eff depends on several parameters, including the number of receptors per cell, the cumulative binding probability changes with R tot and other system parameters ( Figure 1d).
Building on this work, we assume that cells have a fixed number of receptors for detecting ligands, senescent cells secrete ligands at a fixed rate (N E ) and that a cell must bind a certain number of ligands (N D ) to become senescent. We calculate the rate that ligands bind to a cell at distance r from the emitting cell and the rate that senescence is induced in cells at distance r from a senescent cell.
Throughout this paper, we assume that juxtacrine secondary senescent cells do not produce a SASP and cannot, in turn, induce juxtacrine senescence in their neighbours as a blunted SASP response in juxtacrine secondary cells has been reported (Hoare et al., 2016;Teo et al., 2019). We explore different scenarios of senescence spread, such as secondary or tertiary induction of juxtacrine senescence from paracrine senescent cells, as summarised in Figure 1g.
We assume that juxtacrine-mediated senescence induction occurs when cells are in contact, that is, when they are in close spatial proximity.
In a minimal model, we calculate the probability of senescence spread from a single senescent cell, by considering the probability of inducing senescence in an annulus of non-senescent cells around the primary senescent cell, as shown in Figure 1e. By considering the probability of senescence spread to a series of annuli, each with the same width, but an increasing radius, we determine the likelihood that senescence spreads from a single cell.
To model more realistic experimental scenarios we built a stochastic simulation of senescence spreading from a collection of cells over time. We use this to investigate the role of secondary senescent cells, by varying parameters such as their SASP secretion, time delays involved in senescence induction and lesion size.
Detailed descriptions of the minimal model and stochastic simulations can be found in Methods (Section 4).

| Model parameters
Our model of the spread of senescence contains several parameters, such as the constants relating to ligand-receptor interactions, where exact values have not been determined experimentally.
Where this was the case we used values from closely related biological systems, for example, the parameters for EGF and its receptor suggested in Batsilas et al. (2003). These parameters are shown in Table 1.
We assume that a non-senescent cell needs to bind a threshold number of ligands (N D ) in a given time frame for senescence to be induced ( Figure 1f) and that a single senescent cell produces N E ligands per unit time. We calculated the number of SASP proteins produced by a senescent cell per hour ( Figure S7) based on a recent paper by Schafer et al. (2020) that measured the concentrations of SASP proteins produced by radiation-induced senescent cells. For fibroblasts, the data suggest ~275,000 protein molecules per hour across all SASP proteins investigated. However, if we limit the proteins of interest to Interleukin 6 (IL-6), Interleukin 8 (IL-8) and Activin A (secreted proteins thought to be responsible for the induction of senescence; Acosta et al., 2008;Ito et al., 2017;Jochems et al., 2021;Prata et al., 2018;Putavet & de Keizer, 2021), the number produced per hour is much smaller, ~2887.
For simplicity, we did not distinguish between the types of ligands produced. Although realistically, a specific ligand will bind to its complementary receptor, and multiple signalling molecules may be required to induce senescence, we did not include this complexity in the model. Table showing the estimated parameters for the spread of EGF, as described in Batsilas et al. (2003), originally from Wiley et al. (2003). We found that there is both an optimum number of receptors (R tot ) and an optimum forward binding rate ( on ) for the spread of a signal (Figure 2a-e) in the regime where both e and off (the endocytosis and dissociation rates) are negligible, such that ligands which bind to a receptor remain bound. Intuitively, it may seem that more receptors would increase the signal spread, but this is not the case. When ligands do not dissociate from a cell easily, once they are bound they are removed from the system. Therefore, in regimes with large numbers of receptors, ligands are more likely to bind to the cell that emitted them (known as autocrine signalling) or to bind in spatial proximity, so the signal does not spread. Conversely, in regimes with a small number of receptors, the signal diffuses to a much greater distance from the senescent cell, but it is unlikely that enough ligands will bind to a cell to signal effectively (Shvartsman et al., 2001). This is illustrated in Figure 2f.

TA B L E 1
The same logic applies if e > > off such that the probability or internalisation, , tends to 1, meaning that the rate of dissociation is relatively low. If the binding rate is large, many ligands will be sequestered in autocrine trajectories, spreading the signal ineffectively. However, if the binding rate is small it is unlikely that enough ligands will bind to a cell to induce a change. Figure 2e illustrates the effect of increasing the binding rate.
When ligands can dissociate effectively and quickly from cells, the optimal number of receptors (for the signal to spread) is increased ( Figure 2b; black vs. blue curves). In this regime, a single ligand can bind to multiple cells over the course of its lifetime. Thus, we assume only the last binding event is relevant for senescence induction. In this regime, senescence spreads further with a larger number of receptors and a greater binding rate.
A similar effect is observed when increasing the diffusion coefficient ( Figure 2c). With increased diffusivity a ligand can more easily move away from the cell that emitted it, there are fewer autocrine binding events, and so a larger number of receptors per cell can be accommodated without limiting the signal spread.
Compared with the binding rate and receptor number, the cell density has only a small impact on the system's dynamics, with a greater cell density resulting in a greater chance of senescence spread ( Figure 2d).

| Juxtacrine senescence cells limit the spread of senescence
To investigate the effect of juxtacrine senescent cells on the spread of senescence, we first considered a ring of juxtacrine cells around the primary senescent cell in our minimal model. We assume juxtacrine secondary cells do not emit a SASP ( Figure 1g)  (c), but with a 6-day time delay before a cell can induce senescence in its neighbours. (e) The same simulation as in (d), but with immediate induction of juxtacrine senescence and delayed induction of paracrine senescence.
manner. As juxtacrine senescent cells in our model do not produce a SASP, their presence is analogous to removing trees to prevent the spread of a forest fire.

| A time delay in the senescence induction process limits spread and alters the dynamics
After a growing cell binds the required number of ligands to induce senescence, it takes a number of days before the cell becomes fully senescent and it is capable of inducing senescence in its neighbours. To investigate this effect, we introduced a 6-day time delay between senescence induction and a cell becoming fully senescent in the stochastic simulation, as Young et al. (2009) suggest 6 days are required to establish full Ras-induced senescence ( Figure 3b).
We found that introducing a time delay drastically increases the time taken for senescence to spread, as expected ( Figure 3d and Figure S1b). An initial region of paracrine secondary senescent cells is created by the primary lesion, limited in diameter by the reach of the SASP from the primary senescent cells (Figure 3d). Once these cells become fully senescent and produce their own SASP, the signal spreads further, inducing a new wave of paracrine secondary senescent cells.
The time delay also alters the density of juxtacrine secondary senescent cells, reducing their number at larger distances from the primary lesion ( Figure S5c,d). This is due to the assumption that a cell must be fully senescent before it can induce either juxtacrine or paracrine senescence in its neighbours. The delay in juxtacrine induction increases the likelihood that neighbouring cells will become paracrine secondary senescent in the intervening time, increasing the influence of the SASP over senescence spread. The reduction in juxtacrine cells further from the senescent lesion means that as senescence spreads further, the rate of senescence spread increases due to the higher ratio of paracrine secondary senescent cells, which is inconsistent with experimental observations in vitro and in vivo (Acosta et al., 2013;Michaloglou et al., 2005).

| A dynamic SASP is required to avoid exponential increases in cell numbers at late times
To further investigate the effects of a dynamic SASP we simulated the limiting case, in which senescent cells can induce juxtacrine senescence in their neighbours as soon as they begin the journey towards becoming fully senescent. In reality, evidence suggests that juxtacrine induction may actually be possible after a 2-day delay (JAG1, the NOTCH ligand responsible for juxtacrine senescence induction, is upregulated even after 1 or 2 days [Hoare et al., 2016]

| Reduced SASP from secondary senescent cells can lead to controlled senescence spread
We wanted to explore parameter regimes in which the spread of senescence is controlled (local over a human lifetime), but this was difficult with such a large number of unknown system parameters.
We have previously outlined how senescence spread can be slowed with the presence of juxtacrine senescence, a time delay, and juxtacrine signalling, which precedes paracrine signalling. If the spread of senescence can be slowed drastically enough, this will lead to controlled senescence spread, as even over a human lifespan the senescence spread will be local and contained. senescence. One way to achieve this is to ensure that the number of juxtacrine secondary senescent cells is large enough to dilute the SASP. This is also achieved if paracrine secondary senescent cells produce fewer SASP molecules than primary senescent cells, as suggested by experimental evidence (Acosta et al., 2013). In our simulations, we find that in this case, the primary cells in the initial lesion drive most of the senescence spread ( Figure S2a). Some secondary senescent cells are generated over short time scales, but the spread eventually slows down to the point where it remains controlled, in some cases even after 500,000 hours (on the order of a human lifespan). When the lesion size is increased the size of the region of secondary cells will also increase ( Figure S2b).

| Parameters of paracrine senescence spread can be determined through inference
To estimate the parameters of the model, we propose a similar experiment to those described in Acosta et al. (2013), seeding a circle of primary senescent cells and watching the phenotype spread (Section 4.3).
By focusing on the first wave of senescence induction, we can forgo computationally intensive simulations. Due to the time delay before a cell becomes fully senescent, there is a period of time when a cell is unable to pass the senescent phenotype on to neighbouring cells through paracrine signalling. If it takes cells 6 days to become fully senescent, we can stop the experiment before 6 days. Using a marker to detect cells that are becoming senescent, we can determine the region to which senescence has spread through paracrine signalling (Figure 4a).
Focusing on short timescales allows us to build on the minimal model described in Section 4.1 to predict the fraction of senescent cells at a given distance from the primary lesion. By comparing this with the experimentally observed fraction of senescent cells we can determine N E , N D , R tot and k on using Bayesian inference, fixing other model parameters to experimentally determined values (Section 4.3). As a proof of principle for this experiment, we used our model to generate synthetic data and attempted to infer the parameters used to generate the data. As a prior distribution of the parameters, we chose a uniform distribution over the parameter regime given in Table 1 (Figure 4b).
We show that it is possible to infer some of the parameters accurately from such an experiment (Figure 4b). The posterior distribution for R tot is narrow and centred on the true value, and although the posterior for N D is not as narrow it is well centred on the true value. This reflects the intuition we gained from the minimal model, that senescence spread is sensitive to N D and R tot . The posterior distributions for N E and on are much broader and less centred on the true values.

| DISCUSS ION
Currently, it is unclear how senescence can spread in a controlled way from a primary lesion. We systematically investigated senescence spread, using mathematical modelling informed by experimental evidence. We found that the number of receptors per cell and the binding rate of ligands were critical in determining the distance over which senescence spreads. We observed that an intermediate number propose an experiment aimed at determining four critical parameters involved in senescence spread, and we show that Bayesian inference is a viable tool for estimating these parameters from experimental data.
Although determined from a 2D experiment, several of the inferred parameters (N D , N E , R tot , on ) also describe paracrine senescence spread in 3D. We believe our model of senescence spread is well-justified to study the in vitro scenario and of an appropriate complexity to identify which assumptions are important or affect senescence spread. Once parameterized with future in vitro experiments, the model can straightforwardly be expanded to 3D, to describe senescence spread in tissues, in vivo.
While we can quantify parameters associated with the paracrine spread of senescence from a lesion, determining parameters associated with juxtacrine secondary senescent cells is more challenging.
Single-cell transcriptomic data may help to determine the difference between the senescent cell types. However, it is likely that differences between senescent subtypes will vary depending on the cell type and method of senescence induction, suggesting that senescence spread might differ depending on the location in the body and the cause of senescence.
By considering a minimal model of senescence spread we have extracted an understanding of the dynamics of this system. can bind to receptors on the surface of cells, but they reflect from the space between cells (the dish). Thus, they approximated the cellcovered dish as a semi-reflecting boundary layer (which is less accurate between 1 < r ∕ r cell < 2, close to the emitting cell) and simulated the diffusion and binding of a large number of emitted ligands.
We recreated the simulations described by Batsilas et al. and verified that Equation (1) predicts the binding location of ligands with a relative error <5%. For paracrine trajectories, we compared the time at which a simulated ligand first binds to a new cell to the square of the mean distance at which it binds ( Figure S4b,c), showing that the limiting timescale in this process is diffusion (Müller & Schier, 2011). The shape of the cumulative binding probability (P(r)) is shown in Figure 1d, for both small and large numbers of receptors per cell, R tot .
Equation (1) describes paracrine trajectories, not including ligands that bind to the cell that emitted them (autocrine trajectories).
The probability of an autocrine trajectory is given as, where Da is defined in Table 2 (Batsilas et al., 2003).
Equations (1) and (2)  where = e ∕ e + off . The first binding event and internalisation are the same if off = 0 and approach each other when the rate of endocytosis is large and the rate of dissociation is small.
Between the first binding event and internalisation a ligand could bind to many cells. However, this binding may be only for a short period of time if the rate of dissociation is high. In this paper, we consider only the location of the final binding of the ligand to be important for senescence induction. This assumption breaks down if even fleeting binding of ligands to a cell produces a signal that persists intracellulary, and when combined with the signal from many other fleeting ligand binding events also induces senescence. A comparison of the binding and internalisation distances for an example set of parameters is given in Figure S6.
The diffusion of ligands is clearly affected by the height of the media layer, h. However, to sustain cells in vitro, we would require a layer of media with h > 2 mm. Using dimensional analysis and biologically reasonable system parameters, Batsilas et al. (2003) argue that the characteristic trajectory length, D L ∕ , is much smaller than 2 mm, and so h does not need to be accounted for.

| Senescence spread from a single cell to an annulus
Extending the work of Batsilas et al., we find that for each ligand produced, the probability of a ligand binding in an annulus of width 2r cell (illustrated in Figure 1e) at distance r from the emitting cell is, where a = 1.1D L ∕ K eff . The area of this annulus is given by 4 rr cell .
To find the probability that a cell in this annulus will bind the ligand, we multiply this probability by the area of one cell relative to the annulus, to get, If one cell emits N E ligands per unit of time, the expected number of ligands bound to a cell at a distance r, per unit of time, is, assuming that diffusion has reached a steady state, and N E is constant so that the rate of ligands reaching a cell at distance r is also constant.
(2) P au = Da Da + 4 ∕ , (3) P in (r) ≈ r r + 1.1D L ∕ eff , and P in au =  This is justified if the diffusion timescale is much faster than the other timescales in the system.

Da
To find the probability per unit time that a specific cell in an annulus will become senescent in any given second, cell (r), we assume that a cell must bind to a threshold number (N D ) of ligands per unit of time to become senescent. Here, we use the binomial distribution to describe this process, where N E becomes the number of trials, and Equation (5) gives the probability of success.
However, there is more than one cell in an annulus, on average there will be m = 4r r cell . Therefore, spread (r) = m cell (r) is the rate of senescence induction in each annulus.
Following this analysis through, we find the probability of creating at least one senescent cell at distance r from the primary lesion in time t is, which is 1, minus the probability (given by a Poisson process) that no cells become senescent.

| Total senescence spread from a single cell
We now have the probability that one cell at distance r from the primary lesion will become senescent in a time period. We use the Poisson-Binomial distribution to calculate the probability that a senescent cell is created in at least one annulus, as the probability of senescence induction in each annulus is different. Therefore, the probability of senescence induction over all annuli is given by a set of independent Bernoulli trials that are not identically distributed, where each annulus is a Bernoulli trial with p = P r (T ≤ t) , and we calculate the probability that one or more of these trials is successful.
We use this simplified model of the system to determine whether the spread of senescence is uncontrollable, considering the limiting case of a single primary senescent cell. If a single senescent cell spreads senescence to one other cell before the body can remove it (we approximate that this happens in 2 days), with a probability of ~1, we define the spread of senescence as uncontrollable. This uncontrollable spread of senescence is not observed in vitro; therefore, this is an unrealistic parameter regime.

| Senescence spread with juxtacrine cells
We considered two different regimes, with and without juxtacrine secondary senescence. Where juxtacrine signalling induces secondary senescence, we expect that cells surrounding the primary senescent cell will become juxtacrine senescent through contact. Assuming that juxtacrine secondary senescent cells do not emit a SASP, for the senescence to spread, a cell beyond this ring of juxtacrine cells must become senescent within 2 days. If juxtacrine cells instead produce a blunted SASP, this will lead to an underestimate of the spread of senescence .
The distance of the annulus from the primary senescent cell which is considered in each case is given by, The first value, r 1 , is simply the average distance between cells, and r 2 = 2r 1 , accounts for a ring of juxtacrine senescent cells.

| Limitations
This minimal mathematical model provides a tractable way to consider senescence spread from a single cell but is difficult to extend to a system containing many cells. Furthermore, this model cannot be used to investigate the effect of a time delay in senescence induction, and it makes a number of assumptions which may limit the accuracy of the predictions.

| Stochastic simulations
To deepen our understanding of the system and more closely match the reality of senescence spread in vivo, we created a stochastic model for the spread of senescence from multiple initial cells and over many subsequent senescence induction events. The stochastic model lets us introduce the added complexity of paracrine senescent cells which can in turn emit a SASP. We used the simulation to explore the change in dynamics resulting from changing the strength of the SASP produced by these secondary senescent cells. Finally, we introduce a time delay between senescence induction and SASP, more accurately modelling the cellular biology. Cells do not immediately start emitting a SASP once they have bound enough ligands to induce senescence, and instead, take several days to become fully senescent.

| Simulation implementation
We consider three possible cell states: senescent, in the process of becoming senescent, and non-senescent. A cell starts to become paracrine secondary senescent when it binds to more than N D ligands per unit of time (see Section 2.1.1). The simulation is initialised on a square plane, and a 100 × 100 cell radii area is commonly used as larger areas increase the computational time taken for each simulation. Cells are then created so that they randomly populate this area with a specified density (Figure 3a). Primary senescent cells are seeded at time 0. The simulation proceeds by tracking the state of each cell.
In the stochastic simulation, we approximate the number of ligands bound from each emitting senescent cell using the Poisson distribution. This approximation holds when the number of ligands emitted, N E , is large and the probability of binding (equation (5) This approximation dictates the appropriate time frame for the simulation. Considering paracrine signalling using only IL-6, IL-8 and Activin A molecules, we can estimate N E ≈ 2887h −1 .
We can now calculate the probability of a cell becoming senescent in one unit of time from its position relative to the SASPemitting cells. This probability can then be used to give the rate of senescence induction for that cell. Thus, each cell on the grid which is not yet senescent has a rate of senescence induction. By summing these rates we calculate the time at which the next cell becomes paracrine secondary senescent. We use the Gillespie algorithm, meaning that we first determine the time an event occurs based on creation rates, for example, a cell becomes paracrine secondary senescent at t = 1.5 h. We then determined which cell becomes senescent, again based on the rate of senescence induction associated with each cell (Gillespie, 1977).
The situation is complicated further by the time delay: It may be that cells have finished becoming fully senescent before the next paracrine induction event, in which case the SASP that they produce must also be accounted for. This is achieved simply by tracking the number of cells that are in the process of becoming senescent. If a cell finishes becoming senescent before the next reaction (i.e., new senescence induction event) occurs, then the cell is added to the list of fully senescent cells, and the simulation run again from that time point, with the probabilities of each cell becoming senescent re-calculated including any SASP from the newest senescent cell.
This works as the system is Markovian.
There are a number of algorithms that can be used to speed up the delayed reaction calculation (Cai, 2007). However, at any one time, we do not have a large number of cells waiting to become senescent and the bottleneck in the process is the repeated calculation of the probability of inducing senescence in each cell.
The simulation code can be found on GitHub (https://github. com/lkmar tin90/ Senes cence_Spread), and a flowchart of the stochastic simulation implementation is shown in Figure S10.

| Comparison of stochastic simulations to minimal model
We compared our stochastic simulation to our minimal model to verify that the two agree. To do this, we considered the spread of senescence from a single cell, in the case where no juxtacrine secondary senescent cells are created.
We used the stochastic simulation to simulate this scenario 100 times (fixing all parameters but initialising a new population of cells each time) and recorded the proportion of the simulations in which there was no senescence spread. We varied N D and R tot , to ensure that the simulations agree over a range of parameters. We compared the results with our minimal model, which gives us the probability of senescence spread.
From this, we see a good agreement between the two methods at both R tot = 10 5 and R tot = 10 6 , both when cells in the stochastic simulation are created so that they lie on a lattice ( Figure S3a,b), and when cells are randomly distributed ( Figure S3c,d). There are some differences in the regime without juxtacrine senescence induction with R tot = 10 6 , due to the large number of ligands which bind to nearby cells and the assumption we make when constructing the minimal model that cells are uniformly distributed ( Figure S3b, dotted).
In the stochastic simulation, with cells distributed randomly ( Figure S3c), the range of parameters that result in a probability of senescence spread >0 and <1is greater. If a nearby cell happens to be closer than the average distance between cells then there is a higher probability that senescence will spread. This implies that the parameter regime of large R tot is highly dependent on the initialisation of the simulation, in particular, the precise location of the cells. Although the average density in each simulation is the same, the signal spread is possible in some cases and not others. If there are no cells close enough to the primary senescent cell, senescence will not spread.
We investigated whether the optimum number of receptors that we found for senescence spread in the minimal model, as shown in Figure 2, also appears in the stochastic simulation, or whether this is a consequence of the simplifying assumptions made in the minimal model. We see that the random nature of the cell placement in the stochastic simulation means that there is no longer a clear optimal number of receptors to ensure the spread of senescence, however, the general dynamics remain the same, with juxtacrine cells reducing the spread of senescence and the number of receptors altering the distance over which cells communicate ( Figure S3d).

| Limitations
The stochastic simulation is a useful tool to understand senescence spread from a small lesion, but computationally intensive to scale to larger spatial regions. These limitations, as well as the relationships between the models discussed in this paper, are shown in Figure S9.

| Parameter inference
We use Bayesian statistics to infer system parameters from a realistic in vitro experiment, where a circle of primary senescent cells is seeded in the centre of a dish, as described in Section 2.7. We calculate the distance senescence spreads from the circular lesion.
The number of SASP molecules emitted per unit area by a primary senescent lesion is, , and the probability that a single test cell ( Figure 4a) at distance x from a senescent cell binds a ligand emitted by a senescent cell is Integrating for the emission from the area of seeded cells to get the total expected number of ligands binding to the test cell gives Writing x in terms of r, r ′ and , according to the geometry in Figure 4a, gives Here, r seed is the radius of the circle of primary senescent cells. We solve this integral numerically.
Using the method previously described in Section 4.1, we calculate the probability that the cell that binds these ligands becomes senescent, using a Poisson distribution to determine the probability that more than N D ligands are bound per hour. We then calculate the probability that a cell becomes senescent within the duration of the experiment. Finally, we calculate the fraction of cells, which are senescent at a given distance from the seeded senescent region of cells. As this process is probabilistic, the fraction of senescent cells varies on each run of the model. Before comparing the model output with an experiment, we ought to consider whether the assumptions made in the model's construction are reasonable. In particular, in the model, we have not accounted for the time taken for ligands to diffuse away from the cell that emitted them. From Figure S4, we can see that the overall process of signal spread is diffusion-limited. Previous experiments show senescence spread over a distance of the order of 1 mm, indicating that this is a meaningful distance to consider. The time taken for ligands to travel 1 mm is between 50 and 100 seconds for the parameter regimes shown, depending on the diffusion coefficient of the ligand. As long as the time taken for diffusion is small compared with the timescales associated with senescence spread (hours), we can assume that the diffusion occurs instantaneously.
To determine the parameters that lead to the closest agreement between the model and the experimental data, we compare two summary statistics. First, the distance where the fraction of senescent cells crosses 0.5 (S). Second, the fraction of senescent cells at different distances from the lesion ( � ⃗ F). The distances need not be equally spaced, and we use more data points around the region of r = S, where the fraction of senescent cells is <1 and >0. We verified that these summary statistics are sensitive to changes in the parameters ( Figure S8a).
To quantify the agreement between predicted and observed summary statistics, use the distance function: where d 1 = S model − S data and d 2 = ∥ � ⃗ F model − � ⃗ F data ∥.
We use pyABC (Klinger et al., 2018), a likelihood-free inference software using a sequential Monte Carlo scheme, to estimate posterior distributions of model parameters. In this paper, we use synthetic data, created from the same model we use for the inference.
Our inference ran on 60 cores for 13 populations, with 1000 particles, taking 7.5 h in total. While acceptance thresholds for d (epsilon values) have not converged at this stage, indicating that further iterations may improve the inference, the fraction of accepted samples was ≈ 10 −3 , making it computationally wasteful to run for more populations ( Figure S8b). Unit. We thank M. Nicholson for the critical reading of the manuscript. We thank all members of the Chandra lab and Schumacher lab for their input. We thank Acosta et al. (2013) andSchafer et al. (2020) for sharing additional data related to their publications on request.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Code used in this manuscript is available at https://github.com/ lkmar tin90/ Senes cence_Spread.