# A Framework for Studying Transients in Marine Metapopulations††thanks: This
work was supported by the Natural Science and Engineering Research Council of
Canada (NSERC) Discovery Grants (MAL and PvdD).

Peter D. Harrington, Mark A. Lewis22footnotemark: 2 , and P. van den Driessche
Department of Mathematical and Statistical Sciences, University of Alberta,
Edmonton, Canada.
<EMAIL_ADDRESS>mark.lewis@ualberta.ca)Department of Biological
Sciences, University of Alberta, Edmonton, Canada.Department of Mathematics
and Statistics, University of Victoria, Victoria, Canada<EMAIL_ADDRESS>
###### Abstract

Transient dynamics can often differ drastically from the asymptotic dynamics
of systems. In this paper we provide a unifying framework for analysing
transient dynamics in marine metapopulations, from the choice of norms to the
addition of stage structure. We use the $\ell_{1}$ norm, because of its
biological interpretation, to extend the transient metrics of reactivity and
attenuation to marine metapopulations, and use examples to compare these
metrics under the more commonly used $\ell_{2}$ norm. We then connect the
reactivity and attenuation of marine metapopulations to the source-sink
distribution of habitat patches and demonstrate how to meaningfully measure
reactivity when metapopulations are stage-structured.

Keywords: transient dynamics, metapopulation, reactivity, source-sink
dynamics, marine systems

AMS subject classifications: 92D25, 92D40, 34A30, 34C11

##  1 Introduction

Transient dynamics, those that occur over short timescales, can often be
vastly different from the asymptotic or long term dynamics of ecological
systems. However, throughout the history of mathematical biology much of the
work has focused on determining the asymptotic dynamics of biological systems.
While the study of long-term dynamics has given ecologists many tools to
analyze the behaviour of populations, these tools are often not the same as
those required to understand transient dynamics. Recently Hastings et al. [16]
have shown that transient dynamics are much more ubiquitous than previously
assumed and long transients occur in many different ecological systems, from
plankton and coral to voles and grouse. Studying the transient dynamics of an
ecological system can give useful insight into the different processes that
may occur after a disturbance, change in environmental conditions, or change
in human intervention to a system. In some marine systems that are driven by
environmental fluctuations, such as the Dungeness crab, transient dynamics may
in fact be key to understanding how these systems behave [18].

There has also been a recent push to characterize the different types of
systems which display long transient dynamics that differ significantly from
their asymptotic dynamics [32, 16, 15, 14]. Hastings et al. [16] have loosely
categorized four different drivers of long transient dynamics in ecological
systems: ghost attractors and crawl-bys, slow-fast dynamics, high
dimensionality, and stochastic noise. These categories are not always distinct
and certain systems may indeed fall into multiple categories. For example, a
predator prey system may have a crawl-by past a saddle node which drives the
transient dynamics in this system, but this could also be thought of as a
difference in timescales of the predator decline due to lack of prey. For
metapopulations the main driver of transient dynamics is often the high
dimensionality arising due to spatial structure, though often these transient
dynamics are exacerbated by the other drivers as well.

Some of the earliest studies of systems that could generate long transients
were systems with spatial structure [17, 27]. It seems intuitive that spatial
structure or spatial heterogeneity can drive some sort of transient dynamics
in a system. If individuals start in one location in a habitat, especially a
poor habitat, then it will take time before they can spread over the entire
habitat and the long-term population dynamics begin to emerge. What is
surprising is that spatial structure can also give rise to so called long-
lived transients, where the transient dynamics are extensive enough that they
continue on timescales past which we typically measure biological populations
[17].

One method of adding spatial structure to a population is to formulate it as a
metapopulation, where distinct populations live on habitat patches which are
connected via dispersal or migration. Metapopulation models were originally
proposed by Levins [26] to model patch occupancy in habitats consisting of
isolated habitat patches, but these early models used space implicitly rather
than explicitly. Later metapopulation models have included space explicitly by
allowing for differing habitat quality on patches or differing dispersal
between patches [12, 11], though often these models are focused on the
proportion of occupied patches rather than the population size on each patch.
However, many marine metapopulation models as well as epidemiological
metapopulation models explicitly track the number of individuals on each patch
as well as movement or dispersal between patches [27, 2, 9, 3]. In this paper
we model the metapopulation structure following this spatially explicit
framework where individuals are tracked rather than the proportion of occupied
patches.

Another benefit of the metapopulation framework is that habitat patches can be
classified into source patches and sink patches. This classification can occur
in many different ways [9, 39, 25], but commonly a source is a productive
habitat patch and a sink is a poor habitat patch. Early measures of sources
and sinks were mainly focused on connectivity between patches, however more
recently it has been understood that it is the interplay between patch
connectivity and local patch productivity that characterizes patches as
sources or sinks. One of the new and easily tractable metrics that embodies
this relationship comes from the theory of next-generation matrices and the
basic reproduction number, $R_{0}$. This framework, originally developed in
epidemiology, has been used to characterize sources and sinks in populations
of mussels, salmon, and sea lice on salmon farms [25, 22, 13].

While metapopulation theory has previously been used to classify patches as
sources and sinks, other metrics have been used to characterize the transient
dynamics of systems. Reactivity was initially introduced by Neubert and
Caswell [35] to measure the maximum initial growth rate of a system over all
possible perturbations from an equilibrium. If the maximum initial growth rate
is positive, then the system is reactive. Complimenting reactivity is the
amplification envelope, which is the maximum possible amplification at time
$t$ that can be achieved by a perturbation. Later, Townley and Hodgson [49]
introduced attenuation as the opposite metric to measure initial decline of
populations; a system attenuates if the minimum possible growth rate declines
following a perturbation. Reactivity and attenuation are then most interesting
when they are different from the stability of the equilibrium of a system —
when a system attenuates but is unstable, or is reactive but stable — and it
is in these situations that we focus this paper.

It should be noted that reactivity, attenuation and the amplification envelope
are all defined from the linearization of a non-linear system about an
equilibrium. These measures are therefore most useful around hyperbolic
equilibria, where the dynamics of the non-linear system can be well
approximated by the dynamics of the linear system. If an equilibrium is not
hyperbolic then the trajectories in the non-linear system may no longer be
similar to the linearization by which reactivity, attenuation, and the
amplification envelope are defined. Even around a hyperbolic equilibrium the
trajectories of the non-linear and linearsized systems may diverge as they
move away from the equilibrium. Here we use the technique of linearization to
determine reactivity and attenuation as others have before us, but want to
emphasize these caveats as they are often brushed over in the transient
literature.

In this paper we apply these transient measures of growth to a class of
biological metapopulation models where there is no migration between
population patches, only birth on new patches. These are a subset of birth-
jump processes [19] and include models for marine meroplanktonic species,
where larvae can travel through the ocean between population patches but
adults remain confined to a habitat patch. Specific species that exhibit this
structure include sea lice [1], corals and coral reef fish [6, 24], barnacles
[40], Dungeness crabs [4], sea urchins [4], and many other benthic marine
species [7]. This type of system also encompasses many plant species where
seeds are carried between suitable habitat patches [23], and depending on the
census timing could also include insect species where there is one large
dispersal event between habitat patches, such as the spruce budworm [28, 33,
52] and mountain pine beetle [41]. Lastly this class of models also includes
multi-patch or multi-city epidemiological metapopulation models where
infections can spread between patches, for example infected residents of a
city may travel and infect residents of other cities before returning home
[2].

The paper is structured as follows. First we extend the general theory of
transients to marine and other birth-jump metapopulations. We then demonstrate
that even for linear two-patch metapopulations it is possible to have
transient dynamics that are vastly different than the asymptotic dynamics of
these populations and moreover that the time for which these transient
dynamics occur can extend for an arbitrarily long amount of time. Next, we
examine how the structure of the metapopulation can enhance the transient
dynamics of metapopulations. Then we demonstrate how to connect the transient
behaviour of the metapopulation to its source-sink distribution. Finally, we
discuss how to measure the transient dynamics of stage-structured
metapopulations in a useful and biological meaningful manner.

##  2 Extending the general theory of transients to metapopulations

In this section we apply the metrics of reactivity [35] and attenuation [49]
to general systems of single-species metapopulations. We demonstrate that
often reactivity and attenuation in metapopulations have simple closed forms,
especially if we are analysing how population size responds to initial
perturbations. In order to present our work in a general form, we model the
dynamics of a metapopulation of a single species on $n$ patches around the
zero equilibrium with the system:

$x^{\prime}=Ax,$ (1)

where $A=[a_{ij}]$ is a real irreducible matrix of order $n$. This most often
represents the linearization of a non-linear system, which more completely
captures the dynamics of the population but could also represent the full
dynamics of a linear system if density dependence was not important to the
population dynamics.

For the analyses in this paper we focus on biologically realistic single-
species metapopulations where the entries of $x(t)$ are non-negative when
beginning with a non-negative initial condition, $x(0)$. This condition is
equivalent to requiring that $A$ be an essentially non-negative (Metzler
matrix), such that all the off-diagonal entries of $A$ are non-negative (Thm
2.4, [47]). Biologically this means that the presence of individuals on one
patch cannot contribute to the decline of a population on another patch and
that the population on each patch will not become negative.

###  2.1 Reactivity and attenuation using the $\ell_{1}$ norm

To analyze the transient dynamics of this metapopulation let us introduce some
notation from Neubert and Caswell [35]. We say that a population is reactive
if there is an initial condition such that the initial growth rate of the
total population is positive. Using notation from Lutscher and Wang [29],
reactivity is formally defined as:

$\bar{\sigma}=\max_{||x_{0}||\neq
0}\left[\frac{1}{||x||}\frac{d||x||}{dt}\biggr{\rvert}_{t=0}\right].$ (2)

If $\bar{\sigma}>0$ then a system is reactive, and if $\bar{\sigma}\leq 0$
then the system is not reactive. Neubert and Caswell [35] use the $\ell_{2}$
norm to measure the population size, as this allows for the simplification
that $\bar{\sigma}$ is the max eigenvalue of $(A+A^{T})/2$. However, the
$\ell_{2}$ norm lacks a reasonable biological interpretation, and so others
have instead used the $\ell_{1}$ norm to define reactivity [22, 48, 44].
Biologically, the $\ell_{1}$ norm,

$||x||_{1}=\sum_{i=1}^{n}|x_{i}|$

can be interpreted as the total population on all patches of a metapopulation.
Moreover in the metapopulation framework using the $\ell_{1}$ norm is also
convenient to determine reactivity from the population matrix $A$.

Similar to reactivity, we also say that a population attenuates if there is an
initial condition for which the initial growth rate of the total population
declines [49]. This is formally defined as

$\underline{\sigma}=\min_{||x_{0}||\neq
0}\left[\frac{1}{||x||}\frac{d||x||}{dt}\biggr{\rvert}_{t=0}\right].$ (3)

If $\underline{\sigma}<0$ then the system attenuates, and if
$\underline{\sigma}\geq 0$ then the system does not attenuate. Comparing the
definitions of attenuation and reactivity we can see that it is possible for a
system to be both reactive and to attenuate, if there are certain initial
conditions for which $\bar{\sigma}>0$ is achieved and others such that
$\underline{\sigma}<0$. In fact, reactivity and attenuation are most
interesting when they are different from the stability of the system, as this
is when the transient dynamics are different than the long term dynamics. This
is when a system is reactive but stable, so that the total population
initially grows but eventually declines, or when a system attenuates but is
unstable, so that the total population declines but eventually grows. It
should also be noted that the only systems which are not reactive and do not
attenuate are those in which the total population size remains constant for
all time.

The last measures that we define here to use in some later sections are the
amplification envelope and the maximum amplification. The amplification
envelope is the maximum possible amplification at time $t$ from a perturbation
$x_{0}$ and is defined mathematically as:

$\rho(t)=\max_{||x_{0}||\neq 0}\frac{||x(t)||}{||x_{0}||}.$ (4)

The maximum amplification is simply the maximum of the amplification envelope
over all time:

$\rho_{\max}=\max_{t\geq 0}\rho(t)=\max_{\begin{subarray}{c}t\geq 0\\\
||x_{0}||\neq 0\end{subarray}}\frac{||x(t)||}{||x_{0}||}.$ (5)

While reactivity and attenuation quantify the short time response to a
perturbation, the amplification envelope and maximum amplification quantify
how large a perturbation can become and how long growth can last. It is for
these purposes that we use the amplification envelope and maximum
amplification in subsection 3.1 and section 4.

Now before quantifying the reactivity and attenuation of the entire
metapopulation, let us first determine the initial growth rate of the
population if we begin with one individual on patch $j$. We call this initial
growth rate $\lambda_{j}$, and mathematically we define

$\lambda_{j}=\sum_{i=1}^{n}x_{i}^{\prime}(0),$

with $x(0)=e_{j}$, where $e_{j}$ is the vector of length $n$ with $1$ in the
$j$th entry and $0$s elsewhere. In terms of system (1) this simplifies to the
$j$th column sum of $A$,

$\lambda_{j}=\sum_{i=1}^{n}a_{ij}.$

The initial growth rate for a given patch $j$, $\lambda_{j}$, can also be
calculated from the lifecycle graph as the sum of all the outgoing birth rates
from a patch minus the death rate on that patch, where any paths describing
movement of individuals between patches are ignored. We can then connect this
patch specific initial growth rate with the total growth rate, or reactivity,
using the following lemma.

######  Lemma 2.1.

Assuming that population sizes are non-negative on each patch, so that
$x_{j}\geq 0$ for all $j$, then for system (1) under the $\ell_{1}$ norm,

$\bar{\sigma}=\max_{1\leq j\leq n}\lambda_{j}=\max_{1\leq j\leq
n}\sum_{i}a_{ij}.$

###### Proof.

Since $x_{j}\geq 0$ for all $j$, the absolute value signs in (2) can be
dropped and so

