Transit times and mean ages for nonautonomous and autonomous compartmental systems.

We develop a theory for transit times and mean ages for nonautonomous compartmental systems. Using the McKendrick-von Förster equation, we show that the mean ages of mass in a compartmental system satisfy a linear nonautonomous ordinary differential equation that is exponentially stable. We then define a nonautonomous version of transit time as the mean age of mass leaving the compartmental system at a particular time and show that our nonautonomous theory generalises the autonomous case. We apply these results to study a nine-dimensional nonautonomous compartmental system modeling the terrestrial carbon cycle, which is a modification of the Carnegie-Ames-Stanford approach model, and we demonstrate that the nonautonomous versions of transit time and mean age differ significantly from the autonomous quantities when calculated for that model.

Introduction

Compartment models play an important role in the modeling of many biological systems ranging from pharmacokinetics to ecology (Anderson 1983; Godfrey 1983; Jacquez and Simon 1993). Key values in understanding the dynamics of these systems are the transit time: the mean time a particle spends in the compartmental system measured as the mean time from entry into the system to leaving the system (Bolin and Rodhe 1973; Eriksson 1971), and the mean age: the mean age of particles still in the system (Bolin and Rodhe 1973; Eriksson 1971). It is well known that these quantities need not be the same (Bolin and Rodhe 1973; Eriksson 1971; Rothman 2015).

We are motivated by an interest in studying the dynamics of the terrestrial carbon cycle which is typically modeled as a number of discrete pools of carbon in plant biomass, litter and soil organic matter. Many of the best studied models of the dynamics of carbon are linear, which reflects the fact that changes in carbon pools are proportional to the pool size (Bolker et al. 1998). Additionally, most analyses make the further assumption that all parameters describing the dynamics (and the input fluxes) are constant in time, leading to a model in the form of an autonomous linear differential equation. In this autonomous case, it is possible to derive analytic formulae giving expressions for the transit time (Garcia-Meseguer et al. 2003; Manzoni et al. 2009). These formulae for transit time are given in terms of (constant) transfer coefficients among compartments and analogous formulae are available for the mean age of particles in the system.

Many applications of models of terrestrial carbon relate to situations in which constant model parameters are replaced by time-dependent functions. Perhaps the most well-known examples are studies of how terrestrial carbon dynamics respond to climate change. In these, it is often assumed that the specific rates (per unit carbon) of carbon inputs and losses from the system change over time as a function of changes in climate, such as temperature. For example, increases in temperature are normally assumed to increase the rates of soil decomposition (Lloyd and Taylor 1994; Orchard and Cook 1983; Rothman 2015). As a consequence, the compartmental models of interest are nonautonomous, i.e. they depend on time (Luo et al. 2001, 2015; Xia et al. 2012). Nonautonomous compartmental systems are special cases of linear nonautonomous differential equations (Kloeden and Rasmussen 2011), which, in contrast to the linear autonomous case, cannot be solved analytically in general. Yet, both the mean age of particles in the system and the transit time remain of great interest for these time-dependent systems, as both quantities can be potentially measured in the actual systems being modeled (Rothman 2015; Trumbore 2000).

In this paper, we develop a theory for transit times and mean ages of mass in nonautonomous compartmental systems. As noted in one of the first papers to study transit time (Bolin and Rodhe 1973), there is obviously a close connection between age distribution and transit time in compartment models. We will build on this relationship to develop an approach for understanding the definition of transit time. We define a time-dependent version of transit time as the mean age of mass leaving the compartmental system. We use a time-dependent version of the McKendrick–von Förster equation (Brauer and Castillo-Chavez 2012; McKendrick 1926; Thieme 2003), the classic first-order partial differential equation describing age distributions, to prove that the mean age of mass satisfies an (inhomogeneous) linear nonautonomous differential equation. We show that under weak conditions, this equation is exponentially stable. Starting with demographic models highlights another important aspect of our approach. As is well known, solutions of demographic models depend on initial conditions, so quantities like the mean age and transit time also depend on initial conditions, but conventional definitions of these quantities ignore the influence of the initial conditions. For this reason, our nonautonomous approach also provides additional insight for autonomous compartmental systems that are not in equilibrium.

