Failure of complex systems, cascading disasters, and the onset of disease

Complex systems can fail through different routes, often progressing through a series of (rate-limiting) steps and modified by environmental exposures. The onset of disease, cancer in particular, is no different. A simple but very general mathematical framework is described for studying the failure of complex systems, or equivalently, the onset of disease. It includes the Armitage-Doll multi-stage cancer model as a particular case, and has potential to provide new insights into how diseases arise and progress. A method described by E.T. Jaynes is developed to provide an analytical solution for the models, and highlights connections between the convolution of Laplace transforms, sums of random samples, and Schwinger/Feynmann parameterisations. Examples include: exact solutions to the Armitage-Doll model, the sum of Gamma-distributed variables with integer-valued shape parameters, a clonal-growth cancer model, and a model for cascading disasters. The approach is sufficiently general to be used in many contexts, such as engineering, project management, disease progression, and disaster risk for example, allowing the estimation of failure rates in complex systems and projects. The intended result is a mathematical toolkit for the study of failure rates in complex systems and the onset of disease, cancer in particular.


I. INTRODUCTION
Complex systems such as a car can fail through many different routes, often requiring a sequence or combination of events for a component to fail. The same can be true for human disease, cancer in particular [1,2]. For example, cancer can arise through a sequence of steps such as genetic mutations, each of which must occur prior to cancer [3][4][5][6][7]. In addition, there is considerable genetic variation within a single cancer [8], suggesting that similar cancers might arise through a variety of different paths.
Here a simple model is described for the failure of a system through one or more possible routes, each of which consisting of one or more sequential or non-sequential steps. The model is easy to conceptualise and derive, and many specific examples have analytical solutions or approximations, making it ideally suited to the construction of biologically-or physicallymotivated models for the incidence of events such as diseases, disasters, or mechanical failures. A method described by E.T. Jaynes [9] generalises to give an exact analytical formula for the sums of random variables needed to evaluate the sequential model. This is evaluated for specific cases. Moolgavkar's exact solution [10] to the Armitage-Doll multistage cancer model is one example that is derived surprisingly easily, and is easily modified. More generally, the mathematical framework is intended to capture more complex cancer models, that might for example include a clonal expansion prior to cancer detection [4][5][6]. We also consider "cascading disasters" [11], where each disaster can substantially modify the risk of subsequent (possibly different) disasters.

II. FAILURE BY MULTIPLE POSSIBLE ROUTES
Imagine that we can enumerate all possible routes 1 to n by which a failure can occur ( Figure 1). The probability of surviving the ith of these routes after time t is S i (t), and consequently the probability of surviving all of these possible routes to failure S(t) is, or in terms of cumulative hazard functions with S i (t) = e −H i (t) , The system's hazard rate for failure by any of the routes is, and In words, if failure can occur by any of n possible routes, the overall hazard of failure equals the sum of the hazard of failure by all the individual routes.
A few notes on Eq. 2 and its application to cancer modelling. Firstly, if the sth route to failure is much more likely than the others, with H s H j for s = j, then Thirdly, most cancers are sufficiently rare that S ∼ 1. As a consequence, many cancer models (implicity or explicitly) assume S 1 − n s n i=1 H i (t) and f = −dS/dt n s n i=1 h i (t), a limit emphasised in the Appendix of Moolgavkar [10]. Often failure by a particular path will require more than one failure to occur independently. Consider firstly when there are m i steps to failure, and the order of failure is unimportant ( Figure 2). The probability of surviving failure by the ith route, S i (t) is, S i (t) = P (survive any one or more, necessary step for failure) where F ij (t) is the cumulative probability distribution for failure of the jth step on the ith route within time t. Writing S ij (t) = 1 − F ij (t), this can alternately be written as,

