Population growth in discrete time: a renewal equation oriented survey

Abstract

Traditionally, population models distinguish individuals on the basis of their current state. Given a distribution, a discrete time model then specifies (precisely in deterministic models, probabilistically in stochastic models) the population distribution at the next time point. The renewal equation alternative concentrates on newborn individuals and the model specifies the production of offspring as a function of age. This has two advantages: (i) as a rule, there are far fewer birth states than individual states in general, so the dimension is often low; (ii) it relates seamlessly to the next-generation matrix and the basic reproduction number. Here we start from the renewal equation for the births and use results of Feller and Thieme to characterize the asymptotic large time behaviour. Next we explicitly elaborate the relationship between the two bookkeeping schemes. This allows us to transfer the characterization of the large time behaviour to traditional structured-population models.

Keywords:

• Structured-population model
• renewal equation
• growth rate
• basic reproduction number
• next-generation matrix
• reproductive value
• stable distribution

Mathematics Subject Classification:

• 92D25

1. Introduction

In the pioneering paper [7], Jim Cushing and Zhou Yicang

1. defined the net reproductive number (aka the basic reproduction number and in the present paper denoted by ${\mathcal{R}}_0$) in the context of linear discrete time population models, while highlighting its interpretation as the expected lifetime number of offspring,

2. showed that the stability of the extinction state (i.e. the trivial steady state) is governed by the sign of ${\mathcal{R}}_0-1$ (where, in the case of nonlinear models, ${\mathcal{R}}_0$ refers to the linearized model),

3. demonstrated, by way of examples, that often ${\mathcal{R}}_0$ is much easier to determine than the real time growth rate (in the present paper denoted by ρ).

By making systematic use of Perron-Frobenius theory, Chi-Kwong Li and Hans Schneider streamlined the proofs and extended these results in [22]. For far-reaching generalizations we refer to recent papers of Horst Thieme, in particular [28, 29]. And for similar continuous time results in the context of epidemic models to [9, 14, 15, 27]. Also see [6].

Over the years (and in close collaboration with various other authors) we have developed theory for continuous time structured-population models [10–12, 24]. These papers deal with rather general models and may, as a consequence, not be easily accessible. The aim of the present paper is to formulate the discrete time version of the linear theory in a reader-friendly fashion. To this end, we limit our attention to population models in which individuals are distinguished by a number of traits (e.g. species, age, size, spatial location, infection status), but such that the set of conceivable individual states (shortly, i-states) is finite.

In our approach, the renewal equation (RE) for the birth rate takes centre stage. Once one solves the RE constructively, all other relevant quantities are given by explicit expressions in terms of the initial condition and the birth rate. As a consequence, we can deduce the asymptotic large time behaviour of every quantity of interest from the asymptotic large time behaviour of the birth rate. And to determine the latter, we can rely on Feller's celebrated Renewal Theorem for scalar RE and on Thieme's generalization for systems of RE (see Section XIII.10 in [16] and [25]; continuous time results can be found in [4, 8, 17]).

The key feature that makes this approach both attractive and efficient is that, as a rule, there are far fewer individual birth states than individual states in general. In particular, the RE is a scalar equation when all individuals are identical at birth, either in the literal sense (like, for instance, in the case of age as the only structuring trait) or in the stochastic sense, when states-at-birth have a fixed probability distribution (for example, in an epidemiological model where individuals begin their infected life (i.e. are born epidemiologically) as symptomatic with probability p or asymptomatic with probability 1−p, regardless of the i-state of the individual that infected them).

We begin by preparing the ground: in Section 2, we formulate the RE for the birth rate and determine the asymptotic large time behaviour using the theorems of William Feller (for scalar RE) and Horst Thieme (for systems of RE). In Section 3 we first show how classical discrete time structured-population models (SPM) give rise to (typically much smaller) systems of RE for the birth rate and discuss how the two bookkeeping schemes relate to each other in terms of dynamics and next-generation matrices. Next we show how, after having deduced the asymptotic behaviour of solutions of the RE using the results of Section 2, one can determine the asymptotic i-state distribution and the reproductive values in SPM from the ones obtained in the corresponding RE (see Figure 1).

Figure 1. Deducing the asymptotic large time behaviour in (linear discrete time) structured-population models via the corresponding renewal equation for the birth rate.

In Section 6 of [10] and Section 3.2 of [12] it was shown that, for the special situation that ${\mathcal{R}}_0=1$, not only the spectral equivalence $\rho =1$ holds, but that one can also compute the asymptotic distribution (i.e. the right eigenvector) in the real time setting from that in the generation bookkeeping framework. In Section 4 we supplement this result with a corresponding result for the reproductive values (or, in other words, for the left eigenvector).

In [12] it was demonstrated that for a large class of nonlinear models the condition ${\mathcal{R}}_0=1$ arises when looking for steady states. This makes the results about connecting the eigenvectors useful for nonlinear models as well.

In Appendix we present a brief outline of continuous time analogues of various equations and identities derived in Sections 2, 3 and 4.

2. The renewal equation

The first step in the formulation of renewal equations is the identification of states-at-birth. We may interpret this literally in the sense of determining all the i-states in which individuals may begin their lives. For example, think of a model where individuals are structured by age a and size s and where the set of i-states is $\{(a,s):a=\{1,\dots ,k\},s=\{\mathrm{small},\mathrm{large}\}\}.$ Assuming that newborns may be either small or large, there are two states-at-birth: $(1,\mathrm{small})$ and $(1,\mathrm{large}).$ However, when small and large offspring is produced at a fixed ratio $p:(1-p)$ (independently of the i-state of the parent) we can exploit the fact that there is a single state-at-birth in the stochastic sense: all individuals are born with the same probability distribution of their i-state. Another example is an epidemiological model where individuals may begin their infected life (i.e. are born epidemiologically) as either asymptomatic or symptomatic. There are again two states-at-birth in the literal sense, but when new cases arise at a fixed ratio there is but one state-at-birth in the stochastic sense.

Both interpretations of the term state-at-birth are valid. However, while the literal interpretation may be a natural starting point, understanding the term in the stochastic sense allows us to formulate RE systems of minimal size (we elaborate this point in Section 3).

We begin by considering the case where all individuals are identical at birth, i.e. are born with the same probability distribution of their i-state.

2.1. The scalar equation

Imagine that immediately after the yearly breeding season a census is made of the subpopulation of newborn females. Let $b(t)$ denote that number in year t. If we ignore density dependence and adopt a deterministic point of view at the population level (i.e. we assume that the numbers are so large that it makes sense to focus on expected values and to ignore the effects of demographic stochasticity) then (1) $b(t)=\sum _{s=1}^{\mathrm{\infty }}k(s)b(t-s),$(1) where $k(s)$ is the expected number of daughters born to a female in year s after her own birth. We consider the sequence $\{k(t){\}}_{t=1}^{\mathrm{\infty }}$ as the basic model ingredient and make the following assumptions:

(A1)	$k(t)\ge 0$ for all $t\in \mathbb{N}$ with strict inequality for at least one value of t.
(A2)	$\sum _{t=1}^{\mathrm{\infty }}k(t)<\mathrm{\infty }$.

Note that (2) ${\mathcal{R}}_0=\sum _{t=1}^{\mathrm{\infty }}k(t)$(2) is the basic reproduction number, i.e. the expected lifetime number of (female) offspring.

Equation (1(1) $b(t)=\sum _{s=1}^{\mathrm{\infty }}k(s)b(t-s),$(1) ) is linear and translation invariant, so we expect it to have geometric solutions. If we substitute the Ansatz (3) $b(t)=z^t$(3) into (1(1) $b(t)=\sum _{s=1}^{\mathrm{\infty }}k(s)b(t-s),$(1) ) and divide both sides by $z^t$ we obtain the Euler–Lotka characteristic equation (4) $1=\sum _{s=1}^{\mathrm{\infty }}k(s)z^{-s}=:\overline{k}(z)$(4) for the unknown z. Note that $\overline{k}(z)$ is the (one-sided) z-transform of $\{k(t){\}}_{t=0}^{\mathrm{\infty }}$ (with $k(0):=0$) and that the z-transform is the discrete time analogue of the Laplace transform.

For real z, $\overline{k}$ is a strictly decreasing function with limit zero at infinity. So if $\overline{k}$ takes, for a real z, a value larger than one, Equation (4(4) $1=\sum _{s=1}^{\mathrm{\infty }}k(s)z^{-s}=:\overline{k}(z)$(4) ) has precisely one real solution. We denote this solution by ρ and observe that there are precisely three possibilities: $\begin{array}{rl}& 1<\rho <{\mathcal{R}}_0,\\ & 1=\rho ={\mathcal{R}}_0,\\ & {\mathcal{R}}_0<\rho <1\end{array}$whenever $k(s)>0$ for a value $s\ge 2$ (in the degenerate case that only $k(1)>0$ we have $\rho ={\mathcal{R}}_0=k(1)$). To see this, note that

1. $\overline{k}(1)={\mathcal{R}}_0$,

2. if ${\mathcal{R}}_0>1$ then $\overline{k}({\mathcal{R}}_0)<1$,

3. if ${\mathcal{R}}_0<1$ then $\overline{k}({\mathcal{R}}_0)>1$ (possibly infinite).

If z is complex then $|\overline{k}(z)|\le \overline{k}(|z|).$Consequently the solutions of (4(4) $1=\sum _{s=1}^{\mathrm{\infty }}k(s)z^{-s}=:\overline{k}(z)$(4) ) are contained within the closed disc of radius ρ.

We say that $\{k(t){\}}_{t=1}^{\mathrm{\infty }}$ is periodic with period p>1 if $k(s)=0$ except when s is an integer multiple of p. One easily verifies that in that case z is a solution of (4(4) $1=\sum _{s=1}^{\mathrm{\infty }}k(s)z^{-s}=:\overline{k}(z)$(4) ) if z is the product of ρ and a pth root of unity. So for periodic $\{k(t){\}}_{t=1}^{\mathrm{\infty }}$ there are non-real solutions of (4(4) $1=\sum _{s=1}^{\mathrm{\infty }}k(s)z^{-s}=:\overline{k}(z)$(4) ) on the circle of radius ρ in $\mathbb{C}$. We make a further assumption:

(A3)	$\{k(t){\}}_{t=1}^{\mathrm{\infty }}$ is not periodic.

