First passage events in biological systems with non-exponential inter-event times

It is often possible to model the dynamics of biological systems as a series of discrete transitions between a finite set of observable states (or compartments). When the residence times in each state, or inter-event times more generally, are exponentially distributed, then one can write a set of ordinary differential equations, which accurately describe the evolution of mean quantities. Non-exponential inter-event times can also be experimentally observed, but are more difficult to analyse mathematically. In this paper, we focus on the computation of first passage events and their probabilities in biological systems with non-exponential inter-event times. We show, with three case studies from Molecular Immunology, Virology and Epidemiology, that significant errors are introduced when drawing conclusions based on the assumption that inter-event times are exponentially distributed. Our approach allows these errors to be avoided with the use of phase-type distributions that approximate arbitrarily distributed inter-event times.


Competing events
Let us consider that the process X just arrived into state i ∈ S at the current time. When analysing the next event to occur, we focus on the conditioned times T i→ j |i → j and the probabilities of these events taking place, p i j = P(i → j). In general, the probability density of these conditioned times is given by where P(T i→ j < t|i → j) = P(T i→ j ≤ min(t, T i→k : ∀k = j)) P(i → j) and where f T i→ j (t) is the absolute waiting time probability density; that is, the density of T i→ j , and F T i→k (t) represents the absolute cumulative waiting time density of T i→k . The latter integral allows one to combine information about every transition in the absence of the others into a conditioned probability that the event i → j actually occurs before any other transition. Eq. (1), Eq. (2) and Eq. (3) do not seem very useful when the times T i→ j are generally distributed as in Fig. 3b left. However, when all but one absolute waiting times are exponentially or phase-type distributed, as in Fig. 3a or Fig. 3b right, Eq. (1), Eq. (2) and Eq. (3) simplify and allow actual computations. Thus, in the following sections we describe how to apply these equations when, from any state i ∈ S , there are n i directly accessible states j ∈ A(i) = {i 1 , i 2 , . . . , i n i −1 } ∪ {i n i } with corresponding absolute waiting times T i→ j as in Fig. 3a or Fig. 3b right in the main text.

Probability of an exponential transition taking place before a general one
For the sake of simplicity, let us first consider the case where there are only two possible transitions out of state i (n i = 2 arrows leaving i in Fig. 3): transition to a state i 1 , occurring after an exponentially distributed random time T i→i 1 ≡ T Exp(λ i,i 1 ) ∼ Exp(λ i,i 1 ), and to a state i 2 occurring after a generally distributed random time T i→i 2 ≡ T G ∼ G. In this case Eq. (2) above can be written as and Eq. (3) as Eq. (5) is the Laplace-Stieltjes transform of 1 − F T G (t) with argument z = λ i,i 1 and the Laplace transform of the time derivative of Eq. (4) is also the Laplace-Stieltjes transform of 1 − F T G (t) with argument z + λ i,i 1 . Finally, we get and where For the general case in which we have n i − 1 exponentially distributed and one generally distributed waiting times, as in Fig. 3a, it is possible to obtain the following analogous expressions

Probability of a general transition taking place before an exponential one
In a similar fashion, and when considering n i = 2 in Fig. 3a, we can analyse the probability of i → i 2 taking place before i → i 1 .
In particular, it is easy to see that and We note that we can also work in backward fashion. That is, making use of the properties of functional equations 1 , we can invert Eq. (9) to obtain the absolute waiting time distribution as a function of the conditioned one. In particular, we find This result is striking as it allows one to obtain information about the absolute waiting time distribution of an event from the observation of the distribution conditioned that it has been experimentally observed. In fact, by the properties of the Laplace-Stieltjes transform, Eq. (10) can be transformed back to the time domain, leading to the density function Finally, similar formulae can be obtained for a general number n i of states in Fig. 3a. In particular, we get Eq. (2) in the main text, namely

Generalisation to the case of one phase-type waiting time and a general one
The arguments in the previous sections relate to a generally distributed waiting time competing with several exponentially distributed waiting times. It is possible to generalise these arguments for phase-type distributions instead of exponential ones. This can be extremely useful, given that the family of phase-type (PH) distributions is dense in the family of non-negative continuous distributions, and there exist algorithms to approximate any generally distributed random variable by a phase-type random variable 2 . This means that in Fig. 3b left, when trying to consider more than one generally distributed waiting time, we could replace all but one generally distributed waiting times by phase-type distributions, as in Fig. 3b right, and then our approach would also apply.
A phase-type distribution is defined as the time to absorption for a continuous time Markov chain with infinitesimal generator of the form and it can be seen as a generalisation of the exponential distribution (see Chapter 1 in 2 ). The density function of a phase-type distribution PH(α, T) is given by where T 0 ≡ −T · e, e is a column vector of 1's, and α is the vector of initial probabilities. Let us consider, without loss of generality, n i = 2 in Fig. 3b right: from state i the system moves to state i 1 after a time T i→i 1 ≡ T PH that follows a phase-type distribution, and moves to state i 2 after a generally distributed time T i→i 2 ≡ T G . We note that our arguments can be easily generalised to larger values of n i , where n i − 1 phase-type distributed times compete against a generally distributed time T G . The distribution function of this phase-type distribution can be expressed as 3 where N is the number of eigenvalues of T, p i (t) is a polynomial of order L i − 1, and L i is the multiplicity of the eigenvalue −µ i . It is clear that the survival function and the density function of this phase-type random variable are given by where q i (t) = for all i ∈ {1, . . . , N}. One can compute first the conditioned density function so that the conditioned transform is given by The probability P(T G < T PH ) can be obtained as It is clear that , and that Then, the conditioned Laplace-Sieltjes transform can be obtained by following similar arguments than above, and is given by where V (z) is the Laplace transform of the survival function of the general distribution, We finally note that, if the phase-type distribution under analysis corresponds to a diagonalisable matrix T, which is usually the case when the PhaseType R package 4 is used for approximating a general distribution (see case study 2 and the corresponding phase-type distribution in Fig. 9 in the main text), and the density function can be expressed as In this way, expressions above simplify as follows