$\displaystyle\bar{\sigma}$
$\displaystyle=\max_{||x_{0}||=1}\left[\frac{d||x||}{dt}\biggr{\rvert}_{t=0}\right]$
$\displaystyle=\max_{||x_{0}||=1}\left[\frac{d}{dt}\sum_{i=1}^{n}x_{i}\biggr{\rvert}_{t=0}\right]$
$\displaystyle=\max_{||x_{0}||=1}\left[\sum_{i=1}^{n}\frac{d}{dt}x_{i}\biggr{\rvert}_{t=0}\right]$
$\displaystyle=\max_{||x_{0}||=1}\left[1^{T}x^{\prime}\biggr{\rvert}_{t=0}\right].$

Substituting $x^{\prime}=Ax$ from system (1) gives

$\displaystyle\bar{\sigma}$
$\displaystyle=\max_{||x_{0}||=1}\left[1^{T}Ax_{0}\right]$
$\displaystyle=\max_{||x_{0}||=1}\left[\sum_{j=1}^{n}\left(\sum_{i=1}^{n}a_{ij}\right)x_{0j}\right]$

Now let $k$, with $1\leq k\leq n$ maximize $\sum_{i=1}^{n}a_{ij}$, so that
$\sum_{i=1}^{n}a_{ik}=\max_{1\leq j\leq n}\sum_{i=1}^{n}a_{ij}$. Then, under
$||x_{0}||=1$,

$\displaystyle\sum_{j=1}^{n}\left(\sum_{i=1}^{n}a_{ij}\right)x_{0j}$
$\displaystyle\leq\left(\sum_{i=1}^{n}a_{ik}\right)\sum_{j=1}^{n}x_{0j}$
$\displaystyle=\left(\sum_{i=1}^{n}a_{ik}\right),$

because $||x_{0}||=1$, with equality when $x_{0j}=\begin{cases}1&j=k\\\
0&j\neq k\end{cases}$. Therefore

$\displaystyle\bar{\sigma}$ $\displaystyle=\max_{1\leq j\leq
n}\left(\sum_{i=1}^{n}a_{ij}\right)$ $\displaystyle=\max_{1\leq j\leq
n}\lambda_{j}.$

∎

With a similar proof we can connect the patch specific initial growth rate to
attenuation via the following lemma:

######  Lemma 2.2.

Assuming that population sizes are non-negative on each patch, so that
$x_{j}\geq 0$ for all $j$, then for system (1) under the $\ell_{1}$ norm,

$\underline{\sigma}=\min_{1\leq j\leq n}\lambda_{j}=\min_{1\leq j\leq
n}\sum_{i}a_{ij}.$

###  2.2 Comparing reactivity in the $\ell_{1}$ and $\ell_{2}$ norms

In this section we present some examples of systems that are reactive in
$\ell_{1}$ but not in $\ell_{2}$ and vice versa to illuminate the difference
between measuring reactivity in the two norms. It has previously been noted
that reactivity depends on the norm in which it is measured [35, 29], and the
following examples help clarify the underlying biological and mathematical
meaning of the two norms. The $\ell_{1}$ norm,

$||x||_{1}=\sum_{i=1}^{n}|x_{i}|,$

represents the total population size of the metapopulation, whereas the
$\ell_{2}$ norm,

$||x||_{2}=\sqrt{\sum_{i=1}^{n}x_{i}^{2}},$

represents a sort of Euclidean distance of the metapopulation away from the
origin.

Figure 1: The phase plane for system (1) with $A=\begin{bmatrix}-1&0\\\
1&0\end{bmatrix},$ which is reactive in $\ell_{2}$ but not in $\ell_{1}$. The
line $x_{1}+x_{2}=1$ and the circle $x_{1}^{2}+x_{2}^{2}=1$ geometrically
depict $||x||=1$ in the $\ell_{1}$ and $\ell_{2}$ norms respectively. The
derivative vectors for the phase plane are shown in red and two different
initial trajectories are shown in green and blue. The green trajectory is an
example that is reactive in $\ell_{2}$, but not in $\ell_{1}$, and the blue
trajectory is another example that is not reactive in $\ell_{1}$.

### Example 1

First, we present an example that is reactive in $\ell_{2}$ but not in
$\ell_{1}$. Take system (1) with

$A=\begin{bmatrix}-1&0\\\ 1&0\end{bmatrix}.$

This system simply redistributes individuals from patch 1 to patch 2, and the
phase plane is shown in Figure 1. It is not reactive in the $\ell_{1}$ norm
because the total population size is not increasing, but it is reactive in
$\ell_{2}$. This highlights how measuring reactivity in the $\ell_{2}$ norm
can at times defy our biological expectation of what reactivity should mean —
the growth of a population — and reinforces our rationale for using the
$\ell_{1}$ norm to measure reactivity in metapopulations. While the matrix $A$
is reducible and this system is only semi-stable, and thus may be considered a
borderline example, if $a_{22}$ is replaced by a small negative number,
$-\epsilon$, and $a_{12}$ is replaced by a small positive number,
$\epsilon/2$, then for sufficiently small $\epsilon$, $A$ will be irreducible
and the system will now be stable, but will still be reactive in $\ell_{2}$
and not in $\ell_{1}$.

### Example 2

The second example, which is reactive in $\ell_{1}$ but not in $\ell_{2}$ is
system (1) with

$A=\begin{bmatrix}-1&3/2\\\ 1/3&-1\end{bmatrix},$

where the phase plane is shown in Figure 2. Now the system is reactive in
$\ell_{1}$ because if we start with one individual on the second patch (the
dynamics governed by the second row of $A$) the total population grows, but in
such a way that it will not be reactive in $\ell_{2}$. This example
demonstrates that again reactivity in $\ell_{2}$ can defy our biological
expectation of reactivity, but now in the opposite way. Here the total
population grows, yet the system is not reactive in $\ell_{2}$. Note that this
system is equivalent to system (8) with $\epsilon=3$.

Together, the two examples highlight the differences that can occur when
measuring reactivity in different norms and the caution that should be taken
when interpreting reactivity in the $\ell_{2}$ norm biologically. Here we only
present examples which are reactive in $\ell_{2}$ but not in $\ell_{1}$ and
vice versa but it is also possible to find examples of systems which attenuate
in $\ell_{2}$ but not in $\ell_{1}$.

Figure 2: The phase plane for system (1) with $A=\begin{bmatrix}-1&3/2\\\
1/3&-1\end{bmatrix},$ which is reactive in $\ell_{1}$ but not in $\ell_{2}$.
The line $x_{1}+x_{2}=1$ and the circle $x_{1}^{2}+x_{2}^{2}=1$ geometrically
depict $||x||=1$ in the $\ell_{1}$ and $\ell_{2}$ norms respectively. The
derivative vectors for the phase plane are shown in red and two different
initial trajectories are shown in green and blue. The green trajectory is an
example that is reactive in $\ell_{1}$, but not in $\ell_{2}$, and the blue
trajectory is another example that is not reactive in $\ell_{2}$.

###  2.3 The relationship between stability and reactivity/attenuation

Now that we have presented a couple of examples that demonstrate the
difference between reactivity in the $\ell_{1}$ and $\ell_{2}$ norms, we show
that in both norms if a system is stable it attenuates, and if a system is
unstable it is reactive. While we hypothesize that this is true in all norms,
these are necessary properties for a norm to be useful in measuring reactivity
or attenuation and thus we present the proofs in $\ell_{1}$ and $\ell_{2}$ for
completeness.

We begin by presenting the proof for the $\ell_{1}$ norm.

######  Lemma 2.3.

If the $x=0$ equilibrium for

$x^{\prime}=Ax$

is stable and $A$ is a Metzler matrix, then $\min\lambda_{j}<0$, where
$\lambda_{j}=\sum_{i=1}^{n}a_{ij}$, and thus the system attenuates in the
$\ell_{1}$ norm. Likewise if the $x=0$ equilibrium is unstable then
$\max\lambda_{j}>0$ and so the system is reactive in the $\ell_{1}$ norm.

###### Proof.

Let $A$ be a Metzler matrix, then

$\mu(A)=\max\\{\Re(\lambda):\lambda\in\sigma(A)\\},$

and $\mu(A)$ is an eigenvalue of $A$ with a non-negative normalized
eigenvector $v$, with $\sum_{j=1}^{n}v_{j}=1$, such that $Av=\mu(A)v.$ (Thm
A.43, [46]).

Therefore if the initial population is arranged according to the eigenvector
$v$, then the initial growth rate is

$\displaystyle\sum_{i=1}^{n}\frac{dx_{i}}{dt}$
$\displaystyle=1^{T}x^{\prime}|_{x_{0}=v}$ $\displaystyle=1^{T}Av$
$\displaystyle=1^{T}\mu(A)v.$

If the $x=0$ equilibrium for $x^{\prime}=Ax$ is stable, then $\mu(A)<0$, and
therefore

$\sum_{i=1}^{n}\frac{dx_{i}}{dt}=1^{T}\mu(A)v<0.$

Let $\lambda_{k}=\min_{j}\lambda_{j}$, then

$\lambda_{k}\leq\sum_{j=1}^{n}\lambda_{j}v_{j}=1^{T}Av=1^{T}\mu(A)v<0,$

as $v$ is normalized. Thus if the $x=0$ equilibrium for $x^{\prime}=Ax$ is
stable, then $\min\lambda_{j}<0$, and so the system attenuates in the
$\ell_{1}$ norm, and the minimum growth rate occurs when the initial
population is all on a single patch.

Similarly if the $x=0$ equilibrium for $x^{\prime}=Ax$ is unstable, then
$\mu(A)>0$, and therefore

$\sum_{i=1}^{n}\frac{dx_{i}}{dt}=1^{T}\mu(A)v>0.$

Let $\lambda_{k}=\max_{j}\lambda_{j}$, then

$\lambda_{k}\geq\sum_{j=1}^{n}\lambda_{j}v_{j}=1^{T}Av=1^{T}\mu(A)v>0.$

Thus if the $x=0$ equilibrium for $x^{\prime}=Ax$ is unstable, then it is also
reactive in the $\ell_{1}$ norm and the maximum growth rate occurs when the
initial population is all on a single patch.

∎

In the $\ell_{2}$ norm the following lemma gives the corresponding result:

######  Lemma 2.4.

If the $x=0$ equilibrium for $x^{\prime}=Ax$ is stable then

$\min_{||x_{0}||_{2}\neq
0}\left(\frac{1}{||x||_{2}}\frac{d||x||_{2}}{dt}\right)\biggr{\rvert}_{t=0}<0$

and thus the system attenuates in $\ell_{2}$. Similarly if the $x=0$
equilibrium for $x^{\prime}=Ax$ is unstable then

$\max_{||x_{0}||_{2}\neq
0}\left(\frac{1}{||x||_{2}}\frac{d||x||_{2}}{dt}\right)\biggr{\rvert}_{t=0}>0$

and so the system is reactive in $\ell_{2}$.

###### Proof.

In the $\ell_{2}$ norm

$\left(\frac{1}{||x||_{2}}\frac{d||x||_{2}}{dt}\right)\biggr{\rvert}_{t=0}=\frac{x_{0}^{T}H(A)x_{0}}{x_{0}^{T}x_{0}},$
(6)

where $H(A)=(A+A^{T})/2$ is the Hermitian part of $A$ [35]. Furthermore the
maximum of the right hand side of equation (6) over all $x_{0}$ is equal to
$\lambda_{\max}(H(A))$, where $\lambda_{\max}(H(A))$ is the maximum eigenvalue
of $H(A)$ (Thm 4.2.2, [20]). Similarly the minimum of the right hand side of
equation (6) over all $x_{0}$ is equal to $\lambda_{\min}(H(A))$, the minimum
eigenvalue of $H(A)$. Since $A$ is real, $H(A)$ is real and symmetric and so
$\lambda_{\max}(H(A))$ and $\lambda_{\min}(H(A))$ are real.

Therefore $x^{\prime}=Ax$ is reactive in the $\ell_{2}$ norm if
$\lambda_{\max}(H(A))>0$ and attenuates if $\lambda_{\min}(H(A))<0$. To relate
stability to reactivity and attenuation, we use the fact that the real part of
the dominant eigenvalue of $A$ (denoted by $\mu(A)$) is between the maximum
and minimum eigenvalues of the Hermitian part of $A$ (Fact 5, pg 14-2, [50]),

$\lambda_{\min}(H(A))\leq\mu(A)\leq\lambda_{\max}(H(A)).$

If the $x=0$ equilibrium for $x^{\prime}=Ax$ is unstable, then

$\lambda_{\max}(H(A))\geq\mu(A)\geq 0$

and so the system is reactive, and if the $x=0$ equilibrium for
$x^{\prime}=Ax$ is stable, then

$\lambda_{\min}(H(A))\leq\mu(A)\leq 0$

and so the system attenuates.

∎

In this section we have shown how to calculate reactivity and attenuation
using the $\ell_{1}$ norm in metapopulations, proven that if the equilibrium
of a system is unstable/stable then the system must be reactive/attenuate in
both the $\ell_{1}$ and $\ell_{2}$ norms, and demonstrated the difference
between reactivity in the $\ell_{1}$ and $\ell_{2}$ norms using a couple
salient examples. We now return to the motivating feature of this paper —
systems that are reactive and stable or attenuate and are unstable — and in
the following section we provide examples of long lived transients in these
systems.

##  3 Metapopulations with arbitrarily large transient growth or decay

Here we examine two different metapopulations, one of which is reactive and
can exhibit arbitrarily large transient growth, and the other that attenuates
and can decline to arbitrarily small levels. In each case this transient
growth differs from the system’s long term growth trajectory: the
metapopulation that exhibits large growth eventually declines, and the system
that declines eventually grows. Both of these example metapopulations are
linear systems, and therefore the addition of non-linearities to construct
more realistic models could further exacerbate the length of the transient
period. These examples are not meant to imply that there are realistic
biological metapopulations that can grow arbitrarily large before decaying,
but rather to emphasize that the difference between transient dynamics and
asymptotic dynamics can be quite stark even in linear systems.