We apply the theory we have developed to numerically study transit times for a nine-dimensional compartmental system model of the carbon cycle, which is a modified version of the Carnegie–Ames–Stanford approach (CASA) model (Buermann et al. 2007; Potter et al. 1993; Randerson et al. 1996). We compare our nonautonomous quantities to the classical notion of transit time for autonomous systems, where we freeze the nonautonomous system in time to obtain an autonomous system, and we assume that we are in equilibrium. Our simulations illustrate the different and sometimes diverging trajectories of the autonomous and nonautonomous quantities over time. Our results demonstrate the necessity of our theory for the computation of transit times in nonautonomous compartmental systems and in autonomous compartmental systems that are not in equilibrium.

This paper is organized as follows. In Sect. 2, we first review the theory of transit times for autonomous compartmental systems, and we provide a heuristic derivation of the transit time formula. We then define nonautonomous compartmental systems in Sect. 3. In Sect. 4, we prove that under the assumption that the compartmental system is lower block triangular, and the diagonal blocks a diagonally dominant, the nonautonomous compartmental system is exponentially stable. In Sect. 5, we prove that the mean ages satisfy a linear nonautonomous differential equation, and we then use the stability criterion from Sect. 4 to prove exponential stability of the mean age equation. We define the concept of a transit time for nonautonomous compartmental systems in Sect. 6. In Sect. 7, we show that our nonautonomous theory is consistent with the autonomous case, in the sense that we get exactly the well-known autonomous transit time formula when applying the nonautonomous transit time to an autonomous system. Finally, in Sect. 8, we apply the theory to compute transit times for a nonautonomous compartmental model of the carbon cycle, which is a simplified version of the Carnegie–Ames–Stanford approach (CASA) model.

Transit times and mean ages for autonomous compartmental systems

An open (linear) autonomous compartmental system with both inputs and outputs (Anderson 1983; Godfrey 1983; Jacquez and Simon 1993) and with d pools is described by an inhomogeneous linear differential equationwhere is an invertible matrix, , and the entries of the matrix B satisfyThe i-th row of the matrix B describes the dynamics of the mass in pool i: is the rate at which mass moves from pool j to pool i, and is the rate at which mass leaves the pool i which includes transfer to other pools and losses from the system. The flux at which mass enters from outside the system to pool i is given by .

for all ,

for all ,

for all .

We assume that the homogeneous linear system is exponentially stable, i.e. all eigenvalues of B have negative real parts (this is fulfilled e.g. when the matrix B is strictly diagonally dominant). This means that (1) has the exponentially stable equilibrium .

The concept of transit time for compartmental systems describes the mean time a particle spends in the compartmental system before it is released. There is a huge amount of literature on this topic, see e.g. Anderson (1983), Bolin and Rodhe (1973), Eriksson (1971), Garcia-Meseguer et al. (2003), Manzoni et al. (2009), but to our knowledge, the following simple derivation of the transit time formula has not been written down before.

Define as the mean (remaining) transit time in the system for a particle that has entered pool i either from outside the system or from another pool, and note that the transit time in pool i for a particle that has entered pool i either from outside the system or from another pool is given by . Let be the probability that a particle that enters pool i goes next to pool j, and note thatNext, note that the transit times for particles entering any pool i must satisfy the equationreflecting the fact that a particle in pool i spends the average time in pool i, before it either leaves the system or moves with the probability to pool j, after which it spends the mean time before it leaves the system. This reads asand multiplying the i-th row of this equation with yieldsHence .

Let be the fraction of particles that enter the system from outside directly into pool i, i.e.and let . Then the transit time for the whole system is given byNote that this transit time is equal to the turnover time (see Bolin and Rodhe 1973), which follows directly from .

We will show later that if the linear compartmental system (1) is in the equilibrium , then the mean age of the particles in the system is given bywhere , defined bydescribes how mass is distributed when the system is in equilibrium. It is well-known that the mean age M is unequal to the transit time R (Bolin and Rodhe 1973; Rothman 2015), and we will demonstrate this now by means of two very simple compartmental systems.

Example 1