So far we focused on special solutions of (1(1) $b(t)=\sum _{s=1}^{\mathrm{\infty }}k(s)b(t-s),$(1) ) of the form (3(3) $b(t)=z^t$(3) ). Now, in the spirit of [13], we prescribe the history of b up to t = 0 by putting (5) $b(\tau )=\theta (\tau ),\tau =\dots ,-2,-1$(5) for a given function $\theta (\tau )\ge 0$ such that (6) $g(t):=\sum _{s=t+1}^{\mathrm{\infty }}k(s)\theta (t-s)=\sum _{l=1}^{\mathrm{\infty }}k(l+t)\theta (-l)$(6) is finite for $t\in \mathbb{N}$. This allows us to rewrite (1(1) $b(t)=\sum _{s=1}^{\mathrm{\infty }}k(s)b(t-s),$(1) ) in the form (7) $b(t)=\sum _{s=1}^{t}k(s)b(t-s)+g(t),t\in \mathbb{N},$(7) with the convolution product interpreted as zero for t = 0, i.e. (8) $b(0)=g(0).$(8) Let $\overline{b}$ and $\overline{g}$ denote the (one-sided) z-transforms of, respectively, $\{b(t){\}}_{t=0}^{\mathrm{\infty }}$ and $\{g(t){\}}_{t=0}^{\mathrm{\infty }}$, i.e. $\overline{b}(z)=\sum _{s=0}^{\mathrm{\infty }}b(s)z^{-s},\overline{g}(z)=\sum _{s=0}^{\mathrm{\infty }}g(s)z^{-s}.$The z-transform has the useful property that it turns a convolution product into an ordinary product. So in terms of the z-transforms, the convolution equation (7(7) $b(t)=\sum _{s=1}^{t}k(s)b(t-s)+g(t),t\in \mathbb{N},$(7) ) can be written as $\overline{b}=\overline{k}\overline{b}+\overline{g}.$If follows that $\overline{b}=(1-\overline{k}{)}^{-1}\overline{g}$and next, by taking the inverse z-transform, that $b(t)=\frac{1}{2\mathrm{\pi i}}\oint (1-\overline{k}(z){)}^{-1}\overline{g}(z)z^{t-1}\mathrm{d}z$where the closed contour encircles the origin counter-clockwise at a distance great enough to enclose all singularities of the integrand. The observations made above concerning the solutions of (4(4) $1=\sum _{s=1}^{\mathrm{\infty }}k(s)z^{-s}=:\overline{k}(z)$(4) ) now suggest that $b(t)\approx {\rho }^t$ for $t\to \mathrm{\infty }$. Theorem 2.1 in Section XIII.10 of [16] gives a far more informative characterization of the asymptotic large time behaviour of b.

Theorem 2.1

Feller's Renewal Theorem

Consider the renewal Equation (7(7) $b(t)=\sum _{s=1}^{t}k(s)b(t-s)+g(t),t\in \mathbb{N},$(7) ) and assume that (A1), (A2) and (A3) hold. Furthermore, assume that $g(t)\ge 0$ for all $t\in \mathbb{N}$ and that $\sum _{t=0}^{\mathrm{\infty }}g(t)<\mathrm{\infty }$. For ${\mathcal{R}}_0>1$ let $\rho >1$ denote the unique real solution of the Euler-Lotka Equation (4(4) $1=\sum _{s=1}^{\mathrm{\infty }}k(s)z^{-s}=:\overline{k}(z)$(4) ).

(i)	If ${\mathcal{R}}_0<1$ then $b(t)\to 0$ for $t\to \mathrm{\infty }$ and $\sum _{t=0}^{\mathrm{\infty }}b(t)=\frac{\overline{g}(1)}{1-{\mathcal{R}}_0}.$
(ii)	If ${\mathcal{R}}_0=1$ then for $t\to \mathrm{\infty }$ $b(t)\to \frac{\overline{g}(1)}{\sum _{s=1}^{\mathrm{\infty }}\mathrm{sk}(s)}.$
(iii)	If ${\mathcal{R}}_0>1$ then for $t\to \mathrm{\infty }$ ${\rho }^{-t}b(t)\to \frac{\overline{g}(\rho )}{\sum _{s=1}^{\mathrm{\infty }}\mathrm{sk}(s){\rho }^{-s}}.$

Remark 2.2

To facilitate comparison with Feller's formulation in [16], we present the notation-translation table:

Table

Here	Feller [16]
b	v
k	f
g	b

.

More importantly, note that Feller uses the generating function rather than the z-transform. If $\tilde{k}$ denotes the generating function of $\{k(t){\}}_{t=0}^{\mathrm{\infty }}$ then $\tilde{k}(s)=\overline{k}(s^{-1}).$

If (6(6) $g(t):=\sum _{s=t+1}^{\mathrm{\infty }}k(s)\theta (t-s)=\sum _{l=1}^{\mathrm{\infty }}k(l+t)\theta (-l)$(6) ) holds we deduce from the second identity that (9a) $\overline{g}(\rho )=\sum _{l=1}^{\mathrm{\infty }}c(l)\theta (-l)$(9a) with (9b) $c(l):=\sum _{s=0}^{\mathrm{\infty }}{\rho }^{-s}k(l+s).$(9b) The coefficients $c(l)$ thus tell us how the various components of the initial condition (5(5) $b(\tau )=\theta (\tau ),\tau =\dots ,-2,-1$(5) ) contribute to the ultimate geometric population growth with multiplication factor ρ when ${\mathcal{R}}_0>1$ (for the time being, we do not normalize these contributions).

We conclude that Theorem 2.1 gives a complete description of the large time behaviour of solutions of the renewal Equation (7(7) $b(t)=\sum _{s=1}^{t}k(s)b(t-s)+g(t),t\in \mathbb{N},$(7) ) when we prescribe as an initial condition the history of b up to time zero as in (5(5) $b(\tau )=\theta (\tau ),\tau =\dots ,-2,-1$(5) ).

2.2. Systems of renewal equations

Suppose that individuals begin their lives in one of m>1 states-at-birth. The renewal equation now reads (10) $B(t)=\sum _{s=1}^{\mathrm{\infty }}K(s)B(t-s),$(10) where the components $b_j(t)$ of the vector $B(t)\in {\mathbb{R}}^m$ denote the number of females born at time t with state-at-birth j and $K(s)$ is a positive $m\times m$ matrix with elements $\begin{array}{rl}{k}_{\mathrm{ij}}(s)& =\mathit{the}\mathit{expected}\mathit{number}\mathit{of}\mathit{daughters}\mathit{with}\mathit{state}\mathit{-}\mathit{at}\mathit{-}\mathit{birth}\mathit{i}\mathit{produced}\mathit{by}\mathit{a}\\ & \mathit{female}\mathit{that}\mathit{was}\mathit{herself}\mathit{born}\mathit{s}\mathit{time}\mathit{units}\mathit{ago}\mathit{with}\mathit{state}\mathit{-}\mathit{at}\mathit{-}\mathit{birth}\mathit{j}\mathit{.}\end{array}$

Remark 2.3

A matrix A is called positive if $a_{\mathrm{ij}}\ge 0$ for all i and j. In this case, we write $A\ge 0$. We write A>0 when $A\ge 0$ and $A\ne 0$. A matrix A is called strictly positive if $a_{\mathrm{ij}}>0$ for all i and j. In this case, we write $A\gg 0$. In Appendix we use that a matrix A is called positive-off-diagonal if $a_{\mathrm{ij}}\ge 0$ for $i\ne j$ and $a_{\mathrm{ii}}\in \mathbb{R}$ [1].

The sequence of matrices $K=\{K(t){\}}_{t=0}^{\mathrm{\infty }}$ (with $K(0):=0$) forms a kernel. We assume

(A4)	$K(t)\ge 0$ for all $t\in \mathbb{N}$ with $K(t)>0$ for at least one $t\in \mathbb{N}$.
(A5)	$\sum _{t=0}^{\mathrm{\infty }}\Vert K(t)\Vert <\mathrm{\infty }$.

Remark 2.4

In a population dynamical context, it is preferable to equip ${\mathbb{R}}^m$ with the $l_1$-norm, i.e. to define $\Vert X\Vert =\sum _{j=1}^m|x_j|.$ In A5, we then use the corresponding operator norm, $\Vert K(t)\Vert =\underset{1\le j\le m}{max}\sum _{i=1}^m|K(t{)}_{\mathrm{ij}}|.$

If we now look for solutions of the RE (10(10) $B(t)=\sum _{s=1}^{\mathrm{\infty }}K(s)B(t-s),$(10) ) in the form $B(t)=z^t\mathrm{\Phi }$ for some m-vector Φ we observe that the Euler-Lotka equation takes the form (11) $r_{\sigma }(\overline{K}(z))=1,$(11) where $r_{\sigma }$ denotes the spectral radius and where, as before, (12) $\overline{K}(z)=\sum _{s=0}^{\mathrm{\infty }}K(s)z^{-s}.$(12) We call $\overline{K}(1)$ the next-generation matrix and define (13) ${\mathcal{R}}_0:=r_{\sigma }(\overline{K}(1)).$(13) Note (by exploiting the fact that the spectral radius of $\overline{K}(z)$ is a decreasing function of z) that the Euler-Lotka Equation (11(11) $r_{\sigma }(\overline{K}(z))=1,$(11) ) has a unique real solution ρ when ${\mathcal{R}}_0$ is larger than or equal to one. In that case, we necessarily have $1<\rho <{\mathcal{R}}_0$ or $1=\rho ={\mathcal{R}}_0.$ Furthermore, $\overline{K}(1)$ and $\overline{K}(\rho )$ are positive matrices, and therefore for both matrices the spectral radius is a dominant eigenvalue with a positive corresponding eigenvector [1, 2].

The key reference concerning the asymptotic behaviour of the solutions of the renewal equation (10(10) $B(t)=\sum _{s=1}^{\mathrm{\infty }}K(s)B(t-s),$(10) ) is now [25]. In order to present Thieme's renewal theorem, we again prescribe the history of B up to t = 0, (14) $B(\tau )=\mathrm{\Theta }(\tau ),\tau =\dots ,-2,-1$(14) for given positive $\mathrm{\Theta }(\tau )\in {\mathbb{R}}^m$ such that (15) $G(t):=\sum _{s=t+1}^{\mathrm{\infty }}K(s)\mathrm{\Theta }(t-s)=\sum _{l=1}^{\mathrm{\infty }}K(l+t)\mathrm{\Theta }(-l)$(15) defines a vector with finite components for all $t\in \mathbb{N}$. This allows us to write (10(10) $B(t)=\sum _{s=1}^{\mathrm{\infty }}K(s)B(t-s),$(10) ) as (16) $B(t)=\sum _{s=1}^{t}K(s)B(t-s)+G(t),t\in \mathbb{N}$(16) with the convention (17) $B(0)=G(0).$(17) Next, we introduce the resolvent of the kernel K. The resolvent $R=\{R(t){\}}_{t=0}^{\mathrm{\infty }}$ is defined by the equation $R=K\ast R+K=R\ast K+K$where * denotes the convolution product, i.e. (18) $K\ast R(t):=\sum _{s=0}^{t}K(s)R(t-s).$(18)

In our case $K(0)=0$ and consequently $R(0)=0$, meaning that we can write (18(18) $K\ast R(t):=\sum _{s=0}^{t}K(s)R(t-s).$(18) ) as $K\ast R(t)=\sum _{s=1}^{t-1}K(s)R(t-s).$The results in [25] concern a rather general setting in terms of positive linear operators on ordered Banach spaces. When dealing with spaces of functions defined on non-compact domains, it often helps to define positivity in terms of a special element w of the positive cone, with w possibly tending to zero for the argument tending to infinity. Here, however, we only deal with the standard positive cone in ${\mathbb{R}}^m$ and when invoking the results of [25] we simply take the m-vector $w=[1,\dots ,1{]}^T$.

We make further assumptions

(A6)	There exists $t_0\in \mathbb{N}$ such that $R(t_0)\gg 0$.
(A7)	There exist $t_1,t_2\in \mathbb{N}$ such that their greatest common divisor is one and $R(t_i)>0$ (i = 1, 2).

Theorem 2.5

Thieme's Renewal Theorem

