First-passage times and normal tissue complication probabilities in the limit of large populations

The time of a stochastic process first passing through a boundary is important to many diverse applications. However, we can rarely compute the analytical distribution of these first-passage times. We develop an approximation to the first and second moments of a general first-passage time problem in the limit of large, but finite, populations using Kramers–Moyal expansion techniques. We demonstrate these results by application to a stochastic birth-death model for a population of cells in order to develop several approximations to the normal tissue complication probability (NTCP): a problem arising in the radiation treatment of cancers. We specifically allow for interaction between cells, via a nonlinear logistic growth model, and our approximations capture the effects of intrinsic noise on NTCP. We consider examples of NTCP in both a simple model of normal cells and in a model of normal and damaged cells. Our analytical approximation of NTCP could help optimise radiotherapy planning, for example by estimating the probability of complication-free tumour under different treatment protocols.

Here, γ is the rate at which radiation-damaged cells repair their DNA. The fractional change in the population over an infinitesimal timeψ(t)/ψ(t) gives the hazard function h(t). This is found to be given by 5

Evaluation of Approximation 1 for Model 2
In Model 2 we write N t for the number of normal cells at time t and X t for the number of doomed cells. We are interested in the population of total functional cells, S t ≡ N t + X t . Specifically, we are interested in the time S t first passes a boundary L. The master equation can be formulated in terms of S and X: (SI. 4) where P S,X (t) is the probability that random processes S t , X t have the values S, X at time t. The operator E S is the step operator affecting the size of the total population, and E X is the step operator affecting the number of doomed cells, i.e. E S [ f S,X ] = f S+1,X and E X [ f S,X ] = f S,X+1 . We proceed by approximating the master equation via a Kramers-Moyal expansion. First, we introduce re-scaled processes s t = S t /M and x t = X t /M, and then expand the step operators in the limit M 1. We arrive at the Fokker-Planck equation where we have neglected higher-order terms in M −1 . This Fokker-Planck equation can equivalently be written as an SDE: where the drift is given by The diffusion B(s, x,t) is the positive-semidefinite matrix satisfying We proceed by linearising the SDE (SI.6). Let where φ 1 (t) and φ 2 (t) are the deterministic functions of time. Substituting and collecting lowest order terms, we see these functions are given by the ODEs i.e., we recover Eqs. (28b). The random processes ξ 1 t and ξ 2 t describe deviations from this deterministic trajectory, and are of the Ornstein-Uhlenbeck type (SI.10) where A(φ 1 , φ 2 ,t) is given by (SI.11) We note that the argument of B in Eq. (SI.10) is now given by φ 1 and φ 2 , so that the noise is additive rather than multiplicative. We are interested in the variation of the total population size from the deterministic path ξ 1 2 t ; we remark that by construction ξ 1 t = ξ 2 t = 0. The variances and covariance of ξ 1 t and ξ 2 t can be seen to evolve in time as follows 6 d ξ 1 For a given set of parameters, we numerically integrate the five coupled Eqs. (SI.9) and Eqs. (SI.12). This provides the mean and covariance matrix for the bivariate Gaussian distribution of the number of normal and doomed cells as a function of time. For Approximation 1, the time t * is defined by φ 1 (t * ) = ; this is the point in time when the total number of functional cells crosses the threshold for onset of NTC. The variance of the number of functional cells at this time is given by Σ 2 (t * ) = ξ 1

A comparison of Approximations 1 and 2
Fig. SI.1 shows a quantitative comparison of the distributions of first-passage time from the two approximations for Model 1. It indicates that Approximation 2 provides an improvement relative to Approximation 1. Both methods do considerably better than the deterministic approximation in Stocks et al. 5 .
To compare the three approximations we have used the Earth Mover's Distance (EMD), also known as the Kantorovich metric 7 , as a measure of distance between two probability distributions. Intuitively, it is a measure of the amount of 'effort' needed to turn one distribution into the other; it is the amount of probability that needs to be moved weighted by the distance it has to be moved. We choose this rather than, say, the Kullback-Leibler divergence 8 or total variation distance since the distribution of first-passage times from the deterministic approach is a Dirac delta-distribution 5 which results in infinite Kullback-Leibler divergence. The EMD gives a more useful measure of error.

Complication-free tumour control
The objective of radiation therapy is to successfully eliminate cancerous cells while avoiding further complications from damaging normal tissue cells. We have developed analytical approximations for the efficient calculation of NTCPs. Tumour control probabilities-the probability of eliminating all cancer cells-from a stochastic birth-death model have been previously considered by Zaider and Minerbo 9 ; the authors derive a general equation for the probability of the elimination of all tumour cells. Here, we combine these two results for NTCP and TCP to investigate how, in principle, mathematical models can be used to optimise the application of radiation therapy to achieve complication-free tumour control.