(Transit times and mean ages) Consider the two compartmental systemsandIt is easy to see that the transit times and for the two pools satisfies for (4) and for (5). This follows either from using the above explicit formula for the vector r, or by considering the fact that particles can only leave from pool 1 in (4) and from pool 2 in (5). Since the transit time is given in both cases by , and the mean age is a convex combination of and , the transit time will be smaller than the mean age in (4), in contrast to the situation in (5).

Nonautonomous compartmental systems

In contrast to the autonomous case, both the coefficient matrix B and the input vector s of a nonautonomous compartmental system are allowed to depend on time.

Definition 1

(Nonautonomous compartmental system) Let with be a time interval, be a bounded continuous function of invertible matrices and be a bounded continuous function. A (linear) nonautonomous compartmental system with d pools is given by an inhomogeneous linear nonautonomous differential equationwhere we assume that the entries of the matrix B(t) satisfy

for all and ,

for all and ,

for all and .

Let denote the transition operator of the corresponding homogeneous equation , i.e. the function is the solution to fulfilling the initial condition . Then the maximal solution to (6) satisfying the initial condition is given byIn contrast to the autonomous case, nonautonomous compartmental systems of dimension two or higher are not explicitly solvable in general. Solutions can be obtained for systems with no feedbacks between pools, as the following example demonstrates.

Example 2

(Explicitly solvable nonautonomous two-pool model) The nonautonomous compartmental systemwhere and for all , can be solved explicitly as follows: the general solution of the first equation is given byand thus, the second equation reads asand can be solved using (7), since the equation is one-dimensional.

Exponential stability of nonautonomous compartmental systems

In this section, we provide a sufficient condition for global exponential stability of the nonautonomous compartmental system (6). This criterion will concern only the homogeneous part of (6), i.e. the matrix-valued function B, from which stability for the inhomogeneous equation follows. Since the result holds also for linear systems which are not compartmental systems, we formulate it more generally.

Theorem 1

(Sufficient condition for exponential stability) Consider the linear nonautonomous differential equationwith transition operator . Suppose that the function B is of the formfor with bounded functions . Note that . We assume that the linear subsystems , are strictly diagonally dominant, i.e. there exists a such thatThen the linear system (8) is exponentially stable, i.e. there exist constants and such that

for all and ,

for all and ,

for all and .

Proof

Assume first that I is bounded below, and consider the linear systemsfor each . These systems are strictly diagonally dominant, and it follows from (Coppel 1978, Proposition 3, page 55) that there exist and such thatwhere is the transition operator of (11). Next Pötzsche (2016, Theorem 4.1) or Battelli and Palmer (2015, p. 540) yields that the dichotomy spectrum of (8) is bounded above by . This in turn implies the claimed estimate (10). In case , the results from Pötzsche (2016, Theorem 4.1) and Battelli and Palmer (2015, p. 540) are not applicable directly, since they require the system to be defined on a half line, but the result follows by considering the two time intervals and separately (note that we are not interested in the dichotomy spectrum for the entire line, which is not determined by the block diagonal system; we only require an upper bound, which we get from the block diagonal system).

The estimate (10) for the homogeneous system (8) implies that any two solutions of the compartmental system (6) converge to each other exponentially. More precisely, given two solutions of (6), thenwhich follows from the fact that the difference of these two solutions is a solution of the homogeneous system (8), for which the estimate (10) holds. This implies that any solution is forward attracting, and in case the interval I is unbounded below, then there also exists a unique pullback attracting solutionsee Aulbach and Wanner (1996). This solution pullback attracts bounded sets , in the sense ofwhere denotes the maximal solution defined in (7) and denotes the Hausdorff distance. We refer to Kloeden and Rasmussen (2011), Rasmussen (2007) for an introduction to forward and pullback attractors of nonautonomous dynamical systems.

The mean age system

We prove in this section that the mean ages of mass in a nonautonomous compartmental system are solutions of a linear nonautonomous differential equation, which we call the mean age system. We derive this result from the evolution of age distributions, given by the well-known McKendrick–von Förster equation (McKendrick 1926; Brauer and Castillo-Chavez 2012; Thieme 2003), which is a linear first order partial differential equation. We also prove that the mean age system is exponentially stable under additional weak assumptions, by applying the theory developed in Sect. 4.