###  3.1 Arbitrarily large transient growth

First we present a reactive metapopulation that can exhibit arbitrarily large
transient growth, but eventually declines. In this metapopulation individuals
can either give birth to new individuals on the same patch, or give birth to
individuals on the other patch, but there is no migration of individuals
between patches. As mentioned in the introduction, this type of model is
applicable to many marine metapopulations where adults are sedentary but
larvae can disperse, to plant populations where seeds can be carried between
habitat patches, or other populations governed by birth-jump processes. Let
the metapopulation be described by:

$\displaystyle x^{\prime}$ $\displaystyle=(b-d)x+b_{12}y$ (7) $\displaystyle
y^{\prime}$ $\displaystyle=b_{21}x+(b-d)y$

so that $b$ is the on patch birth rate and $d$ is the on patch death rate,
$b_{12}$ is the birth rate of individuals on patch 2 producing new individuals
on patch 1, and $b_{21}$ is the birth rate of individuals on patch 1 producing
new individuals on patch 2. The system is linear, so assuming that
$(b-d)^{2}\neq b_{12}b_{21}$ the only steady state is $x=y=0$.

For the metapopulation to eventually decline, both eigenvalues need to be
negative. For system (7) the eigenvalues are

$\displaystyle\lambda_{1}$ $\displaystyle=(b-d)+\sqrt{b_{12}b_{21}}$
$\displaystyle\lambda_{2}$ $\displaystyle=(b-d)-\sqrt{b_{12}b_{21}}$

and thus we require that $(b-d)<0$ and $(b-d)^{2}>b_{12}b_{21}$. Now in order
for the metapopulation to be reactive we need either $b_{12}>-(b-d)$ or
$b_{21}>-(b-d)$. Here we choose $b_{21}>-(b-d)$, so that if we start with one
individual on patch 1, i.e. $x(0)=1$, $y(0)=0$, the metapopulation initially
grows.

To prove that the metapopulation can grow arbitrarily large, we show that the
limit as some parameter approaches 0 of $\max_{t}(x(t)+y(t))$ is unbounded. In
terms of the maximum amplification defined by (5), this is equivalent to
showing that the limit of the maximum amplification in the $\ell_{1}$ norm,
$\rho_{\max}$, becomes unbounded. This is because the initial condition chosen
above is such that $x(t)+y(t)=\rho(t)$ in the $\ell_{1}$ norm . To take the
limit, we must first reduce the parameters in our system until we are left
with a single parameter which we can let approach 0, while still maintaining
the inequalities above that govern the stability and reactivity of the system.
Let $b-d=-1$, $b_{12}=\epsilon/2$, and $b_{21}=1/\epsilon$, where $\epsilon$
is a small positive parameter that approaches 0. Our reduced system can now be
written as:

$\displaystyle x^{\prime}$ $\displaystyle=-x+\frac{{\epsilon}}{2}y$ (8)
$\displaystyle y^{\prime}$ $\displaystyle=\frac{1}{{{\epsilon}}}x-y$
$\displaystyle x(0)$ $\displaystyle=1$ $\displaystyle y(0)$ $\displaystyle=0.$

This system is stable and the digraph for this system is shown in Figure 3.
This system is reactive in $\ell_{1}$ and $\ell_{2}$ for small $\epsilon$ and
the solution is:

$\displaystyle x(t)$
$\displaystyle=\frac{1}{2}\left(e^{-(1-\frac{1}{\sqrt{2}})t}+e^{-(1+\frac{1}{\sqrt{2}})t}\right)$
$\displaystyle y(t)$
$\displaystyle=\frac{1}{\sqrt{2}\epsilon}\left(e^{-(1-\frac{1}{\sqrt{2}})t}-e^{-(1+\frac{1}{\sqrt{2}})t}\right).$

$x$$y$$\epsilon^{-1}$$\epsilon/2$$-1$$-1$ Figure 3: Digraph for system (8).
The directed edges represent the birth rate of individuals on the outgoing
patch producing new individuals on the incoming patch. The self loops are the
birth rate minus the death rate on a patch.

We want to show that

$\lim_{\epsilon\rightarrow
0}\max_{t}\left(x(t)+y(t)\right)=\lim_{\epsilon\rightarrow
0}\rho_{\max}=\infty.$

Normally to calculate the maximum we would take the derivative of
$(x(t)+y(t))$, set it equal to 0, solve for $t$, and then evaluate
$(x(t)+y(t))$ at this value of $t$. However it turns out this is rather
complicated, so we will simplify this process by first noting that $x(t)>0$
for all $t$. Therefore

$\max_{t}\left(x(t)+y(t)\right)>\max_{t}y(t).$

Now we only have to perform the above process on $y(t)$, rather than
$(x(t)+y(t))$.

Setting $y^{\prime}(t)=0$ and solving for $t$, we find that the time that the
maximum of $y(t)$ is achieved, $t_{\rm max}$, along with the corresponding
maximum in $y$, $y(t_{\rm max})$, are:

$\displaystyle t_{\rm max}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\log(1+\sqrt{2})-\log(-1+\sqrt{2})\right)$
$\displaystyle y(t_{\rm max})$
$\displaystyle=\frac{\sqrt{2}}{\epsilon}(1+\sqrt{2})^{(-\frac{1}{2}-\frac{1}{\sqrt{2}})}(-1+\sqrt{2})^{(\frac{1}{\sqrt{2}}-\frac{1}{2})}.$

We can clearly see that $\lim_{\epsilon\rightarrow 0}y(t_{\rm max})=\infty$
and thus also $\lim_{\epsilon\rightarrow
0}\max_{t}x(t)+y(t)=\allowbreak\lim_{\epsilon\rightarrow
0}\rho_{\max}\allowbreak=\infty$.

Therefore even in a two patch metapopulation that is asymptotically stable,
there is always a parameter combination for which the total population, and
thus also the maximum amplification in the $\ell_{1}$ norm, $\rho_{\max}$, can
initially grow arbitrarily large before they decay. Again, this is not meant
to imply that there are realistic biological metapopulations that can grow
arbitrarily large before decaying, but to emphasize how different the
transient and asymptotic dynamics of a system can be.

###  3.2 Transient decay to arbitrarily small levels

Similar to the previous section, we now present an example of a metapopulation
that attenuates and can decay to an arbitrarily small population size before
eventually growing. This metapopulation is unstable mathematically because it
eventually grows, but biologically the metapopulation could first go extinct
if the total population size decays below one individual before it eventually
increases. We again use a metapopulation where individuals can either give
birth to new individuals on their patch or on the other patch, but cannot
migrate between patches. The difference between this metapopulation and the
example used in the previous section, is that now the on patch birth and death
rates differ between patches, but the between patch birth rates are the same.
Let the metapopulation be described by:

$\displaystyle x^{\prime}$ $\displaystyle=(b_{1}-d_{1})x+\epsilon y$
$\displaystyle y^{\prime}$ $\displaystyle=\epsilon x+(b_{2}-d_{2})y$

where $b_{1}$ and $d_{1}$ are the birth and death rates on patch 1, $b_{2}$
and $d_{2}$ are the birth and death rates on patch 2, and $\epsilon$ is the
interpatch birth rates for both patches.

$x$$y$$\epsilon$$\epsilon$$-1$$\epsilon$ Figure 4: Digraph for system (9). The
directed edges represent the birth rate of individuals on the outgoing patch
producing new individuals on the incoming patch. The self loops are the birth
rate minus the death rate on a patch.

In order for the metapopulation to eventually grow, we assume that the birth
rate is greater than the death rate on one of the patches. We choose this to
be patch 2, thus we require $b_{2}-d_{2}>0$. We also want our population to
initially decline when starting on patch 1, for this to occur we assume
$b_{1}-d_{1}+\epsilon<0$. To prove that the metapopulation can decay to an
arbitrarily small population size we reduce the system to have a single
parameter and then show that the limit as the parameter approaches 0 of
$\min_{t}(x(t)+y(t))=0$. Let $b_{1}-d_{1}=-1$ and $b_{2}-d_{2}=\epsilon$, then
our system can be written in terms of a single positive parameter, $\epsilon$,
as:

$\displaystyle x^{\prime}$ $\displaystyle=-x+\epsilon y$ (9) $\displaystyle
y^{\prime}$ $\displaystyle=\epsilon x+\epsilon y$ $\displaystyle x(0)$
$\displaystyle=1$ $\displaystyle y(0)$ $\displaystyle=0.$

This system is unstable and the corresponding digraph is shown in Figure 4. It
attenuates in both the $\ell_{1}$ and $\ell_{2}$ norms for small $\epsilon$.

It is possible to show that the minimum population size can grow arbitrarily
small in a manner similar to the previous section, though the calculations are
somewhat more complicated. Instead in this section, we perform an asymptotic
expansion in terms of $\epsilon$ to demonstrate the limiting behaviour of
system (9). Let $x(t)=x_{0}(t)+\epsilon x_{1}(t)+O(\epsilon^{2})$ and
$y(t)=y_{0}(t)+\epsilon y_{1}(t)+O(\epsilon^{2})$. Then the zero order system
is:

$\displaystyle x_{0}^{\prime}(t)$ $\displaystyle=-x_{0}(t)$ $\displaystyle
y_{0}^{\prime}(t)$ $\displaystyle=0$ $\displaystyle x_{0}(0)$
$\displaystyle=1$ $\displaystyle y_{0}(0)$ $\displaystyle=0,$

that has the solution $x_{0}(t)=e^{-t}$ and $y_{0}=0$. We can proceed in a
similar manner to solve the first order terms, and then our solution up to
order $\epsilon$ is given by:

$\displaystyle x(t)$ $\displaystyle=e^{-t}+O(\epsilon^{2})$ $\displaystyle
y(t)$ $\displaystyle=\epsilon(1-e^{-t})+O(\epsilon^{2}).$

This solution is valid for small $t$, and is therefore our inner
approximation. To find our outer approximation for large $t$, we rescale
$t=\tau/\epsilon$ and arrive at the system:

$\displaystyle\epsilon X^{\prime}$ $\displaystyle=-X+\epsilon Y$
$\displaystyle\epsilon Y^{\prime}$ $\displaystyle=\epsilon X+\epsilon Y.$

We can again solve the zero order and first order equations and arrive at the
following solution with two undetermined coefficients:

$\displaystyle X(\tau)$ $\displaystyle=\epsilon Ce^{\tau}+O(\epsilon^{2})$
$\displaystyle Y(\tau)$ $\displaystyle=Ce^{\tau}+\epsilon(C\tau
e^{\tau}+(C+K)e^{\tau})+O(\epsilon^{2}).$

To solve our undetermined coefficients we require that
$\lim_{t\rightarrow\infty}x(t)=\lim_{\tau\rightarrow 0}X(\tau)$, and
$\lim_{t\rightarrow\infty}y(t)=\lim_{\tau\rightarrow 0}Y(\tau)$. From
$x(\infty)=X(0^{+})$, we find $C=0$. Substituting $C=0$ into
$y(\infty)=Y(0^{+})$ to solve for $K$ we find $K=1$. Adding our inner and
outer solutions together and subtracting the overlap ($x(\infty)=X(0^{+})=0$
and $y(\infty)=Y(0^{+})=\epsilon$) we find

$\displaystyle x(t)$ $\displaystyle=e^{-t}+O(\epsilon^{2})$ $\displaystyle
y(t)$ $\displaystyle=\epsilon(e^{\epsilon t}-e^{-t})+O(\epsilon^{2}),$

thus our total population size behaves as

$x(t)+y(t)=e^{-t}+\epsilon(e^{\epsilon t}-e^{-t})+O(\epsilon^{2}).$ (10)

We can see from equation (10) and Figure 5 that for very small $\epsilon$, the
total population size behaves similarly to $e^{-t}$ before eventually growing.
Thus for an minimum population threshold, we can always find an $\epsilon$
small enough, such that the solution crosses the threshold before the
population grows. We can prove this by solving the full system and taking the
limit of the minimum, though we leave this out as the calculations become
rather messy.

Figure 5: Asymptotic approximation of the total population size (red) compared
to the true total population size (blue) for system (9), with $\epsilon=0.01$.
The asymptotic approximation is given by equation (10).

Here we have shown that there are metapopulations for which the transient
population can grow arbitrarily large or small, no matter the asymptotic
stability of the system. In the next section we demonstrate how increasing the
patch number can increase transient growth in advective metapopulations.

##  4 Increasing patch number increases transient timescale

In this section we examine the role that strong advection, coupled with a
large number of patches, can have on the transient dynamics of
metapopulations. In aquatic systems, habitat patches may be quite productive,
but strong advection can sweep most larvae to the next patch, leading to large
transient growth on downstream patches before the population eventually
disappears from the last patch. This phenomenon can occur in metapopulations
situated in rivers, ocean channels, or reef systems where reefs are arranged
along a coastline with a directional current.

$x_{1}$$x_{2}$$b_{2}$$\epsilon$$b-d$$b-d$$\dots$$x_{j}$$\dots$$x_{n-1}$$x_{n}$$b_{2}$$\epsilon$$b-d$$b-d$$b-d$$b_{2}$$\epsilon$$b_{2}$$\epsilon$$b_{2}$$\epsilon$$b_{2}$$\epsilon$
Figure 6: Digraph for system (11). The directed edges represent the birth rate
of individuals on the outgoing patch producing new individuals on the incoming
patch. The self loops are the birth rate minus the death rate on a patch.

Consider a metapopulation on $n$ patches where the dynamics are described by
the following system of equations:

$\displaystyle x^{\prime}$ $\displaystyle=Ax\quad$ (11) $\displaystyle x(0)$
$\displaystyle=e_{1}$ $\displaystyle A$
$\displaystyle=\begin{bmatrix}b-d&\epsilon&0&\dots&0\\\
b_{2}&b-d&\epsilon&\ddots&\vdots\\\ 0&b_{2}&b-d&\ddots&\\\
\vdots&\ddots&\ddots&\ddots&\epsilon\\\ 0&\dots&0&b_{2}&b-d\end{bmatrix},$

where $b$ is the on patch birth rate for each patch, $d$ is the death rate on
each patch, $b_{2}$ is the birth rate of patch $j-1$ on patch $j$, $\epsilon$
is the birth rate from patch $j+1$ to patch $j$, and $e_{1}$ is a vector with
1 in the first entry and 0s elsewhere. The digraph for this system is shown in
Figure 6.

The instantaneous measures of growth, $\lambda_{j}$, and the reactivity,
$\bar{\sigma}$ in the $\ell_{1}$ norm, are therefore

$\displaystyle\lambda_{1}$ $\displaystyle=b-d+b_{2}$
$\displaystyle\lambda_{j}$ $\displaystyle=b-d+b_{2}+\epsilon\quad
j=2,\dots,n-1$ $\displaystyle\lambda_{n}$ $\displaystyle=b-d+\epsilon$
$\displaystyle\bar{\sigma}$ $\displaystyle=b-d+b_{2}+\epsilon.$

Let $b-d+\epsilon<0$ and $b-d+b_{2}>0$, then the system is reactive
($\bar{\sigma}>0$), and this maximum initial growth rate is achieved if the
initial individual is on any patch except for patch 1 or $n$, though if the
individual starts on patch 1 the initial growth rate is still positive. In
system (11) $A$ is a tridiagonal Toeplitz matrix, so it has eigenvalues [38]

$\displaystyle\lambda_{h}$
$\displaystyle=b-d+2\sqrt{b_{2}\epsilon}\cos\left(\frac{h\pi}{n+1}\right)$
$\displaystyle h=1,\dots,n,$

and corresponding right eigenvectors, $v_{h}$, where the $k$th entry is given
by

$\displaystyle v_{h,k}$
$\displaystyle=(b_{2}/\epsilon)^{k/2}\sin\left(\frac{hk\pi}{n+1}\right)$
$\displaystyle k=1,\dots,n;h=1,\dots,n.$

The solution to system (11) can therefore be written as

$x(t)=We^{Jt}W^{-1}e_{1},$

where $W$ is a matrix containing the eigenvectors, $v_{h}$, and $J$ is a
diagonal matrix with the eigenvalues, $\lambda_{h}$, on the diagonal. For all
but very small $t$, this is equivalent to the amplification envelope in the
$\ell_{1}$ norm, $\rho(t)$, defined by equation (4). Through examination of
the eigenvalues, this system is stable if $\epsilon$ is small enough such that

$b-d+2\sqrt{b_{2}\epsilon}<0.$

Parameters that satisfy the inequalities that determine stability and
reactivity in the $\ell_{1}$ norm can be found in the caption of Figure 7. In
this case the maximum total population size, and also maximum amplification,
are given by

$x_{\max}=\rho_{\max}=\max_{t\geq 0}1^{T}We^{Jt}W^{-1}e_{1},$

with the corresponding time $t_{\max}$, which is the value of $t$ for which
the maximum occurs. The last measure of transience that is useful in this
system is the total transient time, $t_{\rm total}$, which we define as the
time it takes for the population size to decline below one, after initially
starting with one individual, or

$t_{\rm total}=\min\\{t>0:1^{T}We^{Jt}W^{-1}e_{1}\leq 1\\}.$

So how does the number of patches affect the magnitude and length of
transients? In Figure 7, which compares a 5 and 15 patch system, we can see
that increasing the patch number increases both the magnitude of growth and
the duration.

(a)  (b)  (c)

Figure 7: The population sizes on each patch for the advective system (11),
starting with one individual on patch 1. In a) the population sizes are shown
on a log scale for a metapopulation of 5 patches, and in b) and c) the
population sizes are shown for a metapopulation of 5 and 15 patches
respectively. On the untransformed scale only the population size on the last
patch can be seen as it is far larger than on any of the other patches,
whereas on the log scale the population sizes of all patches can be seen.
Parameters for this simulation are $b=0.09655$, $\epsilon=0.000001$,
$b_{2}=2$, and $d=0.1$, chosen so that system (11) is reactive but stable.

Here it can be difficult to see the duration of transience exhibited by all
patches on a regular scale, but on the log scale we can see that all patches
except for patch 1 experience a large period of transient growth, before they
decay below 1 individual (dashed line). Patch 1 does not experience a large
period of growth because the internal growth rate, $b-d$, is negative and the
birth rate from patch 2 to patch 1, $\epsilon$, is too small to overcome this
negative internal growth rate.

What cannot be seen from Figure 7 is the dependence of transient growth on
system parameters. We find that decreasing $b_{2}$ in system (11) results in a
large decrease in the maximum population size (and maximum amplification),
$x_{\rm max}$ $(\rho_{max}$), and the total transient time, $t_{\rm total}$,
but only a small decrease in the time at which the maximum population size is
achieved, $t_{\max}$. Decreasing $b$ however, results in a large decrease in
$x_{\rm max}$ ($\rho_{\max}$), $t_{\rm max}$, and $t_{\rm total}$.

The relationship between increased transient time and number of patches can
also be found for a linear metapopulation where all patches have negative
initial growth rates, $\lambda_{j}$, except for the last patch which has a
positive initial growth rate. In this case the total population size decays
for a long time before it eventually grows, and the time that it decays
depends on the number of patches.

We can see then that for a linear metapopulation, the length of the linear
array can accentuate the transient growth that is possible in the system and
that this is especially true for advective systems where there is some sort of
directed birth in one direction in the array. Systems with this type of
advective flow include marine metapopulations located in channels near the
mouth of rivers, or long coral reefs that are captured inside of a dominant
coastal current flow.

Having presented some illuminating metapopulation examples that demonstrate
the magnitude that transient dynamics can differ from asymptotic dynamics, we
now turn back to the general theory of transients in metapopulations and
connect it to the source-sink distribution of habitat patches.

##  5 Connecting the source-sink dynamics to the transient dynamics

In this section we demonstrate how to connect the transient measures of
initial population growth to the source-sink distribution of the
metapopulation. For the transient measure of the patch specific contribution
to the initial growth of the total population we use $\lambda_{j}$, previously
defined in section 2. To measure the source-sink distribution we use the the
classification derived from the next-generation matrix, $K$. This measure has
previously been used to calculate the source-sink distribution in heterogenous
environments [25, 31, 21, 13]. In order to calculate the next-generation
matrix for system (1) we decompose $A=F-V$, where $F$ is a non-negative matrix
with positive entries that describe the birth of new individuals in the
metapopulation, and $V$ is a non-singular $M$ matrix with entries that
describe the transfer of individuals between compartments or in this case
habitat patches, and also includes the death of individuals. Because $V$ is a
non-singular $M$ matrix, $V^{-1}$ is non-negative. The next-generation matrix,
$K=[k_{ij}]$, can then be calculated as $K=FV^{-1}$. This next generation
matrix is then commonly used to calculate the basic reproduction number
$R_{0}$, which is the average number of new individuals produced by one
initial individual, as $R_{0}=\rho(K)$, where $\rho$ denotes the spectral
radius. However we can also define $R_{j}$ as the number of new individuals
produced on all patches from one initial individual starting on patch $j$,
that can then be calculated as:

$R_{j}=\sum_{i=1}^{n}k_{ij}.$

Using this measure we classify patch $j$ as a source if $R_{j}>1$, as then one
individual on patch $j$ would produce more than one individual in the total
metapopulation. Likewise we classify patch $j$ as a sink if $R_{j}<1$.

###  5.1 Expressing $R_{0}$ as a weighted sum of $R_{j}$

Before examining the connection between the source-sink distribution $R_{j}$
and the initial growth $\lambda_{j}$, we first highlight a connection betwen
$R_{j}$ and $R_{0}$. It turns out, as shown in the following Lemma, that
$R_{0}$ can be calculated as a weighted sum of each $R_{j}$, and surprisingly
this relationship between the spectral radius and the column sums of a matrix
does not require any further assumptions on the matrix structure, though if
the matrix is not non-negative, the components of the right eigenvector need
not be real. Here $1$ is the (column) vector with each entry equal to 1, and
$e_{j}$ is the vector with the only non-zero entry being 1 in the $j$th row.

######  Lemma 5.1.

Let $v=[v_{i}]$ be the right eigenvector associated with the dominant
eigenvalue of the next-generation matrix, $\mathcal{R}_{0}=\rho(K)$,
normalized so $\sum_{1\leq i\leq n}v_{i}=1$. Then the basic reproduction
number $\mathcal{R}_{0}=\sum_{1\leq j\leq n}R_{j}v_{j}$, where
$R_{j}=1^{T}Ke_{j}=\sum_{i=1}^{n}k_{ij}.$

###### Proof.

First, we can rewrite $\mathcal{R}_{0}$ as

$\displaystyle\mathcal{R}_{0}$ $\displaystyle=\mathcal{R}_{0}1^{T}v$
$\displaystyle=1^{T}\mathcal{R}_{0}v,$

because the eigenvector has been normalized to sum to 1. Then, as
$\mathcal{R}_{0}$ is the eigenvalue of $K$ associated with $v$ and the column
sums of $K$ are $R_{j}$,

$\displaystyle\mathcal{R}_{0}$ $\displaystyle=1^{T}\mathcal{R}_{0}v$
$\displaystyle=1^{T}Kv$
$\displaystyle=\begin{bmatrix}R_{1}&R_{2}&\dots&R_{n}\end{bmatrix}v$
$\displaystyle=\sum_{j=1}^{n}R_{j}v_{j}$

∎

The entries $v_{j}$ of the right eigenvector can be interpreted as the
probability that a new individual begins on patch $j$ [8]. Therefore
$\mathcal{R}_{0}$ can be interpreted as the sum over all patches, of the
probability that an individual is born on patch $j$, multiplied by the number
of new individuals it will produce on all other patches over its lifetime.

Similarly, if we define $\lambda_{0}$ to be the dominant eigenvalue of $A$,
with the associated normalized eigenvector $u$, then

$\displaystyle\lambda_{0}$ $\displaystyle=1^{T}Au$
$\displaystyle=\sum_{j=1}^{n}\sum_{i=1}^{n}a_{ij}u_{j}$
$\displaystyle=\sum_{j=1}^{n}\lambda_{j}u_{j},$

where it should be noted that $\lambda_{j}$ is the $j$th column sum of $A$,
rather than an eigenvalue of $A$.

###  5.2 Connecting the source-sink distribution, $R_{j}$, to the initial
growth rate, $\lambda_{j}$

Now that we have decomposed the dominant eigenvalues, $R_{0}$ and
$\lambda_{0}$, into weighted sums of the columns of $K$ and $A$ respectively,
we proceed to connect the source-sink distribution of a particular patch,
$R_{j}$, to the initial growth from an individual on that patch,
$\lambda_{j}$. To do so there are some restrictions that we need to impose on
our metapopulation system and this is where we limit our study to marine or
birth-jump metapopulation models where juveniles or seeds can disperse between
patches while adults remain confined to habitat patches. The mathematical
restriction defined by this class of models comes from the decomposition of
$A$ into $F-V$. Here $V$ contains all entries that describe the transfer of
individuals between compartments or patches. For the results presented in this
paper, we require that $V$ has the following reducible form:

$V=\begin{bmatrix}V_{11}&0\\\ V_{21}&D\end{bmatrix},$

where $V_{11}$ is $k\times k$, $D=\operatorname{diag}(d_{k+1},\dots,d_{n})$
with $d_{k+1},\dots,d_{n}$ all positive, $0\leq k\leq n-1$, and $V$ is a non-
singular $M$ matrix.

With this structure, individuals on patches $j=k+1,\dots,n$ cannot migrate
between patches, but can still give birth to new individuals on any patch. If
$V$ is completely diagonal, then there is no migration between any patches,
only birth on other patches. This is the case for models of plants with seed
dispersal, or simplified marine metapopulation models if the juvenile stage is
not explicitly modelled. Under this structure, we first present a proofs
connecting our instantaneous and generational growth measures, $\lambda_{j}$
and $R_{j}$, before presenting a two patch example.

######  Lemma 5.2.

Let $A=F-V$ for system (1), where $F$ is non-negative, and $V$ is a non-
singular $M$ matrix with the following form:

$V=\begin{bmatrix}V_{11}&0\\\ V_{21}&D\end{bmatrix},$

where $V_{11}$ is $k\times k$, $D=\operatorname{diag}(d_{k+1},\dots,d_{n})$
with $d_{k+1},\dots,d_{n}$ all positive, and $0\leq k\leq n-1$. For $k+1\leq
j\leq n$, $\lambda_{j}$ is positive if and only if $R_{j}>1$.

###### Proof.

First, we can write $\lambda_{j}$ as

$\displaystyle\lambda_{j}$ $\displaystyle=\sum_{i=1}^{n}a_{ij}$
$\displaystyle=1^{T}Ae_{j}.$

Then decomposing $A$ into $F-V$, and inserting $V^{-1}V$

$\displaystyle\lambda_{j}$ $\displaystyle=1^{T}(F-V)e_{j}$
$\displaystyle=1^{T}(F-V)V^{-1}Ve_{j}$ $\displaystyle=1^{T}(FV^{-1}-I)Ve_{j}.$

For $k+1\leq j\leq n$, $V$ is diagonal, so

$Ve_{j}=d_{j}e_{j}.$

Therefore

$\displaystyle\lambda_{j}$ $\displaystyle=1^{T}(FV^{-1}-I)d_{j}e_{j}$
$\displaystyle=(R_{j}-1)d_{j}.$

Now $d_{j}>0$, and thus $\lambda_{j}>0$ if and only if $R_{j}>1$. ∎