Consider the renewal equation (16(16) $B(t)=\sum _{s=1}^{t}K(s)B(t-s)+G(t),t\in \mathbb{N}$(16) ) and assume that (A4), (A5), (A6) and (A7) hold. Furthermore, assume that there exists M>0 such that $\Vert G(t)\Vert \le M$ for all $t\in \mathbb{N}$. Assume that ${\mathcal{R}}_0\ge 1$. The Euler–Lotka equation (11(11) $r_{\sigma }(\overline{K}(z))=1,$(11) ) has a unique real solution ρ with either $1<\rho <{\mathcal{R}}_0$ or $\rho =1={\mathcal{R}}_0$. Furthermore let Φ, with $\Vert \mathrm{\Phi }\Vert =1$ and Ψ, with $\mathrm{\Psi \Phi }=1$, denote, respectively, the right and the left normalized positive eigenvector of $\overline{K}(\rho )$ corresponding to eigenvalue 1.

Then for $t\to \mathrm{\infty }$ (19a) ${\rho }^{-t}B(t)\to c\mathrm{\Phi }$(19a) with (19b) $c=\mathrm{\Psi }\frac{\overline{G}(\rho )}{c_0}\text{and}{c}_0=\mathrm{\Psi }\left(\sum _{s=1}^{\mathrm{\infty }}s{\rho }^{-s}K(s)\mathrm{\Phi }\right).$(19b)

Remark 2.6

1. The following table aims to facilitate comparison with Thieme's results in [25]:

   Table

   Here	Thieme [25]
   B	u
   K	A
   G	$\overline{u}$

2. Note that one can normalize Ψ in any way one wants, since in (19b(19b) $c=\mathrm{\Psi }\frac{\overline{G}(\rho )}{c_0}\text{and}{c}_0=\mathrm{\Psi }\left(\sum _{s=1}^{\mathrm{\infty }}s{\rho }^{-s}K(s)\mathrm{\Phi }\right).$(19b) ) Ψ is a factor in both the numerator and the denominator. Our choice, $\mathrm{\Psi \Phi }=1$, is handy in the context of many applications.

3. Structured-population models

When density dependence is ignored and the set of i-states is finite, a discrete time structured-population model leads to the system of linear recursion relations (20) $X(t)=(F+T)X(t-1).$(20) Here the components $x_j(t)$ of the vector $X(t)\in {\mathbb{R}}^n$ correspond to the density of females with i-state j at time t, while F and T are positive $n\times n$ matrices that describe, respectively, reproduction and survival $\mathrm{\&}$ development, cf. [3, 5]. More precisely, $F=(f_{\mathrm{ij}}{)}_{i,j=1}^n$ and $T=(t_{\mathrm{ij}}{)}_{i,j=1}^n$ where $\begin{array}{rl}{f}_{\mathrm{ij}}& =\mathit{the}\mathit{expected}\mathit{number}\mathit{of}\mathit{daughters}\mathit{with}\mathit{i}\mathit{-}\mathit{state}\mathit{i}\mathit{,}\\ & \mathit{produced}\mathit{by}\mathit{a}\mathit{female}\mathit{with}\mathit{i}\mathit{-}\mathit{state}\mathit{j}\end{array}$and $\begin{array}{rl}{t}_{\mathrm{ij}}& =\mathit{the}\mathit{probability}\mathit{that}\mathit{an}\mathit{individual}\mathit{with}\mathit{i}\mathit{-}\mathit{state}\mathit{j}\\ & \mathit{is}\mathit{alive}\mathit{and}\mathit{has}\mathit{i}\mathit{-}\mathit{state}\mathit{i}\mathit{one}\mathit{time}\mathit{unit}\mathit{later}\mathit{.}\end{array}$Since no individual is immortal we assume that $r_{\sigma }(T)<1$.

The aim of this section is twofold: (i) we show how the recurrence (20(20) $X(t)=(F+T)X(t-1).$(20) ) and the renewal equation (10(10) $B(t)=\sum _{s=1}^{\mathrm{\infty }}K(s)B(t-s),$(10) ) relate to each other and (ii) we discuss how the next-generation matrix and the basic reproduction number corresponding to (20(20) $X(t)=(F+T)X(t-1).$(20) ) are connected with (13(13) ${\mathcal{R}}_0:=r_{\sigma }(\overline{K}(1)).$(13) ). We begin with the former.

Suppose there are m states-at-birth in the literal sense ($1\le m\le n$). There exist a positive $n\times m$ matrix V (normalized such that all m columns of V have $l_1$-norm equal to one) and a positive $m\times n$ matrix U such that (21) $F=\mathrm{VU}.$(21) In particular, if the states-at-birth are $j_1,\dots ,j_m\in \{1,\dots ,n\}$ we define V as the matrix with columns $e_{j_1},\dots ,e_{j_m}$ and U as the fertility matrix obtained by taking the rows $j_1,\dots ,j_m$ from F.

Directly from the interpretation we discover the first relation, namely, (22) $B(t)=\mathrm{UX}(t-1).$(22) Using the generation expansion we then deduce from (20(20) $X(t)=(F+T)X(t-1).$(20) ), (21(21) $F=\mathrm{VU}.$(21) ) and (22(22) $B(t)=\mathrm{UX}(t-1).$(22) ) that $\begin{array}{rl}X(t)& =\mathrm{VB}(t)+\mathrm{TX}(t-1)\\ & =\mathrm{VB}(t)+\mathrm{TVB}(t-1)+T^{2}X(t-2)\\ & =\dots ,\end{array}$leading in the limit to the second relation (23) $X(t)=\sum _{s=0}^{\mathrm{\infty }}{T}^s\mathrm{VB}(t-s).$(23) The two Equations (22(22) $B(t)=\mathrm{UX}(t-1).$(22) ) and (23(23) $X(t)=\sum _{s=0}^{\mathrm{\infty }}{T}^s\mathrm{VB}(t-s).$(23) ) can both be understood on the basis of their interpretation and together they provide a complete description of the dynamics. By substituting (23(23) $X(t)=\sum _{s=0}^{\mathrm{\infty }}{T}^s\mathrm{VB}(t-s).$(23) ) into (22(22) $B(t)=\mathrm{UX}(t-1).$(22) ) one obtains (10(10) $B(t)=\sum _{s=1}^{\mathrm{\infty }}K(s)B(t-s),$(10) ) with (24) $K(s)=U{T}^{s-1}V.$(24) Likewise one recovers (20(20) $X(t)=(F+T)X(t-1).$(20) ) from (22(22) $B(t)=\mathrm{UX}(t-1).$(22) ) and (23(23) $X(t)=\sum _{s=0}^{\mathrm{\infty }}{T}^s\mathrm{VB}(t-s).$(23) ) by splitting off the $s=0$ term in (23(23) $X(t)=\sum _{s=0}^{\mathrm{\infty }}{T}^s\mathrm{VB}(t-s).$(23) ) to obtain (25) $X(t)=\mathrm{VB}(t)+\mathrm{TX}(t-1)$(25) and next use (22(22) $B(t)=\mathrm{UX}(t-1).$(22) ) and (21(21) $F=\mathrm{VU}.$(21) ).

The two Equations (22(22) $B(t)=\mathrm{UX}(t-1).$(22) ) and (23(23) $X(t)=\sum _{s=0}^{\mathrm{\infty }}{T}^s\mathrm{VB}(t-s).$(23) ) relate the time course of B to that of X and vice versa, when defined for all t, not just for $t\ge 0$. To deal with the initial value problem for X, we can rewrite (20(20) $X(t)=(F+T)X(t-1).$(20) ) in the variation-of-constants form (26) $X(t)=T^{t}X(0)+\sum _{s=1}^t{T}^{s-1}\mathrm{FX}(t-s).$(26) By acting on both sides of this identity with F we obtain (27) $Y(t)=F{T}^{t}X(0)+\sum _{s=1}^{t}F{T}^{s-1}Y(t-s)$(27) where $Y(t)=\mathrm{FX}(t).$Once we solve the RE (27(27) $Y(t)=F{T}^{t}X(0)+\sum _{s=1}^{t}F{T}^{s-1}Y(t-s)$(27) ), we can rewrite (26(26) $X(t)=T^{t}X(0)+\sum _{s=1}^t{T}^{s-1}\mathrm{FX}(t-s).$(26) ) as the explicit expression (28) $X(t)=T^{t}X(0)+\sum _{s=1}^t{T}^{s-1}Y(t-s).$(28) Here Y takes values in ${\mathbb{R}}^n$, but whenever F has the form (21(21) $F=\mathrm{VU}.$(21) ), we know that (29) $Y(t)=\mathrm{VB}(t)$(29) for a function B taking values in ${\mathbb{R}}^m$. We can then rewrite (27(27) $Y(t)=F{T}^{t}X(0)+\sum _{s=1}^{t}F{T}^{s-1}Y(t-s)$(27) ) in the form (30) $B(t)=U{T}^{t}X(0)+\sum _{s=1}^{t}U{T}^{s-1}\mathrm{VB}(t-s),$(30) that is, we get (16(16) $B(t)=\sum _{s=1}^{t}K(s)B(t-s)+G(t),t\in \mathbb{N}$(16) ) with (31a) $\begin{array}{rl}K(s)& =U{T}^{s-1}V,\end{array}$(31a) (31b) $\begin{array}{rl}G(s)& =U{T}^{s}X(0).\end{array}$(31b) Note that when the rank of U is not maximal, there are fewer states-at-birth in the stochastic sense than there are in the literal sense. In such a case, a further reduction in RE system size is possible with a different choice of U and V (see one such example below).

Hence, if the conditions of Theorem 2.5 (or Theorem 2.1 in the case where (30(30) $B(t)=U{T}^{t}X(0)+\sum _{s=1}^{t}U{T}^{s-1}\mathrm{VB}(t-s),$(30) ) amounts to a scalar equation) are satisfied, we can first deduce the asymptotic behaviour of $B(t)$ from that theorem and then use (28(28) $X(t)=T^{t}X(0)+\sum _{s=1}^t{T}^{s-1}Y(t-s).$(28) ), with Y given by (29(29) $Y(t)=\mathrm{VB}(t)$(29) ), to determine the large time behaviour of the population state $X(t)$. We work out the details of this procedure in the two subsections that follow. But first, we discuss how the next-generation matrix and the basic reproduction number corresponding to (20(20) $X(t)=(F+T)X(t-1).$(20) ) are connected with (13(13) ${\mathcal{R}}_0:=r_{\sigma }(\overline{K}(1)).$(13) ).

For $t\in \mathbb{N}$, the $(i,j)$th element of $T^t$ is the probability that an individual with i-state j is alive t units of time later and has then i-state i. If $r_{\sigma }(T)<1$ then $(I-T{)}^{-1}$ exists and the $(i,j)$th element of $(I-T{)}^{-1}=\sum _{t=0}^{\mathrm{\infty }}{T}^t$equals the expected time an individual starting with i-state j will spend in i-state i.

The matrix (32) ${\mathcal{L}}_{\mathrm{L}}=F(I-T{)}^{-1}$(32) is the next-generation matrix (NGM) with large domain, with the $(i,j)$th element of ${\mathcal{L}}_{\mathrm{L}}$ specifying the expected number of future offspring with i-state i produced by an individual presently having i-state j. Often one defines (33) ${\mathcal{R}}_0=r_{\sigma }({\mathcal{L}}_L).$(33) The elements of the NGM with large domain (32(32) ${\mathcal{L}}_{\mathrm{L}}=F(I-T{)}^{-1}$(32) ) specify the expected numbers of future offspring of individuals for all i-states. But, with the factorization of F as described directly after (21(21) $F=\mathrm{VU}.$(21) ), we can focus on the true NGM (34) $\mathcal{L}=U(I-T{)}^{-1}V.$(34) We invite the reader to verify that the $(i,j)$th element of $\mathcal{L}$ gives the expected lifetime number of offspring with state-at-birth i produced by a newborn individual with state-at-birth j.