Model definitions
We consider a model, which we call Model 3, that contains both normal cells N and cancerous cells C . The two populations are assumed to be spatially separated from each other. The normal cells are as described in Model 1. Cancerous cells, on the other hand, undergo mitosis with a constant rate b 2 9 ; numerical evidence suggests that the resulting exponential growth characterise tumours of small sizes well 10 . Cancerous cells are also subject to a natural death with a rate d 2 and to death from a 3/6 source of radiation with hazard function h 2 (t). The model can be summarised by the following reactions: (SI.13) Although both types of cell are subject to the same source of radiation, the hazard functions h 1 (t) and h 2 (t) can differ. This is because each cell type differs in its susceptibility to radiation and in their ability to repair damaged DNA. We again consider the case of brachytherapy; the hazard function is as in Eq. (29), where the parameters α 1,2 , β 1,2 , and γ 1,2 depend on the cell type. We also assume that, due to the presumed spatial separation of normal tissue and cancerous cells, the treatment can be targeted such that each cell type absorbs a different fraction of the total dose rate. This is incorporated into the hazard function by replacing the initial dose rate R 0 with an effective dose rate θ 1,2 R 0 . The parameters describing the initial dose rate R 0 and the decay rate λ are characteristics of the radioactive implant and are thus common to the hazard function of both cell types. As before, we initialise the population of normal cells in its stationary state. We let there be initially C 0 cancer cells.

Tumour control probability, normal-tissue complication probability, and probability of complication-free control
The probability as a function of time of eliminating all cancer cells, TCP(t), is mathematically also a first-passage time problem. Zaider and Minerbo 9 developed an analytical description for TCP for the linear dynamics of cancerous cells described above. This was achieved using a generating-function. This approach is feasible due to two features of the problem: (i) the model is linear (i.e., cells do not interact with each other), and (ii) the boundary of interest for TCP is at zero (i.e., extinction of tumour cells). Their result for TCP(t) is 9 where C(t) is the deterministic path for the expected number of cancerous cells, given by (SI.15) While Eq. (SI.15) cannot be solved analytically in most cases, it can be integrated numerically for an efficient calculation of TCP(t). Complication-free tumour control (CFC) refers to the elimination of all cancer cells while maintaining enough normally functioning tissue for an organ to operate without complications 11 . The probability of CFC as a function of time is therefore given by 11 (SI.16) We remark that Eq. (SI.16) implies an equal weighting of the importance of tumour control and NTCs. In the most extreme cases, for example where NTCs relate to organ failure this is justified. In other cases, for example when NTC refers to increased urinal frequency, a complication may be preferable to a potentially life-threatening tumour. In such cases, Eq. (SI.16) can be modified by appropriately weighting the two probabilities to maximise a 'quality of life' measure in accordance with clinical experience 12 . Fig. SI.2 (a) shows the probabilities 1 − NTCP(t) and TCP(t) for Model 3 and for a specific choice of parameters (see Table SI.1). These quantities are obtained by Approximation 2 for NTCP, and Eq. (SI.14) for TCP. Similarly, Fig. SI.2 (b) shows CFC(t) and compares the results from our approximation to those of numerical integration of the master equation. For this choice of parameters we find a non-trivial time (∼ 20 days) which maximises the probability of CFC. In the case of a temporary brachytherapy implant, this would indicate the optimum moment for removal.
In order to achieve CFC, we assume we are able to control the initial dose rate R 0 (i.e., the size of the radioactive seed) and the time at which the implant is removed. Fig. SI.2 (c) shows the probability of CFC for different values of time and initial dose, again efficiently generated using Approximation 2 for NTCP and Eq. (SI.14) for TCP. With the exception of the population sizes, the parameters we choose here were previously used to model the treatment of prostate carcinoma 5 consistent with experimentally collected parameters 13 . In this context NTC refers to acute radiation proctitis 14 . For these parameters, the optimal strategy involves an initial dose of size 1.7 G y and removal at a time over 50 days. Using this initial dose, the probability of CFC(t) does not decrease at large times, providing a large window for the removal of the implant or allowing the use of a permanent implant. This is not the case for all parameters; the optimum strategy may require the timely removal of the implant. An example of this is shown in Fig. SI.2 (d), which shows CFC(t) for parameters where the cancer cells have a three-fold higher growth rate than normal cells. The probability of CFC is peaked when implanting a high dose of radiation for a short time. For this case, we see the band where CFC is likely is narrow, indicating that such a treatment may be very sensitive to the time of removal of the implant.   Fig. SI.2, along with λ = 0.0117 day −1 . The parameters in the upper two rows were previously used to model brachytherapy as a treatment for prostate cancer, where the normal tissue complication refers to rectal proctitis 5 . The parameters in the bottom row are hypothetical, used to show that a change in the optimum treatment strategy may result upon variation of parameters.