######  Corollary 5.3.

In the notation of Lemma 5.2, if $V$ is diagonal, then $\lambda_{j}>0$ if and
only if $R_{j}>1$ for $j=1,\dots,n$.

######  Corollary 5.4.

Under the same conditions as Lemma 5.2, in the $\ell_{1}$ norm
$\bar{\sigma}>0$ if
$\max_{k+1\leq j\leq n}R_{j}>1$, and $\underline{\sigma}<0$ if $\min_{k+1\leq
i\leq n}R_{i}<1$.

###### Proof.

Under the conditions in Lemma 5.2, we know that $R_{j}-1$ has the same sign as
$\lambda_{j}$ for $k+1\leq j\leq n$. Therefore if $\max_{k+1\leq j\leq
n}R_{j}>1$ then $\bar{\sigma}=\max_{1\leq j\leq n}\lambda_{j}>0$, i.e. the
system is reactive. Similarly if $\min_{k+1\leq j\leq n}R_{j}<1$ then
$\underline{\sigma}=\min_{1\leq j\leq n}\lambda_{j}<0$, i.e. the population
attenuates. ∎

######  Corollary 5.5.

Under the same conditions of Lemma 5.2, only with $V$ diagonal, then in the
$\ell_{1}$ norm $\bar{\sigma}>0$ if and only if $\max_{1\leq j\leq n}R_{j}>1$
and $\underline{\sigma}<0$ if and only if $\min_{1\leq j\leq n}R_{j}$.

###### Proof.

From Corollary 5.3, $\lambda_{j}>0$ if and only if $R_{j}>1$ for each patch
$j$. Therefore if $\bar{\sigma}=\max_{1\leq j\leq n}\lambda_{j}>0$, then
$\max_{1\leq j\leq n}R_{j}>1$, and likewise if $\max_{1\leq j\leq n}R_{j}>1$,
then $\bar{\sigma}>0$. The same argument holds for $\min_{1\leq j\leq
n}\lambda_{j}$ and $\min_{1\leq j\leq n}R_{j}$. ∎

Now that we have presented theory connecting the initial growth rate,
$\lambda_{j}$, and the source-sink distribution, $R_{j}$, we present an
example to illustrate how to calculate these measures and how Lemma 5.2 and
the following corollaries can be used to connect them.

### Example 3

Here we present an example of a metapopulation consisting of two habitat
patches, patch 1 and patch 2. New individuals can be born on either patch, but
no individuals can migrate between patches. This system represents a
simplification of the adult dynamics of many marine meroplanktonic
metapopulations, where dispersal between patches occurs at the larval stage,
rather than the sedentary adult stage. This system could also represent plant
metapopulations that spread through seed dispersal, if the habitat landscape
is patchy. The metapopulation dynamics can be represented with the following
set of ODEs:

$\displaystyle x_{1}^{\prime}$
$\displaystyle=b_{11}x_{1}+b_{12}x_{2}-d_{1}x_{1}$ (12) $\displaystyle
x_{2}^{\prime}$ $\displaystyle=b_{21}x_{1}+b_{22}x_{2}-d_{2}x_{2},$

where $b_{ij}$ is the birth rate for births from patch $j$ to patch $i$, and
$d_{i}$ is the death rate on patch $i$. The lifecycle graph for this system is
shown in Figure 8.

$x_{1}$$x_{2}$$b_{21}$$b_{12}$$b_{11}-d_{1}$$b_{22}-d_{2}$ Figure 8: Digraph
for system (12). The directed edges represent the birth rate of individuals on
the outgoing patch producing new individuals on the incoming patch. The self
loops are the birth rate minus the death rate on a patch.

We then decompose $A=F-V$ and construct the next generation matrix,
$K=FV^{-1}$:

$\displaystyle A$ $\displaystyle=\begin{bmatrix}b_{11}-d_{1}&b_{12}\\\
b_{21}&b_{22}-d_{2}\end{bmatrix},\quad F=\begin{bmatrix}b_{11}&b_{12}\\\
b_{21}&b_{22}\end{bmatrix},\quad V=\begin{bmatrix}d_{1}&0\\\
0&d_{2}\end{bmatrix},$ $\displaystyle K$
$\displaystyle=FV^{-1}=\begin{bmatrix}b_{11}/d_{1}&b_{12}/d_{2}\\\
b_{21}/d_{1}&b_{22}/d_{2}\end{bmatrix}.$

For an initial individual starting on patch 1, the expected lifetime is
$1/d_{1}$, and the rate that the individual is producing new individuals on
both patches is $b_{11}+b_{21}$. Therefore $R_{1}=(b_{11}+b_{21})/d_{1}$ is
the total number of individuals born onto both patch 1 and patch 2 over one
generation. It is clear that $\lambda_{1}=b_{11}+b_{21}-d_{1}>0$ if and only
if $R_{1}=\frac{b_{11}}{d_{1}}+\frac{b_{21}}{d_{1}}>1$, in accordance with
Corollary 5.3. Similarly $\lambda_{2}=b_{12}+b_{22}-d_{2}>0$ if and only if
$R_{2}=\frac{b_{12}}{d_{2}}+\frac{b_{22}}{d_{2}}>1$. The system is therefore
reactive if $\max(R_{1},R_{2})>1$ (Corollary 5.5).

At first glance it seems obvious that if an individual starts on a source
patch the population should have a positive initial growth rate or if the
population starts on a sink patch it should have a negative initial growth
rate, and we have shown from Corollary 5.3 and Example 3 that this is indeed
the case for marine metapopulations. What is perhaps surprising is that this
is not the case for general metapopulations when adults can migrate between
habitat patches, and thus when the conditions of Lemma 5.2 and Corollary 5.3
are not met. In the general case it is possible to start with an individual on
a source patch, but for the population to initially decline and likewise to
start on a sink patch but for the population to initially grow. An example of
such a metapopulation is shown in Appendix A.

Here in this section we have shown that for marine metapopulations and other
metapopulations where the population dynamics are defined by birth jump
processes, there is a one-to-one relationship between the source-sink
distribution of patches and the initial growth rate when starting with one
adult on a patch. That is, the initial population growth rate is positive if
we start with one adult on patch $j$ if and only if patch $j$ is a source, and
the initial growth rate is negative if and only if patch $j$ is a sink. This
is a useful relationship biologically as there are several marine
metapopulations where patches have already been classified into sources and
sinks, and thus the transient dynamics for these systems can now be better
understood.

##  6 Stage structure

In this section we add stage structure to demonstrate some of the nuances in
analysing transients in stage structured metapopulation models. The main issue
with analysing reactivity and attenuation in models with stage structure is
due to the fact that adults often give birth to many more juveniles than will
survive to become adults, and that juveniles cannot normally give birth to new
juveniles. This presents a few complications.

The first complication is the fact that if we want to analyse the initial
growth or decay of a population, starting with an individual in a patch, it
now depends if the individual is a juvenile or an adult. If we start with a
juvenile, then there is no way that the total population, or even the patch
population, can grow, given that the juvenile has to first survive to the
adult stage to give birth to new juveniles. Thus we want to start with one
adult on a patch.

However, if we start with an adult in a patch, and it gives birth to new
juveniles, how do we count these new juveniles? If we are considering a marine
metapopulation do we count every larvae as a new individual? If so, every
marine metapopulation would exhibit transience, as each adult often produces
thousands of larvae. This then begs the question: in a stage structured
metapopulation, can we scale the juvenile population so that transient
measures of population growth, such as reactivity and attenuation, are useful
for stage structured models and measure the biologically relevant quantities?

To motivate the necessity of an honest scaling we highlight a discrete time
example of transients in Dungeness crabs from Caswell and Neubert [5].
Dungeness crabs give birth to an enormous number of larvae, many of which do
not survive to settle and become juveniles after 1 year. In this case the
discrete time model requires a census time to measure new crabs after 1 year.
If the census is taken pre breeding, then the system exhibits little
reactivity, as many of the larvae that where initially born have not survived
to become 1 year old juveniles. However, it the census is taken post breeding,
then all of the eggs or larvae are counted and the initial amplification is
increased by $10^{5}$. The models considered in this paper are continuous time
and do not face this exact problem, but it is easy to see that the addition of
a larvae stage in a marine stage structured metapopulation has large effects
on the reactivity of the system.

Returning then to our stage structured model with only juvenile and adult
stages, how should the juvenile stage be scaled so that an initial growth in
juveniles also corresponds in some sense to growth in the total population?
Ideally, we would scale the juvenile population so that each juvenile is
scaled by the probability that it will become an adult. If we scale our
population in this way then the measures of reactivity and attenuation regain
their original meaning. If the maximum initial growth rate of our population,
now scaled to be in terms of adults, is positive then our system is reactive,
and if the minimum is negative then it attenuates.

A biologically relevant measure of reactivity in a stage structured model must
then be focused on the initial growth rate of the population, calculated so
that the growth rates of juveniles are scaled by their contribution to the
adult population. Under this scaling if any adult on any habitat patch
produces many juveniles, but less than one become viable adults, then such a
metapopulation is not reactive. Whereas if there is a patch such that one
adult produces many juveniles and more than one survive to adulthood then the
metapopulation is reactive, because the stage structured population, where
juveniles are scaled according to their contribution to the adult population,
is growing.

In the following sections we formally define such a scaling using a weighted
$\ell_{1}$ norm, and contrast it with the unweighted $\ell_{1}$ norm that we
have previously been using to calculated reactivity in metapopulations without
stage structure.

###  6.1 Unweighted $\ell^{1}$

We first present the unweighted $\ell^{1}$ measure of the initial growth rate
to demonstrate the mathematical framework that we use to examine reactivity in
a stage-structured population with juveniles and adults. Now that there are
both juveniles and adults we want to measure reactivity and attenuation as the
total initial growth rate of the population, measured using either the
weighted or unweighted norm, when we start with one adult on a patch.

Consider a population on $n$ patches with a juvenile and adult stage. Let the
population dynamics be described by

$x^{\prime}=Ax,$ (13)

where $A$ is a $2n\times 2n$ matrix. Arrange $A$ so that all ODEs describing
the change in the adult populations are in rows $n+1$ to $2n$. Assume that
$A=F-V$, where $F$ is non-negative and $V$ is a non-singular M matrix so that
$V^{-1}$ is non-negative. We are interested in metapopulations where adults
can give birth to juveniles, but juveniles cannot give birth to new juveniles,
so $F$ and $V$ can be written in block form as follows:

$F=\begin{bmatrix}0&F_{12}\\\ 0&0\end{bmatrix},\quad
V=\begin{bmatrix}V_{11}&0\\\ V_{21}&V_{22}\end{bmatrix}.$ (14)

With this decomposition $F_{12}$ contains all the new juvenile births from
each adult patch, $V_{11}$ is a diagonal matrix that contains the rates of
juvenile mortality on each patch as well as the maturation from juveniles to
adults, $V_{22}$ is a diagonal matrix that contains the rates of adult
mortality on each patch, and $V_{21}$ contains the negative of the rates of
maturation/migration from juveniles to adults.

We define $\tilde{\lambda}_{j}$ to be the initial population growth rate,
starting with one adult on patch $j$, measured using the $\ell_{1}$ norm. This
can be defined mathematically for $1\leq j\leq n$ as

$\tilde{\lambda}_{j}=\underbrace{\sum_{i=1}^{n}x_{i}^{\prime}(0)}_{\text{juvenile}}+\underbrace{\sum_{i=n+1}^{2n}x_{i}^{\prime}(0)}_{\text{adult}},\quad
x(0)=e_{j+n}.$

In terms of $F$ and $V$

$\tilde{\lambda}_{j}=\sum_{i=1}^{n}f_{12ij}-\sum_{i=1}^{n}v_{22ij},$

where $f_{12ij}$ and $v_{22ij}$ are the $(i,j)$ entries of $F_{12}$ and
$V_{22}$ respectively.

We use the tilde to differentiate the initial growth rate in the stage
structured population, where we specifically begin with one adult on a patch,
from the initial growth rate in a population without stage structure, where
there is no difference in the type of individual that we start with. Having
presented the mathematical framework that we use to measure reactivity in a
stage structured population using the unweighted $\ell_{1}$ norm, we now use a
weighted $\ell_{1}$ norm that better captures the biological meaning of
reactivity.

###  6.2 Weighted $\ell^{1}$ for each patch

In order to measure reactivity in a biologically meaningful fashion, we
introduce a new measure of the initial population growth rate,
$\tilde{\lambda}_{j}^{p}$. This initial population growth rate is calculated
using a weighted $\ell_{1}$ norm so that the adult population is still
measured using the regular $\ell_{1}$ norm, but the juvenile population on
each patch is scaled by the probability that the juveniles survive to
adulthood; the patch specific nature of the weighing is why we denote the
initial growth rate $\tilde{\lambda}_{j}^{p}$. In this fashion
$\tilde{\lambda}_{j}^{p}$ measures the initial growth rate of the total
population if every member of the metapopulation was weighted according to how
much they will contribute to the adult population. Adults are therefore not
weighted, and juveniles are weighted by the probability that they survive to
adulthood. This weighting recaptures the biological meaning of reactivity,
where a system will only be reactive if the adult population will grow, and a
system will not be reactive if there is only an initial spike in the juvenile
population.

We use the same framework as the previous section to mathematically calculate
$\tilde{\lambda}_{j}^{p}$, where we decompose $A=F-V$ and $F$ and $V$ are
shown in block form in equation (14). Then we weight the juvenile population
growth on each patch $i$ by a factor $s_{i}$, where $s_{i}$ is the probability
that a juvenile from patch $i$ eventually becomes an adult. From the block
form $V$, $s_{i}$ can be calculated as

$s_{i}=\sum_{k=1}^{n}(-V_{21}V_{11}^{-1})_{ki}.$ (15)