Note that when the rank of U is smaller than m (which happens whenever there are fewer states-at-birth in the stochastic sense than there are in the literal sense), a further reduction to the NGM with small domain ${\mathcal{L}}_{\mathrm{S}}$ is possible by altering U and V (see one example below and also [14] for the derivation of NGM with small and large domains in the context of epidemiological models in continuous time).

One easily verifies that

1. the NGM with large domain (32(32) ${\mathcal{L}}_{\mathrm{L}}=F(I-T{)}^{-1}$(32) ) and the NGM (34(34) $\mathcal{L}=U(I-T{)}^{-1}V.$(34) ) have the same non-zero eigenvalues and therefore (33(33) ${\mathcal{R}}_0=r_{\sigma }({\mathcal{L}}_L).$(33) ) amounts to (35) ${\mathcal{R}}_0=r_{\sigma }(\mathcal{L}).$(35) (And when further reduction to ${\mathcal{L}}_{\mathrm{S}}$ is possible then also ${\mathcal{R}}_0=r_{\sigma }({\mathcal{L}}_{\mathrm{S}}).$)

2. If Φ is a right eigenvector of $\mathcal{L}$ corresponding to a non-zero eigenvalue λ then $\tilde{\mathrm{\Phi }}=V\mathrm{\Phi }$ is a right eigenvector of ${\mathcal{L}}_{\mathrm{L}}$ corresponding to eigenvalue λ.

3. If $\tilde{\mathrm{\Psi }}$ is a left eigenvector of ${\mathcal{L}}_{\mathrm{L}}$ corresponding to a non-zero eigenvalue λ then $\mathrm{\Psi }=\tilde{\mathrm{\Psi }}V$ is a left eigenvector of $\mathcal{L}$ corresponding to eigenvalue λ.

Furthermore, when $r_{\sigma }(T)<1$ we deduce from (12(12) $\overline{K}(z)=\sum _{s=0}^{\mathrm{\infty }}K(s)z^{-s}.$(12) ) and (31a(31a) $\begin{array}{rl}K(s)& =U{T}^{s-1}V,\end{array}$(31a) ) that (36) $\overline{K}(1)=U(I-T{)}^{-1}V.$(36) That is, when U and V are defined as described after (21(21) $F=\mathrm{VU}.$(21) ) then $\overline{K}(1)$ is the NGM and we recover (13(13) ${\mathcal{R}}_0:=r_{\sigma }(\overline{K}(1)).$(13) ).

3.1. One state-at-birth

The simplest situation arises when there is only one state-at-birth, either in the deterministic sense where all individuals are identical at birth (e.g. when age is the only structuring variable) or in the stochastic sense where states-at-birth have a given probability distribution. In that case F has one dimensional range.

Then F = V U, where $V\in {\mathbb{R}}^n$ (with $\Vert V\Vert =1$) spans the range of F and U corresponds to the row n-vector (i.e. a $1\times n$ matrix) of fertility rates. In the notation of Section 2.1 we then define $b(t)=\mathrm{UX}(t-1)$and consider the scalar RE (7(7) $b(t)=\sum _{s=1}^{t}k(s)b(t-s)+g(t),t\in \mathbb{N},$(7) ) with (37a) $\begin{array}{rl}k(t)& =U{T}^{t-1}V\end{array}$(37a) (37b) $\begin{array}{rl}g(t)& =U{T}^{t}X(0).\end{array}$(37b) Let's assume that $\{k(t){\}}_{t=1}^{\mathrm{\infty }}$ and $\{g(t){\}}_{t=1}^{\mathrm{\infty }}$ thus defined satisfy the conditions in Theorem 2.1.

When $r_{\sigma }(T)<1$ we have (38) $\overline{k}(z)=U(\mathrm{zI}-T{)}^{-1}V$(38) for $|z|\ge 1$ and in particular (39) ${\mathcal{R}}_0=U(I-T{)}^{-1}V.$(39) If ${\mathcal{R}}_0<1$ then $b(t)\to 0$ as $t\to \mathrm{\infty }$ and consequently $X(t)\to 0$ as $t\to \mathrm{\infty }$.

Let's focus on the case ${\mathcal{R}}_0\ge 1$. We can then compute the real time growth rate ρ as the unique solution of the Euler–Lotka equation $\overline{k}(z)=1,$with necessarily either $\rho =1={\mathcal{R}}_0$ or $1<\rho <{\mathcal{R}}_0$. Since (38(38) $\overline{k}(z)=U(\mathrm{zI}-T{)}^{-1}V$(38) ) can be written as (40) $\overline{k}(z)=z^{-1}U(I-z^{-1}T{)}^{-1}V$(40) we observe that the population growth rate ρ is the value of z for which both the generation and the real time process become stationary when we multiply all fertilities and all transition probabilities by $z^{-1}$. When $\rho >1$ the latter can be interpreted as introducing an additional, state-independent, death probability per time step $1-{\rho }^{-1}$.

Then for $t\to \mathrm{\infty }$ $b(t)=c{\rho }^t+o({\rho }^t)$for some constant c and from (28(28) $X(t)=T^{t}X(0)+\sum _{s=1}^t{T}^{s-1}Y(t-s).$(28) ) and (29(29) $Y(t)=\mathrm{VB}(t)$(29) ) it now follows that $\begin{array}{rl}X(t)& =c{\rho }^t\sum _{s=1}^t{\rho }^{-s}{T}^{s-1}V+o({\rho }^t)\\ & =c{\rho }^t(\mathrm{\rho I}-T{)}^{-1}V+o({\rho }^t).\end{array}$That is, the solution $X(t)$ of (20(20) $X(t)=(F+T)X(t-1).$(20) ) with $X(0)>0$ grows, for large t, geometrically with rate ρ while converging to the stable distribution (41) ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (\mathrm{\rho I}-T{)}^{-1}V\Vert }(\mathrm{\rho I}-T{)}^{-1}V.$(41) Here the index r refers to ‘real time’ (as opposed to ‘generation’). By rewriting (41(41) ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (\mathrm{\rho I}-T{)}^{-1}V\Vert }(\mathrm{\rho I}-T{)}^{-1}V.$(41) ) as (42) ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (I-{\rho }^{-1}T{)}^{-1}V\Vert }(I-{\rho }^{-1}T{)}^{-1}V$(42) we can (with the $l_1$ norm) interpret the denominator as the life expectancy of a newborn individual when we introduce into the system an additional, state-independent, death probability per time step $1-{\rho }^{-1}$ (or, in other words, discount again and again the next year by a factor ${\rho }^{-1}$).

In the spirit of (9b(9b) $c(l):=\sum _{s=0}^{\mathrm{\infty }}{\rho }^{-s}k(l+s).$(9b) ) we next ask: how do the components of $X(0)$ contribute to future population sizes, i.e. to the constant c? From (37b(37b) $\begin{array}{rl}g(t)& =U{T}^{t}X(0).\end{array}$(37b) ) it follows that (43) $\overline{g}(z)=\mathrm{zU}(\mathrm{zI}-T{)}^{-1}X(0).$(43) Now recall Theorem 2.1(iii) and conclude that the j-th component of the row vector (44) ${\mathrm{\Psi }}_{\mathrm{r}}=\mathrm{cU}(\mathrm{\rho I}-T{)}^{-1}$(44) specifies the (relative) contribution of individuals with i-state j to future population size. Or, in the now generally accepted terminology introduced by Fisher [18], specifies the reproductive value of state j.

It is straightforward to check that

1. F + T has dominant eigenvalue ρ with the stable distribution ${\mathrm{\Phi }}_{\mathrm{r}}$ as the corresponding right eigenvector and the vector of reproductive values ${\mathrm{\Psi }}_{\mathrm{r}}$ as the corresponding left eigenvector.

2. ${\mathrm{\Psi }}_{\mathrm{r}}$ is also a left eigenvector of the scaled NGM with large domain ${\mathcal{L}}_{\mathrm{L}}^{(\rho )}:=F(\mathrm{\rho I}-T{)}^{-1}={\rho }^{-1}F(I-{\rho }^{-1}T{)}^{-1}$corresponding to eigenvalue 1. That is, it is also (modulo normalization) the vector of generation-based reproductive values when we discount fertilities and transition probabilities by ${\rho }^{-1}$.

We normalize the vector of reproductive values by requiring (45) ${\mathrm{\Psi }}_{\mathrm{r}}{\mathrm{\Phi }}_{\mathrm{r}}=1$(45) and obtain (46) ${\mathrm{\Psi }}_{\mathrm{r}}=\frac{\Vert (\mathrm{\rho I}-T{)}^{-1}V\Vert }{U(\mathrm{\rho I}-T{)}^{-2}V}U(\mathrm{\rho I}-T{)}^{-1}.$(46) If we now observe that $(I-{\rho }^{-1}T{)}^{-2}=I+2{\rho }^{-1}T+3{\rho }^{-2}{T}^2+\dots ,$we can rewrite ${\mathrm{\Psi }}_{\mathrm{r}}$ as (47) ${\mathrm{\Psi }}_{\mathrm{r}}=\frac{\Vert (I-{\rho }^{-1}T{)}^{-1}V\Vert }{U(I-{\rho }^{-1}T{)}^{-2}V}U(I-{\rho }^{-1}T{)}^{-1}$(47) and interpret the denominator as ${\rho }^{-1}$ times the expected discounted age of the parent of a newborn individual.

Example 3.1

Consider the Leslie model with three age classes, i.e. $T=\left[\begin{array}{ccc}0& 0& 0\\ t_{21}& 0& 0\\ 0& t_{32}& 0\end{array}\right]\text{and}F=\left[\begin{array}{ccc}{f}_1& f_2& f_3\\ 0& 0& 0\\ 0& 0& 0\end{array}\right].$(see Figure 2 for the schematic representation of the model). All individuals are born into the first age class so F = V U for $V=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]\text{and}U=\left[\begin{array}{ccc}{f}_1& f_2& f_3\end{array}\right].$A straightforward calculation reveals that $(\mathrm{zI}-T{)}^{-1}=\left[\begin{array}{ccc}{z}^{-1}& 0& 0\\ t_{21}{z}^{-2}& z^{-1}& 0\\ t_{32}{t}_{21}{z}^{-3}& t_{32}{z}^{-2}& z^{-1}\end{array}\right],$which yields the Euler–Lotka equation $\overline{k}(z)=f_1{z}^{-1}+f_2{t}_{21}{z}^{-2}+f_3{t}_{32}{t}_{21}{z}^{-3}=1.$Hence $\begin{array}{rl}1& =f_1{\rho }^{-1}+f_2{t}_{21}{\rho }^{-2}+f_3{t}_{32}{t}_{21}{\rho }^{-3},\\ {\mathcal{R}}_0& =f_1+f_2{t}_{21}+f_3{t}_{32}{t}_{21}.\end{array}$Using (42(42) ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (I-{\rho }^{-1}T{)}^{-1}V\Vert }(I-{\rho }^{-1}T{)}^{-1}V$(42) ), we obtain the asymptotic age distribution ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{1+t_{21}{\rho }^{-1}+t_{32}{t}_{21}{\rho }^{-2}}\left[\begin{array}{c}1\\ t_{21}{\rho }^{-1}\\ t_{32}{t}_{21}{\rho }^{-2}\end{array}\right],$while (47(47) ${\mathrm{\Psi }}_{\mathrm{r}}=\frac{\Vert (I-{\rho }^{-1}T{)}^{-1}V\Vert }{U(I-{\rho }^{-1}T{)}^{-2}V}U(I-{\rho }^{-1}T{)}^{-1}$(47) ) gives the row vector of reproductive values $\begin{array}{rl}{\mathrm{\Psi }}_{\mathrm{r}}& =c\left[\begin{array}{ccc}{f}_1{\rho }^{-1}+f_2{t}_{21}{\rho }^{-2}+f_3{t}_{32}{t}_{21}{\rho }^{-3},& f_3{t}_{32}{\rho }^{-2}+f_2{\rho }^{-1},& f_3{\rho }^{-1}\end{array}\right]\\ & =c\left[\begin{array}{ccc}1,& f_3{t}_{32}{\rho }^{-2}+f_2{\rho }^{-1},& f_3{\rho }^{-1}\end{array}\right]\end{array}$where the constant c is determined such that ${\mathrm{\Psi }}_{\mathrm{r}}{\mathrm{\Phi }}_{\mathrm{r}}=1$.

