Linear Autonomous Compartmental Models as Continuous-Time Markov Chains: Transit-Time and Age Distributions

Linear compartmental models are commonly used in different areas of science, particularly in modeling the cycles of carbon and other biogeochemical elements. The representation of these models as linear autonomous compartmental systems allows different model structures and parameterizations to be compared. In particular, measures such as system age and transit time are useful model diagnostics. However, compact mathematical expressions describing their probability distributions remain to be derived. This paper transfers the theory of open linear autonomous compartmental systems to the theory of absorbing continuous-time Markov chains and concludes that the underlying structure of all open linear autonomous compartmental systems is the phase-type distribution. This probability distribution generalizes the exponential distribution from its application to one-compartment systems to multiple-compartment systems. Furthermore, this paper shows that important system diagnostics have natural probabilistic counterparts. For example, in steady state the system’s transit time coincides with the absorption time of a related Markov chain, whereas the system age and compartment ages correspond with backward recurrence times of an appropriate renewal process. These relations yield simple explicit formulas for the system diagnostics that are applied to one linear and one nonlinear carbon-cycle model in steady state. Earlier results for transit-time and system-age densities of simple systems are found to be special cases of probability density functions of phase-type. The new explicit formulas make costly long-term simulations to obtain and analyze the age structure of open linear autonomous compartmental systems in steady state unnecessary.