To see how this corresponds to the probability of survival of a juvenile on
patch $i$, consider the different block components of $F$ and $V$. The matrix
$V_{11}^{-1}$ is diagonal, with the $(j,j)$ entry representing the average
residence time of a juvenile born onto patch $j$. The matrix $-V_{21}$
contains the rates of maturation/migration of juveniles becoming adults on
different patches, so the $(i,j)$ entry is the rate of maturation/migration of
a juvenile on patch $j$ becoming an adult on patch $i$. This means that when
we multiply $-V_{21}$ by $V_{11}^{-1}$ we are multiplying each of these rates
by the residence times of the juveniles in the appropriate patches. In this
way, the $(i,j)$ entry of $-V_{21}V_{11}^{-1}$ is then the probability that a
juvenile leaving patch $j$ arrives on patch $i$. Therefore the $j$th column
sum of $-V_{21}V_{11}^{-1}$ is the probability that a juvenile starting on
patch $j$ becomes an adult on any other patch.

The initial growth rate using the weighted norm, $\tilde{\lambda}_{j}^{p}$,
can then be calculated as the sum of the juvenile growth rates, each
multiplied by the patch specific survival $s_{i}$, and the adult growth rates.
Mathematically, this is defined as:

$\displaystyle\tilde{\lambda}_{j}^{p}$
$\displaystyle=\underbrace{\sum_{i=1}^{n}s_{i}x_{i}^{\prime}(0)}_{\text{juvenile}}+\underbrace{\sum_{i=n+1}^{2n}x_{i}^{\prime}(0)}_{\text{adult}},\quad
x(0)=e_{j+n}$
$\displaystyle=\sum_{i=1}^{n}s_{i}f_{12ij}-\sum_{i=1}^{n}v_{22ij}.$

In order to demonstrate that the initial growth rate calculated using the
weighted $\ell_{1}$ norm, $\tilde{\lambda}_{j}^{p}$, indeed measures the
growth rate of the population if all individuals are weighted according to
their contribution to the adult population, we show that
$\tilde{\lambda}_{j}^{p}$ is equivalent to scaling the juvenile population on
each patch by the probability of survival to adulthood, and then measuring the
initial growth rate using the unweighted $\ell_{1}$ norm, defined previously
as $\tilde{\lambda}_{j}$.

######  Lemma 6.1.

If each juvenile population on patch $j$ in system (13) is rescaled by the
patch specific survival probability,
$s_{i}=\sum_{k=1}^{n}(-V_{21}V_{11}^{-1})_{ki}$, then the initial growth rate
using the unweighted $\ell^{1}$ norm, $\tilde{\lambda}_{j}$, is equal to the
patch specific weighted initial growth rate, $\tilde{\lambda}_{j}^{p}$ for the
unscaled system.

###### Proof.

Rescale the juvenile population on patch $i$ by the patch specific survival
probability $s_{i}$ given in equation (15). In terms of system (13) this means
that $x_{i}^{*}=s_{i}x_{i}$ for $i=1,\dots,n$, $x_{i}^{*}=x_{i}$ for
$i=n+1,\dots,2n$. Rewriting the system of equations

$\displaystyle x^{*\prime}$ $\displaystyle=A^{*}x^{*}$ $\displaystyle
A^{*}=F^{*}-V^{*},\quad F^{*}$ $\displaystyle=\begin{bmatrix}0&SF_{12}\\\
0&0\end{bmatrix},\quad V^{*}=\begin{bmatrix}V_{11}&0\\\
V_{21}S^{-1}&V_{22}\end{bmatrix},$

where $S=\operatorname{diag}(s_{1},\dots,s_{n})$. The unweighted initial
growth rate for the scaled system is then

$\displaystyle\tilde{\lambda}_{j}^{*}$
$\displaystyle=\sum_{i=1}^{n}(SF_{12})_{ij}-\sum_{i=1}^{n}v_{22ij}$
$\displaystyle=\sum_{i=1}^{n}s_{i}f_{12ij}-\sum_{i=1}^{n}v_{22ij}$
$\displaystyle=\tilde{\lambda}_{j}^{p}.$

Thus the unweighted initial growth rate for the scaled system is equal to the
patch-weighted initial growth rate of the unscaled system,
$\tilde{\lambda}_{j}^{p}$.

∎

We believe that it is more intuitive to measure reactivity in a stage-
structured system using a weighted norm, rather than scaling the juvenile
population and using the unweighted $\ell_{1}$ norm, but for other systems
this may not be the case. Recently Mari et al. [30] have developed a new
measure of reactivity called generalized reactivity, or g-reactivity, so that
the reactivity of any specific combination of state variables in a system can
be measured, and we demonstrate how to place our work in this context. The
general framework of g-reactivity allows the reactivity of only the predator
to be measured in a predator-prey system, or a single patch in a
metapopulation model. Moreover in a stage-structured model, g-reactivity can
be used to allow for a differential contribution of the juvenile and adult
populations to the reactivity of the system, and so we can compare the
calculation of g-reactivity to our calculation using the weighted $\ell_{1}$
norm. To calculate the g-reactivity of a system $x^{\prime}=Ax$, a linear
transformation is introduced, $y=Cx$, where $C$ is a matrix that defines the
required contribution of each of the state variables, and then reactivity is
calculated for $y$ using equation (2). For system (13), g-reactivity is
equivalent to reactivity in our weighted norm,
$\max_{j}\tilde{\lambda}_{j}^{p}$, if $C$ is a $2n\times 2n$ identity matrix,
but with the first $n$ diagonal entries replaced with $s_{1},\dots,s_{n}$.

Returning to our measure of initial growth rate using a patch weighted norm,
we present a couple of examples below to illustrate the calculation of
$\tilde{\lambda}_{j}^{p}$ in different systems.

### Example 4

Consider a two patch system where juveniles are born onto all patches but only
mature into adults on the patch where they were born:

$\displaystyle j_{1}^{\prime}$
$\displaystyle=b_{11}a_{1}+b_{12}a_{2}-m_{1}j_{1}-d_{j1}j_{1}$ (16)
$\displaystyle j_{2}^{\prime}$
$\displaystyle=b_{21}a_{1}+b_{22}a_{2}-m_{2}j_{2}-d_{j2}j_{2}$ $\displaystyle
a_{1}^{\prime}$ $\displaystyle=m_{1}j_{1}-d_{a}a_{1}$ $\displaystyle
a_{2}^{\prime}$ $\displaystyle=m_{2}j_{2}-d_{a}a_{2}.$

Here $j_{i}$ is the number of juveniles on patch $i$, $a_{i}$ is the number of
adults on patch $i$, $b_{ij}$ is the birth rate of juveniles on patch $i$ from
adults on patch $j$, $m_{i}$ is the maturation rate of juveniles on patch $i$
into adults on patch $i$, $d_{ji}$ is the death rate of juveniles on patch
$i$, and $d_{a}$ is the death rate of adults, which is the same on both
patches. The lifecycle graph for this system is shown in Figure (9).

$j_{1}$$j_{2}$$a_{1}$$a_{2}$$b_{11}$$m_{1}$$b_{21}$$b_{12}$$m_{2}$$b_{22}$$d_{j1}$$d_{a1}$$d_{j2}$$d_{a2}$
Figure 9: Digraph for system (16). Here $b_{ij}$ is the birth rate of
juveniles on patch $i$ from adults on patch $j$, $m_{i}$ is the maturation
rate of juveniles on patch $i$ to adults on patch $i$, $d_{ji}$ is the death
rate of juveniles on patch $i$ and $d_{ai}$ is the death rate of adults on
patch $i$.

In this case

$\displaystyle F_{12}=\begin{bmatrix}b_{11}&b_{12}\\\
b_{21}&b_{22}\end{bmatrix},\quad V_{11}$
$\displaystyle=\begin{bmatrix}m_{1}+d_{j1}&0\\\
0&m_{2}+d_{j2}\end{bmatrix},\quad V_{21}=\begin{bmatrix}-m_{1}&0\\\
0&-m_{2}\end{bmatrix},$ $\displaystyle V_{22}=\begin{bmatrix}d_{a}&0\\\
0&d_{a}\end{bmatrix},\quad-V_{21}V_{11}^{-1}$
$\displaystyle=\begin{bmatrix}\frac{m_{1}}{m_{1}+d_{j1}}&0\\\
0&\frac{m_{2}}{m_{2}+d_{j2}}\end{bmatrix},$
$\displaystyle\tilde{\lambda}_{1}^{p}$
$\displaystyle=s_{1}b_{11}+s_{2}b_{21}-d_{a}$
$\displaystyle\tilde{\lambda}_{2}^{p}$
$\displaystyle=s_{1}b_{12}+s_{2}b_{22}-d_{a}$ $\displaystyle s_{1}$
$\displaystyle=\frac{m_{1}}{m_{1}+d_{j1}}$ $\displaystyle s_{2}$
$\displaystyle=\frac{m_{2}}{m_{2}+d_{j2}}.$

If we look at $\tilde{\lambda}_{1}^{p}$, we see that $s_{1}$ is the
probability that a juvenile born onto patch 1 survives to become an adult and
it is multiplying $b_{11}$, the birth rate of juveniles onto patch 1 from
adults on patch 1. Therefore the first component of $\tilde{\lambda}_{1}^{p}$
represents the rate of birth of new juveniles onto patch 1 from one adult on
patch 1, but scaled by the probability that these juveniles survive to become
adults. Likewise the second component, $s_{2}b_{21}$, is the rate of birth of
new juveniles onto patch 2 from one adult on patch 1, scaled by the
probability that those juveniles also become adults. Thus
$\tilde{\lambda}_{1}^{p}$ is the initial growth rate of the total population,
scaled in terms of the contribution to the adult population, when the
population begins with one adult on patch 1.

### Example 5

We also consider a system in which juveniles are born onto the same patch as
adults, but can then migrate between patches as they mature into adults:

$\displaystyle j_{1}^{\prime}$
$\displaystyle=b_{11}a_{1}-m_{11}j_{1}-m_{21}j_{1}-d_{j1}j_{1}$ (17)
$\displaystyle j_{2}^{\prime}$
$\displaystyle=b_{22}a_{2}-m_{22}j_{2}-m_{12}j_{2}-d_{j2}j_{2}$ $\displaystyle
a_{1}^{\prime}$ $\displaystyle=m_{11}j_{1}+m_{12}j_{2}-d_{a}a_{1}$
$\displaystyle a_{2}^{\prime}$
$\displaystyle=m_{22}j_{2}+m_{21}j_{1}-d_{a}a_{2}.$

from which we calculate

$\displaystyle F_{12}=\begin{bmatrix}b_{11}&0\\\ 0&b_{22}\end{bmatrix},\quad$
$\displaystyle V_{11}=\begin{bmatrix}m_{11}+m_{21}+d_{j1}&0\\\
0&m_{12}+m_{22}+d_{j2}\end{bmatrix},\quad
V_{21}=\begin{bmatrix}-m_{11}&-m_{12}\\\ -m_{21}&-m_{22}\end{bmatrix},$
$\displaystyle V_{22}=\begin{bmatrix}d_{a}&0\\\ 0&d_{a}\end{bmatrix},\quad-$
$\displaystyle
V_{21}V_{11}^{-1}=\begin{bmatrix}\frac{m_{11}}{m_{11}+m_{21}+d_{j1}}&\frac{m_{12}}{m_{12}+m_{22}+d_{j2}}\\\
\frac{m_{21}}{m_{11}+m_{21}+d_{j1}}&\frac{m_{22}}{m_{12}+m_{22}+d_{j2}}\end{bmatrix},$
$\displaystyle\tilde{\lambda}_{1}^{p}$ $\displaystyle=s_{1}b_{11}-d_{a}$
$\displaystyle\tilde{\lambda}_{2}^{p}$ $\displaystyle=s_{2}b_{22}-d_{a}$
$\displaystyle s_{1}$
$\displaystyle=\frac{m_{11}+m_{21}}{m_{11}+m_{21}+d_{j1}}$ $\displaystyle
s_{2}$ $\displaystyle=\frac{m_{12}+m_{22}}{m_{12}+m_{22}+d_{j2}}.$

$j_{1}$$j_{2}$$a_{1}$$a_{2}$$b_{11}$$m_{11}$$m_{12}$$m_{21}$$m_{22}$$b_{22}$$d_{j1}$$d_{a1}$$d_{j2}$$d_{a2}$
Figure 10: Digraph for system (17). Here $b_{ii}$ is the birth rate of
juveniles on patch $i$ from adults on patch $i$, $m_{ij}$ is the maturation
rate of juveniles on patch $j$ to adults on patch $i$, $d_{ji}$ is the death
rate of juveniles on patch $i$ and $d_{ai}$ is the death rate of adults on
patch $i$.

Now if we examine the first component of $\tilde{\lambda}_{1}^{p}$,
$s_{1}=(m_{11}+m_{21})/(m_{11}+m_{21}+d_{j1})$, we can see that because
juveniles from patch 1 can now migrate (as they mature) to both patches,
$s_{1}$ is the probability that juveniles from patch 1 become adults on either
patch. Likewise $s_{2}$ is the probability that juveniles from patch 2 become
adults on either patch.

Here we have shown that if we use a weighted $\ell_{1}$ norm to scale the
initial growth rate so that the juvenile population is scaled by the patch
specific probability that juveniles become adults, then our scaled initial
growth rate, $\tilde{\lambda}_{j}^{p}$ matches the biological intuition that
we would like when measuring initial growth of the population. It is positive
if the population, scaled so that every individual is measured by its
contribution to the adult population, is growing, and negative if the
population is decreasing. Measures of reactivity and attenuation then also
represent their intuitive biological properties, and we are no longer in the
situation (as if the initial growth rate was unscaled) that most marine
metapopulations are reactive.