Figure 2. A schematic representation of the Leslie model in Example 3.1.

For the specific case with $t_{21}=0.7,t_{32}=0.8,f_1=0,f_2=1=f_3$ we find that $\rho =1.1$ and ${\mathcal{R}}_0=1.26$ (see Figure 3(a)). The asymptotic age distribution is ${\mathrm{\Phi }}_{\mathrm{r}}=\left[\begin{array}{c}0.476\\ 0.303\\ 0.221\end{array}\right],$while the normalized vector of reproductive values is ${\mathrm{\Psi }}_{\mathrm{r}}=\left[\begin{array}{ccc}0.867,& 1.363,& 0.788\end{array}\right]$(see panels (b) and (c) of Figure 3). Note that if we define $N(t,a)=x_a(t)$ we can write $\begin{array}{rl}N(t+1,a+1)& =t_{(a+1)a}N(t,a)\\ N(t+1,1)& =\sum _{a=1}^3{f}_{a}N(t,a)\end{array}$which clearly exposes the Leslie matrix model as the analogue of the PDE formulation of a continuous age and time model. The RE (7(7) $b(t)=\sum _{s=1}^{t}k(s)b(t-s)+g(t),t\in \mathbb{N},$(7) ) is, of course, the analogue of Lotka's RE in the continuous age and time setting. Note also that $b(t)=N(t,1).$

Figure 3. The plots of (a) $z\mapsto \overline{k}(z)$, (b) the asymptotic age distribution and (c) the reproductive values for the (concrete example of the) Leslie model in Example 3.1.

3.2. More than one state-at-birth

Even when individuals may be born in different i-states, the number of individual states-at-birth is typically significantly lower than the total number of i-states, thus making the much smaller system sizes of RE an attractive alternative to (20(20) $X(t)=(F+T)X(t-1).$(20) ).

Suppose there are m states-at-birth with 1<m<n (this can now be understood literally, or, if aiming for the minimal RE system size, in a stochastic sense). We then write F = V U for some positive $n\times m$ matrix V (normalized such that all m columns of V have $l_1$-norm equal to one) and a positive $m\times n$ matrix U and consider the RE (16(16) $B(t)=\sum _{s=1}^{t}K(s)B(t-s)+G(t),t\in \mathbb{N}$(16) ) with $\{K(s){\}}_{s=1}^{\mathrm{\infty }}$ and $\{G(s){\}}_{s=1}^{\mathrm{\infty }}$ as in (31a(31a) $\begin{array}{rl}K(s)& =U{T}^{s-1}V,\end{array}$(31a) ). Next, require that the assumptions of Theorem 2.5 hold.

If $r_{\sigma }(T)<1$ the asymptotic dynamics of (30(30) $B(t)=U{T}^{t}X(0)+\sum _{s=1}^{t}U{T}^{s-1}\mathrm{VB}(t-s),$(30) ) is completely determined by the second term. Furthermore, (48) $\overline{K}(z)=U(\mathrm{zI}-T{)}^{-1}V$(48) for $|z|\ge 1$. In particular, $\overline{K}(1)$ is the next-generation matrix (possibly with small domain) and (49) ${\mathcal{R}}_0=r_{\sigma }(U(I-T{)}^{-1}V).$(49) When ${\mathcal{R}}_0<1$ we have $B(t)\to 0$ as $t\to \mathrm{\infty }$ and consequently $X(t)\to 0$ as well.

Suppose now that ${\mathcal{R}}_0\ge 1$. The real time growth rate ρ is determined as the unique solution of the Euler-Lotka equation $r_{\sigma }(\overline{K}(z))=1,$where again either $\rho =1={\mathcal{R}}_0$ or $1<\rho <{\mathcal{R}}_0$. Note that rewriting (48(48) $\overline{K}(z)=U(\mathrm{zI}-T{)}^{-1}V$(48) ) as (50) $\overline{K}(z)=z^{-1}U(I-z^{-1}T{)}^{-1}V,$(50) we can once more observe that the population growth rate ρ is the value of z for which both the generation and the real time process become stationary when we multiply all fertilities and all transition probabilities by $z^{-1}$ (and when $\rho >1$ the latter can be interpreted as introducing an additional, state-independent, death probability per time step $1-{\rho }^{-1}$). Note also that the matrix $\overline{K}(\rho )$ is the NGM (or NGM with small domain) when we discount fertilities and transition rates by a factor ${\rho }^{-1}$, $\overline{K}(\rho )={\rho }^{-1}U(I-{\rho }^{-1}T{)}^{-1}V.$The m-vector $B(t)$ will for large t grow geometrically with rate ρ while converging to the distribution Φ, where Φ is the normalized positive right eigenvector of $\overline{K}(\rho )$ corresponding to eigenvalue 1. We can then deduce from (28(28) $X(t)=T^{t}X(0)+\sum _{s=1}^t{T}^{s-1}Y(t-s).$(28) ) and (29(29) $Y(t)=\mathrm{VB}(t)$(29) ) that, for ${\mathcal{R}}_0>1$ and a non-trivial initial condition $X(0)>0$, the vector $X(t)$ grows geometrically with rate ρ while converging to the asymptotic distribution (51a) $\begin{array}{rl}{\mathrm{\Phi }}_{\mathrm{r}}& =\frac{1}{\Vert (\mathrm{\rho I}-T{)}^{-1}V\mathrm{\Phi }\Vert }(\mathrm{\rho I}-T{)}^{-1}V\mathrm{\Phi }\end{array}$(51a) (51b) $\begin{array}{rl}& =\frac{1}{\Vert (I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }\Vert }(I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }.\end{array}$(51b) When using the $l_1$ norm, we can interpret the denominator in (51b(51b) $\begin{array}{rl}& =\frac{1}{\Vert (I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }\Vert }(I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }.\end{array}$(51b) ) as the expected duration of life of a newborn when newborns are sampled from the asymptotic distribution Φ and with additional, state-independent, death probability $1-{\rho }^{-1}$ introduced into the system (or alternatively, when we discount again and again the next year by ${\rho }^{-1}$).

To see how the components of $X(0)$ contribute to future population growth (i.e. to the constant c in (19b(19b) $c=\mathrm{\Psi }\frac{\overline{G}(\rho )}{c_0}\text{and}{c}_0=\mathrm{\Psi }\left(\sum _{s=1}^{\mathrm{\infty }}s{\rho }^{-s}K(s)\mathrm{\Phi }\right).$(19b) )) we compute (52) $\overline{G}(\rho )=\mathrm{\rho U}(\mathrm{\rho I}-T{)}^{-1}.$(52) Using (19b(19b) $c=\mathrm{\Psi }\frac{\overline{G}(\rho )}{c_0}\text{and}{c}_0=\mathrm{\Psi }\left(\sum _{s=1}^{\mathrm{\infty }}s{\rho }^{-s}K(s)\mathrm{\Phi }\right).$(19b) ) we now conclude that, if Ψ is a positive left eigenvector of $\overline{K}(\rho )$ corresponding to eigenvalue 1 then the j-th element of the row n-vector (53) ${\mathrm{\Psi }}_{\mathrm{r}}=c\mathrm{\Psi U}(\mathrm{\rho I}-T{)}^{-1}$(53) specifies the (relative) contribution of individuals with i-state j to future population growth. Again, we normalize this vector by requiring ${\mathrm{\Psi }}_{\mathrm{r}}{\mathrm{\Phi }}_{\mathrm{r}}=1$. Thus we obtain (54a) $\begin{array}{rl}{\mathrm{\Psi }}_{\mathrm{r}}& =\frac{\Vert (\mathrm{\rho I}-T{)}^{-1}V\mathrm{\Phi }\Vert }{\mathrm{\Psi }U(\mathrm{\rho I}-T{)}^{-2}V\mathrm{\Phi }}\mathrm{\Psi }U(\mathrm{\rho I}-T{)}^{-1}\end{array}$(54a) (54b) $\begin{array}{rl}& =\frac{\Vert (I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }\Vert }{\mathrm{\Psi }U(I-{\rho }^{-1}T{)}^{-2}V\mathrm{\Phi }}\mathrm{\Psi }U(I-{\rho }^{-1}T{)}^{-1}.\end{array}$(54b) With the same reasoning as in the previous subsection, we can now interpret the j-th component of $U(I-{\rho }^{-1}T{)}^{-2}V\mathrm{\Phi }$ in the denominator of (54b(54b) $\begin{array}{rl}& =\frac{\Vert (I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }\Vert }{\mathrm{\Psi }U(I-{\rho }^{-1}T{)}^{-2}V\mathrm{\Phi }}\mathrm{\Psi }U(I-{\rho }^{-1}T{)}^{-1}.\end{array}$(54b) ) as ${\rho }^{-1}$ times the expected age of the parent of a newborn individual with i-state j when we sample parents according to Φ and take into account an additional, state-independent, death probability $1-{\rho }^{-1}$.

Again, one easily verifies that

1. F + T has dominant eigenvalue ρ with the stable distribution ${\mathrm{\Phi }}_{\mathrm{r}}$ as the corresponding right eigenvector and the vector of reproductive values ${\mathrm{\Psi }}_{\mathrm{r}}$ as the corresponding left eigenvector.

2. ${\mathrm{\Psi }}_{\mathrm{r}}$ is also a left eigenvector of the scaled NGM with large domain ${\mathcal{L}}_{\mathrm{L}}^{(\rho )}:=F(\mathrm{\rho I}-T{)}^{-1}={\rho }^{-1}F(I-{\rho }^{-1}T{)}^{-1}$corresponding to eigenvalue 1. That is, it is also the vector of (generation-based) reproductive values when we discount fertilities and transition probabilities by ${\rho }^{-1}$.

Note that the conditions $A6$ and $A7$ imposed on the resolvent can be interpreted biologically. Indeed, by definition the kernel $K=\{K(t){\}}_{t=1}^{\mathrm{\infty }}$ contains information about the expected numbers of daughters produced by one female at various ages. The matrix $K^{2\ast }=K\ast K$ yields the expected numbers of granddaughters (i.e. second generation offspring) and the j-th convolution of K with itself, $K^{j\ast }$, yields the expected number of j-th generation offspring. Summing over all generations of offspring we obtain the clan kernel, $K^{(c)}:=\sum _{j=1}^{\mathrm{\infty }}{K}^{j\ast }.$Since clan members of every female are either her daughters or clan members of one of her daughters, or, alternatively, either her daughters or daughters of a member of the clan, we must have $K^{(c)}=K+K\ast K^{(c)}=K+K^{(c)}\ast K.$That is, the clan kernel $K^{(c)}$ is the resolvent of K. The assumptions $A6$ and $A7$ can therefore be interpreted as follows:

(A6)	Consider a newborn individual. At some time $t_0$ after her birth, clan members with any state-at-birth will be born, irrespective of the birth state of the focus individual.
(A7)	Consider two newborn individuals. There exist $t_1,t_2$ with greatest common divisor one, such that it is possible to choose the birth states of the focus individuals so that at times $t_1$ and $t_2$ after their birth, clan members of at least one of them will be born.

As a summary of the results, we provide at the end of this Section an algorithm for deducing the large time behaviour of (discrete time) structured-population models (SPM) via the renewal equation formulation (and for an application see the example below).

Example 3.2

Consider a population in which individuals are characterized by both age a and some indicator of size s (for example, both age and size are important determinants of population dynamics of Rhododendron maximum shoots, see [3, 23]) and let $\{(a,s):a\in \{1,2,3\},s\in \{1,2\}\}$be the set of i-states. We enumerate the i-states in the following way:

Table

(a,s)	(1,1)	(2,1)	(3,1)	(1,2)	(2,2)	(3,2)
label	1	2	3	4	5	6

.

Let (55) $T=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ t_{21}& 0& 0& 0& 0& 0\\ 0& t_{32}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ t_{51}& 0& 0& t_{54}& 0& 0\\ 0& t_{62}& 0& 0& t_{65}& 0\end{array}\right]\text{and}F=\left[\begin{array}{cccccc}0& f_{12}& f_{13}& 0& f_{15}& f_{16}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& f_{45}& f_{46}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$(55) describe the survival $\mathrm{\&}$ state transitions and fertility rates, respectively (see Figure 4 for a schematic representation of the model). There are two states at birth, i.e. $(1,1)$ and $(1,2)$, corresponding to labels 1 and 4, respectively. Then F = V U for (56) $V=\left[\begin{array}{cc}1& 0\\ 0& 0\\ 0& 0\\ 0& 1\\ 0& 0\\ 0& 0\end{array}\right]\text{and}U=\left[\begin{array}{cccccc}0& f_{12}& f_{13}& 0& f_{15}& f_{16}\\ 0& 0& 0& 0& f_{45}& f_{46}\end{array}\right].$(56) The next-generation matrix is $\begin{array}{rl}\mathcal{L}& =U(I-T{)}^{-1}V\\ & =\left[\begin{array}{cc}{f}_{12}{t}_{21}+f_{13}{t}_{32}{t}_{21}+f_{15}{t}_{51}+f_{16}(t_{62}{t}_{21}+t_{65}{t}_{51})& f_{15}{t}_{54}+f_{16}{t}_{65}{t}_{54}\\ f_{45}{t}_{51}+f_{46}(t_{62}{t}_{21}+t_{65}{t}_{51})& f_{45}{t}_{54}+f_{46}{t}_{65}{t}_{54}\end{array}\right]\end{array}$and ${\mathcal{R}}_0$ is the dominant eigenvalue of $\mathcal{L}.$

Let us first consider the case where all the $t^{\prime }$s and the $f^{\prime }$s are strictly positive. Using (31a(31a) $\begin{array}{rl}K(s)& =U{T}^{s-1}V,\end{array}$(31a) ) we find $\begin{array}{rl}K(1)& =0,\\ K(2)& =\left[\begin{array}{cc}{f}_{12}{t}_{21}+f_{15}{f}_{51}& f_{15}{t}_{54}\\ f_{45}{t}_{51}& f_{45}{t}_{54}\end{array}\right]\gg 0\\ K(3)& =\left[\begin{array}{cc}{f}_{13}{t}_{32}{t}_{21}+f_{16}(t_{62}{t}_{21}+t_{65}{t}_{51})& f_{16}{t}_{65}{t}_{54}\\ f_{46}(t_{62}{t}_{21}+t_{65}{t}_{51})& f_{46}{t}_{65}{t}_{54}\end{array}\right]\gg 0\\ K(t)& =0\mathrm{for}t\ge 4\end{array}$and using (31b(31b) $\begin{array}{rl}G(s)& =U{T}^{s}X(0).\end{array}$(31b) ) that $G(t)=0$ for $t\ge 3$. Therefore $\begin{array}{rl}R(1)& =K(1)=0,\\ R(2)& =K(2)\gg 0\\ R(3)& =K(3)\gg 0\\ R(t)& \ge K(2)R(t-2)\gg 0\mathrm{for}t\ge 4.\end{array}$The positivity of the resolvent for $t\ge 2$ can also be deduced by interpretation. Indeed, if all the t's and the f's are strictly positive then individuals with either birth state will produce offspring of either size at ages 2 and 3. For $t\ge 2$ there will therefore be born clan members with any state-at-birth, irrespective of the birth state of the mother.

Figure 4. A schematic representation of (a) the state transitions and (b) fertilities in the model described in Example 3.2.

All the assumptions of Theorem 2.5 are therefore satisfied. We have $\overline{K}(z)=\sum _{s=1}^{\mathrm{\infty }}K(s)z^{-s}=K(2)z^{-2}+K(3)z^{-3}.$Let's assume that ${\mathcal{R}}_0>1$. We can then compute the growth rate $\rho >1$ as the unique real solution z of the Euler-Lotka equation $r_{\sigma }(\overline{K}(z))=1.$As a concrete example consider $T=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0.3& 0& 0& 0& 0& 0\\ 0& 0.3& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0.3& 0& 0& 0.5& 0& 0\\ 0& 0.3& 0& 0& 0.5& 0\end{array}\right],F=\left[\begin{array}{cccccc}0& 1& 1& 0& 1& 1\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right].$Then $\overline{K}(z)=\left[\begin{array}{cc}0.6z^{-2}+0.33z^{-3}& 0.5z^{-2}+0.25z^{-3}\\ 0.3z^{-2}+0.24z^{-3}& 0.5z^{-2}+0.25z^{-3}\end{array}\right].$We can find the solution of the Euler-Lotka equation numerically and get $\rho =1.18$, while ${\mathcal{R}}_0=1.48$ (see Figure 5(a)). The asymptotic distribution of states-at-birth is found as the ( normalized) positive eigenvector of $\overline{K}(\rho )$ corresponding to eigenvalue 1. We obtain $\mathrm{\Phi }=\left[\begin{array}{c}0.578\\ 0.422\end{array}\right]$and then from (51b(51b) $\begin{array}{rl}& =\frac{1}{\Vert (I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }\Vert }(I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }.\end{array}$(51b) ) conclude that the asymptotic population distribution is given by ${\mathrm{\Phi }}_{\mathrm{r}}=\left[\begin{array}{c}0.343\\ 0.087\\ 0.022\\ 0.251\\ 0.193\\ 0.104\end{array}\right].$We furthermore find that the left eigenvector of $\overline{K}(\rho )$ corresponding to eigenvalue 1 (normalized such that $\mathrm{\Psi \Phi }=1$) is $\mathrm{\Psi }=[0.986,1.02],$from which we conclude using (54b(54b) $\begin{array}{rl}& =\frac{\Vert (I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }\Vert }{\mathrm{\Psi }U(I-{\rho }^{-1}T{)}^{-2}V\mathrm{\Phi }}\mathrm{\Psi }U(I-{\rho }^{-1}T{)}^{-1}.\end{array}$(54b) ) that the normalized vector of reproductive values is ${\mathrm{\Psi }}_{\mathrm{r}}=\left[\begin{array}{cccccc}0.713& 1.067& 0.603& 0.737& 1.745& 1.226\end{array}\right]$(see also panels (b) and (c) of Figure 5).

Figure 5. The plots of (a) $z\mapsto \mathrm{spectral}\mathrm{radius}\mathrm{of}\overline{K}(z)$, (b) the asymptotic i-state distribution and (c) the reproductive values for (the specific example of) the age-size structured-population model in Example 3.2.

Note that the condition of strict positivity of t's and f's can be relaxed somewhat. For example, if we assume that (i) individuals of size 1 only reproduce at age 3 (and then only produce offspring of size 1) and (ii) individuals of size 2 reproduce at age 2 (then producing size 1 individuals) and 3 (when they produce offspring of either size) then it is clear that $R(t)$ is strictly positive for $t\ge 3$.

In some situations, a further reduction of the NGM and the RE system is possible. Such is the case when we assume that (i) small individuals do not reproduce (i.e. $f_{12}=0=f_{13}$) and that (ii) large individuals produce size 1 and size 2 offspring at a given ratio $p:(1-p)$, regardless of the age of the parent. That is, $f_{15}=p{f}_5,f_{16}=p{f}_6$ and $f_{45}=(1-p)f_5,f_{46}=(1-p)f_6$ for some 0<p<1 and $f_5,f_6>0$, leading to $F=\left[\begin{array}{cccccc}0& 0& 0& 0& p{f}_5& p{f}_6\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& (1-p)f_5& (1-p)f_6\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right].$In this case, F (or U in (56(56) $V=\left[\begin{array}{cc}1& 0\\ 0& 0\\ 0& 0\\ 0& 1\\ 0& 0\\ 0& 0\end{array}\right]\text{and}U=\left[\begin{array}{cccccc}0& f_{12}& f_{13}& 0& f_{15}& f_{16}\\ 0& 0& 0& 0& f_{45}& f_{46}\end{array}\right].$(56) )) has rank 1. There are still two states-at-birth in the literal sense. However, if we interpret the term state-at-birth in the stochastic sense then all individuals are identical at birth in the sense that they are all born with i-state distribution $V=\left[\begin{array}{c}p\\ 0\\ 0\\ 1-p\\ 0\\ 0\end{array}\right].$Then F = V U for $U=\left[\begin{array}{cccccc}0& 0& 0& 0& f_5& f_6\end{array}\right]$and the NGM with the small domain amounts to the scalar ${\mathcal{L}}_{\mathrm{S}}=U(I-T{)}^{-1}V=p(f_5{t}_{51}+f_6(t_{62}{t}_{21}+t_{65}{t}_{51}))+(1-p)(f_5{t}_{54}+f_6{t}_{65}{t}_{54})={\mathcal{R}}_0.$In this special case, the basic reproduction number ${\mathcal{R}}_0$ can also easily be deduced from interpretation (and we invite the reader to do so).

We now have $\begin{array}{rl}k(1)& =0,\\ k(2)& =p{f}_5{t}_{51}+(1-p)f_5{t}_{54}>0\\ k(3)& =p{f}_6(t_{62}{t}_{21}+t_{65}{t}_{51})+(1-p)f_6{t}_{65}{t}_{54}\\ k(t)& =0\mathrm{for}t\ge 4\end{array}$and the scaler Euler-Lotka equation takes the form (57) $\overline{k}(z)=\sum _{s=1}^{\mathrm{\infty }}k(s)=k(2)z^{-2}+k(3)z^{-3}=1.$(57) If ${\mathcal{R}}_0>1$ we can then first find the growth rate ρ as the unique real solution of the Euler–Lotka equation (57(57) $\overline{k}(z)=\sum _{s=1}^{\mathrm{\infty }}k(s)=k(2)z^{-2}+k(3)z^{-3}=1.$(57) ) and then use (42(42) ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (I-{\rho }^{-1}T{)}^{-1}V\Vert }(I-{\rho }^{-1}T{)}^{-1}V$(42) ) and (47(47) ${\mathrm{\Psi }}_{\mathrm{r}}=\frac{\Vert (I-{\rho }^{-1}T{)}^{-1}V\Vert }{U(I-{\rho }^{-1}T{)}^{-2}V}U(I-{\rho }^{-1}T{)}^{-1}$(47) ) to determine the asymptotic i-state distribution and the asymptotic reproductive values, respectively.

