Finding Reproduction Numbers for Epidemic Models & Predator-Prey Models of Arbitrary Finite Dimension Using The Generalized Linear Chain Trick
aa r X i v : . [ q - b i o . P E ] A ug Finding Reproduction Numbers forEpidemic Models & Predator–Prey Models ofArbitrary Finite Dimension UsingThe Generalized Linear Chain Trick
Preprint
Hurtado, Paul J.
University of Nevada, RenoORCID: 0000-0002-8499-5986 [email protected]
Richards, Cameron
University of Nevada, RenoORCID: 0000-0002-1620-9998August 18, 2020
Abstract
Reproduction numbers, like the basic reproduction number R , play an important role inthe analysis and application of dynamic models, including contagion models and ecologicalpopulation models. One difficulty in deriving these quantities is that they must be computed on amodel-by-model basis, since it is typically not practical to obtain a general reproduction numberexpressions applicable to a family of related models, especially if some are of different dimensions.For example, this is the case for SIR-type infectious disease models derived using the linearchain trick (LCT). Here we show how the next generation operator approach for computingreproduction numbers can be used in conjunction with the generalized linear chain trick (GLCT)to find reproduction numbers expressions that apply generally to such families of models ofvarying dimensions. We further show how the GLCT enables modelers to draw insights fromthese results by using theory and intuition from continuous time Markov chains (CTMCs) andtheir absorption time distributions (i.e., the phase-type probability distributions). To do this,we first review the GLCT and other connections between mean-field ODE model assumptions,CTMCs, and phase-type distributions. We then apply this technique to a generalized familyof SEIRS models of arbitrary finite dimension, and to a generalized family of predator-prey(Rosenzweig-MacArthur) models of arbitrary finite dimension. These results highlight the utilityof the GLCT for the derivation and analysis of mean field ODE models, especially when usedin conjunction with theory from CTMCs and their associated phase-type distributions. K eywords basic reproductive ratio; linear chain trick; gamma chain trick; phase-type distribution;Coxian distribution; Erlang distribution; SIR; consumer-resource reprint - August 18, 2020 Contents R ) for SEIR and SEIRS Models . . . . . . . 112.4 Predator-Prey Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 R for the SEIRS model with Erlang latent and infectious periods . . . . . 18 reprint - August 18, 2020 The basic reproduction number R is perhaps the most well known threshold quantity derivedfrom epidemic disease models (Brauer and Castillo-Chavez, 2011; Brauer, van den Driessche, andWu, 2008; Diekmann and Heesterbeek, 2000; Diekmann, Heesterbeek, and Metz, 1990; Heesterbeek,2002; Heffernan, Smith, and Wahl, 2005; Hethcote, 2000; Hyman and Li, 2000, 2005; Kermack andMcKendrick, 1927, 1932, 1933, 1991a,b,c; Roberts and Heesterbeek, 2012; van den Driessche andWatmough, 2002). Along with its time-varying counterpart, the effective reproduction number R t (also called the replacement number; Hethcote, 1976, 2000), it plays a central role in how we under-stand and attempt to control the spread of infectious diseases, since these quantities summarize howthe overall transmission process reflects the combined impact of various biological processes, e.g.,host susceptibility, infectiousness, recovery, between-host contact processes, and control measuressuch as vaccination and quarantine. Such reproduction numbers can also be derived using modelsof other kinds of contagion and other biological populations, e.g., multispecies ecological models(Hilker and Schmitz, 2008; Hurtado, Hall, and Ellner, 2014; Roberts and Heesterbeek, 2012) andin cellular-level models of cancer (e.g., Eftimie et al., 2011) and viral infection (e.g., Hews et al.,2009).Reproduction number expressions are, however, only as good as the model assumptions from whichthey are derived. For example, the importance of incorporating non-exponentially distributedlatent and infectious periods in SEIR-type models (see section 2.3 below) is well known (Wearing,Rohani, and Keeling, 2005). The standard linear chain trick (LCT; Hurtado and Kirosingh, 2019;Smith, 2010) has been widely used for decades as a way of replacing the assumption of exponentiallydistributed dwell times in ODE models with Erlang distributions (i.e., gamma distributions withinteger shape parameters). This yields a system of ODEs whose dimension depends on the Erlangshape parameter value(s), which determine(s) the number of new sub-states (and their governingequations) introduced into the model. Deriving a new model using the LCT, therefore, actuallydefines a countably large family of new models, each of a different (finite) dimension. One challengewith this approach is that it has proven difficult to then derive a single general formula for thebasic reproduction numbers for such models, as a function of the Erlang shape parameter(s), andinstead these are typically derived on a case-by-case basis using a single, fixed shape parametervalue (although, see Bonzi et al. (2010) and Hyman and Li (2005), for examples). The generalizedlinear chain trick (GLCT; Hurtado and Kirosingh, 2019) can help overcome this challenge.To illustrate how, here we introduce a generalized SEIRS model that encompasses multiple SEIRS-type models as special cases, and we show how the next generation operator approach (Diekmann,Heesterbeek, and Roberts, 2009; Diekmann and Heesterbeek, 2000; Diekmann, Heesterbeek, andMetz, 1990; Heesterbeek, 2002; Roberts and Heesterbeek, 2012; van den Driessche and Watmough,2002) can be used to calculate a general expression for the basic reproduction number R that holdsfor all instances of this general model, regardless of the dimension of those models so long as theyremain finite. We also show how this same technique can be used to find a general expression for thepredator reproduction number in a generalized Rosenzweig-MacArthur predator-prey (consumer-resource) model.The key to the success of this approach is that we conduct this analysis on a matrix-vector formof the model, as in Hyman and Li (2005), but we formulate the model using the GLCT andinterpret the resulting reproduction number expressions through the lens of continuous time Markovchains (CTMCs) and phase-type distributions, i.e., the broad family of hitting time distributionsfor CTMCs; these include exponential, generalized Erlang, and Coxian distributions (Bladt and3 reprint - August 18, 2020 Nielsen, 2017a,b; Horv´ath and Telek, 2017; Reinecke, Bodrog, and Danilkina, 2012).In the sections below, we first review the GLCT, LCT, and other connections between mean fieldODE models, CTMCs, and phase-type distributions. We then review a standard SEIRS model, R ,and introduce our generalized SEIRS and predator-prey models before deriving general reproductionnumber expressions for each. In this section, we give a brief overview of the phase-type family of univariate, matrix exponentialprobability distributions, and their connection to continuous time Markov chain (CTMCs) , basedon Hurtado and Kirosingh (2019) and Hurtado and Richards (2020a,b). Some familiar examplesof distributions in the phase-type family are the exponential distribution, Erlang distribution (i.e.,gamma distributions with integer shape parameters k ), hypoexponential (or generalized Erlang)distribution, hyper-Erlang (finite Erlang mixture) distribution, and the Coxian distribution. Vari-ous statistical tools exist for fitting these phase-type distributions to data (e.g., Horv´ath and Telek,2017, 2020). For more details, see Altiok (1985), Bladt and Nielsen (2017a), Horv´ath, Scarpa, andTelek (2016), Horv´ath et al. (2012), Reinecke, Bodrog, and Danilkina (2012), and Reinecke, Krau,and Wolter (2012).Phase-type distributions describe the time it takes to first reach an absorbing state in a CTMC.For a given phase-type distribution, the corresponding CTMC, with k transient states and the k + 1state an absorbing state, has the ( k + 1) × ( k + 1) transition rate matrix " A a0 (1)where A is k × k and a is a length k column vector. As detailed in eqs. (4), t he correspondingphase-type distribution is parameterized by A and a length k vector α , which is the CTMC inintialdistribution vector for the k transient states (and thus the initial state is the absorbing state withprobability α ∗ = 1 − P ki =1 α i ).Since the rows in a transition rate matrix must sum to zero, the vector of loss rates ( a ) from eachtransient state to the absorbing state can be written in terms of A as a = − A 1 (2)(here, and below, denotes an appropriately long column vector of ones). The negative valuesalong the diagonal1 of A are the loss rates from the corresponding transient states to any other For readers only familiar with discrete time Markov chains, CTMCs are very similar, except that only statetransitions out of a given state are considered (so transition probability p ii = 0), and rather than state transitionsoccurring after a fixed time step, they occur after an exponentially distributed duration of time (i.e., an exponentiallydistributed dwell time ), where the time spent in the i th state is exponential with its own rate r i >
0. CTMCs can beparameterized with a transition probability matrix, just like a discrete time Markov chain, and also a rate vector of r i values. But, these quantities are instead more commonly combined into a transition rate matrix with − r i alongthe diagonal, and off-diagonal entries in row i column j of the form r i p ij where p ij is the transition probability fromthe i th state to the j th state, and thus r i p ij is the rate at which individuals move from the i th state into the j th state. reprint - August 18, 2020 state, and the non-negative off diagonal entries1 of A are the rates of influx into each transientstate from the other transient states (Bladt and Nielsen, 2017a; Hurtado and Kirosingh, 2019;Hurtado and Richards, 2020a,b).One important family of phase-type distributions are the Erlang distributions. These are the gammadistributions with integer shape parameters, which can be thought of as the sum of finitely many i.i.d. exponential distributions. Specifically, the sum of k i.i.d. exponential random variables withrate r is Erlang distributed with rate r and shape k . Erlang distributions can also be parameterizedin terms of their mean τ and coefficient of variation c v (the standard deviation σ over the mean τ ), where τ = kr , σ = kr , c v = 1 √ k , and thus, k = τ σ = 1 c v , and r = τσ = 1 c v τ = kτ . (3)The generalized Erlang distributions (also called hypoexponential distributions) are equivalent tothe sums of k independent exponential distributions, each with potentially distinct rates r i .Another important family of phase-type distributions are the Coxian distributions. These are theabsorption time distributions for CTMCs where each of k transient states has its own rate r i (similarto generalized Erlang distributions) but for each state there is some probability p i of entering thenext state in the ‘‘chain’’, or with probability 1 − p i transitioning straight to the absorbing state (seeAppendix B for an example). Phase-type distributions can be classified as acyclic (transient statescannot be revisited once left) and cyclic (one or more transient states can be revisited multipletimes due to cycles in the transition rate matrix), and all acyclic phase-type distributions haveCoxian representations (Cumani, 1982; O’Cinneide, 1991, 1993). We make use of the following properties of phase-type distributions Markov chains.The phase-type family of distributions is closed under various operations including addition (con-volution), minimum, maximum, and finite mixtures (Bladt and Nielsen, 2017a; Neuts, 1981). Theydo not have unique parameterizations, generally speaking. Phase-type distributions have a prob-ability density function, cumulative density function (CDF), Laplace-Stieltjes Transform of theCDF, and j th moment generating function given, respectively, by the following formulas adaptedfrom Reinecke, Bodrog, and Danilkina (2012): f ( t ) = α T e A t ( − A1 ) (4a) F ( t ) = 1 − α T e A t (4b) L ( s ) = α ∗ + α T ( s I − A ) − ( − A 1 ) (4c) E ( T j ) = j ! α T ( − A ) − j . (4d)As above, α ∗ is the probability of starting the CTMC in the absorbing state (typically α ∗ = 0), isa column vector of ones, and I is the identity matrix, where each has appropriate dimensions given A . If a phase-type distribution with parameters α and A has α ∗ >
0, then it can be thought ofas the zero-inflated mixture distribution of the phase-type distribution with parameters α / (1 − α ∗ )and A , and a Dirac delta distribution (point mass w.p. 1) at 0 with respective mixing probabilities5 reprint - August 18, 2020 − α ∗ and α ∗ . Let U = − A − (this is called the Green matrix, and is analogous to the fundamentalmatrix in discrete time Markov chains). The entries u ij represent the expected time spent in the j th state prior to reaching the absorbing state, given that the initial state was the i th state (Bladtand Nielsen, 2017a). Expected Rewards.
The following is a property of the reward process associated with a givenCTMC. Specifically, suppose that g i is the reward rate associated with the i th transient state inan absorbing CTMC, where the reward amount accrued while spending a duration of time T i inthe i th state is g i T i . Let Y be the total reward accrued prior to reaching the absorbing state, andlet g be the column vector of reward rates for the transient states. Then the expected reward is(see Appendix A for a proof) E ( Y ) = α T ( − A − ) g . (5) Minimum of Phase-Type Random Variables.
As mentioned above, the minimum of twoindependent phase-type random variables, parameterized by α i , A i ( i = 1 , α min = α ⊗ α (6a) A min = A ⊕ A . (6b) CTMC Absorption Probabilities.
It is useful to know how individuals are distributed acrossmultiple absorbing states upon leaving the transient states in an absorbing CTMC. Standard CTMCtheory can be used to compute the probabilities of individuals going to a particular absorbing state.In the context of the GLCT (detailed below) this can be a useful tool for thinking through modelingscenarios in which individuals leave a single state but are distributed across multiple subsequentstates (e.g., transitioning from an infected state to being either recovered or deceased). Morespecifically, these probabilities can be calculated using the transition probability matrix P (for theembedded jump process) implicit in a transition rate matrix like eq. (1).Matrix P can be obtained from A as follows, where we let D v denote a diagonal matrix with entriesof vector v along its diagonal, and we let vector diag( M ) denote the diagonal entries of matrix M .The vector of exponential dwell-time rates for the transient states r = − diag( A ). Define matrix Q by dividing each row of A + D r by the corresponding rate r i for that row, which yields therelationship D r − A = Q − I (7)thus, Q = I + D r − A and − A = D r (cid:0) I − Q (cid:1) . (8)If there are one or more absorbing states, then P = " Q R0 I (9) This reward rate can be thought of as a mean rate of reward accrual in a more general renewal reward processcontext. reprint - August 18, 2020 where R has as many rows as there are transient states, and as many columns as there are absorbingstates, and I is an identity matrix. R is needed to compute these absorption probabilities, but it is only uniquely determined by A ifthere is exactly 1 absorbing state. Otherwise, R can be constructed as follows. The loss rate vector a = − A1 in eqs. (1) and (2) can be divided (elementwise) by the rate vector r which yields vector q = D r − ( − A1 ) . (10)The entries q i are the probabilities of entering an absorbing state upon leaving the i th transientstate. Let C be a matrix that is the same dimension as R , and the i th row and j th column entryof C is the probability that an individual enters the j th absorbing state given that it left the i th transient state. Multiplying the i th row of C by q i yields the i th row of R , thus R = D q C . (11)The desired absorption probabilities, for hitting the j th absorbing state given the initial state wasthe i th transient state, are the ij entries in the following matrix (Bladt and Nielsen, 2017a; Resnick,2002): (cid:0) I − Q (cid:1) − R = (cid:0) − D r − A (cid:1) − R = ( − A − ) D r R (12)For additional properties of phase-type distributions, see Altiok (1985), Bladt and Nielsen (2017a),Horv´ath, Scarpa, and Telek (2016), Horv´ath et al. (2012), Reinecke, Bodrog, and Danilkina (2012),and Reinecke, Krau, and Wolter (2012) and references cited in Hurtado and Kirosingh (2019). The following overview of the Generalized Linear Chain Trick (GLCT) is adapted from the more de-tailed presentation of the GLCT found in Hurtado and Kirosingh (2019) and Hurtado and Richards(2020a,b). As we show below, ODE models derived using the GLCT can be used to prove generalresults for models that could otherwise be obtained using the standard Linear Chain Trick (LCT).This is partly the result of the equations being in a specific matrix-vector form – which can beanalyzed without the constraint of a fixed model dimension – but also because the analysis ofsuch models can give rise to other phase-type distributions or related quantities, and it is useful torecognize them as such.The GLCT enables modelers to rigorously formulate new ODE models relatively quickly by by-passing the need to explicitly derive those ODEs using mean field integral equations and theirderivatives (Hurtado and Kirosingh, 2019; Hurtado and Richards, 2020b). This can be done inroughly one of two ways: using the GLCT to generalize an existing model, or to derive a new ODEmodel from first principles. These approaches can be briefly summarized as follows.First, just as the LCT is used to introduce Erlang distributed dwell times into an existing ODEmodel, the GLCT can be used to take an existing (e.g., ODE or DDE) model and modify itsassumptions to introduce phase-type distributed dwell times resulting in a new ODE model. Thiscan most easily be done by first applying the standard LCT, then writing that new set of ODEsin matrix-vector form a la
Theorem 1 below to apply the GLCT for phase-type distributions. Forexamples, see Hurtado and Richards (2020b). 7 reprint - August 18, 2020
Second, new ODE models can be derived from first principles by using the GLCT to organizeunderlying (stochastic) assumptions about how individuals transition among different states, andthen quickly produce a mean field ODE model. This is done by framing assumptions in the contextof CTMCs or (non-homogeneous) Poisson processes, as detailed in Hurtado and Kirosingh (2019).Theorem 1 below is a re-statement of the GLCT for phase-type distributions (Corollary 2) inHurtado and Kirosingh (2019), as it appears in Hurtado and Richards (2020b). See Hurtadoand Kirosingh (2019) for additional results that further clarify the link between stochastic modelassumptions and mean field ODE model equations, including a more general statement of theGLCT.
Theorem 1 ( GLCT for Phase-Type Distributions ) . Assume individuals enter a state (call itstate X) at rate I ( t ) ∈ R and that the distribution of time spent in state X follows a continuousphase-type distribution given by the length k initial probability vector α and the k × k matrix A .Then partitioning X into k sub-states X i , and denoting the corresponding amount of individuals instate X i at time t by x i ( t ) , the mean field equations for these sub-states x i are given by ddt x ( t ) = α I ( t ) + A T x ( t ) (13) where the rate of individuals leaving each of these sub-states of X is given by the vector ( − A 1 ) ◦ x ,where ◦ is the Hadamard (element-wise) product of the two vectors, and thus the total rate ofindividuals leaving state X is given by the sum of those terms, i.e., ( − A 1 ) T x = − T A T x . The standard linear chain trick (LCT) is a special case of Theorem 1. For completeness, it isprovided below as stated in Hurtado and Richards (2020b). We have chosen to give the LCT forgeneralized Erlang (hypoexponential) distributions (i.e., the sum of k independent exponentiallydistributed random variables, each with potentially different rates r i ) as this only changes a fewsubscripts in the mean field equations from the Erlang-based form of the LCT. Corollary 1 ( Linear Chain Trick for Erlang and Hypoexponential Distributions ) . Consider the GLCT above (Theorem 1). Assume that the dwell-time distribution is a generalizedErlang (hypoexponential) distribution with rates r = [ r , r , . . . , r k ] T , where r i > , or an Erlangdistribution with rate r (all rates r i = r ) and shape k (or if written in terms of shape k and mean τ = k/r , use r = k/τ ). Then the corresponding mean field equations are dx dt = I ( t ) − r x dx dt = r x − r x ...dx k dt = r k − x k − − r k x k . (14) Proof.
The phase-type distribution formulation of a generalized Erlang distribution with rates8 reprint - August 18, 2020 r i > i = 1 , . . . , k is given by eqs. (15). Substituting these into eq. (13) yields equations (14). α = ... and A = − r r · · · − r r . . . ... . . . . . . . . . . . . . . . − r k − r k − · · · − r k . (15) (cid:4) In addition, the following intuition – regarding individuals transitioning from one state to oneof multiple states – is useful for further understanding the links between underlying stochasticmodel assumptions and the corresponding mean field ODE model structure. First, recall that theminimum of multiple independent exponentially distributed random variables (each with rate r i > i = 1 , . . . , n ) is itself an exponentially distributed with rate r = r + · · · + r n . Consider thefollowing pair of simple scenarios, which yield the same mean field equations according to Theorem2 below.Suppose the time and individual spends in a given state is assumed to be the minimum of multipleexponentially distributed random variables, where each corresponds to a certain event time andthe individual leaves their current state once the first of those events occurs. Further assume thatif it is the i th of these events that occurs first, then the individual transitions to the i th recipientstate with probability 1. It follows that the mean field ODE terms that correspond to the scenariojust described are equivalent to the mean field ODE terms obtained by instead assuming that thedwell time in the focal state is exponentially distributed with rate r = r + · · · + r n , and that uponleaving that state individuals are distributed across the n recipient states with probabilities r i /r .The following theorem (a special case of Theorem 7 in Hurtado and Kirosingh (2019)) gives a moreformal, and more general, statement of the mean field equivalence of the two scenarios describedabove. Theorem 2 ( Mean field equivalence of proportional outputs & competing exponentialevent times ) . Suppose state X has a dwell time given by random variable T = min i T i , whereeach T i is exponentially distributed with rate r i , i = 1 , . . . , n and individuals transition to one of m states Y ℓ , ℓ = 1 , . . . , m , with probability p iℓ ( T ) when T = T i . The corresponding mean field modelis equivalent to having instead assumed that the X dwell time is exponentially distributed with rate r = P ni =1 r i , and the transition probability vector for leaving X and entering one of the states Y ℓ is given by p ℓ = P ni =1 p iℓ r i /r . Example:
Consider the well known SIR model. The net loss rate from the infected state is oftenwritten − ( γ + µ + ν ) I where the recovery rate is γ , the baseline mortality rate is µ , and the disease-induced mortality rate is ν . Accordingly, the fraction of those individuals leaving the infected stateand entering the recovered state is the product of that total loss rate and recovery fraction: γγ + µ + ν and thus the rate of individuals entering the recovered state is( γ + µ + ν ) I γγ + µ + ν = γ I. More generally, in the context of nonhomogeneous Poisson processes, the minimum of two 1 st event timescorresponding to two independent Poisson processes, with rates r ( t ) and r ( t ), can be thought of as the 1 st eventtime under a Poisson process with rate r ( t ) = r ( t ) + r ( t ). See Hurtado and Kirosingh (2019) for details. reprint - August 18, 2020 Compare the above to the SEIR model terms to those in eqs. (17).
The well known SEIR model with mass action transmission, a fixed population size, and no birthsor deaths is given by eqs. (16) as stated in Hurtado and Richards (2020a). A more general SEIRmodel was introduced in Hurtado and Richards (2020a), which has phase-type distributed latentand infectious periods. Below, we further extend this model to include births and deaths, as wellas waning immunity, where that duration of immunity is also phase-type distributed. But first, wereview these simpler models and their basic reproduction numbers.The state variables below have all been scaled by the total population size N so that S + E + I + R = 1,where S is the susceptible fraction of the population, E the fraction with latent infections (i.e., notyet symptomatic or contagious), I the fraction who are infected and contagious, and R is therecovered fraction of the population. dSdt = − β S I (16a) dEdt = β S I − r E E (16b) dIdt = r E E − r I I (16c) dRdt = r I I (16d)The model eqs. (16) can be interpreted as the mean field model for a continuous time stochasticSEIR model in which the time individuals spend in states E and I follow exponential distributionswith respective rates r E and r I , or equivalently, with respective means τ E = 1 /r E and τ I = 1 /r I .The SEIRS model given by eqs. (17) extend the above model to include births, deaths, and waningimmunity. To provide a simple illustrative example with a single locally asymptotically stabledisease free equilibrium, we assume a constant birth rate Λ and per-capita mortality rates µ whichdo not differ among the different infectious states, thus dN/dt = Λ − µ N . Upon recovery, immunitylasts for a period of time before individuals return to the susceptible class at rate ǫ . Again, wecan interpret these parameters as the mortality rate µ or mean lifetime τ N = 1 /µ , and waningimmunity loss rate ǫ or the mean duration of immunity τ R = 1 /ǫ . dSdt = Λ − µ S − β S I + ǫ R (17a) dEdt = β S I − r E E − µ E (17b) dIdt = r E E − r I I − µ I (17c) dRdt = r I I − µ R − ǫ R (17d)The SEIRS model with exponentially distributed latent and infectious periods (eqs. (17)) can be10 reprint - August 18, 2020 further extended using the LCT to instead assume Erlang distributed latent and infectious periods ,with the same means τ E and τ I as the models above, but with coefficients of variation cv E and cv I where we choose these so that the square of the inverse of the cv is an integer, i.e., where k E = 1 / √ cv E , k I = 1 / √ cv I are integers .Equations for the resulting SIERS model with Erlang distributed latent and infectious periods are dSdt = Λ − β I S − µ S + ǫ R (18a) dE dt = β S I − k E τ E E − µ E (18b) dE j dt = k E τ E E j − − k E τ E E j − µ E j , j = 2 , . . . , k E − dE k E dt = k E τ E E k E − − k E τ E E k E − µ E k E (18d) dI dt = k E τ E E k E − k I τ I I − µ I (18e) dI i dt = k I τ I I i − − k I τ I I i − µ I i , i = 2 , . . . , k I − dI k I dt = k I τ I I k I − − k I τ I I k I − µ I k I (18g) dRdt = k I τ I I k I − µ R − ǫ R. (18h)Let us now review the derivation and interpretation of the basic reproduction number R for thethree models above. R ) for SEIR and SEIRS Models Using the method in van den Driessche and Watmough (2002), the expression for R for eqs. (16)is R = β S r I = β S |{z} rate of newinfectionsmean infectiousperiod z}|{ r I . (19)Observe that we can factor R into the product of a (per infected individual) rate of new infectionstimes the mean duration of infectiousness. We therefore interpret R as the expected number ofnew infectious cases per infectious individual added to a population at the DFE. Similarly, R forthe SEIRS model eqs. (17), with births, deaths, and waning immunity, is A simple generalization of this model would be to allow different rates for each sub-state of E and I . Then thedwell times would follow the sum of independent but non-identically distributed exponentials each with their ownrates r i . The distribution of such sums are known as the hypoexponential or generalized Erlang distributions. Alternatively, one could smoothly specify the cv using the mixture distribution approximation in the appendixof Hurtado and Kirosingh (2019) reprint - August 18, 2020 R = β S r E ( r E + µ )( r I + µ ) = β S |{z} rate of newinfectionslatent periodsurvival probability z }| { r E r E + µ r I + µ , | {z } mean infectiousperiod (20)which can be similarly factored and interpreted but includes an additional term for the expectedfraction of individuals who survive the latent period to become infectious.To find R for eqs. (18), the SEIRS model with Erlang distributed latent and infectious periods, weagain use the next generation operator approach van den Driessche and Watmough (2002), but weassume fixed values of k E and k I . This yields R expressions for each particular case considered,which each have a similar post hoc interpretation of terms as in eq. (20). While it is possibleto inspect a few specific cases and conjecture a general expression for R in such situations, it isoften the case that no general expression for R is formally obtained (e.g., Johnson, Landguth, andStone, 2016). In section 3 we show how this can be done by writing the model in a form suggestedby the GLCT (Theorem (1)).Consider the R expression obtained from eqs. (18), where k E , k I ≥
1, rates r E,k E = k E /τ E and r I = 1 /τ I . For example, if k E = 3 and k I = 4, then R = β S ( r E, ) ( r E, + µ ) (cid:18) r I, + µ + r I, ( r I, + µ ) + r I, ( r I, + µ ) + r I, ( r I, + µ ) (cid:19) . (21)This might prompt a modeler to a conjecture that, for arbitrary positive integers k E and k I , R = β S ( r E,k E ) k E ( r E,k E + µ ) k E (cid:18) r I,k I + µ + r I,k I ( r I,k I + µ ) + · · · + r I,k I k I − ( r I,k I + µ ) k I (cid:19) . (22)However this is only a conjecture without any additional supporting analyses, and such termscan sometimes be challenging to interpret with confidence: As with the R expressions for the twosimpler models above, we can interpret this expression post hoc as the per-infectious-individual rateof new infections, times the probability of surviving all sub-states of E, times the expected periodof infectiousness (although it is not immediately obvious that this is the proper interpretation).In section 3, we introduce a novel phase-type distributed SEIRS model (a further extension of eqs.(20)) using the GLCT, and illustrate how the natural matrix-vector formulation of the model thatarises from such a derivation can be used to find a very general form of R . In doing so, as detailedin section 3.3, we confirm that the above conjecture and R interpretation regarding eqs. (22)holds true for arbitrary positive integers k E and k I . In Hurtado and Richards (2020a), an extension of the Rosenzweig-MacArthur model was introducedwhich assumes an immature/mature stage structure in the predator population. Here we consider aslight extension of that model (the aforementioned model is the η = 0 case of the equations below)given by dNdt = r N (cid:18) − NK (cid:19) − a ( P m + η P imm ) k + N N (23a)12 reprint - August 18, 2020 d x dt = χ a Nk + N P m α x + A x T x (23b) d y dt = − T A x T x | {z } scalar α y + A y T y . (23c)The prey population ( N ) follows a logistic growth model in the absence of predators, and areremoved by predators according to a Holling type-II functional response. For simplicity, the im-mature predators that eventually reach the mature stage are tracked by k state variables in thevector x = [ x , . . . , x k ] T ( P imm = P x i ), following a maturation time that is phase-type distributedwith parameter vector α x and matrix A x . Upon maturing, predators then survive as adults fora period of time that is also phase-type distributed (tracked by m state variables in the vector y = [ y , . . . , y m ] T , where P m = P y i ) defined by parameter vector α y and matrix A y (Hurtadoand Richards, 2020a).In the next section, we show how reproduction numbers for this generalized predator-prey modelcan be obtained for arbitrary phase-type distributions, and then do the same for the similarlygeneralized SEIRS model introduced below. Mathematically, the method for finding the basic reproduction number in van den Driessche andWatmough (2002) can be reinterpreted to find the predator reproduction number in a model likeeqs. (23). The result of applying this technique to find the predator reproduction number R pred for this model are summarized in the following theorem. Theorem 3.
The predator-free (prey-only) equilibrium state for eqs. (23) is locally asymptoticallystable if R pred < , but unstable if R pred > , where R pred is given by R pred = χ aNk + N | {z } birth rate mean time predatorsspend as adults z }| {(cid:0) α y T ( − A y ) − (cid:1) (1 − α ∗ ) | {z } maturation prob. (24) where the product of the first two terms is the expected number of new immature predators createdby a single mature predator over the average predator’s reproductive lifetime, when introduced toa prey population at the predator-free (prey-only) equilibrium, and the last term is the fraction ofpredators that survive to reach maturity. Here − α x ∗ is the sum of entries in α x (see eqs. (4) ).Proof. Using the next generation operator approach detailed in van den Driessche and Watmough13 reprint - August 18, 2020 (2002) to find R pred , we first rewrite the system as ˙ x ˙ y ˙ N = α x χ a Nk + N P m | {z } F − − D diag( A x ) x − D diag( A y ) y a Nk + N ( P m + η P imm ) | {z } V − − ( A x − D diag( A x ) ) T x α y ( − T A x T x ) + ( A y − D diag( A y ) ) T y r N (cid:18) − NK (cid:19) | {z } V + The Jacobians of F and V evaluated at the prey-only equilibrium yield the matrices F = " | x |×| x | α x | y | T χ aNk + N | y |×| x | | y |×| y | and V = " − A x T α y | x | T A x T − A y T . (25)The subscripts on the matrices indicates their dimensions (e.g., | x | × | y | indicates there are as manyrows as sub-states in vector x , and as many columns as the length of vector y ). Below we alsouse the notation I to indicate an appropriately sized identity matrix, and k to indicate a columnvector of ones, that is length k .The general form for the inverse of a block matrix like eq. (36) is given by " A C D − = " A − − D − C A − D − (26)which gives V − = " ( − A x T ) − − ( − A y T ) − α y | x | T A x T ( − A x T ) − ( − A y T ) − . (27)Recall that R pred is equal to the spectral radius of FV − . Observe that FV − = " α x | y | T χ aNk + N ( − A y T ) − α y | x | T α x | y | T χ aNk + N ( − A y T ) − (28)is an upper triangular block matrix with a matrix of zeros as one of its diagonal blocks, so anynonzero eigenvalues will come from the upper left block of the matrix. This yields R pred = ρ (cid:18) α x T χ aNk + N ( − A y T ) − α y | x | T (cid:19) , (29)where ρ ( M ) is the spectral radius of matrix M . Also note that χ aNk + N is a scalar value and that T ( − A y T ) − α y = α y T ( − A y ) − is the mean survival time of predators after they reach maturity according to eqs. (4). Factoringthese scalars out of the spectral radius computation yields R pred = χ aNk + N (cid:0) α y T ( − A y ) − (cid:1) (cid:2) ρ (cid:0) α x | x | T (cid:1)(cid:3) . (30)14 reprint - August 18, 2020 Since α x is a column vector, and | x | T is a row vector of the same length, then the spectral radiusof square matrix α x | x | T is the dot product of α x and | x | , which equals the sum of the entries in α x . Recalling the definition of α ∗ from eqs. (4), define α x ∗ = 1 − ρ ( α x | x | T ), then it follows that R pred = χ aNk + N | {z } birth rate (cid:0) α y T ( − A y ) − (cid:1)| {z } mean time predatorsspend as adults maturation prob. z }| { (1 − α x ∗ ) . (cid:4) Observe that the above reproduction number expression holds for any choice of phase-type distri-bution assumptions (which determines the dimension of model eqs. (23)). Here we have used asomewhat simplistic set of assumptions regarding the predator birth rate and survival of offspringto maturity, to more clearly illustrate the process of computing the reproduction number using thenext generation operator approach, and to more clearly illustrate how the expression that is ob-tained can be interpreted using properties of phase-type distributions. We next consider a slightlymore complex example, to further illustrate the implementation of this approach.
Consider the SEIRS model eqs. (18). We apply the GLCT to generalize this model as follows(cf. the steps used to derive the simpler SEIR model in Hurtado and Richards (2020a)). Assumethat the latent period (time spent in state E) follows a phase-type distribution parameterized by alength k E vector α E and k E × k E matrix A E , the infectious period (time spent in state I) followsa phase-type distribution parameterized by a length k I vector α I and k I × k I matrix A I , andthe duration of immunity after recovery (time spent in state R) follows a phase-type distributiondefined by a length k R vector α R and k R × k R matrix A R .Let x = [ E , . . . , E k E ] T , y = [ I , . . . , I k I ] T , and z = [ R , . . . , R k R ] T be the column vectors of sub-states of E, I, and R, respectively, where E = P E i , I = P I i , and R = P R i . Also, let the forceof infection be given by λ ( t ) = β I + · · · + β k I I k I and let β = [ β , . . . , β k I ] T .The mean field ODEs for this generalized SEIRS model are dSdt = Λ( S, x , y , z ) − µ S − λ ( t ) S + (cid:0) − T A R T z (cid:1)| {z } scalar (31a) d x dt = α E λ ( t ) S + A E T x − µ x (31b) d y dt = α I (cid:0) − T A E T x (cid:1)| {z } scalars + A I T y − µ y (31c) d z dt = α R z }| {(cid:0) − T A I T y (cid:1) + A R T z − µ z . (31d) It is known that an n × n matrix that can be written as a length n column vector times a length n row vectorhas n − reprint - August 18, 2020 Here we have also assumed a more general birth rate Λ( S, x , y , z ) which we assume is a non-negative smooth function that yields a locally asymptotically stable disease free equilibrium (DFE)( S , , ,
0) in the absence of the pathogen. The terms labeled as scalars are the sums of terms ineach loss rate vector from the different exposed, infectious, and recovered (immune) states (cf. a in eq. (1) and see Theorem 1).Using the next generation operator approach to find R (van den Driessche and Watmough, 2002)yields the following result, where G E = A E − µ I k E × k E and G I = A I − µ I k I × k I (we will discussthe interpretation of these matrices below). Theorem 4.
The DFE for eqs. (31) is locally asymptotically stable if R < , but unstable if R > , where R is given by R = R P E → I (32) where R is the expected number of new infections created by a single infectious individualintroduced to the population at the DFE, given by R = α I T ( − G I − )( β S ) , (33) and P E → I is the probability of surviving the exposed state (E) and transitioning into the infectiousstate (I), which is given by the dot product P E → I = α E · (( − G E − )( − A E k E )) . (34) Proof.
The right hand side of eqs. (31) can be written as follows, yielding the given definitions for F and V which satisfy the requirements of Theorem 2 in van den Driessche and Watmough (2002). ˙ x ˙ y ˙ z ˙ S = ~α E λ ( t ) S | {z } F − − D diag( A E ) x + µ x − D diag( A I ) y + µ y − D diag( A R ) z + µ z λ ( t ) S + µS | {z } V − − ( A E − D diag( A E ) ) T x α I ( − T A E T x ) + ( A I − D diag( A I ) ) T y α R ( − T A I T y ) + ( A R − D diag( A R ) ) T z Λ( S, x , y , z ) + − T A R T z | {z } V + , The upper left block ( F ) of the Jacobian for F evaluated at the DFE is F = " k E × k E α E β T S k I × k E k I × k I = " α E β T S . (35)The upper left block V of the Jacobian for V = V − − V + evaluated at the DFE is V = " µ I − A E T α I T k E A E T µ I − A I T = " − G E T α I T k E A E T − G I T (36)16 reprint - August 18, 2020 Using eq. (26) yields V − = " ( − G E − ) T ( − G I − ) T ( − α I T k E A E T )( − G E − ) T ( − G I − ) T . (37)The basic reproduction number is the spectral radius of FV − , i.e., R = ρ ( FV − ), where FV − = " α E β T S ( − G I − ) T ( − α I T k E A E T )( − G E − ) T α E β T S ( − G I − ) T . (38)Observe that, as in the previous example, FV − is an upper triangular block matrix with a diagonalthat is all zeros except for the top left block. Therefore, R = ρ (cid:0) α E β T S ( − G I − ) T ( − α I T k E A E T )( − G E − ) T (cid:1) = ρ (cid:0) α E (cid:0) α I T ( − G I − )( β S ) (cid:1) T ( − A E k E ) T ( − G E − ) T (cid:1) (39)The above expression for R can be simplified as follows.First, recall the expected reward equation eq. (5). The first part of the above expression, α I T ( − G I − )( β S ) (40)is exactly an expected reward accrued prior to reaching an absorbing state, for a reward Markovprocess associated with a phase-type distribution parameterized by vector α I and matrix G I = A I − µ I k I × k I with reward rate vector β S . This phase-type distribution describes the duration oftime spent in the infected state, and is formally the phase-type distribution that is the minimum (seeeq. (6)) of an exponential distribution with rate µ and a phase-type distribution with parameters α I and A I . Therefore, eq. (40) is the expected number of new infectious created by an averageinfectious individual over the duration of the infectious period. Since this quantity is a 1 × R = α I T ( − G I − )( β S ) , (41)then it follows that R = R ρ (cid:18) α E ( − A E k E ) T ( − G E − ) T (cid:19) . (42)Let Q E denote the transient block of the transition probability matrix for the embedded jumpprocess associated with the phase-type distribution parameterized by α E and G E , and denote thevector of dwell time rates for the E sub-states by r E = − diag( G E ). Then by eq. (8) R = R ρ (cid:18) α E ( − A E k E ) T (( D r E ( I − Q E )) − ) T (cid:19) = R ρ (cid:18) α E (( I − Q E ) − D r E − ( − A E k E )) T (cid:19) (43)Let r i denote the rate for the i th sub-state of E when µ = 0, i.e., it is the negative of the loss ratein the i th diagonal entry of A E . Then the i th entry of the column vector (see eq. (10), but notethe definition of r E ) 17 reprint - August 18, 2020 D r E − ( − A E k E ) (44)is the probability of leaving state E from the i th sub-state of E, times the rate r i divided by the i th entry in r E , which is r i + µ . In other words, it’s the probability of leaving state E from the i th sub-state of E, times the probability of entering state I (given by r i / ( r i + µ )) as opposed toentering the deceased state (not tracked in the model), i.e., it is the column of matrix R in eq. (12)corresponding to the infected (versus deceased) state. Thus, according to eq. (12), the vector U E = (( I − Q E ) − D r E − ( − A E k E ) (45)is the column vector in which the i th entry is the probability of surviving the exposed state (E)and entering the infectious state (I), given that the initial state was the i th transient state.Therefore, it follows that R = R ρ (cid:18) α E ( U E ) T (cid:19) = R α E · U E (46)where the scalar product ρ (cid:0) α E ( U E ) T (cid:1) = α E · U E is the sum of each initial probability, α i , timesthe exposed-state survival probability, U E i . Thus, P E → I = α E · U E is the overall probability ofsurviving the exposed state E and transitioning to the infectious state I, and therefore R = (cid:0) α I T ( − G I − )( β S ) (cid:1)| {z } R ,new α E · (( − G E − )( − A E k E )) | {z } P E → I . (47) (cid:4) R for the SEIRS model with Erlang latent and infectious periods Observe that the SEIRS model with Erlang latent and infectious periods, and exponentially dis-tributed duration of immunity (eqs. (18)) is a special case of eqs. (31). Therefore, the followingcorollary to Theorem 4 proves the previous conjecture about eqs. (18) and eq. (22).
Corollary 2.
The basic reproduction number R for eqs. (18) – the SEIRS model with mortalityrate µ , exponentially distributed duration of immunity, and Erlang latent and infectious periodswith respective means τ E and τ I and arbitrary integer-valued shape parameters k E , k I ≥ – isgiven by R = β S (cid:18) r E r E + µ (cid:19) k E (cid:18) r I + µ + r I ( r I + µ ) + r I ( r I + µ ) + · · · + r k I − I ( r I + µ ) k I (cid:19) . (48) where r E = k E /τ E and r I = k I /τ I . The DFE is locally asymptotically statble if R < andunstable if R > .The first term, β S , is the expected rate of new infections per unit time per infected individualadded to a population at the DFE. The second term is the probability of surviving the exposedclass (E) and entering the infected class (I). The third term is the expected vale of the random reprint - August 18, 2020 variable T I – the duration of time spent in the infected class – which follows a Coxian distributiondefined as the minimum of an exponential distribution with rate µ and the Erlang infectious perioddistribution.Proof. Consider eqs. (31). Let the birthrate Λ( S, x , y , R ) = Λ be constant, and the force ofinfection λ ( t ) = β P k I i =0 I i ( t ). As in Hurtado and Richards (2020a), the Erlang distributed latentperiod is phase-type distributed with α E = ... and A E = − r E r E · · · − r E r E . . . ... . . . . . . . . . . . . . . . − r E r E · · · − r E . (49)The Erlang distributed infectious period can be similarly parameterized using α I = ... and A I = − r I r I · · · − r I r I . . . ... . . . . . . . . . . . . . . . − r I r I · · · − r I (50)For the exponentially distributed duration of immunity, α R = [1] and A R = [ − ǫ ].Let G I denote the matrix given by A I with additional − µ terms along the diagonal ( G I = A I − µ I and similarly G E = A E − µ I ). Then, by Theorem 6 in Appendix B, the parameter pairs α E , G E and α I , G I define two Coxian distributions describing the effective latent period and infectiousperiod distributions after accounting for losses from deaths. Observe that since all β i = β then thevector β = β . Then by Theorem 4 it follows that R is given by the product of R = α I T ( − G I − )( β S )= βS α I T ( − G I − ) (51)and P E → I = α E · (( − G E − )( − A E k E ))= α E · (( − G E − )([0 · · · T r E )) . (52)If we let T I denote the time spent in state I, thenE( T I ) = α I T ( − G I − ) = (cid:18) r I + µ + r I ( r I + µ ) + r I ( r I + µ ) + · · · + r k I − I ( r I + µ ) k I (cid:19) . (53)19 reprint - August 18, 2020 Since − G E − is of the form (see Appendix B) − G E − = 1 r E + µ r E r E + µ (cid:0) r E r E + µ (cid:1) · · · (cid:0) r E r E + µ (cid:1) k E − r E r E + µ . . . (cid:0) r E r E + µ (cid:1) k E − ... . . . . . . . . . . . . . . . r E r E + µ · · · (54)then ( − G E − ) ... r E = (cid:0) r E r E + µ (cid:1) k E (cid:0) r E r E + µ (cid:1) k E − ... (cid:0) r E r E + µ (cid:1) r E r E + µ (55)and thus, P E → I = α E · ( − G E − ) ... r E = ... · (cid:0) r E r E + µ (cid:1) k E (cid:0) r E r E + µ (cid:1) k E − ... (cid:0) r E r E + µ (cid:1) r E r E + µ = (cid:18) r E r E + µ (cid:19) k E . (56)It follows that R = β S (cid:18) r E r E + µ (cid:19) k E (cid:18) r I + µ + r I ( r I + µ ) + r I ( r I + µ ) + · · · + r k I − I ( r I + µ ) k I (cid:19) . (57) (cid:4) Here we have shown how reproduction numbers for ODE models of arbitrary finite dimension canbe derived and interpreted using the generalized linear chain trick (GLCT; Hurtado and Kirosingh,2019; Hurtado and Richards, 2020b) in conjunction with the standard linear chain trick (LCT),continuous time Markov chain (CTMC) theory, the associated theory of phase-type (i.e., CTMCabsorption time) distributions (Bladt and Nielsen, 2017b; Reinecke, Bodrog, and Danilkina, 2012),and the next generation operator method (Diekmann, Heesterbeek, and Roberts, 2009; Diekmann,Heesterbeek, and Metz, 1990; van den Driessche and Watmough, 2002). This approach can yieldgeneral expressions for basic reproduction numbers ( R ) in epidemic models, like the SEIR modelwith Erlang latent and infectious periods, but can also be used to find similar threshold quantitieslike population reproduction numbers in multispecies ecological models.20 reprint - August 18, 2020 The success of this approach relies on two related features of mean field ODE models that were(or could have been) derived using the GLCT. The first is that formulating the model (or familyof models) in a matrix-vector form consistent with the GLCT yields equations that are agnosticof the actual number of individual ODE equations. That fact can then be exploited during theapplication of the next generation operator approach for computing reproduction numbers. Thesecond feature of such models is that the resulting general reproduction number expressions can besimplified and interpreted through the lens of phase-type distributions and Markov chain theory,allowing for the kinds of insights that are often important in applications.We have also introduced a novel SEIRS family of contagion models with phase-type distributedlatent periods, infectious periods, and durations of immunity. This general model framework couldbe further generalized (e.g., to consider other assumptions about the functional form of the force ofinfection, or less simplistic birth and death processes) but nonetheless there are many existing SEIR-type models used in applications that are special cases of this more general model. These include,e.g., models currently being applied to the ongoing SARS-CoV-2 pandemic. To illustrate the utilityof this approach for mathematical analyses of such models, we used these results to prove the formof a general R expression for an SEIR model of arbitrarily many (finite) dimensions derived usingthe standard LCT (i.e., it assumes Erlang distributed latent and infectious periods).Not only can the phase-type form of such models be an useful analytically, but it can be used tospeed up computations in some cases, as shown in Hurtado and Richards (2020a). Also, since thematrix and vector parameters of a phase-type distribution can be estimated from data (Horv´ath andTelek, 2017, 2020), the GLCT can be used to build approximate empirical dwell time distributionsinto ODE models (e.g., see Kim, Byun, and Jung (2019) for an example where Coxian distributedlatent periods are derived from data and used to formulate a similar SEIR model to the moregeneral one introduced above).It is relatively uncommon to see the next generation operator approach used to compute populationreproduction numbers in multispecies population models, which may be a reflection of the somewhatlow-dimensional nature of these models (or the single-species components of such models). It mightalso be a consequence of these quantities being derived during a more direct equilibrium stabilityanalysis of different steady states in such models. However, as such models are extended to includeindividual heterogeneity in the form of different discrete individual types or sub-populations, orthere is a discrete stage structure added to these populations a la the LCT or GLCT, the approachoutlined above may prove to be a useful addition to the modelers’ toolkit.It is also worth noting that there are other threshold quantities, such as type reproduction numbers(Heesterbeek et al., 2015; Heesterbeek and Roberts, 2007; Roberts and Heesterbeek, 2003; Shuai,Heesterbeek, and van den Driessche, 2012), and methods for finding these quantities might alsobe able to exploit the matrix-vector structure of GLCT-derived models in ways that yield generalresults comparable to the analyses and results presented here. Furthermore, R can be obtainedby other means, which sometimes yield different R expressions (but with the same thresholdbehavior near 1), especially for age structured models with differential sensitivity and multi-host(e.g., vector-born) diseases (Heffernan, Smith, and Wahl, 2005; Hyman and Li, 2000, 2005; Yang,2014). Nonetheless, analyses that employ other methods to find reproduction numbers may stillbenefit from using the matrix-vector form of models consistent with the GLCT, as well as theapplication of CTMC and phase-type distribution theory in the interpretation of the results.Mean field ODE models are used throughout the sciences, and those model equations reflect muchof the structure associated with their corresponding (often unwritten) individual-based, stochastic21 reprint - August 18, 2020 model analogs. The GLCT provides one way to more clearly see these deeper connections betweenmean field ODEs and their stochastic counterparts. Here we have attempted to highlight howthe GLCT can be a useful tool for developing and analyzing mean field ODE models by leveragingexisting theory, statistical tools, and intuition from CTMCs and phase-type distributions. We hopeothers find these techniques to be useful for the development, analysis, and application of meanfield dynamical systems models. Acknowledgments:
The authors thank Dr. Deena Schmidt, the Mathematical Biology lab groupat UNR, Dr. Marisa Eisenberg, and Dr. Michael Cortez for conversations and comments thatimproved this manuscript.
Funding:
This work was supported by a grant awarded to PJH by the Sloan Scholars MentoringNetwork of the Social Science Research Council with funds provided by the Alfred P. Sloan Foun-dation; and this material is based upon work supported by the National Science Foundation underGrant No. DEB-1929522.
Disclosure statement:
The authors declare that they have no conflict of interest.
Appendix A Expected Reward For Phase-Type Distributions
Theorem 5.
Suppose g i is the reward rate for time spent in the i th state ( T i ) such that thecumulative reward is g i T i . Let Y be the reward accrued up to hitting the absorbing state in a CTMCwith n transient states and a single absorbing state, with corresponding phase-type distribution givenby α and A . Then E ( Y ) = α T ( − A − ) g . Proof.
Following Theorem 3.1.16 in Bladt and Nielsen (2017a), E ( Y ) = n X i =1 α i E (cid:18) n X j =1 T j r j (cid:12)(cid:12)(cid:12)(cid:12) initial state = i th state (cid:19) = n X i =1 α i n X j =1 u ij r j = α T Ug = α T ( − A − ) g . (cid:4) Appendix B Minimum of an Exponential and (Generalized) Er-lang Distribution
Theorem 6.
Let random variable T be the minimum of an exponentially distributed random vari-able with rate µ , and a generalized Erlang distributed (hypoexponential) random variable (i.e., thesum of k independent exponentially distributed random variables) with parameters λ , . . . , λ k . Then T obeys a Coxian distribution parameterized by reprint - August 18, 2020 α = ... and G = − λ − µ λ · · · − λ − µ λ . . . ... . . . . . . . . . . . . . . . − λ k − − µ λ k − · · · − λ k − µ . (B1) and the mean value of T is given byE ( T ) = 1 λ + µ + 1 λ + µ (cid:18) λ λ + µ (cid:19) + 1 λ + µ (cid:18) λ λ + µ λ λ + µ (cid:19) + · · · + 1 λ k + µ k − Y i =1 λ i λ i + µ (B2) Proof.
The generalized Erlang (hypoexponential) and exponential distributions can each be writtenas phase-type distributions, with parameterization α = ... and A = − λ λ · · · − λ λ . . . ... . . . . . . . . . . . . . . . − λ k − λ k − · · · − λ k (B3)where A is a k × k matrix, and α = [1] and A = [ − µ ] . (B4)By eq. (6), this minimum random variable T is also phase-type distributed, and more specifically,it has parameters α ⊗ α = α ⊗ [1] = α (B5)and A ⊕ A = A ⊗ [1] + I k × k ⊗ [ − µ ] = G . (B6)By the structure of G , one can be more specific and say that T has a Coxian distribution. For thisparticular case, direct inspection confirms that − G − = λ + µ λ + µ λ λ + µ λ + µ λ λ + µ λ λ + µ · · · λ k + µ Q k − i =1 λ i λ i + µ λ + µ λ + µ λ λ + µ . . . λ k + µ Q k − i =2 λ i λ i + µ ... . . . . . . . . . . . . . . . λ k − + µ λ k + µ λ k − λ k − + µ · · · λ k + µ (B7)23 reprint - August 18, 2020 and thus, since α is a 1 followed by k − T ) = α ( − G − ) , is the sum of the entries in the first row of − G − . (cid:4) reprint - August 18, 2020 References
Altiok, Tayfur (1985). ‘‘On the Phase-Type Approximations of General Distributions’’. In:
IIETransactions doi : .Bladt, Mogens and Bo Friis Nielsen (2017a). Matrix-Exponential Distributions in Applied Proba-bility . Springer US. doi : .– (2017b). ‘‘Phase-Type Distributions’’. In: Matrix-Exponential Distributions in Applied Probabil-ity . Springer US, pp. 125–197. doi : .Bonzi, B. et al. (2010). ‘‘Stability of differential susceptibility and infectivity epidemic models’’. In: Journal of Mathematical Biology doi : .Brauer, Fred and Carlos Castillo-Chavez (2011). Mathematical Models in Population Biology andEpidemiology . 2nd (2012). Texts in Applied Mathematics (Book 40). Springer-Verlag.Brauer, Fred, Pauline van den Driessche, and Jianhong Wu, eds. (2008).
Mathematical Epidemiol-ogy . Lecture Notes in Mathematics: Mathematical Biosciences Subseries. Springer-Verlag BerlinHeidelberg. doi : .Cumani, Aldo (1982). ‘‘On the canonical representation of homogeneous markov processes mod-elling failure - time distributions’’. In: Microelectronics Reliability doi : .Diekmann, O., J. A. P. Heesterbeek, and M. G. Roberts (2009). ‘‘The construction of next-generation matrices for compartmental epidemic models’’. In: Journal of The Royal SocietyInterface doi : .Diekmann, O. and J.A.P. Heesterbeek (2000). Mathematical epidemiology of infectious diseases:Model building, analysis and interpretation . Wiley Series in Mathematical and ComputationalBiology. John Wiley & Sons, LTD, New York.Diekmann, O., J.A.P. Heesterbeek, and J.A.J. Metz (1990). ‘‘On the definition and the computationof the basic reproduction ratio R in models for infectious diseases in heterogeneous populations’’.In: Journal of Mathematical Biology doi : .Eftimie, Raluca et al. (2011). ‘‘Multi-Stability and Multi-Instability Phenomena in a Mathemat-ical Model of Tumor-Immune-Virus Interactions’’. In: Bulletin of Mathematical Biology doi : .Heesterbeek, H. et al. (2015). ‘‘Modeling infectious disease dynamics in the complex landscape ofglobal health’’. In: Science doi : .Heesterbeek, J.A.P. (2002). ‘‘A Brief History of R and a Recipe for its Calculation’’. In: ActaBiotheoretica doi : .Heesterbeek, J.A.P. and M.G. Roberts (2007). ‘‘The type-reproduction number T in mod-els for infectious disease control’’. In: Mathematical Biosciences doi : .Heffernan, J.M, R.J Smith, and L.M Wahl (2005). ‘‘Perspectives on the basic reproductive ratio’’.In: Journal of The Royal Society Interface doi : .25 reprint - August 18, 2020 Hethcote, Herbert W. (1976). ‘‘Qualitative analyses of communicable disease models’’. In:
Mathe-matical Biosciences doi : .– (2000). ‘‘The Mathematics of Infectious Diseases’’. In: SIAM Review doi : .Hews, Sarah et al. (2009). ‘‘Rich dynamics of a hepatitis B viral infection model with lo-gistic hepatocyte growth’’. In: Journal of Mathematical Biology doi : .Hilker, Frank M. and Kirsten Schmitz (2008). ‘‘Disease-induced stabilization of preda-tor–prey oscillations’’. In: Journal of Theoretical Biology doi : .Horv´ath, Andr´as, Marco Scarpa, and Mikl´os Telek (2016). ‘‘Phase Type and Matrix Expo-nential Distributions in Stochastic Modeling’’. In: Principles of Performance and ReliabilityModeling and Evaluation: Essays in Honor of Kishor Trivedi on his 70th Birthday . Ed. byLance Fiondella and Antonio Puliafito. Cham: Springer International Publishing, pp. 3–25. doi : .Horv´ath, G´abor and Mikl´os Telek (2017). ‘‘BuTools 2: A rich toolbox for Markovian performanceevaluation’’. English. In: ValueTools 2016 - 10th EAI International Conference on PerformanceEvaluation Methodologies and Tools . Association for Computing Machinery, pp. 137–142. doi : .– (2020). BuTools V2.0 . http://webspn.hit.bme.hu/~telek/tools/butools/doc/ . Accessed:2020-05-15.Horv´ath, G´abor et al. (2012). ‘‘Efficient Generation of PH-Distributed Random Variates’’. In: Pro-ceedings of the 19th international conference on Analytical and Stochastic Modeling Techniquesand Applications . Ed. by Khalid Al-Begain, Dieter Fiems, and Jean-Marc Vincent. Berlin, Hei-delberg: Springer Berlin Heidelberg, pp. 271–285. doi : .Hurtado, Paul J., Spencer R. Hall, and Stephen P. Ellner (2014). ‘‘Infectious disease in consumerpopulations: Dynamic consequences of resource-mediated transmission and infectiousness’’. In: Theoretical Ecology doi : .Hurtado, Paul J. and Adam S. Kirosingh (2019). ‘‘Generalizations of the ‘Linear Chain Trick’:incorporating more flexible dwell time distributions into mean field ODE models’’. In: Journalof Mathematical Biology doi : .Hurtado, Paul J. and Cameron Richards (2020a). Building Mean Field State Transition ModelsUsing The Generalized Linear Chain Trick and Continuous Time Markov Chain Theory . arXiv: .– (2020b).
Time Is Of The Essence: Incorporating Phase-Type Distributed Delays And DwellTimes Into ODE Models . arXiv: .Hyman, James M. and Jia Li (2000). ‘‘An intuitive formulation for the reproductive number for thespread of diseases in heterogeneous populations’’. In:
Mathematical Biosciences doi : . 26 reprint - August 18, 2020 Hyman, James M. and Jia Li (2005). ‘‘The reproductive number for an HIV model with differentialinfectivity and staged progression’’. In:
Linear Algebra and its Applications doi : .Johnson, Tammi L., Erin L. Landguth, and Emily F. Stone (2016). ‘‘Modeling Relapsing DiseaseDynamics in a Host-Vector Community’’. In: PLOS Neglected Tropical Diseases doi : .Kermack, W. O. and A. G. McKendrick (1927). ‘‘A Contribution to the Mathematical Theory ofEpidemics’’. In: Proceedings of the Royal Society of London. Series A, Containing Papers of aMathematical and Physical Character
Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathe-matical and Physical Character
Proceedings of the Royal Society of London. Series A, ContainingPapers of a Mathematical and Physical Character
Bulletin of Mathemat-ical Biology doi : .– (1991b). ‘‘Contributions to the mathematical theory of epidemics—II. The problem of endemic-ity’’. In: Bulletin of Mathematical Biology doi : .– (1991c). ‘‘Contributions to the mathematical theory of epidemics—III. Further studies ofthe problem of endemicity’’. In: Bulletin of Mathematical Biology doi : .Kim, Sungchan, Jong Hyuk Byun, and Il Hyo Jung (2019). ‘‘Global stability of an SEIR epidemicmodel where empirical distribution of incubation period is approximated by Coxian distribution’’.In: Advances in Difference Equations doi : .Neuts, Marcel F. (1981). Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Ap-proach (Updated Edition) . Baltimore: The Johns Hopkins University Press.O’Cinneide, Colm Art (1991). ‘‘Phase-type distributions and invariant polytopes’’. In:
Advances inApplied Probability doi : .– (1993). ‘‘Triangular order of triangular phase-type distributions ∗ ’’. In: Communications in Statis-tics. Stochastic Models doi : .Reinecke, Philipp, Levente Bodrog, and Alexandra Danilkina (2012). ‘‘Phase-Type Dis-tributions’’. In: Resilience Assessment and Evaluation of Computing Systems . Ed. byKatinka Wolter et al. Berlin, Heidelberg: Springer Berlin Heidelberg, pp. 85–113. doi : .Reinecke, Philipp, Tilman Krau, and Katinka Wolter (2012). ‘‘Cluster-based fitting of phase-typedistributions to empirical data’’. In: Computers & Mathematics with Applications doi : .Resnick, Sidney I. (2002). Adventures in Stochastic Processes . Birkhuser Boston. doi : . 27 reprint - August 18, 2020 Roberts, M. G. and J. A. P. Heesterbeek (2003). ‘‘A new method for estimating the effort requiredto control an infectious disease’’. In:
Proceedings of the Royal Society of London. Series B:Biological Sciences doi : .– (2012). ‘‘Characterizing the next-generation matrix and basic reproduction number in eco-logical epidemiology’’. In: Journal of Mathematical Biology doi : .Shuai, Zhisheng, J. A. P. Heesterbeek, and P. van den Driessche (2012). ‘‘Extending the typereproduction number to infectious disease control targeting contacts between types’’. In: Journalof Mathematical Biology doi : .Smith, Hal (2010). An introduction to delay differential equations with applications to the lifesciences . Vol. 57. Springer Science & Business Media.van den Driessche, P. and James Watmough (2002). ‘‘Reproduction numbers and sub-thresholdendemic equilibria for compartmental models of disease transmission’’. In:
Mathematical Bio-sciences doi : .Wearing, Helen J, Pejman Rohani, and Matt J Keeling (2005). ‘‘Appropriate Models for the Man-agement of Infectious Diseases’’. In: PLOS Medicine doi : .Yang, Hyun Mo (2014). ‘‘The basic reproduction number obtained from Jacobian and next gener-ation matrices – A case study of dengue transmission modelling’’. In: Biosystems doi :10.1016/j.biosystems.2014.10.002