It is also possible to create a weighted norm where the juvenile populations
on each patch are weighted by the same probability of survival, rather than by
patch-specific probabilities $s_{i}$. In some cases it may be useful to scale
all patches by the same survival probability, though under this weighted norm
reactivity no longer corresponds exactly to the intuitive biological meaning
mentioned previously.

##  7 Discussion

Transient dynamics often differ drastically from the asymptotic dynamics of a
system and require different analytical tools. In this paper we have presented
a unifying framework for analysing transient dynamics in marine and other
birth-jump metapopulations, from the choice of norms to the incorporation of
stage structure. We began by using the $\ell_{1}$ norm to define reactivity
and attenuation in single species metapopulations and used examples to compare
reactivity in the $\ell_{1}$ norm with reactivity in the more commonly used
$\ell_{2}$ norm. We presented a couple of two patch models that gave rise to
long transients: one stable system that exhibits a long period of growth
before eventual decay, and one unstable system that exhibits a long period of
decay before growth. In marine metapopulations, where patches are connected
via larval dispersal, we examined how strong advective flow, coupled with a
large number of patches, can lead to large transient growth. We then connected
the initial growth rate of the metapopulation to the source-sink distribution
of patches and lastly we demonstrated how to measure reactivity meaningfully
in stage-structured marine metapopulations.

We are by no means the first to analyze the transient dynamics of systems, and
in fact there has been an increase in the study of transient dynamics over the
last couple of decades. In a pair of recent papers, several authors have
identified certain mechanisms as the main causes of long transients in
ecological systems [16, 32]. These identified mechanisms that cause the long
transients present in the examples in this paper are slow-fast systems, crawl-
bys, and high dimensionality. Slow-fast systems cause long transients when the
system rapidly converges to a slow manifold, then moves slowly towards or away
from an equilibrium, depending on the stability of the system. This occurs in
both examples in section 3. The second example in section 3 is also an
instance of a crawl-by where the initial condition is near a saddle
equilibrium but the movement away from the equilibrium occurs over a long
timescale. Lastly in section 4 we explicitly demonstrated how increasing the
dimension of a system, by increasing the patch number in a linear
metapopulation, leads to longer transients.

Our work also reinforces the fact that reactivity is a property specific to
the norm under which it is measured. This has been mentioned in the first
paper on reactivity by Neubert and Caswell[35], who also recognize that it is
always possible to find a norm such that a stable system is never reactive. It
has also been noted by Lutscher and Wang [29], who mention that reactivity
must be analysed in the dimensional version of a system rather than the non-
dimensionalized version. The reactivity may be different between the two
systems but the dimensional system is where the measure of reactivity is
biologically meaningful. When analysing reactivity in metapopulations this
fact is significant in two ways: first by using the $\ell_{1}$ norm rather
than the $\ell_{2}$ norm to measure reactivity we can explicitly measure the
growth rate of a population, and second by using a weighted $\ell_{1}$ norm we
prevent the juvenile population from disproportionately affecting the
reactivity of the system.

Differentially weighting certain classes of a population to calculate
reactivity has been mentioned in passing by Verdy and Caswell [51], and more
extensively by Mari et al. [30] who developed a new measure of reactivity
called general reactivity, or g-reactivity. This is a method of only measuring
the reactivity of the components of interest in a population, e.g. predators
in a predator-prey model, and can also be used more generally to scale the
contribution of different components of the population. Our method of using a
weighted $\ell_{1}$ norm for stage-structured models has an equivalent
formulation using g-reactivity that is discussed in section 6.2, though Mari
et al. [30] use the $\ell_{2}$ norm to measure reactivity, rather than the
$\ell_{1}$ norm, and are thus using a different measure of population growth.

While we believe the $\ell_{1}$ norm is the most biologically relevant norm to
measure reactivity, we are among the first to use it to analyse reactivity in
continuous time models. Townley et al. [48] show how to calculate reactivity
for stage-structured models in continuous time using the $\ell_{1}$ norm, but
in following papers proceed to analyse reactivity in the $\ell_{1}$ norm only
in discrete time systems [49, 43, 44, 45]. Most authors measure reactivity
with the $\ell_{2}$ norm, presumably due to the nice mathematical property
that reactivity in the $\ell_{2}$ norm can be measured simply as
$\lambda(H(A))$, where $H(A)=(A+A^{T})/2$ [35, 5, 36, 37, 51, 29, 42]. But
while mathematically tractable, the biological meaning of Euclidean distance
($\ell_{2}$) is less clear than population size ($\ell_{1}$) and as shown in
section 2.2, there are times when reactivity in $\ell_{2}$ does not correspond
to an increase in population size.

In addition to being biologically meaningful, reactivity in the $\ell_{1}$
norm is also simple to calculate for single species metapopulations due to the
fact that the populations on each patch are non-negative, and thus the full
solution remains in the non-negative cone. Mathematically, a metapopulation
governed by $x^{\prime}=Ax$ will be reactive in the $\ell_{1}$ norm if
$\max_{j}\sum_{i=1}^{n}a_{ij}>0$. Geometrically, a system is reactive in the
$\ell_{1}$ norm if the dot product of the derivative vector of any initial
condition and the outward normal vector of the plane $x_{1}+x_{2}+\dots
x_{n}=1$ is positive. If we instead want to consider the reactivity in a
system where this assumption is relaxed, say a predator-prey metapopulation,
then the mathematical calculation using the $\ell_{1}$ norm is no longer so
simple, but our geometric intuition is similar. Instead of verifying if the
dot product of the derivative vector of any initial condition and the outward
normal vector of the plane $x_{1}+x_{2}+\dots x_{n}=1$ is positive, see for
example Figure 1, we would instead need to take the dot product with the
outward normal vector to the hypercube $|x_{1}|+|x_{2}|+\dots|x_{n}|=1$. If we
continue to use $\max_{j}\sum_{i=1}^{n}a_{ij}>0$ to calculate reactivity for
systems that do not remain in the non-negative cone then it is possible for
systems to not be reactive and yet also be unstable, thus an interesting area
for future work would be to mathematically formulate reactivity in terms of
the matrix $A$ for these metapopulations.

No matter the norm in which reactivity and attenuation are measured, they are
defined in terms of the linearization of a non-linear system around an
equilibrium. As mentioned in the Introduction, reactivity and attenuation are
therefore most relevant around hyperbolic equilibria, where the dynamics of
the non-linear system are well approximated by the linearized system. Another
caveat is that it is possible for an equilibrium of a non-linear system to not
be reactive, but for a perturbation of the non-linear system to still cause a
large excursion away from the stable equilibrium before eventually returning.
Excitable systems, such as the FitzHugh-Nagumo system, have stable equilibria
with attracting regions, but small perturbations still trigger large
excitations [10, 34]. These systems may not be reactive from the linearized
definition of reactivity, but can still exhibit similar behaviour to
reactivity in the non-linear system, given a sufficient perturbation.

The final extension that we would like to highlight is the relationship
between reactivity of continuous time models and reactivity of their discrete
counterparts. Many marine metapopulations are modelled in discrete time due to
yearly breeding cycles, however some are modelled in discrete time due to ease
of simulation. For these models, where the time step is on the order of hours
or days, we can connect the reactivity of the continuous time system with the
discrete time system using a Taylor expansion. The continuous time system,
$x^{\prime}=Ax$, has the solution $x(t)=e^{At}x_{0}$ that could be sampled at
discrete time steps $\tau$ to create the discrete time system
$x(t+\tau)=Bx(t),$ where $B=e^{A\tau}$.

The continuous time system $x^{\prime}=Ax$ is reactive in $\ell_{1}$ if $A$
has a positive column sum (Lemma 2.1). In discrete time the system
$x(t+\tau)=Bx(t)$ is reactive in $\ell_{1}$ if $B$ has a column sum that is
greater than 1 [48]. Assuming $\tau$ is a small timestep then we can
approximate $B=e^{A\tau}=I+A\tau+O(\tau^{2})$. Thus we can see that if the
system is reactive in continuous time, i.e. there is a positive column sum of
$A$, then we can find a sufficiently short time step $\tau$ such that the
discrete time system is also reactive, i.e. there is a column sum of $B$
greater than 1. However for a predetermined time-step $\tau$ there are
continuous time systems $x^{\prime}=Ax$ that are reactive but for which their
discrete counterparts $x(t+\tau)=e^{A\tau}x(t)$ are not reactive. One such
example is system (8) with $\epsilon=0.9$ and $\tau=1$.

Lastly, we hope our work can be used to better understand the transient
dynamics in marine metapopulations for which habitat patches have already been
classified as sources and sinks. For these systems the transient dynamics that
may occur following a disturbance depend directly on the new distribution of
the population. If the remaining population is distributed amongst sink
patches then it initially declines, even if it eventually recovers. Likewise
if the population is distributed amongst source patches then it initially
grows, though this growth may not necessarily occur on the source patch
itself. In addition, the relationship between transient dynamics and the
source-sink distribution of marine metapopulations may also be useful when
examining the dynamics that can occur following the protection of new marine
environments, such as newly implemented Marine Protected Areas.

##  Appendix A Two patch example with migration

Consider a different metapopulation on two patches where individuals are born
only onto their patch, but can now also migrate between patches. This is the
case for many terrestrial species that live on patchy landscapes, where
individuals can migrate between habitat patches. The dynamics of this
metapopulation is then described with the following set of ODEs:

$\displaystyle x_{1}^{\prime}$
$\displaystyle=b_{1}x_{1}-m_{21}x_{1}+m_{12}x_{2}-d_{1}x_{1}$ (18)
$\displaystyle x_{2}^{\prime}$
$\displaystyle=b_{2}x_{2}-m_{12}x_{2}+m_{21}x_{1}-d_{2}x_{2},$

where $b_{i}$ is the birth rate on patch $i$, $m_{ij}$ is the migration rate
from patch $j$ to patch $i$, and $d_{i}$ is the death rate on patch $i$. The
lifecycle graph for this system is shown in Figure 11.

$x_{1}$$x_{2}$$m_{21}$$m_{12}$$b_{11}-d_{1}$$b_{22}-d_{2}$ Figure 11: Digraph
for system (18). The directed edges represent the movement of individuals from
the outgoing patch to the incoming patch. The self loops are the birth rate
minus the death rate on a patch.

We again decompose $A=F-V$ and construct our next generation matrix,
$K=FV^{-1}$, though now $V$ is not diagonal nor in the same form as required
for Lemma 5.2.

$\displaystyle A$ $\displaystyle=\begin{bmatrix}b_{1}-m_{21}-d_{1}&m_{12}\\\
m_{21}&b_{2}-m_{12}-d_{2}\end{bmatrix}$ $\displaystyle F$
$\displaystyle=\begin{bmatrix}b_{1}&0\\\ 0&b_{2}\end{bmatrix}$ $\displaystyle
V$ $\displaystyle=\begin{bmatrix}d_{1}+m_{21}&-m_{12}\\\
-m_{21}&d_{2}+m_{12}\end{bmatrix}$ $\displaystyle K=FV^{-1}$
$\displaystyle=\begin{bmatrix}\frac{b_{1}(d_{2}+m_{12})}{d_{1}d_{2}+d_{1}m_{12}+d_{2}m_{21}}&\frac{b_{1}m_{12}}{d_{1}d_{2}+d_{1}m_{12}+d_{2}m_{21}}\\\
\frac{b_{2}m_{21}}{d_{1}d_{2}+d_{1}m_{12}+d_{2}m_{21}}&\frac{b_{2}(d_{1}+m_{21})}{d_{1}d_{2}+d_{1}m_{12}+d_{2}m_{21}}\end{bmatrix}$

The entries of $K$ may seem counter intuitive, but they represent the infinite
sum of a geometric series. Consider the first entry, $k_{11}$. If we are
tracking the total number of new individuals produced by one individual
starting on patch 1, then this individual can either produce new offspring in
patch 1 immediately, or it can migrate to patch 2, then back to patch 1 and
produce offspring, or migrate again and produce more offspring. The entry
$k_{11}$ is then the birth rate in patch $1$ multiplied by the residence time
in patch $1$, multiplied by a geometric series where the ratio is the
probability of surviving the migration from patch 1 to patch 2 and then back
to patch 1. Mathematically

$k_{11}=\frac{b_{1}}{m_{21}+d_{1}}\left[\sum_{i=0}^{\infty}\left(\frac{m_{21}}{m_{21}+d_{1}}\frac{m_{12}}{m_{12}+d_{2}}\right)^{i}\right]$

Now consider $\lambda_{1}$ and $R_{1}$, the measures of transient growth for
patch 1:

$\displaystyle\lambda_{1}$ $\displaystyle=b_{1}-d_{1}$ $\displaystyle R_{1}$
$\displaystyle=\frac{b_{1}d_{2}+b_{1}m_{12}+b_{2}m_{21}}{d_{1}d_{2}+d_{1}m_{12}+d_{2}m_{21}}.$

We can see that

$\lim_{d_{2}\rightarrow\infty}R_{1}=\frac{b_{1}}{d_{1}+m_{21}}.$

Therefore even if $\lambda_{1}>0$, and so $b_{1}>d_{1}$, as $m_{21}$ becomes
large, $R_{1}<1$.

In the other direction, if $d_{1}=2b_{1}$, then

$\lim_{b_{1}\rightarrow 0}R_{1}=\frac{b_{2}}{d_{2}}.$

Therefore even if $R_{1}>1$, we can still have $\lambda_{1}=b_{1}-d_{1}<0$. We
present this example to demonstrate that if the assumptions of Lemma 5.1 are
not met, there is no longer a one-to-one relationship between $\lambda_{j}$
and $R_{j}$.

Now we might also consider moving the off diagonal entries of $V$ into $F$ so
that $V$ becomes diagonal. This is similar to considering migrating
individuals as new individuals entering a patch. In this case