Introduction
Compartmental models are widely used to describe a range of biological and physical processes that rely on mass conservation principles (Anderson 1983;Jacquez and Simon 1993). In particular, most models that describe the carbon cycle can be generalized as compartmental models (Luo and Weng 2011;Sierra and Müller 2015;Rasmussen et al. 2016). These models are very important for predicting interactions between ecosystems and the global climate. However, the predictions of these models diverge widely (Friedlingstein et al. 2006(Friedlingstein et al. , 2014 since models have different structures and parameter values. To compare diverse models, key quantities such as system age and transit time can be considered. These quantities provide relevant information about the time scales at which material is cycled in compartmental systems and facilitate comparisons among different model structures and parameterizations. Formulas for their means were recently provided by Rasmussen et al. (2016) for linear nonautonomous compartmental systems. For autonomous systems, that is, systems with constant coefficients, the existence of explicit solutions raises hope for explicit formulas not only for the means, but also for the densities of system age and transit time.
A first attempt in this direction was started by Nir and Lewis (1975) who established formulas for transit-time and age densities in dependence on a not explicitly known system response function. This response function was the basis for Thompson and Randerson (1999) to compute the desired densities numerically by long-term simulations in two carbon-cycle models. Impulsive inputs and how systems respond to them were later investigated by Manzoni et al. (2009) to obtain explicitly the system response function for a set of carbon-cycle models of very simple structure. Despite the simple structure of these systems, the derivation of the according density formulas is involved, because the system's response needs to be transformed from the phase domain to the Laplace domain and back. In this manuscript, the dynamics of linear autonomous compartmental systems are considered as a stochastic process to show that the impulse response function is the probability density function of a phase-type distribution. To this end, open linear autonomous compartmental systems are related to absorbing continuous-time Markov chains.
This manuscript is organized as follows. In Sect. 2, linear autonomous compartmental systems are introduced, and the idea of looking at them from the perspective of one single particle is framed. The transit time of this particle through the system coincides with the absorption time of a continuous-time Markov chain. A short presentation of the basic ideas of continuous-time Markov chains follows in Sect. 3. The probability distribution of the absorption time is shown to be an example of a phase-type distribution, and its cumulative distribution function, probability density function, and moments are computed. In Sect. 4, concepts from renewal theory are used to construct a regenerative process, which is the basis for the computation of the system-and compartment-age densities. The relationship of Markov chains and lin-ear autonomous compartmental systems is established in Sect. 5, where the previous probabilistic results are used to obtain simple explicit formulas for the densities of the system age, compartment age, and transit time. In Sect. 6, the derived formulas are applied to two autonomous terrestrial carbon-cycle models. Both models are considered in steady state, which allows the treatment of a linear model (Emanuel et al. 1981), as well as a nonlinear one (Wang et al. 2014) within the proposed framework. Relations to earlier work and some other aspects are discussed in Sect. 7, and in Sect. 8 conclusions are presented. "Appendix A" contains the proof of the main theorem from Sect. 4, and "Appendix B" covers some examples for systems with simple structure.

Linear Autonomous Compartmental Systems
Compartmental differential equations are useful tools to describe flows of material between units called compartments under the constraint of mass conservation. Following Jacquez and Simon (1993), a compartment is an amount of some material that is kinetically homogeneous. This means that the material of a compartment is at all times homogeneous; any material entering the compartment is instantaneously mixed with the material already there. Hence, compartments are always well-mixed.
Consider the d-dimensional linear system of differential equations d dt This equation system represents how mass flows at time t = 0 into a set of compartments, flows and cycles among the compartments, and eventually leaves the system. Therefore, the vector x 0 of initial contents and the constant input vector u are required to have nonnegative entries only. The law of conservation of mass imposes restrictions on the d × d-matrix B.
If now B is a compartmental matrix, (1) is called a linear autonomous compartmental system. The term autonomous refers to the fact that both B and u are independent of time t. A detailed treatment of such systems can be found in Anderson (1983) and Jacquez and Simon (1993). For t ≥ 0 the vector x(t) ∈ R d represents the contents of the different compartments at time t. The off-diagonal entries b i j of B are called fractional transfer coefficients. They are the rates at which mass moves from compartment j to compartment i. For j = 1, . . . , d, the nonnegative value z j = − d i=1 b i j is the rate at which mass leaves the system from compartment j. If at least one of these output or release rates is greater than zero, system (1) is called open, otherwise it is called closed. This paper is concerned with open systems only, since closed systems with positive input accumulate mass indefinitely.

The One-Particle Perspective
Looking at the behavior of the entire system, the initial content x 0 distorts what can be seen happening to the system content. There are two ways to get rid of this disturbing influence of the initial content. One way is to consider the system after it has run for an infinite time such that all the initial content has left. Another way is to look at one single particle that just arrives at the system through new input u.
Since all compartments are well-mixed, this particle's way through the system is not influenced by the presence or absence of any other particles. It enters the system at a compartment according to u and then, at each time step, whether it stays or moves on is decided on the basis of its current position and its schedule. If the decision is to move on, then it can move to another compartment or leave the system, depending only on the connections of the current compartment. The particle follows a schedule and a map given by the matrix B. Its diagonal entries govern how long the particle stays in a certain compartment, and the off-diagonal entries provide the connections to other compartments. By leaving the system, the particle finishes a cycle and starts a new one by reentering the system. At each cycle, the sequence of compartments to which the particle belongs at successive time steps constitutes a stochastic process called discrete-time Markov chain. Letting the size of the time steps tend to zero, the particle's future becomes continuously uncertain. The path of the particle traveling through the system is then represented by a continuous-time Markov chain (Norris 1997).
In Sect. 2.2, the structure and properties of continuous-time Markov chains are introduced, and the reader should always have in mind the traveling particle. When the Markov chain changes its state from j to i, the particle is considered to move from compartment j to compartment i. When the Markov chain is absorbed, the particle leaves the system. The time that elapses from the moment the particle enters the system to the moment of its exit is called the particle's transit time. To get a grasp on it, from the beginning of the cycle only a look into the future is necessary.
While the particle is still in the system, the time that has passed since the particle's entry is called its system age. If the age of the particle is considered at a random time, a look into the past of the particle is needed, not knowing when it entered the system the last time. Consequently, it has to be assumed that the particle has ran through the system already infinitely often. Otherwise the existence of a maximum age of the particle would be implied, but any choice of this maximum value would be arbitrary and ill-founded. Therefore, a particular regenerative process will be considered: a sequence of absorbing continuous-time Markov chains.

From One Particle to All Particles in the System
Assume that the system at time t ≥ 0 contains n particles with system ages represented by n independent and identically distributed random variables. Hence, the system age A k of particle k is assumed to have the cumulative distribution function F A for all k = 1, 2, . . . , n. The share of particles with system age less than or equal to y ≥ 0 equals Here, 1 S (s) denotes the indicator function. It equals 1 if s ∈ S and it equals 0 otherwise. The function F n is called empirical distribution function of A 1 , A 2 , . . . , A n , and from the Glivenko-Cantelli theorem (Dudley 1999) it is known that it converges almost surely uniformly to F A if the number n of particles goes to infinity. Consequently, if the system contains an infinite amount of particles, the share of the total system content x(t) that has system age less than or equal to y is where f A is the common probability density function of the particles' ages.

Continuous-Time Markov Chains
Markov chains are the most important examples of random processes. Given their simple structure and high diversity, they are applied to many different scientific problems. In particular, they are the simplest mathematical models for random phenomena evolving in time. Their characteristic property is that they retain no memory on the states of the system in the past. Only the current state of the process can influence where it goes next. If the process can assume only a finite or countable set of states, it is called a Markov chain. Discrete-time Markov chains are usually defined on a set of integers, whereas continuous-time Markov chains live on a subset of the real line. In the following, the basic theory of continuous-time Markov chains is introduced along the line of Norris (1997), from which also most unproven results are taken. A vector λ ∈ R d is called a distribution if all its components are nonnegative and sum to one. Furthermore, a matrix P = ( p i j ) i, j∈J is called stochastic on a finite set J if every column ( p i j ) i∈J is a distribution. Note that in standard literature on Markov chains row sums are considered. Here, column sums are used instead, because thereby the connection to compartmental matrices will become more obvious later. The definition of a Markov chain can best be formalized in terms of its corresponding transition rate matrix. A transition rate matrix is a special compartmental matrix where all columns sum to zero.
Definition 2 Let J be a finite state space. A transition rate matrix on J is a matrix Q = (q i j ) i, j∈J satisfying the conditions (i) 0 ≤ q i j for all i = j, (ii) 0 ≤ −q j j < ∞ for all j, (iii) i∈J q i j = 0 for all j.
Each off-diagonal entry q i j can be interpreted as the rate of going from state j to state i and the diagonal entries q j j as the rate of leaving state j.
Definition 3 Let Q be a transition rate matrix and λ a distribution on a finite state space J . A stochastic process X = (X t ) t≥0 is a continuous-time Markov chain on J with initial distribution λ and transition rate matrix Q if (i) P(X 0 = j) = λ j for all j ∈ J ; (ii) for all n = 0, 1, 2, . . ., all times 0 ≤ t 0 ≤ t 1 ≤ · · · ≤ t n+1 , and j 0 , j 1 , . . . , j n+1 ∈ J , where for any quadratic matrix M the expression e M denotes the matrix exponential Since Q does not depend on time, X is called homogeneous. Property (ii) states that the future evolution of a Markov process depends only on its current state and not on its history. This is called Markov property. Let X be a continuous-time Markov chain on J with initial distribution λ and transition rate matrix Q. Then, for i, j ∈ J , the probability of being in state i at time t having started in state j is equal to and the probability of being in state i at time t (conditioned on the initial distribution λ) is

Absorbing Continuous-Time Markov Chains
Consider a continuous-time Markov chain X = (X t ) t≥0 with a special structure. Its finite state space J is assumed to be equal to {1, 2, . . . , d, d + 1} for some natural number d ≥ 1, and its transition rate matrix has the shape Here, 0 is the d-dimensional column vector containing only zeros. Let S := {1, 2, . . . , d} ⊆ J . The d × d-matrix B = (b i j ) i, j∈S is assumed to meet the requirements (i) and (ii) of a transition rate matrix, but instead of property (iii) of Definition 2 it fulfills only the weaker condition i∈S b i j ≤ 0 for all j ∈ S.
Since Q is required to be a transition rate matrix, the vector z ∈ R d must contain the missing parts to make the columns sum to zero. Consequently, z j = − i∈S b i j , which in matrix notation is written as where 1 T denotes the d-dimensional row vector comprising ones. This means that the z j are nonnegative and denote the transition rates from j to d +1. The (d +1)st column of Q contains only zeros. Consequently, the process X cannot change its state anymore once it has reached state d + 1. For that reason, d + 1 is called the absorbing state of X . The trivial case in which the process starts in its absorbing state is excluded by considering only initial distributions λ with λ d+1 = 0. Then the new initial distribution of X is defined as A standard linear algebra argument shows that where the asterisk * is a place holder for a d-dimensional row vector. This means for i, j ∈ S that

The Absorption Time
The absorption time is a random variable that tells the time when the absorbing state is reached. Another name for it is hitting time of the absorbing state. It is an interesting question whether the absorbing state will always be reached no matter in which state the process starts.
The following lemma (Neuts 1981, Lemma 2.2.1) provides an answer to this question.

Lemma 1 If B is nonsingular, then the absorption time T is finite with probability one.
Definition 4 A Markov chain X on J = {1, 2, . . . , d, d + 1} whose transition rate matrix Q has the structure of Eq.
(2) and that will eventually be absorbed to d + 1 with probability one is called absorbing. The matrix B is called its transition rate matrix, S = {1, 2, . . . , d} its state space, and d + 1 its absorbing state. If λ denotes the initial distribution of X , then β as defined in Eq. (4) is called the initial distribution of the absorbing chain.
The following results related to absorbing continuous-time Markov chains will be used in the forthcoming.
Remark 1 For an absorbing continuous-time Markov chain eventual absorption is certain, independent of the initial state. Hence, Remark 2 Since Q is a transition rate matrix, e t Q is stochastic. This makes all of its entries and hence all entries of e t B nonnegative. From Eq. (6) it follows that also all entries of −B −1 are nonnegative.
From now on, let B be the nonsingular transition rate matrix of an absorbing continuous-time Markov chain X = (X t ) t≥0 on S = {1, 2, . . . , d} with initial distribution β and absorbing state d + 1. As it turns out, the absorption time T follows a phase-type distribution that depends on the initial distribution β and on the transition rate matrix B.

Probability Distributions of Phase-Type
Phase-type distributions constitute a highly versatile class of probability distributions and are closely related to the solutions of systems of linear differential equations with constant coefficients. As mixtures of exponential distributions, they generalize the Erlang, hypoexponential, and hyperexponential distribution. These are widely used in queuing theory and the field of renewal processes. An introduction to phase-type distributions and the unifying matrix formalism used in this paper can be found in Neuts (1981).
At time t, the cumulative distribution function F T (t) = P(T ≤ t) of the absorption time T is equal to the probability of not being in any of the states j ∈ S. Consequently, using Eq. (5) leads to A probability distribution according to this density is called phase-type with initial distribution β and transition rate matrix B. Neuts (1981) introduced the notation T ∼ PH(β, B).

Moments of the Phase-Type Distribution
Let T ∼ PH(β, B) with nonsingular B and z T = −1 T B. Lemma 2 can be used to calculate the nth moment of the phase-type distribution. Repeated integration by parts leads to In particular, the first moment of the phase-type distribution is where v := d j=1 |v j | denotes the norm of a vector v ∈ R d . The last equality holds since both the matrix −B −1 and the vector β are nonnegative. For −B −1 this follows from Remark 2 and for β because it is a distribution.
Note that the variance, that is the second central moment, of the phase-type distribution can be obtained from σ 2

A Closure Property of the Phase-Type Distributions
The class of phase-type distributions has several closure properties, one of which is given by the following lemma (Neuts 1981, Theorem 2.2.3).

Lemma 3 If T ∼ PH(β, B), then a random variable with cumulative distribution function
and transition rate matrix B.

Occupation Time
Before its absorption, X takes on different states j ∈ S. Its occupation time of state j is defined as the time the process spends in state j. It is given as the nonnegative random variable Using E[1 X t = j ] = P(X t = j) and Eq. (6), its expected value can be computed to In particular, (−B −1 ) i j is the mean time the Markov chain is in state i before absorption, under the condition that it started in state j.

The Last State Before Absorption
Let E denote the state from which X jumps to the absorbing state. From the homogeneity of X follows where O j is the occupation time of state j.

A Regenerative Process
Assume that an absorbing continuous-time Markov chain is restarted immediately every time it hits the absorbing state. The goal of this section is to determine the distribution of the age of this new process at a random time τ drawn from the positive half-line (technical details on τ are presented in "Appendix A"). Age means here the time that has passed since the last restart. First, the new situation is set up, then the distribution of the age at a random time is derived, and as a last step the age distribution at a random time under the condition of being in a fixed state is computed. Important tools are elements from renewal theory such as renewal and regenerative processes. A comprehensive treatise of this topic can be found in Asmussen (2003).

Definition of the Regenerative Process
Let X be an absorbing continuous-time Markov chain on a finite state space J with transition rate matrix B and initial distribution β. Now, the process X is stopped at its absorption time T , an immediate restart is executed, and this procedure is repeated over and over again. This gives an infinite sequence (X k ) k=1,2,... of independent and identically distributed cycle processes, where each cycle process behaves like the process X up to absorption. This sequence constitutes a regenerative process Z = The process that counts the number of restarts is called a renewal process.
A regenerative process is called alternating, if it has only the two possible states 1 and 0. Let Y = (Y t ) t≥0 be an arbitrary alternating regenerative process with the same cycle lengths as Z , that Such alternating processes together with the following theorem build the main tool for the computation of the age distributions.

Theorem 1
The probability that an alternating process Y with the same cycle lengths as Z is on at a random time τ ≥ 0 equals the ratio of expected on-time during the first cycle to the average cycle length. This leads to Despite the intuitive nature of this result, its proof requires quite some technical effort and it is postponed to "Appendix A".

The Age of the Regenerative Process
Consider the age process (A t ) t≥0 defined such that the random variable A t describes the time that has passed by t since the last restart. It is also called backward recurrence time of the renewal process. It can be considered the age of the regenerative process Z at time t. For y ≥ 0 and a randomly chosen time τ , the probability F A τ (y) = P(A τ ≤ y) is to be computed.
To obtain that probability, an alternating process Y is constructed that is on at time t if and only if the time elapsed since the last restart is less than or equal to y. Obviously, Write A instead of A τ and apply Theorem 1 to the alternating process Y to obtain the cumulative distribution function of the age of the regenerative process Z . Hence, this age is again phase-type distributed by Lemma 3, now with initial probability vector η.

The Age of the Regenerative Process in a Fixed State
Fix j ∈ S and let a j be the random variable that describes the age of the process Z given that it is in state j. For y ≥ 0 this means The conditional probability on the right hand side can be expressed as The numerator and the denominator can be computed using two alternating processes Y andỸ , respectively. The process Y is defined to be on at time t if and only if A t ≤ y and Z t = j, whereasỸ is defined to be on if and only if Z is in state j. Invoke Theorem 1 to Y andỸ to get and respectively. By plugging Eq. (12) for the numerator and Eq. (13) for the denominator in Eq. (11), the cumulative distribution function of a j is given by This means that the probability density function of the age of Z in state j is given by

Compartmental Systems and Markov Chains
Consider again the linear autonomous compartmental system (1) and the (d Owing to the properties of B and z = (z j ) j=1,2,...,d , the extended matrix Q meets Definition 2. Consequently, it is the transition rate matrix of a continuous-time Markov chain X = (X t ) t≥0 on the state space J = {1, 2, . . . , d, d + 1}. Assume that mass u ∈ R d comes into system (1) at a time t 0 . Since the system is linear, the way how the mass will be distributed through it can be modeled by a system without input and with initial value u given by Furthermore, the fractional transfer coefficients of this system are time independent, so the entire system can be shifted to the left and considered to have started at time t 0 = 0. The content of compartment j at time t ≥ 0 is then given byx j (t) = e t B u j . Let and β := (λ 1 , λ 2 , . . . , λ d ) T . Then the probability of the continuous-time Markov chain X with initial distribution β being in state j ∈ S := {1, 2, . . . , d} at time t is by an application of Eq. (5) given by P(X t = j) is the proportion of the initially present mass of system (15) that is in compartment j at time t. Consequently, the continuous-time Markov chain X describes the stochastic travel of a single particle through the compartmental system. When the traveling particle leaves the compartmental system, the process X jumps to the absorbing state d + 1.

Global Asymptotic Stability
Assume the compartmental matrix B to be nonsingular. Then the constant vector x * = −B −1 u is a steady-state solution or equilibrium of system (1), that is, d x * /dt = 0.
Any solution x of the linear system (1) has the form (Anderson 1983) If t → ∞, an application of Eq. (6) and Remark 1 leads to which means that all solutions of system (1) converge to x * , independently of their initial value x 0 . This makes the equilibrium x * = −B −1 u globally asymptotically stable and by the classic Lyapunov theorem (Engel and Nagel 2000, Theorem 2.10) all eigenvalues of B have negative real part. It is well known that system (1) is globally asymptotically stable, if B is strictly diagonally dominant, that is d i=1 b i j < 0 for all j ∈ S. However, here B is not strictly diagonally dominant (it is diagonally dominant, though), but it is a nonsingular compartmental matrix. This means that B is the transition rate matrix of an absorbing continuous-time Markov chain. For the absorbing state to be reached eventually, at least one rate z j = − d i=1 b i j of leaving the system from compartment j must be strictly greater than zero. This makes system (1) an open system (Jacquez and Simon 1993).

Steady-State Compartment Contents
For j ∈ S the steady-state content of compartment j is Consequently, the steady-state compartment content x * j is proportional to the expected occupation time of state j by the absorbing continuous-time Markov chain X .

Release from the System
For j ∈ S, the release of mass from compartment j to the environment is denoted by r j . As a function of time t, it can be computed as product of the rate z j of mass leaving compartment j toward the environment and the mass x j (t) contained in compartment j at time t. For a system in steady state, x j (t) = x * j remains constant and consequently r j = z j x * j remains constant as well. Let now r = (r 1 , r 2 , . . . , r d ) T ∈ R d be the vector of mass release or output from the system in steady state. Probabilistically, one would expect r j to be connected to the probability of the absorbing continuous-time Markov chain X to be absorbed through state j. From Eq. (9) the probability of j being the last state before absorption is P(E = j) = −z j B −1 β j . Use β = u/ u and x * = −B −1 u to get This means that the release through compartment j is proportional to the probability of X being absorbed through state j. Since absorption is certain, j∈S P(E = j) = 1, hence r = j∈S r j = u , reflecting that in steady state total system output equals total system input.

Age Distribution of the System in Steady State
Suppose that the total system content in steady state has an unknown age density f A . In Sect. 4.2, it was already considered a particle that travels over and over again through the system, and its age at a random time was computed by methods from renewal theory. The according regenerative process expresses the need for an infinite history. Consequently, it reflects the behavior of the compartmental system in steady state, when all initial mass has left the system. The age random variable A was shown to be phase-type distributed with transition rate matrix B and initial distribution Thus, the density is given by Eq. (7) with η replacing β. This gives the density The mean age of the system can be obtained from Eq. (8) as

Age Distribution of the Compartments in Steady State
Now, suppose that compartment j ∈ S has an unknown age density f a j . In Sect. 4.3 the age of the regenerative process Z was calculated under the condition that its state equals j. This is the age that a traveling particle has when it is in compartment j. From Eq. (14) and x * = −B −1 u, it follows that Integration by parts gives then for the mean age in compartment j ∈ S the equation Defining X * := diag(x * 1 , . . . , x * d ), this gives a vector valued function of age densities in the compartments f a (y) = (X * ) −1 e y B u, y ≥ 0, and its expectation is the vector of mean ages in the compartments. It is called mean age vector. For n ≥ 1 let a n denote the vector (a n 1 , . . . , a n d ) T . Multiple integration by parts leads to a formula for the vector whose components contain the nth moments of the compartment ages. It is given by Note that a computation of the average of the components of the mean age vector weighted by the respective steady-state compartment contents leads to the same equation for the mean system age as calculating the expected value of PH(η, B) as it was done to obtain Eq. (18).

Transit Time in Steady State
For compartmental systems two types of transit times can be considered (Nir and Lewis 1975). The forward transit time (FTT) at time t a is the time t a particle needs to travel through the system after it arrives at time t a . The backward transit time (BTT) at time t e gives the age y a particle has at the moment it leaves the system, that is, the time it needed to travel through the system given that it exits at time t e . For an autonomous system in steady state one would expect the two types of transit time to coincide.
Obviously, the FTT is the absorption time of the absorbing continuous-time Markov chain X . As soon as the traveling particle leaves the system, X hits its absorbing state. Hence, the FTT is phase-type distributed with initial distribution β and transition rate matrix B. This makes its density function f FTT (t) = z T e t B β, t ≥ 0, and its mean For the BTT the above constructed regenerative process can be considered. The BTT is the age of Z at the time of a restart. Hence, it follows the same distribution as the cycle length and is identically distributed to the absorption time. So, from the probabilistic point of view, FTT and BTT are obviously the same. But the density of the BTT can also be calculated as a weighted average of ages of particles that are leaving the system by Intuition does not betray in this case: FTT and BTT coincide in steady state. Furthermore, the transit time depends neither on the arrival time t a nor on the exit time t e .

A Linear Autonomous Compartmental System in Steady State
As an example, consider the global carbon-cycle model introduced by Emanuel et al. (1981), for which Thompson and Randerson (1999) numerically calculated age and transit-time distributions using the impulse response function approach. This model is therefore a good test-case for the derivations presented above. It comprises five compartments: non-woody tree parts x 1 , woody tree parts x 2 , ground vegetation x 3 , detritus/decomposers x 4 , and active soil carbon x 5 . The model is given by d dt where the input vector is given by u = 77 0 36 0 0 T Pg C year −1 and the compart- The numbers are chosen exactly as in Thompson and Randerson (1999). For comparison of the results, also refer to that paper. The input vector is expressed in units of petagrams of carbon per year (Pg C year −1 ), and the fractional transfer coefficients in units of year −1 . Because B is a lower triangular matrix with diagonal entries different from zero, it is nonsingular and all the earlier obtained formulas can be applied to system (21). The steady-state vector of carbon contents is x * = −B −1 u = 37.00 452.00 69.00 81.00 1121.00 T Pg C, The expected value E[T ] = 15.58 year is identical to the value found by Thompson and Randerson (1999), and the standard deviation of the transit time is σ T = 45.01 year. Furthermore, the system-age density is given by Its expectation E[A] = 72.83 year is very similar to that reported by Thompson and Randerson (1999)  From these density functions the system's and the compartments' contents can be plotted with respect to their age (Fig. 2). This gives useful information about the range of ages for each compartment and how they contribute to the system-age distribution. In comparison to the results of Thompson and Randerson (1999), the approach presented here not only provides mean values for ages and transit times, but also exact formulas for their probability distribution. In their approach, these authors obtained results that depended on the simulation time and, therefore, included numerical errors, something that can be easily avoided with the approach presented here.

A Nonlinear Autonomous Compartmental System in Steady State
Consider the nonlinear autonomous compartmental system d dt where B : R d → R d×d is a matrix-valued mapping. In this setup the fractional transfer coefficients are not constant but depend on the system's content. Consequently, the transition rates of the traveling particle are not constant. Assume now that system (22) is in a steady state x * . From d x * /dt = 0 follows that the compartment contents x j do not change and the mapping B turns into a matrix B with constant coefficients. Hence, assuming the nonlinear autonomous compartmental system (22) is in a steady state, it can be treated as a linear autonomous compartmental system and the entire theory developed for those systems can be applied to it.
As an example, consider the nonlinear two-compartment carbon-cycle model described by Wang et al. (2014). It is given by d dt

Fig. 2
Carbon content with respect to age in the model of Emanuel et al. (1981). Dotted lines indicate the mean age μ; the standard deviation is denoted by σ . The first panel is for the entire system, whereas the other panels correspond to the different compartments Here, B depends on x = C s C b T through λ's dependence on x, which is given by Steady-state formulas for the compartment contents can be computed as From Wang et al. (2014) take the parameter values F npp = 345 g C m −2 year −1 , μ b = 4.38 year −1 , ε = 0.39, and K s = 53,954.83 g C m −2 . Since the description of V s is missing in the original publication, let it be equal to 59.13 year −1 to approximately  Wang et al. (2014). Vertical dashed lines represent the mean μ; the standard deviation is denoted by σ . All panels show graphs for two different values of carbon use efficiency meet the given steady-state compartment contents C * s = 12, 650 g C m −2 and C * b = 50.36 g C m −2 .
With the given parameters, the steady-state transition rate matrix B = B(x * ) and the input vector u are given by The derived formulas are then used to calculate mean transit time and mean ages together with the according densities for different values of the model's parameters (Fig. 3). In particular, the effect of different values of the parameter ε on the steadystate ages and transit times can be explored. This parameter controls the proportion of carbon that is transferred from the substrate C s to the microbial biomass compartment C b , and it is commonly referred to as the carbon use efficiency. Interestingly, if the carbon use efficiency ε increases, the mean transit time and mean ages of the model at steady state decrease (Fig. 4), a behavior that at first glance appears counterintuitive. It can be explained by two opposing effects. On the one hand, an increase of carbon use efficiency keeps a higher fraction of carbon in the system due to lower respiration. This has an increasing effect on the transit time. On the other hand, a higher carbon use efficiency ε implies a lower steady-state content of compartment C s and a higher one of compartment C b . Consequently, from Eq. (23), an increasing value of is obtained. This value is the process rate of the C s compartment in steady state. The higher it is, the faster the particles travel through the system. The latter effect wins here and a decrease in the transit time can be observed. The graph of the mean system age for this model with ε = 0.39 (Fig. 4) lies directly on the one of the mean transit time. The huge difference in the compartments' steadystate contents causes very little difference in the initial distributions β and η for the traveling particle. This results in very similar distributions of transit time and system age.

Discussion
Simple, explicit, and general formulas for the densities, higher order moments, and first moments for the transit time, system age, and age vector of open linear autonomous compartmental systems were derived. They can be found in different places in this manuscript. For convenience, Table 1 provides a quick overview. Furthermore, a . , x * d is the diagonal matrix comprising the components of the steady-state vector small Python package for linear autonomous compartmental models, called LAPM, is provided on https://github.com/MPIBGC-TEE/LAPM. It deals with symbolic and numerical computations of the formulas given in Table 1.

Relation to Earlier Results
Over many years compartmental models with fixed coefficients have been studied (Anderson 1983;Bolin and Rodhe 1973;Eriksson 1971) by mainly investigating the mean age of particles in the system and the mean transit time of particles. The topic became of increasing interest again with the need for modeling the terrestrial carbon cycle.
The general situation of linear nonautonomous compartmental models was investigated by Rasmussen et al. (2016). For the special case of autonomous systems they provide the formulas (19), (18), and (20) for the mean age vector, mean system age, and mean transit time, respectively. The probabilistic approach used here showed that these formulas are just expected values of corresponding random variables with according probability distributions f a (y) = −(X * ) −1 e y B u, A ∼ PH(η, B), and T ∼ PH(β, B). In particular, the terms mean system age and mean transit time refer to the expected absorption times of related continuous-time Markov chains. A general formula for the moments of the distribution of the transit time in a compartmental system without external inputs is given by Hearon (1972). He also states a formula for the backward transit-time density in this situation. Manzoni et al. (2009) (also Priestley 1982Thompson and Randerson 1999) calculated a so-called transfer function or impulse response function Ψ that describes the relative contribution of the input at time t − T to the output at time t. They furthermore give a survivor functionÃ that describes the probability of the particle to remain in the system for at least a time t. In the present setup this means which leads by to the conclusion that the impulse response function Ψ is nothing else than the probability density function f T of a phase-type distribution PH(β, B) given by Eq. (7). They also formally provided an age density distribution which is exactly the probability density function f A of the system age given by Eq. (17). This becomes obvious from Eq. (10), which perfectly relates the distributions of transit time and system age.
For special cases of simple compartmental systems, they further provide formulas for the densities and means of transit time and system age. These formulas were obtained by simulating an impulsive input and considering the impulse response in the Laplace space. It turns out that those formulas are special cases of the general shape of probability density functions of phase-type. Examples are presented in "Appendix B".

Singular Compartmental Matrices
In this paper, only nonsingular compartmental or transition rate matrices were considered, because they ensure that every particle that enters the system will eventually leave it. It can be shown (Neuts 1981) that the transit time of a continuous-time Markov chain with singular transition rate matrix does not lead to a certainly finite transit time. In this case, the system contains traps from which particles are not able to get out once they have entered them. So, some particles will remain in the system forever. Although this seems unrealistic, there are successful carbon-cycle models with this structure.
A famous example is the RothC model (Jenkinson and Rayner 1977). It contains an inert compartment that has neither inputs nor outputs, but it has a positive constant content. To apply the stochastic framework of the present paper to that model, it is necessary to cut this compartment out of the system and to look at the steady state of the remaining smaller system. The transit time of newly arriving particles will not be affected. However, the age distribution of the system will depend on time since the particles in the inert compartment become older and older. A detailed analysis of how to treat systems with traps can be found in Jacquez and Simon (1993) and references therein.

Stochastic Versus Deterministic Approach
As noted by Purdue (1979), in modeling ecological systems, there are powerful arguments for the use of stochastic processes. Even if nature is completely deterministic, ecological systems are too complex for a complete theoretical understanding or descriptive tools. This lack of complete knowledge can be accounted for by using probabilistic methods. One of the first authors to integrate stochastic behavior in biological models was Bartholomay (1958). Purdue (1979) gives an early review on stochastic compartmental models with ecological applications. He showed how compartmental models can be considered as a deterministic representation of an underlying stochastic process by establishing the link to continuous-time Markov chains. While Eisenfeld (1979) elaborates on this link and applies it to a liver model, the present paper connects the continuous-time Markov chain approach to the theory of phase-type distributions (Neuts 1981) and focuses on transit time and ages. The deterministic theory describes the ideal or average behavior of a system, whereas stochastic models can examine the deviations from such ideal or average behavior. The stochastic approach allows to think of key quantities as random variables. Knowledge of their probability distributions provides a whole new set of tools not only to study the mean, but also higher order moments of such random variables. A high variance of the transit time, for example, expresses that it is not the usual behavior of particles to travel through the system as long as indicated by the mean transit time. There is rather a large variety of different transit times depending on the particular path taken by the particle. Long tails of the distribution mean that there are in fact particles that need very long to travel through the system and hence there are very old masses in the system. With the probability density functions in hand, it is also possible to think of quantiles and confidence intervals of transit time and system age, something not possible with the deterministic approach.
Linear autonomous compartmental carbon-cycle models relate to continuous-time Markov chains with homogeneous transition probabilities. Turning to nonautonomous models (Rasmussen et al. 2016), nonhomogeneous Markov chains would have to be considered in the probabilistic approach. In this situation the Kolmogorov equations govern the evolution of the transition probabilities. There is plenty of theory available on continuous-time Markov chains (Norris 1997;Ross 2010), and much of it might be applied to carbon-cycle models.

Conclusions
Open linear autonomous compartmental systems were modeled by absorbing continuous-time Markov chains and it was shown that the behavior of those systems is governed by the phase-type distribution. It applies for the transit time as well as for the system age, only with different initial distributions. This knowledge provides simple general formulas for the densities, means, and higher order moments of transit time, compartment ages, and system age. Furthermore, it explains the underlying structure of the intrinsic connection between system age and transit time (Bolin and Rodhe 1973). This approach also revealed the potential and opportunities of expressing deterministic compartmental models as stochastic processes. For example, important system diagnostics in reservoir theory such as system ages and transit times have analogous counterparts in probabilistic terms. The theory of stochastic processes can further help to address important questions in reservoir theory with applications in other fields even beyond global biogeochemical cycles.
Since a random variable cannot be drawn uniformly from the positive half-line, τ is drawn uniformly from the interval [0, L] for some L > 0. To illustrate τ 's dependence on L, it is denoted by τ L . The aim is to find P(Y τ L = 1) and let L tend to infinity. The existence of a random variable Y τ according to the limiting distribution is then guaranteed by Skorokhods representation theorem (Billingsley 1968).
Note that the random variable Y τ is stochastic in two ways. The value of Y at a deterministic time is stochastic and furthermore the time τ itself is stochastic. To overcome this difficulty, the probability P(Y τ L = 1) is conditioned on τ L = t. By the law of total probability it can be rewritten as Let p(t) := P(Y t = 1). If the existence of lim t→∞ p(t) were known, concluding would be allowed.
The problem of calculating lim t→∞ p(t) is left. To solve it, a major result from renewal theory is needed. Recall that f T denotes the common probability density function of the cycle durations.
Theorem 2 (Key renewal theorem) Suppose that the function z in the so-called general renewal equation is directly Riemann integrable. Then The proof and the rather technical definition of direct Riemann integrability can be found in Asmussen (2003, Sect. V.4).
To apply the key renewal theorem to p, it is necessary to show that p satisfies the general renewal equation (25) for some directly Riemann integrable function z. To this end, a standard technique from renewal theory is applied by conditioning on the event that the first restart occurs at time T 1 = x. Recall that T 1 has the probability density function f T . Then Now, this integral is split into a period before and a period after t to get The second integral describes the probability of Y being on during the first cycle, since t ≤ T 1 , and is equal to z(t) := P(Y 1 t = 1). The conditional probability P(Y t = 1 | T 1 = x) under the first integral equals P(Y t−x = 1) = p(t − x), since the process restarts at the end T 1 = x of the first cycle. Hence, satisfies the general renewal equation (25). Furthermore, the function z can be shown to be directly Riemann integrable. Therefore, the key renewal theorem states that The integral can be rewritten to and hence describes the expected time of Y being on during the first cycle. The equality of the limits in Eq. (24) finishes the proof.