Deducing the asymptotic behaviour in (linear discrete time) SPM from the corresponding RE

Starting with a structured-population model $X(t+1)=(F+T)X(t),X(t)\in {\mathbb{R}}^n,$the recipe is as follows.

Step 1.	Identify states-at-birth and write the fertility matrix in the form F = V U for some normalised $n\times m$ birth-state matrix V and the $m\times n$ fertility matrix U. If V has rank 1 proceed with the (a) part of the algorithm, otherwise follow the (b) part. Remark. For practical purposes, it may be handy to order the i-states such that the birth states come first. Then in the case of m states-at-birth (in the literal sense) V is the matrix with columns $e_1,\dots ,e_m$ and U is obtained by taking the rows $1,\dots ,m$ from F.
Step 2.	Derive $k(t)$ and $g(t)$ using (37a(37a) $\begin{array}{rl}k(t)& =U{T}^{t-1}V\end{array}$(37a) ) and check that the assumptions of Theorem 2.1 hold. Furthermore, check that the spectral radius of T is below one 1. Derive $K(t)$ and $G(t)$ using (31a(31a) $\begin{array}{rl}K(s)& =U{T}^{s-1}V,\end{array}$(31a) ) and check that the assumptions of Theorem 2.5 hold. Furthermore, check that the spectral radius of T is below one1.
Step 3.	Compute ${\mathcal{R}}_0$ in (39(39) ${\mathcal{R}}_0=U(I-T{)}^{-1}V.$(39) ). If ${\mathcal{R}}_0<1$ conclude that the population dies out. If ${\mathcal{R}}_0\ge 1$ compute the growth rate $\rho \ge 1$ as the unique real solution of the E-L Equation (4(4) $1=\sum _{s=1}^{\mathrm{\infty }}k(s)z^{-s}=:\overline{k}(z)$(4) ) with $\overline{k}$ in (38(38) $\overline{k}(z)=U(\mathrm{zI}-T{)}^{-1}V$(38) ) (with $\rho =1$ whenever ${\mathcal{R}}_0=1$). Conclude that: the population grows geometrically with rate ρ and the asymptotic i-state distribution is given by ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (\mathrm{\rho I}-T{)}^{-1}V\Vert }(\mathrm{\rho I}-T{)}^{-1}V,$ the (real time) reproductive values are collected in ${\mathrm{\Psi }}_{\mathrm{r}}=\frac{\Vert (\mathrm{\rho I}-T{)}^{-1}V\Vert }{U(\mathrm{\rho I}-T{)}^{-2}V}U(\mathrm{\rho I}-T{)}^{-1}.$ Compute ${\mathcal{R}}_0$ in (49(49) ${\mathcal{R}}_0=r_{\sigma }(U(I-T{)}^{-1}V).$(49) ). If ${\mathcal{R}}_0<1$ conclude that the population dies out. If ${\mathcal{R}}_0\ge 1$ compute the growth rate $\rho \ge 1$ as the unique real solution of the E-L Equation (11(11) $r_{\sigma }(\overline{K}(z))=1,$(11) ) with $\overline{K}$ in (48(48) $\overline{K}(z)=U(\mathrm{zI}-T{)}^{-1}V$(48) ) (with $\rho =1$ whenever ${\mathcal{R}}_0=1$). Determine the asymptotic state-at-birth distribution Φ and the asymptotic (generation-based) reproductive values Ψ as, respectively, the right and the left positive normalised eigenvectors of $\overline{K}(\rho )$ corresponding to eigenvalue 1. Conclude that: the population grows geometrically with rate ρ and the asymptotic i-state distribution is given by ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (\mathrm{\rho I}-T{)}^{-1}V\mathrm{\Phi }\Vert }(\mathrm{\rho I}-T{)}^{-1}V\mathrm{\Phi },$ the (real time) reproductive values are collected in ${\mathrm{\Psi }}_{\mathrm{r}}=\frac{\Vert (\mathrm{\rho I}-T{)}^{-1}V\mathrm{\Phi }\Vert }{\mathrm{\Psi }U(\mathrm{\rho I}-T{)}^{-2}V\mathrm{\Phi }}\mathrm{\Psi }U(\mathrm{\rho I}-T{)}^{-1}.$

¹ A sufficient condition is that all column sums of T are strictly below 1.

4. The case ${\mathcal{R}}_0=1=\rho$

We now show directly that in the special case when ${\mathcal{R}}_0=1=\rho$ we can compute

1. the asymptotic distribution (i.e. the right eigenvector) and

2. the reproductive values (i.e. the left eigenvector)

in the real time setting from those obtained in the generation-bookkeeping framework. Note that the former result has been demonstrated in a more general setting in [10, 12].

Suppose that there are m states-at-birth in the literal sense and as before write F = V U for some normalized $n\times m$ birth-state matrix V and the $m\times n$ fertility matrix U. Then $\mathcal{L}=U(I-T{)}^{-1}V$ is the next-generation matrix.

Let ${\mathrm{\Phi }}_g$ and ${\mathrm{\Phi }}_{\mathrm{r}}$ denote the normalized positive right eigenvector in, respectively, generation and real time framework, corresponding to eigenvalue 1. That is (58a) $U(I-T{)}^{-1}V{\mathrm{\Phi }}_g={\mathrm{\Phi }}_g,\Vert {\mathrm{\Phi }}_g\Vert =1$(58a) and (58b) $(\mathrm{VU}+T){\mathrm{\Phi }}_{\mathrm{r}}={\mathrm{\Phi }}_{\mathrm{r}},\Vert {\mathrm{\Phi }}_{\mathrm{r}}\Vert =1.$(58b) Now apply V from the left on both sides of (58a(58a) $U(I-T{)}^{-1}V{\mathrm{\Phi }}_g={\mathrm{\Phi }}_g,\Vert {\mathrm{\Phi }}_g\Vert =1$(58a) ) and write $\mathrm{VU}(I-T{)}^{-1}V{\mathrm{\Phi }}_g=(I-T)(I-T{)}^{-1}V{\mathrm{\Phi }}_g.$If we now denote $\mathrm{\Xi }=(I-T{)}^{-1}V{\mathrm{\Phi }}_g$ then clearly $(\mathrm{VU}+T)\mathrm{\Xi }=\mathrm{\Xi }$. We thus find that the asymptotic distributions in real time and in the generation setting are related by (59) ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (I-T{)}^{-1}V{\mathrm{\Phi }}_g\Vert }(I-T{)}^{-1}V{\mathrm{\Phi }}_g,$(59) as indeed follows from (51b(51b) $\begin{array}{rl}& =\frac{1}{\Vert (I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }\Vert }(I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }.\end{array}$(51b) ) when we observe that with ${\mathcal{R}}_0=1=\rho$, $\overline{K}(\rho )$ is the NGM and $\mathrm{\Phi }={\mathrm{\Phi }}_g$. When using the $l_1$ norm, the denominator in (59(59) ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (I-T{)}^{-1}V{\mathrm{\Phi }}_g\Vert }(I-T{)}^{-1}V{\mathrm{\Phi }}_g,$(59) ) is the expected duration of life of a newborn when we sample newborns according to ${\mathrm{\Phi }}_g$. As observed before, $V{\mathrm{\Phi }}_g$ is the asymptotic distribution in the generation-bookkeeping framework when we consider the NGM with large domain ${\mathcal{L}}_{\mathrm{L}}=F(I-T{)}^{-1}$.

Now let ${\mathrm{\Psi }}_g$ and ${\mathrm{\Psi }}_{\mathrm{r}}$ denote the normalized positive left eigenvector in, respectively, generation and real time framework, corresponding to eigenvalue 1, i.e. (60a) ${\mathrm{\Psi }}_{g}U(I-T{)}^{-1}V={\mathrm{\Psi }}_g,{\mathrm{\Psi }}_g{\mathrm{\Phi }}_g=1$(60a) and (60b) ${\mathrm{\Psi }}_{\mathrm{r}}(\mathrm{VU}+T)={\mathrm{\Psi }}_{\mathrm{r}},{\mathrm{\Psi }}_{\mathrm{r}}{\mathrm{\Phi }}_{\mathrm{r}}=1.$(60b) Acting with U from the right on both sides of (60a(60a) ${\mathrm{\Psi }}_{g}U(I-T{)}^{-1}V={\mathrm{\Psi }}_g,{\mathrm{\Psi }}_g{\mathrm{\Phi }}_g=1$(60a) ) we find that ${\mathrm{\Psi }}_{g}U(I-T{)}^{-1}\mathrm{VU}={\mathrm{\Psi }}_{g}U(I-T{)}^{-1}(I-T).$

Then for $\mathrm{\Xi }={\mathrm{\Psi }}_{g}U(I-T{)}^{-1}$ we have $\mathrm{\Xi }(\mathrm{VU}+T)=\mathrm{\Xi }$. We thus find that the reproductive values in the real time setting relate to the ones in the generation-bookkeeping by ${\mathrm{\Psi }}_{\mathrm{r}}=c{\mathrm{\Psi }}_{g}U(I-T{)}^{-1},$where the constant c is chosen such that with ${\mathrm{\Phi }}_{\mathrm{r}}$ in (59(59) ${\mathrm{\Phi }}_{\mathrm{r}}=\frac{1}{\Vert (I-T{)}^{-1}V{\mathrm{\Phi }}_g\Vert }(I-T{)}^{-1}V{\mathrm{\Phi }}_g,$(59) ) we have ${\mathrm{\Psi }}_{\mathrm{r}}{\mathrm{\Phi }}_{\mathrm{r}}=1$. We obtain (61) ${\mathrm{\Psi }}_{\mathrm{r}}=\frac{\Vert (I-T{)}^{-1}V{\mathrm{\Phi }}_g\Vert }{{\mathrm{\Psi }}_{g}U(I-T{)}^{-2}V{\mathrm{\Phi }}_g}{\mathrm{\Psi }}_{g}U(I-T{)}^{-1}.$(61) Alternatively, we can deduce this from (54b(54b) $\begin{array}{rl}& =\frac{\Vert (I-{\rho }^{-1}T{)}^{-1}V\mathrm{\Phi }\Vert }{\mathrm{\Psi }U(I-{\rho }^{-1}T{)}^{-2}V\mathrm{\Phi }}\mathrm{\Psi }U(I-{\rho }^{-1}T{)}^{-1}.\end{array}$(54b) ) by taking into account that with ${\mathcal{R}}_0=1=\rho$ we have $\mathrm{\Psi }={\mathrm{\Psi }}_g$. The j-th component of $U(I-T{)}^{-2}V{\mathrm{\Phi }}_g$ in the denominator of (61(61) ${\mathrm{\Psi }}_{\mathrm{r}}=\frac{\Vert (I-T{)}^{-1}V{\mathrm{\Phi }}_g\Vert }{{\mathrm{\Psi }}_{g}U(I-T{)}^{-2}V{\mathrm{\Phi }}_g}{\mathrm{\Psi }}_{g}U(I-T{)}^{-1}.$(61) ) is the expected age of the parent of a newborn individual with i-state j when we sample parents according to ${\mathrm{\Phi }}_g$. As observed before, ${\mathrm{\Psi }}_r$ is also the vector of generation-based reproductive values when we consider the NGM with large domain ${\mathcal{L}}_{\mathrm{L}}=F(I-T{)}^{-1}$.