$\displaystyle F$ $\displaystyle=\begin{bmatrix}b_{1}&m_{12}\\\
m_{21}&b_{2}\end{bmatrix}$ $\displaystyle V$
$\displaystyle=\begin{bmatrix}d_{1}+m_{21}&0\\\ 0&d_{2}+m_{12}\end{bmatrix}$
$\displaystyle K=FV^{-1}$
$\displaystyle=\begin{bmatrix}\frac{b_{1}}{d_{1}+m_{21}}&\frac{m_{12}}{d_{2}+m_{12}}\\\
\frac{m_{21}}{d_{1}+m_{21}}&\frac{b_{2}}{d_{2}+m_{12}}\end{bmatrix}$

Here $\lambda_{1}=b_{1}-d_{1}>0$ is equivalent to $b_{1}>d_{1}$, which is then
equivalent to $R_{1}=(b_{1}+m_{21})/(d_{1}+m_{21})>1$. However, in this case
$R_{1}$ no longer tracks the total number of new individuals produced on all
patches over one generation. This new $R_{1}$ could perhaps be interpreted as
the total number of new individuals produced on patch 1 by a single individual
on patch 1 before that individual dies or migrates, plus the probability that
the individual migrates to patch 2 before it dies. However, this will no
longer be a biologically useful measure of a source or a sink.

## Acknowledgments

The authors would like to thank the members of the Lewis Lab for many helpful
discussions and suggestions. PDH gratefully acknowledges an NSERC-CGSM
scholarship, Queen Elizabeth II scholarship, and Alberta Graduate Excellence
Scholarship; MAL gratefully acknowledges the Canada Research Chair program and
an NSERC Discovery grant; and PvdD gratefully acknowledges an NSERC Discovery
grant.

## References

  * [1] T. P. Adams, R. Proud, and K. D. Black, Connected networks of sea lice populations: dynamics and implications for control, Aquac. Environ. Interact., 6 (2015), pp. 273–284, https://doi.org/10.3354/aei00133.
  * [2] J. Arino and P. van den Driessche, A multi-city epidemic model, Math. Popul. Stud., 10 (2003), pp. 175–193, https://doi.org/10.1080/08898480306720.
  * [3] P. R. Armsworth, Recruitment limitation, population regulation, and larval connectivity in reef fish metapopulations, Ecology, 83, pp. 1092–1104, https://doi.org/10.1890/0012-9658(2002)083[1092:RLPRAL]2.0.CO;2.
  * [4] L. W. Botsford, C. L. Moloney, A. Hastings, J. L. Largier, T. M. Powell, K. Higgins, and J. F. Quinn, The influence of spatially and temporally varying oceanographic conditions on meroplanktonic metapopulations, Deep. Res. Part II, 41 (1994), pp. 107–145, https://doi.org/10.1016/0967-0645(94)90064-7.
  * [5] H. Caswell and M. G. Neubert, Reactivity and transient dynamics of discrete-time ecological systems, J. Differ. Equations Appl., 11 (2005), pp. 295–310, https://doi.org/10.1080/10236190412331335382.
  * [6] R. K. Cowen, C. B. Paris, and A. Srinivasan, Scaling of connectivity in marine populations, Science, 311 (2006), pp. 522–527, https://doi.org/10.1126/science.1122039.
  * [7] R. K. Cowen and S. Sponaugle, Larval dispersal and marine population connectivity, Ann. Rev. Mar. Sci., 1 (2009), pp. 443–466, https://doi.org/10.1146/annurev.marine.010908.163757.
  * [8] J. M. Cushing and O. Diekmann, The many guises of $R_{0}$ (a didactic note), J. Theor. Biol., 404 (2016), pp. 295–302, https://doi.org/10.1016/j.jtbi.2016.06.017.
  * [9] W. F. Figueira and L. B. Crowder, Defining patch contribution in source-sink metapopulations: the importance of including dispersal and its relevance to marine systems, Popul. Ecol., 48 (2006), pp. 215–224, https://doi.org/10.1007/s10144-006-0265-0.
  * [10] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophys. J., 1 (1961), pp. 445–466, https://doi.org/10.1016/S0006-3495(61)86902-6.
  * [11] M. Gyllenberg, A. Hastings, and I. Hanski, 5 - Structured metapopulation models, in Metapopulation Biology, I. Hanski and M. E. Gilpin, eds., Academic Press, San Diego, 1997, pp. 93–122, https://doi.org/10.1016/B978-012323445-2/50008-0.
  * [12] I. Hanski and C. D. Thomas, Metapopulation dynamics and conservation: a spatially explicit model applied to butterflies, Biol. Conserv., 68 (1994), pp. 167–180, https://doi.org/10.1016/0006-3207(94)90348-4.
  * [13] P. D. Harrington and M. A. Lewis, A next-generation approach to calculate source–sink dynamics in marine metapopulations, Bull. Math. Biol., 82 (2020), pp. 1–44, https://doi.org/10.1007/s11538-019-00674-1.
  * [14] A. Hastings, Transient dynamics and persistence of ecological systems, Ecol. Lett., 4 (2001), pp. 215–220, https://doi.org/10.1046/j.1461-0248.2001.00220.x.
  * [15] A. Hastings, Transients: the key to long-term ecological understanding?, Trends Ecol. Evol., 19 (2004), pp. 39–45, https://doi.org/10.1016/j.tree.2003.09.007.
  * [16] A. Hastings, K. C. Abbott, K. Cuddington, T. Francis, G. Gellner, Y. C. Lai, A. Morozov, S. Petrovskii, K. Scranton, and M. L. Zeeman, Transient phenomena in ecology, Science, 361 (2018), https://doi.org/10.1126/science.aat6412.
  * [17] A. Hastings and K. Higgins, Persistence of transients in spatially structured ecological models, Science, 263 (1994), pp. 1133–1136, https://doi.org/10.1126/science.263.5150.1133.
  * [18] K. Higgins, A. Hastings, J. N. Sarvela, and L. W. Botsford, Stochastic dynamics and deterministic skeletons: population behavior of Dungeness crab, Science, 276 (1997), pp. 1431–1434, https://doi.org/10.1126/science.276.5317.1431.
  * [19] T. Hillen, B. Greese, J. Martin, and G. de Vries, Birth-jump processes and application to forest fire spotting, J. Biol. Dyn., 9 (2015), pp. 104–127, https://doi.org/10.1080/17513758.2014.950184.
  * [20] R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, New York, NY, USA, 2nd ed., 2012.
  * [21] Q. Huang, Y. Jin, and M. A. Lewis, $R_{0}$ analysis of a benthic-drift model for a stream population, SIAM J. Appl. Dyn. Syst., 15 (2016), pp. 287–321, https://doi.org/10.1137/15M1014486.
  * [22] Q. Huang and M. A. Lewis, Homing fidelity and reproductive rate for migratory populations, Theor. Ecol., 8 (2015), pp. 187–205, https://doi.org/10.1007/s12080-014-0243-7.
  * [23] B. C. Husband and S. C. H. Barrett, A metapopulation perspective in plant population biology, J. Ecol., 84 (1996), pp. 461–469, https://doi.org/10.2307/2261207.
  * [24] G. Jones, G. Almany, G. Russ, P. Sale, R. Steneck, M. van Oppen, and B. Willis, Larval retention and connectivity among populations of corals and reef fishes: history, advances and challenges, Coral Reefs, 28 (2009), pp. 307–325, https://doi.org/10.1007/s00338-009-0469-9.
  * [25] M. Krkošek and M. A. Lewis, An $R_{0}$ theory for source-sink dynamics with application to Dreissena competition, Theor. Ecol., 3 (2010), pp. 25–43, https://doi.org/10.1007/s12080-009-0051-7.
  * [26] R. Levins, Some demographic and genetic consequences of environmental heterogeneity for biological control, Bull. Entomol. Soc. Am., 15 (1969), pp. 237–250, https://doi.org/10.1088/2058-9565/aaba1a.
  * [27] A. L. Lloyd and R. M. May, Spatial heterogeneity in epidemic models, J. Theor. Biol., 179 (1996), pp. 1–11, https://doi.org/10.1006/jtbi.1996.0042.
  * [28] D. Ludwig, D. D. Jones, and C. S. Holling, Qualitative analysis of insect outbreak systems: the spruce budworm and forest, J. Anim. Ecol., 47 (1978), pp. 315–332, https://doi.org/10.2307/3939.
  * [29] F. Lutscher and X. Wang, Reactivity of communities at equilibrium and periodic orbits, J. Theor. Biol., 493 (2020), p. 110240, https://doi.org/10.1016/j.jtbi.2020.110240.
  * [30] L. Mari, R. Casagrandi, A. Rinaldo, and M. Gatto, A generalized definition of reactivity for ecological systems and the problem of transient species dynamics, Methods Ecol. Evol., 8 (2017), pp. 1574–1584, https://doi.org/10.1111/2041-210X.12805.
  * [31] H. W. Mckenzie, Y. Jin, J. Jacobsen, and M. A. Lewis, $R_{0}$ analysis of a spatiotemporal model for a stream population, SIAM J. Appl. Dyn. Syst., 11 (2012), pp. 567–596, https://doi.org/10.1137/100802189.
  * [32] A. Morozov, K. Abbott, K. Cuddington, T. Francis, G. Gellner, A. Hastings, Y. C. Lai, S. Petrovskii, K. Scranton, and M. L. Zeeman, Long transients in ecology: theory and applications, Phys. Life Rev., 32 (2020), pp. 1–40, https://doi.org/10.1016/j.plrev.2019.09.004.
  * [33] R. F. Morris, The dynamics of epidemic spruce budworm populations, Mem. Entomol. Soc. Canada, 95 (1963), pp. 1–12, https://doi.org/10.4039/entm9531fv.
  * [34] J. Nagumo, S. Arimoto, and S. Yoshizawa, An active pulse transmission line simulating nerve axon, Proc. IRE, 50 (1962), pp. 2061–2070, https://doi.org/10.1109/JRPROC.1962.288235.
  * [35] M. G. Neubert and H. Caswell, Alternatives to resilience for measuring the responses of ecological systems to perturbations, Ecology, 78 (1997), pp. 653–665, https://doi.org/10.1890/0012-9658(1997)078[0653:ATRFMT]2.0.CO;2.
  * [36] M. G. Neubert, H. Caswell, and J. D. Murray, Transient dynamics and pattern formation: reactivity is necessary for Turing instabilities, Math. Biosci., 175 (2002), pp. 1–11, https://doi.org/10.1016/S0025-5564(01)00087-6.
  * [37] M. G. Neubert, T. Klanjscek, and H. Caswell, Reactivity and transient dynamics of predator-prey and food web models, Ecol. Modell., 179 (2004), pp. 29–38, https://doi.org/10.1016/j.ecolmodel.2004.05.001.
  * [38] S. Noschese, L. Pasquini, and L. Reichel, Tridiagonal Toeplitz matrices: properties and novel applications, Numer. Linear Algebra Appl., 20 (2013), pp. 302–326, https://doi.org/10.1002/nla.1811.
  * [39] H. R. Pulliam, Sources, sinks and population regulation, Am. Nat., 132 (1988), pp. 652–661, https://doi.org/10.1086/284880.
  * [40] J. Roughgarden, S. Gaines, and H. Possingham, Recruitment dynamics in complex life cycles, Science, 241 (1988), pp. 1460–1466, https://doi.org/10.1126/science.11538249.
  * [41] L. Safranyik and A. L. Carroll, The biology and epidemiology of the mountain pine beetle in lodgepole pine forests., Canadian Forest Service, Victoria, Canada, 2007, pp. 3–66.
  * [42] R. E. Snyder, What makes ecological systems reactive?, Theor. Popul. Biol., 77 (2010), pp. 243–249, https://doi.org/10.1016/j.tpb.2010.03.004.
  * [43] I. Stott, M. Franco, D. Carslake, S. Townley, and D. Hodgson, Boom or bust? A comparative analysis of transient population dynamics in plants, J. Ecol., 98 (2010), pp. 302–311, https://doi.org/10.1111/j.1365-2745.2009.01632.x.
  * [44] I. Stott, S. Townley, and D. J. Hodgson, A framework for studying transient dynamics of population projection matrix models, Ecol. Lett., 14 (2011), pp. 959–970, https://doi.org/10.1111/j.1461-0248.2011.01659.x.
  * [45] I. Stott, S. Townley, and D. J. Hodgson, A framework for studying transient dynamics of population projection matrix models, Ecol. Lett., 14 (2011), pp. 959–970, https://doi.org/10.1111/j.1461-0248.2011.01659.x.
  * [46] H. R. Thieme, Mathematics in Population Biology, Princeton University Press, 2003.
  * [47] H. R. Thieme, Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity, SIAM J. Appl. Math., 70 (2009), pp. 188–211, https://doi.org/10.1137/080732870.
  * [48] S. Townley, D. Carslake, O. Kellie-Smith, D. Mccarthy, and D. Hodgson, Predicting transient amplification in perturbed ecological systems, J. Appl. Ecol., 44 (2007), pp. 1243–1251, https://doi.org/10.1111/j.1365-2664.2007.01333.x.
  * [49] S. Townley and D. J. Hodgson, Erratum et addendum: Transient amplification and attenuation in stage-structured population dynamics, J. Appl. Ecol., 45 (2008), pp. 1836–1839, https://doi.org/10.1111/j.1365-2664.2008.01562.x.
  * [50] M. Tsatsomeros, Matrix equalities and inequalities, in Handbook of Linear Algebra, L. Hogben, ed., CRC Press, 2006, ch. 14.
  * [51] A. Verdy and H. Caswell, Sensitivity analysis of reactive ecological dynamics, Bull. Math. Biol., 70 (2008), pp. 1634–1659, https://doi.org/10.1007/s11538-008-9312-7.
  * [52] D. W. Williams and A. M. Liebhold, Spatial synchrony of spruce budworm outbreaks in eastern North America, Ecology, 81 (2000), pp. 2753–2766, https://doi.org/10.1890/0012-9658(2000)081[2753:SSOSBO]2.0.CO;2.