IV. RELATION TO RECENT MULTI-STAGE CANCER MODELS
It may be helpful to explain how Eqs. 1 and 4 are used in recently described multi-stage cancer models [13][14][15]. If we take a rate of mutations µ j per cell division for each of the rate-limiting mutational steps 1 to j, and d i divisions of cell i, then the probability of a stem cell surviving without the jth rate limiting mutation is Similarly, the probability of a given stem cell having mutation j is This is the solution of Zhang et al. [14] to the recursive formula of Wu et al. [13] (see Appendix of Zhang et al. [14] for details). Using Eq. 4, the survival of the ith stem cell is described by, Now assuming all n stem cells are equivalent and have equal rates µ i = µ j for all i, j, and consider only one path to cancer with m mutational steps, then, The probability of cancer within m divisions, often referred to as "theoretical lifetime intrinsic cancer risk", is, This is the equation derived by Calabrese and Shibata [15], and that Zhang found as the solution to the model of Wu et al [13,14].
Therefore, in addition to the models of Wu and Calabrese being equivalent cancer models needing m mutational steps, the models also assume that the order of the steps is not important. This differs from the original Armitage-Doll model that considered a sequential set of rate-limiting steps, and was exactly solved by Moolgavkar [10]. Eqs. 8 and 9 are equivalent to assuming: (i) equivalent stem cells, (ii) a single path to cancer, (iii) equivalent divisions per stem cell, and, (iv) equivalent mutation rates for all steps.
Despite the differences in modelling assumptions for Eq. 9 and the Armitage-Doll model, their predictions can be quantitatively similar. To see this, use the Armitage-Doll approximation of µd 1, to expand, If cell divisions are approximately uniform in time, then we can replace µd with µt, with µ now a rate per unit time. Then expanding exp(−µt) 1 − µt, gives, The incidence rate h = f /S is then h n s µ m t m−1 , the same as the original (approximate) Armitage-Doll solution [1]. This approximate solution is expected to become inaccurate at sufficiently long times.
An equivalent expression to Eq. 8 was known to Armitage, Doll, and Pike since at least 1965 [16], as was its limiting behaviour for large n. The authors [16] emphasised that many different forms for the F i (t i ) could produce approximately the same observed F (t), especially for large n, with the behaviour of F (t) being dominated by the small t behaviour of F i (t). As a result, for sufficiently small times power-law behaviour for F (t) is likely, and if longer times were observable then an extreme value distribution would be expected [3,16,17]. However the power-law approximation can fail for important cases with extra rate-limiting steps such as a clonal expansion [4][5][6]. It seems likely that a model that includes clonal expansion and cancer detection is needed for cancer modelling, but the power law approximation could be used for all but the penultimate step, for example. A general methodology that includes Some failures require a sequence of independent events to occur, each following the one before ( Figure 3). A well-known example is the Armitage-Doll multistage cancer model, that requires a sequence of m mutations (failures), that each occur with a different constant rate. The probability density for failure time is the pdf for a sum of the m independent times t j to failure at each step in the sequence, each of which may have a different probability density function f j (t j ). A general method for evaluating the probability density is outlined below, adapting a method described by Jaynes (page 569).
is read as "A and B and C", and expand using the product rule P (A, B) = P (A|B)P (B), To evaluate the integrals, take the Laplace transform with respect to t, to give, This factorises as, Giving a general analytical solution as, where L −1 is the inverse Laplace transform. Eq. 15 is similar to the relationship between whose derivation is analogous to Eq. 17 but with integrals replaced by sums. The survival and hazard functions for f (t) can be obtained from Eq. 16 in the usual way. For example, that can be used in combination with Eq. 1. A number of valuable results are easy to evaluate using Eq. 16, as is illustrated in Section VI.
A useful related result is, that can be inferred from Eq. 16 with m = 2, by replacing f 1 (t 1 ) with f ( n−1 j=1 t j ) and f 2 (t 2 ) with f n (t n ). Eq. 20 can be solved using the convolution theorem for Laplace transforms, that gives, which is sometimes easier to evaluate than two Laplace transforms and their inverse. In general, solutions can be presented in terms of multiple convolutions if it is preferable to do so. Eqs. 19 and 21 are particularly useful for combining a known solution for the sum of (n − 1) samples such as for cancer initiation, with a differently distributed nth sample, such as the waiting time to detect a growing cancer. A final remark applies to the sum of random variables whose domain extends from −∞ to ∞, as opposed to the range 0 to ∞ considered so far. In that case an analogous calculation using a Fourier transform with respect to t in Eq. 13 leads to analogous results in terms of Fourier transforms, with Eq. 22 is mentioned for completeness, but is not used here.
A general solution to Eq. 16 can be given in terms of definite integrals, with, This can sometimes be easier to evaluate or approximate than Eq. 16. A derivation is given in Appendix A1. Eq. 23 allows a generalised Schwinger/Feynmann parameterisation [18] to be derived. Writing g j (s) = L [f j (t j )] and taking the Laplace transform of both sides of Eq.
23, gives, which includes some well known Schwinger/Feynmann parameterisations as special cases.
This is discussed further in Appendix A1.