5. Concluding remarks

The specification of a linear physiologically structured population model involves two rules, one for reproduction and one for development/maturation/movement and survival. Once these ingredients are specified, one can constructively define next-population-state operators. For constant environments, i.e. autonomous dynamics, a first question is: will the population ultimately grow exponentially (in which case one has to ponder the issue of density dependence) or decline and go extinct? As emphasized by Jim Cushing and Zhou Yicang [7], an efficient way of answering the question is to adopt a generation perspective, by focusing on expected lifetime offspring production by newborn individuals, and to characterize the appropriate average as the dominant eigenvalue of the next-generation matrix. On the other hand, people familiar with the Perron+-Frobenius theory of positive dynamical systems (as presented in, for instance, Berman and Plemmons [2], Bátkai et al. [1]) will be inclined to compute the spectral bound of the generator of the real time dynamics. Here, in the spirit of the elegant Li-Schneider paper [22], we have uncovered how these two approaches relate to each other. Our main contribution has been to highlight the connecting role of the Renewal Equation and to revive the powerful results of Feller and Thieme that describe the asymptotic large time behaviour of its solutions. As a one sentence summary of the present paper we offer: use the ingredients to formulate the RE, apply Feller/Thieme, and everything you might be interested in can next be computed explicitly.

Acknowledgments

Hans Metz benefited from the support from the ‘Chaire Modélisation Mathématique et Biodiversité of Veolia Environnement-École Polytechnique-Museum National d'Histoire Naturelle-Fondation X’.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Barbara Boldin acknowledges the support of the Slovenian Research Agency [I0-0035, research program P1-0285 and research projects J1-9186, J1-1715].

 

Appendix. Concise summary of the continuous time formalism

As the analogue of (20(20) $X(t)=(F+T)X(t-1).$(20) ) we consider the linear ODE system (A1) $\frac{\mathrm{d}X}{\mathrm{d}t}=(F+T)X,$(A1) where X takes values in ${\mathbb{R}}^n$, F is a positive $n\times n$ matrix describing reproduction and T is a positive-off-diagonal $n\times n$ matrix describing survival and development. We assume that the spectral bound $s(T)$ is negative, allowing us to write (A2) $-T^{-1}={\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{aT}}\mathrm{da}.$(A2) The variation-of-constants version of (A1(A1) $\frac{\mathrm{d}X}{\mathrm{d}t}=(F+T)X,$(A1) ) reads (A3) $X(t)={\mathrm{e}}^{\mathrm{tT}}X(0)+{\int }_0^t{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{FX}(t-\tau )\mathrm{d}\tau .$(A3) Putting (A4) $Y(t)=\mathrm{FX}(t)$(A4) and applying F to both sides of (A3(A3) $X(t)={\mathrm{e}}^{\mathrm{tT}}X(0)+{\int }_0^t{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{FX}(t-\tau )\mathrm{d}\tau .$(A3) ) we obtain the RE (A5) $Y(t)=F{\mathrm{e}}^{\mathrm{tT}}X(0)+{\int }_0^{t}F{\mathrm{e}}^{\mathrm{\tau T}}Y(t-\tau )\mathrm{d}\tau$(A5) and, once we have constructed Y by solving this RE, we can interpret (A3(A3) $X(t)={\mathrm{e}}^{\mathrm{tT}}X(0)+{\int }_0^t{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{FX}(t-\tau )\mathrm{d}\tau .$(A3) ) as an explicit expression for X: (A6) $X(t)={\mathrm{e}}^{\mathrm{tT}}X(0)+{\int }_0^t{\mathrm{e}}^{\mathrm{\tau T}}Y(t-\tau )\mathrm{d}\tau .$(A6) Whenever F has the form (21(21) $F=\mathrm{VU}.$(21) ), we can introduce $B(t)=\mathrm{UX}(t)$. Then $Y(t)=\mathrm{VB}(t)$ and we can rewrite (A5(A5) $Y(t)=F{\mathrm{e}}^{\mathrm{tT}}X(0)+{\int }_0^{t}F{\mathrm{e}}^{\mathrm{\tau T}}Y(t-\tau )\mathrm{d}\tau$(A5) ) as (A7) $B(t)=U{\mathrm{e}}^{\mathrm{tT}}X(0)+{\int }_0^{t}U{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{VB}(t-\tau )\mathrm{d}\tau .$(A7) The translation invariant version of (A7(A7) $B(t)=U{\mathrm{e}}^{\mathrm{tT}}X(0)+{\int }_0^{t}U{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{VB}(t-\tau )\mathrm{d}\tau .$(A7) ) is (A8) $B(t)={\int }_0^{\mathrm{\infty }}U{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{VB}(t-\tau )\mathrm{d}\tau .$(A8) Solutions of (A1(A1) $\frac{\mathrm{d}X}{\mathrm{d}t}=(F+T)X,$(A1) ) and (A8(A8) $B(t)={\int }_0^{\mathrm{\infty }}U{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{VB}(t-\tau )\mathrm{d}\tau .$(A8) ) defined for all times, i.e. for $t\in \mathbb{R}$, are related to each other by (A9) $\begin{array}{rl}B(t)& =\mathrm{UX}(t),\\ X(t)& ={\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{VB}(t-\tau )\mathrm{d}\tau .\end{array}$(A9) Clearly (A1(A1) $\frac{\mathrm{d}X}{\mathrm{d}t}=(F+T)X,$(A1) ) has a solution of the form (A10) $X(t)={\mathrm{e}}^{\mathrm{rt}}{\mathrm{\Phi }}_{\mathrm{r}}$(A10) if and only if r is an eigenvalue of F + T and ${\mathrm{\Phi }}_{\mathrm{r}}$ is a corresponding eigenvector, i.e. (A11) $(F+T){\mathrm{\Phi }}_{\mathrm{r}}=r{\mathrm{\Phi }}_{\mathrm{r}}.$(A11) By substitution we see that (A8(A8) $B(t)={\int }_0^{\mathrm{\infty }}U{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{VB}(t-\tau )\mathrm{d}\tau .$(A8) ) has a solution of the form (A12) $B(t)={\mathrm{e}}^{\mathrm{rt}}\mathrm{\Phi }$(A12) if and only if $U(\mathrm{rI}-T{)}^{-1}V$ has eigenvalue 1 and Φ is a corresponding eigenvector, i.e. (A13) $\mathrm{\Phi }=U(\mathrm{rI}-T{)}^{-1}V\mathrm{\Phi }.$(A13) According to (A9(A9) $\begin{array}{rl}B(t)& =\mathrm{UX}(t),\\ X(t)& ={\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{VB}(t-\tau )\mathrm{d}\tau .\end{array}$(A9) ) we have (A14) $\begin{array}{rl}\mathrm{\Phi }& =U{\mathrm{\Phi }}_r,\\ {\mathrm{\Phi }}_r& =(\mathrm{rI}-T{)}^{-1}V\mathrm{\Phi }.\end{array}$(A14) We recall how this can be used in practice:

Step 1.	Find the real number r for which the spectral radius of $U(\mathrm{rI}-T{)}^{-1}V$ equals one, by exploiting that the spectral radius is a monotone decreasing function of r (see [19–21]).
Step 2.	Find Φ as the dominant positive eigenvector of $U(\mathrm{rI}-T{)}^{-1}V$.
Step 3.	Compute ${\mathrm{\Phi }}_r$ by using (A14(A14) $\begin{array}{rl}\mathrm{\Phi }& =U{\mathrm{\Phi }}_r,\\ {\mathrm{\Phi }}_r& =(\mathrm{rI}-T{)}^{-1}V\mathrm{\Phi }.\end{array}$(A14) ); if desired, renormalize to obtain the multiple of ${\mathrm{\Phi }}_r$ that has norm one.

For the corresponding left eigenvectors we have (A15) $\begin{array}{rl}{\mathrm{\Psi }}_r(F+T)& =r{\mathrm{\Psi }}_r,\\ \mathrm{\Psi }& =\mathrm{\Psi F}(\mathrm{rI}-T{)}^{-1}.\end{array}$(A15) Rewriting the first of these as follows: ${\mathrm{\Psi }}_{r}F={\mathrm{\Psi }}_r(\mathrm{rI}-T{)}^{-1}⟺{\mathrm{\Psi }}_{r}F(\mathrm{rI}-T{)}^{-1}={\mathrm{\Psi }}_r$we see that, modulo normalization, ${\mathrm{\Psi }}_r$ and Ψ are identical. Note that the NGM with large domain is now given by (A16) ${\mathcal{L}}_{\mathrm{L}}:=-F{T}^{-1}$(A16) and that, with ${\mathcal{R}}_0$ defined as the spectral radius of $\mathcal{L}$, the relation (A17) $\mathrm{sign}({\mathcal{R}}_0-1)=\mathrm{sign}r$(A17) holds. In the special case ${\mathcal{R}}_0=1,r=0$, Ψ and Φ are, respectively, the left and the right eigenvector of $\mathcal{L}$ corresponding to its dominant eigenvalue.

The formalism easily extends to an i-state space that is a continuum and to measures see [10, 19, 20, 26]. A convenient starting point is the following generalization of (A9(A9) $\begin{array}{rl}B(t)& =\mathrm{UX}(t),\\ X(t)& ={\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{\tau T}}\mathrm{VB}(t-\tau )\mathrm{d}\tau .\end{array}$(A9) ): (A18) $\begin{array}{rl}b(t,\omega )& ={\int }_{\mathrm{\Omega }}\beta (\eta ,\omega )x(t,\mathrm{d}\eta ),\\ x(t,\omega )& ={\int }_0^{\mathrm{\infty }}{\int }_{\mathrm{\Omega }}u(a,\xi ,\omega )b(t-a,\mathrm{d}\xi )\mathrm{da}.\end{array}$(A18) By substituting the second of these identities into the first we obtain the RE (A19) $b(t,\omega )={\int }_0^{\mathrm{\infty }}{\int }_{\mathrm{\Omega }}K(a,\xi ,\omega )b(t-a,\mathrm{d}\xi )\mathrm{da}$(A19) with (A20) $K(a,\xi ,\omega ):={\int }_{\mathrm{\Omega }}\beta (\eta ,\omega )u(a,\xi ,\mathrm{d}\eta ).$(A20) Clearly, (A19(A19) $b(t,\omega )={\int }_0^{\mathrm{\infty }}{\int }_{\mathrm{\Omega }}K(a,\xi ,\omega )b(t-a,\mathrm{d}\xi )\mathrm{da}$(A19) ) has a solution of the form (A21) $b(t,\omega )={\mathrm{e}}^{\mathrm{rt}}\mathrm{\Phi }(\omega )$(A21) if and only if (A22) $\mathrm{\Phi }(\omega )={\int }_{\mathrm{\Omega }}\left({\int }_0^{\mathrm{\infty }}K(a,\xi ,\omega ){\mathrm{e}}^{-\mathrm{ra}}\mathrm{da}\right)\mathrm{\Phi }(\mathrm{d}\xi ).$(A22) We refer to [19] for a detailed analysis.