The mean age system is pivotal for the analysis of transit times for nonautonomous compartmental systems, since in order to compute the average time the mass spends in the system, we do not need to look at the full age distribution of ages, but only at the mean ages.

Let be the density function on age a for the mass in pool i at time t, where the age is the time since the mass entered the system. Note that the following formulation is valid in principle even if all rates are age-dependent, i.e.  also depends on a, but we will not treat this situation here. The McKendrick–von Förster equation is given bywith boundary conditionNote that one also needs to specify initial conditions .

A componentwise solution can be written as follows. Note that if there are no loops, the solution is explicit, otherwise it is only implicit. The formula for isNote that an analogous formula for exists.

We are particularly interested in the transit time of (6) at a particular time t, which corresponds to the mean age of mass leaving the system at time t. For this purpose, we do not need the full age distribution determined by (13), since the situation is fully described by the mean age of mass in pool i, denoted as . The following theorem says that the evolution of the mean ages is determined by an ordinary differential equation.

Theorem 2

(Mean age system) Consider the nonautonomous compartmental system (6) with a fixed solution of positive entries. Let be the density function on age a for the mass in pool i at time t (note that ), and define the mean age of mass in pool i byThen the mean ages solve the ordinary differential equationwith

By using (integration by parts), it follows thatThis finishes the proof.

Combining the Eqs. (6) and (15) yieldswhich is a 2d-dimensional ordinary differential equation of skew product type, i.e. the x-equation does not depend on , but the equation for depends on x. Note that (16) is a nonlinear equation, but given a solution of (6), the age equation (15) is an inhomogeneous linear nonautonomous differential equation, which reads aswherewith for all .

We will show now that under additional weak assumptions, the mean age equation is exponentially stable.

Theorem 3

(Exponential stability of the mean age system) Consider the nonautonomous compartmental system (6) with a fixed solution of positive entries that are bounded and bounded away from zero, and suppose that (6) satisfies the assumptions of Theorem 1 with . In addition, assume thatThen the mean age system (15) is exponentially stable. More precisely, there exist and such that the transition operator of the homogenous equation satisfies the estimate

for all and , and

for all and , there exists a such that for all .

We show now that the three conditions (i)–(iii) of Theorem 1 are satisfied with replaced by . Note first that the matrix A(t, x(t)) has the same block decomposition as the matrix B(t), which is described in (9).

Condition (i) of Theorem 1 follows from (a) (in case of ) or (b) (in case ; note that the sum of the entries in the i-th row of the matrix A(t, x(t)) equals to , and (b) guarantees that the diagonal entry is negative even though might be zero). Condition (ii) of Theorem 1 follows from the fact that the original system (6) is a compartmental system, and the solution x(t) of (6) has positive entries. Finally, condition (iii) of Theorem 1 follows from fact that the sum of the i-th row of the matrix A(t, x(t)) equals to , and the positive contribution of at least , with i and j chosen as in (b), will not be considered in the sum in condition (iii) of Theorem 1 and for this reason contributes negatively to this sum.

A natural choice for the solution in the above theorem is the exponentially stable solution defined in (12) if the interval I is unbounded below. If the interval I is unbounded below, this will be the only bounded solution of the system, i.e. the norm of all other solutions converges to in the limit , so the solution (12) is the only solution to which the theorem can be applied. However, if the interval I is bounded below, then all solutions of the nonautonomous compartmental system (6) are bounded and exponentially stable, and they are also bounded away from zero due to assumption (a) of Theorem 3.

Nonautonomous transit times

We define transit time as the mean age of mass leaving the system at a particular time t. Note that in our nonautonomous context, this quantity depends on the actual time t. We also provide a formula that corresponds to the mean age of mass currently residing in the compartmental system.

Definition 2

(Nonautonomous transit time and mean age) Consider the skew product system (16) consisting of the nonautonomous compartmental system (6) and the mean age system (15). The transit time of a solution , of this system is then defined asand then mean age of this solution is defined by

The transit time is the mean age of carbon leaving the system at time t, where as the mean age is the mean age of carbon in the system at time t.