VI. MODELLING SEQUENTIAL EVENTS -EXAMPLES
In the following examples we consider the time t = m i=1 t i for a sequence of events, with possibly different distributions f i (t i ) for the time between the (i − 1)th and ith event.
Some of the results are well-known but not usually presented this way, others are new or poorly known. We will use the Laplace transforms (and their inverses), of, Sums of gamma distributed samples (equal rates): A sum of gamma distributed variables with equal rate parameters can be calculated from, The Armitage- as was used in the original Armitage-Doll paper. Note that an equivalent distribution can be produced by a different combination of (fewer) hazard functions with h i ∼ t p i i provided t P i (p i +1) = t m−1 . This flexibility is removed in Moolgavkar's exact solution to the model [10], that is described next.
For example, if n = 3 then, Using partial fractions, we can write, Allowing the inverse Laplace transforms to be easily evaluated, giving, Note that the result is independent of the order of sequential events, but unlike the approximate solution to the Armitage Doll model [1], the exact solution allows less variability in the underlying models that can produce it. Also note that the leading order terms of an expansion in t cancel exactly, to give identical leading-order behaviour as for a power-law approximation (with p = 0) A general solution can be formed using a Schwinger/Feynmann parameterisation [18] of,  Replacing µ i with s + µ i in Eq. 32, then we can write Eq. 28 as, (which is simpler, but equivalent in effect, to repeatedly using the convolution formula).
Completing the integrals will generate Moolgavkar's solution for a given value of m. For example, taking m = 3 and integrating once gives, Integrating a second time, and simplifying, gives Eq. 31. The relationship between Schwinger/Feynmann parameterisations, Laplace transforms, and the convolution theorem is discussed further in Appendix A.
Moolgavkar [10] used induction to provide an explicit formula for f (t), with, where, For small times the terms in a Taylor expansion of Eq. 35 cancel exactly, so that f (t) For example, if the shape parameters p i of Gamma distributed samples are integer-valued, then we can write, where the product of differential operators can be taken outside the product of Laplace transforms because ∂/∂µ i (1/(s + µ j )) is zero for i = j. Using Eq. 37 we can replace the product of Laplace transforms with a sum, then only keep derivatives of non-zero terms, giving, Finally, take the inverse Laplace transform L −1 [1/(s + µ i )] = e −µ i t , to give, a general solution for sums of Gamma distributed samples with integer-valued shape parameters p i (and arbitrary rate parameters µ i ).
For example, if p i = 2 for all i, then evaluation of Eq. 40 can be accomplished by considering the ith term, and noting that for j = i there is only one factor involving µ j in χ i (Eq. 36). All the derivatives with respect to j = i then make Taking the final derivative with respect to µ i and cancelling factors of (-1) gives, for the sum of Gamma distributions with shape parameters p i = 2 and arbitrary rate parameters, and χ i (m) as defined in Eq. 36.
Sums of samples with different distributions: An advantage of the method described above, is that it is often easy to calculate pdfs for sums of differently distributed samples. For the first example, consider two samples from the same (or very similar) exponential distribution, and a third from a different exponential distribution. The result can be obtained by writing µ 3 = µ 2 + in Eq. 31, and letting → 0. Considering the terms involving exponents of µ 2 and µ 3 , Using Eq. 31 and letting → 0, gives, More generally, this argument shows that a sum of exponentially distributed samples with different rates, smoothly approximates a gamma distribution as the rates become increasingly similar. As shown in Eq. 25, the sum of exponentially distributed samples with exactly the same rate, has a gamma distribution Failure involving a combination of sequential and non-sequential steps: If a path to failure involves a combination of sequential and non-sequential steps, then the necessary set of sequential steps can be considered as one of the non-sequential steps, with overall survival given by Eq. 1 and the survival for any sequential set of steps calculated from Eq. 18 ( Figure 4). Eq. 5), and non-sequential steps (e.g. (n,1) to (n,m n ) using Eq. 16), that may be dependent on each other (e.g. Eq. 50). For the purposes of modelling, a sequence of dependent or multiple routes can be regarded as a single step (e.g. (n − 1,j)).
Clonal-expansion cancer models: Some cancer models have a clonal expansion of cells as a rate-limiting step prior to a cancer's detection [4][5][6]. For example, Michor et al. This gives a survival function for cancer detection of, where, a, c, are rate constants, and N is the total number of cells prior to cancer initiation. Noting that t 0 x(y)dy = log(e ct + (N − 1))/c → t, as t → ∞ and x(t) → 1, then the tail of the survival curve falls exponentially towards zero with time.
Alternatively, we might expect the likelihood of cancer being diagnosed to continue to increase with time since the cancer is initiated. For example, a hazard function that is linear in time would give a Weibull distribution with S(t) = e −at 2 . It is unlikely that either this or the logistic model would be an equally good description for the detection of all cancers, although they may both be an improvement on a model without either. Both models have a single peak, but differ in the tail of their distribution, that falls as ∼ e −act for the logistic model and ∼ e −at 2 for the Weibull model. Qualitatively, we might expect a delay between cancer initiation and the possibility of diagnosis, and diagnosis to occur almost inevitably within a reasonable time-period. Therefore a Weibull or Gamma distributed time to diagnosis seems reasonable for many cancers, with the shorter tail of the Weibull distribution making it more suitable approximation for cancers whose diagnosis is almost inevitable. (The possibility of misdiagnosis or death by another cause is not considered here.) For example, noting that Moolgavkar's solution is a linear combination of exponential distributions, to combine it with a Weibull distribution for cancer detection f 1 (t 1 ) = −d/dt 1 (e −bt 2 1 /2 ), we can consider a single exponential term at a time. Taking f 2 (t 2 ) = ae −at 2 , and using the convolution formula Eq. 21, we get, where we integrated by parts to get the last line. This may be written as, with erf(x) = 2 √ π x 0 e −z 2 dz. Similarly for a Gamma distribution with f 1 = b p t p−1 e −bt /Γ(p) and an exponential, f 2 (t 2 ) = ae −at , then assuming b > a, where γ(p, t(b − a)) is the normalised lower incomplete Gamma function, which is available in most computational mathematics and statistics packages. If a > b then f 1 and f 2 must be exchanged and the result is most easily evaluated numerically.
For failure by a particular route ABC we need the probability for the sequence of events, A&(B&C), then (B&C)|A, then C|(AB). We can calculate this using Eq. 16, for example giving, from which we can construct S 1 (t) = ∞ t f ABC (y)dy. Although in principle every term in e.g. Eqs. 49 and 50 need evaluating, there will be situations where results simplify. For example, if one route is much more probable than another -e.g. if it is approximately true that landslides only occur after deforestation, that may be due to fire, then we only need to evaluate the probability distribution for that route.
As another example, if all the f i are exponentially distributed with different rates, then f ABC will be described by Moolgavkar for some a i > 0 and p i > 0. Then F (t) = 1 − S(t), f (t) = −dS/dt, and h(t) f (t), can be approximated by a sum of power series in time. If one route is much more likely than the others then both f (t) and h(t) can be approximated as a single power of time, with the approximation best at early times, and a cross-over to different power-law behaviour at later times.

VIII. SUMMARY
The purpose of this article is to provide a simple mathematical framework to describe existing multi-stage cancer models, that is easily adaptable to model events such as failure of complex systems, cascading disasters, and the onset of disease. of Beta functions to factorise and re-sum the resulting expression. In mathematical notation, . Now noting that the Beta function has,  Using Eq. A4 to replace 1/Γ ( m i=1 (n i + 1)) in Eq. A1, and grouping terms,     For example, taking g j (s) = 1/(s + a j ) p j and noting that L −1 {1/(s + a j ) p j } = t p j −1 e −a j t /Γ(p j ), then we get,