Note that, in general, and are different, see Example 1 for the autonomous case. In the following example, we show that transit times and mean ages are the same for one-dimensional compartmental systems.

Example 3

(Transit time and mean ages for one-dimensional compartmental systems) Let be an interval, and consider the one-dimensional nonautonomous compartmental systemwhere and are bounded continuous functions. Fix a positive solution of this system. Note that the solution is given explicitly bywhere and are initial time and condition. Then the mean age equation is given byand also this equation can be solved explicitly using (17). Note that, in this one-dimensional context, the formulae for transit time and mean age from Definition 2 are given by exactly the solution to this equation:

Consistency with the autonomous case

In this section, we derive simple expressions for the transit time and mean age from Definition 2 in the special case of an autonomous compartmental system. The expression for the autonomous transit time coincides with the heuristically obtained formula (2), and we confirm the expression for the mean ages stated in (3).

Consider an autonomous compartmental systemwith an invertible matrix and . We assume that the homogeneous system satisfies the assumptions of Theorems 1 and 3. Note that (18) has the exponentially stable equilibrium .

Lemma 1

Consider the autonomous differential equation (18). Then the mean age equation (15) for the equilibrium reads asand has the exponentially stable equilibriumwhere .

Note that the mean age equation (15) is given bySince , we get for all , so the mean age equation (15) gets simplified toHence the attractive equilibrium of (15) is given bywhich finishes the proof of this lemma.

Let be the fraction of particles that enter the system from outside directly into pool i, i.e.and let . Moreover, define bywhich describes how mass is distributed when the system is in equilibrium. Note that .

Proposition 1

(Autonomous transit times and mean ages) Consider the autonomous compartmental system (18). The transit time with respect to the equilibrium solution does not depend on time and is given byand the mean age of mass is given by

Using Definition 2, we havefor the transit time andfor the mean age. Note that both quantities do not depend on t, and this finished the proof of this proposition.

Note that derivation of the autonomous quantities for transit time R and mean age M in Proposition 1 required the autonomous compartmental system (18) to be in equilibrium, and the classical approach to transit times, as outlined in Sect. 2, is not applicable for autonomous systems not in equilibrium. It is very important to note that Definition 2 is useful for autonomous systems also, since it is applicable to systems that are not in equilibrium. For such autonomous systems, transit times and mean ages will depend on time in general, and although they converge to R and M in the limit , they might be very different to R and M.

Mean ages and transit times for the CASA model

Here we illustrate predicted changes in the mean age of carbon leaving and remaining in the system for a terrestrial carbon model under a climate change scenario. We consider a modification of the CASA model as used in Buermann et al. (2007) globally without resolving the spatial details of carbon pools using nine pools representing the global terrestrial carbon (e.g. three pools for plant biomass, or litter or soil organic matter). This caused the model to be precisely of the form of (6). Climate change was simulated by increasing atmospheric over time, which affected both B(t) and s(t) in (6). Increased directly increases carbon inputs s(t) through carbon dioxide fertilization. They also directly increase mean global temperatures. This increases the carbon loss rates from some of the carbon pools, changing components of B(t), and also has an effect on s(t). Thus increased alters the input and loss rates of components of the terrestrial carbon cycle, making both the sign and magnitude of the net change in carbon storage dependent the sensitivity of carbon inputs and loss rates.

We simulated changes in atmospheric usingwhere is the atmospheric carbon dioxide concentration in parts per million and t is years since the year 1850. This represents a plausible time course of atmospheric from year under a zero-mitigation, business as usual global change scenario (Raupach et al. 2011) (illustrated in Fig. 1a).

Fig. 1

Forcing functions and solution of the simplified CASA model. a Nonautonomous dynamics are driven by changes in atmospheric over time as given by . b The increased alters total carbon inputs per unit time via . c Increasing also increases temperatures which increases litter and soil carbon decomposition rates via . d the resulting solution of total terrestrial carbon over time. Parameters for this model are as given in the text but also with

The effect of on mean global temperatures is modelled aswhere is the mean land surface temperature in 1850, and is the sensitivity of global temperatures to . We chose an upper extreme of based on the literature because the resulting simulation emphasises well the interplay between increased carbon input rates and carbon loss rates Scheffer et al. (2006). Changes in carbon input rates are simulated usingwith , where is the proportion of carbon input going to the different carbon pools, is the proportion of gross primary production that remains after respiration and is the sensitivity of s(t) to and , given bywhere is the ratio of the intracellular to , and is given bysee Polglase and Wang (1992). The solution of (21) with changes in and as described above is illustrated in Fig. 1b. The matrix controlling the rates of carbon transfer and loss from the system is given byindicating that it is the loss rates of pools that change with time. The coefficients are listed in the legend to Fig. 1, andwhere is the scaling of decomposition rates at degrees Celsius. Equation (22) is illustrated in Fig. 1c.

To define the model initial conditions we assume that for all and that x(0) has reached the positive equilibrium solution of the resulting system of autonomous equations. The model is then simulated forward from this initial condition using (19) as the forcing function. Under this simulated scenario, total land carbon increases then decreases over time as shown in Fig. 1d. This would represent an initial net uptake of carbon from the atmosphere due to carbon dioxide fertilization followed ultimately by a net carbon loss from the land back to the atmosphere due to global warming (Fig. 2 shows how this carbon change over time is distributed amongst the different components of x(t)).

Fig. 2

Breakdown of the contributions of the different vegetation carbon pools to the change in the overall terrestrial carbon storage dynamics illustrated in Fig. 1d. Pools and are carbon in leaves, roots and wood, respectively; pools to are carbon in different forms of litter, and pools to are carbon in different forms of soil

Calculations of the transit time and mean age of carbon in the system (according to Definition 2), for the nine-pool model for the climate change simulation described above, show an order of magnitude difference in the absolute values of and (Fig. 3). This indicates that the average age of carbon stored on land is much older than the average age of carbon leaving the land. Note that at (corresponding to the year 1850), we assumed that , where M is the mean age of the equilibrium solution at according to Proposition 1.

Fig. 3

Mean transit time , and mean age , compared with the instantaneous quantities R and M

Perhaps surprisingly, the monotonic forcing of B(t) and s(t) translates into non-monotonic effects on and . A detailed mathematical investigation of this phenomenon is outside the scope of the present study.

The nonautonomous properties and show contrasting trajectories to the instantaneous properties R and M (which we computed according to Proposition 1, but note that, since the system is nonautonomous, the assumptions of this proposition are not fulfilled). For example the latter properties change monotonically over time. This must be because the long term outcome of an increase in the input rate of young carbon and an increase in the output rate of old carbon is a decrease in the age of carbon both leaving and remaining in the system. Over the course of the simulation the numerical values of the autonomous and nonautonomous properties become visibly different (Fig. 3). This is because it will take a long time for the values of and to approach R and M due to the small loss rate of the ninth soil pool.

Conclusions

Models for terrestrial carbon cycling have led to renewed interest in the properties of compartment models. Key quantities that have been studied over many years in compartment models with parameters fixed in time (Eriksson 1971; Bolin and Rodhe 1973; Anderson 1983) are the mean age of particles in the system and the transit time of particles leaving the system. Formulae for these quantities that give the mean age and transit time in terms of parameters of the system in the long time limit have led to insights, but cannot be applied to the case of changing parameters.

As parameters change, for example in a model of carbon cycling due to climate change, it is not correct to calculate the mean age or transit time from the instantaneous parameter values. Using the theory of nonautonomous differential equations as a tool, and beginning with time dependent age structured models, we are able to define and derive formulae for the transit time and mean age for particles in the case of temporally changing parameters. These definitions lead to quantities that reduce to the analogous formulae for the autonomous (constant parameter) case when parameters do not change in time. However, the formulae for the nonautonomous case also highlight the fact that even in the constant parameter case the transit time and mean age do depend on initial conditions; some of the standard formulae do not include this dependence.

The difference between a transit time or mean age that is computed based on the parameters at a given instant and the better approach of taking into account the history of the system can be substantial as we illustrate using a variant of the CASA model. Thus, the approach we develop here is not just of mathematical interest but is of substantial practical importance as well.