Stability analysis of the coexistence equilibrium of a balanced metapopulation model

We now derive the metapopulation model used in this paper. We start by deriving a general metapopulation model that is based on the seminal work of Levin⁵⁰. Assuming that the inter-patch migrations are detailed-balanced, we make use of the formulation in Eq. (8) to derive a balanced metapopulation model. We then show that the balanced model admits a unique coexistence equilibrium that is asymptotically stable if the dispersal network is heterogeneous, whereas the same equilibrium is neutrally stable in the case of a homogeneous network.

General metapopulation model

Mathematical models based on traditional metapopulation theory usually assume that the metapopulation is made up of many neighboring spatially homogeneous habitat patches connected via dispersal. Consider an interconnected network of m discrete patches each being inhabited by the same n species. In addition, assume that species can migrate from one patch to some or all of the other patches. The rate of migration of each species between two patches is directly proportional to the proportion of the particular species in the originating patch, with a (nonnegative) constant of proportionality being the same across species. This constant of proportionality will be referred to as the rate constant associated with the migration. It is assumed that if there is migration between two given patches, then it is bidirectional, i.e., the rate constant of migration from j to k is strictly positive if and only the same holds for the migration from k to j. Just like in the case of a reversible single-species chemical reaction network, inter-patch migrations may be described by a weighted symmetric directed graph (G_2=(V_2,E_2)) where (V_2={1,ldots ,m}) is the set of patches (vertices) and an edge ((j,k)in E_2) means that every species can migrate from patch j to patch k. Finally, it is also assumed that the graph (G_2) corresponding to the inter-patch migration is connected, i.e., there is a path between every two distinct vertices of the graph.

The flow of species between the patches can be summarized in a weighted (mtimes m) adjacency matrix ({mathbf {A}}) with entry (A_{jk}) being equal to the rate constant of migration of species from the (j{ {text {th}}}) to the (k{ {text {th}}}) patch. The diagonal elements of ({mathbf {A}}) are hence equal to 0. Due to the bidirectional nature of migration, it holds that (A_{jk}>0 Leftrightarrow A_{kj}>0) and (A_{jk}=0 Leftrightarrow A_{kj}=0), for any (jne k). Let (Delta =text {diag}(delta _1,ldots ,delta _m)) denote the m-dimensional diagonal matrix whose (j{ {text {th}}}) entry is given by

Define ({mathbf {L}}:=Delta -{mathbf {A}}^top ). Note that

Let ({mathbf {x}}in S^{mn}), with (x_{i,j}) the proportion of species i in patch j across the entire metapopulation, then the net migration rate (psi _{i,j}) of species i from other patches to patch j is given by

Let us denote (Psi _i:=left( psi _{i,1},psi _{i,2},ldots ,psi _{i,m}right) ^{top }) and ({mathbf {r}}_i:=left( x_{i,1},x_{i,2},ldots ,x_{i,m}right) ^{top }), then

Within each patch, the proportions of species are affected by other patches only via migration. Let (phi _{i,j}) denote the rate of change of the proportion of species i in patch j in the absence of migration. Since the dominance relationships among the species (described by a tournament matrix ({mathbf {T}})) are assumed to be the same for all patches and since the habitat patches are spatially homogeneous, the expression for (phi _{i,j}) is given by the right-hand side of System (1):

where ({mathbf {p}}_j:=left( x_{1,j},x_{2,j},ldots , x_{n,j}right) ^{top }), (i=1,ldots ,n) and (j=1,ldots ,m). Assuming migration among the patches, the proportion of a species within a patch is influenced by two factors: the first is the interaction with other species within the patch and the second is the migration of that particular species to or from other patches. Thus, the metapopulation model describing the dynamics of the n species in the m-patch network is described by the system of mn differential equations;

This system evolves on the unit simplex (S^{mn}).

Proposition 2

The unit simplex (S^{mn}) is positively invariant for System (11).

Proof

To show the invariance of the unit simplex (S^{mn}) under the flow of System (11), it suffices to show that each of the faces of the simplex cannot be crossed, i.e., the vector field points inward from the faces of (S^{mn}).

On the one hand, if (x_{i,j}=0) for some i, j, then

which implies that (x_{i,j}=0) cannot be crossed from positive to negative. In an ecological context, this condition simply states the obvious fact that an extinct species is in no danger of declining. On the other hand, if (x_{i,j}=1) for some i, j, then obviously (x_{l,k}=0) for any (lne i) or (kne j) and

Hence, the vector field associated with System (11) points inward from the faces of (S^{mn}). So, (S^{mn}) is positively invariant under the flow of System (11). (square )

Note that Proposition 2 does not exclude the solution trajectories of System (11) from approaching the boundary equilibria of the system as (trightarrow infty ). We call metapopulation model (11) persistent if for every ({mathbf {x}}_0in S^{mn}_{+}), the (omega )-limit set (omega ({mathbf {x}}_0)) does not intersect the boundary of (S^{mn}). In other words, a metapopulation model is persistent if the initial existence of all the species implies that none of the species goes extinct with the passage of time.

Balanced homogeneous and heterogeneous metapopulation models

We say that the inter-patch migration of a metapopulation model is detailed balanced if the overall migration rate of any species between any two patches is zero for a certain positive set of proportions of that species in the different patches. From the theory of detailed-balanced reaction networks described in “Detailed-balanced single species mass action reaction networks” section, it follows that a detailed-balanced inter-patch migration network corresponds to a detailed-balanced single species mass action reaction network. Let B denote the incidence matrix corresponding to the directed graph (G_2) describing the inter-patch migrations and let r denote the number of edges in (G_2). Comparing Eqs. (8) and (9), it follows that if the inter-patch migration is detailed balanced, then there exist diagonal matrices ({mathcal {K}}in {mathbb {R}}^{rtimes r}) and ({mathbf {Z}}^*in {mathbb {R}}^{mtimes m}) with positive diagonal entries such that ((mathbb {1}^m)^{top }{mathbf {Z}}^*mathbb {1}^m=1) and

Let ({mathbf {Z}}^*=text { diag}({mathbf {z}}^*)). Equation (9) can now be rewritten as

Henceforth in this manuscript, we restrict our analysis to metapopulation models of type (11) for which the interactions within each patch correspond to a tournament with a completely mixed optimal strategy and whose inter-patch migration is detailed balanced. Such metapopulation models will be referred to as balanced metapopulation models.

We have seen earlier in “Species interactions and tournament matrices” section that if the interactions within every patch correspond to a tournament with a completely mixed optimal strategy, then the corresponding mean-field model admits a unique coexistence equilibrium ({mathbf {y}}^*in S^{n}_{+}) with ({mathbf {T}}{mathbf {y}}^*={mathbf {0}}^n). Thus, for a balanced metapopulation model, System (10) can be rewritten as

where ({mathbf {Y}}^*:=) diag(({mathbf {y}}^*)). Consequently, from Eqs. (11)–(13), it follows that the dynamics of a balanced metapopulation model containing n species and m patches can be described by mn differential equations

If all the elements of ({mathbf {z}}^*) in the above equation are equal, i.e., if (z_j^*=frac{1}{m}) for (j=1,ldots ,m), then we say that the balanced metapopulation model is homogeneous, otherwise we call it heterogeneous. Whether a balanced metapopulation model is homogeneous or not can be checked from the adjacency matrix ({mathbf {A}}) corresponding to its inter-patch migration graph (G_2). If ({mathbf {A}}) is symmetric, then the model is homogeneous, otherwise it is heterogeneous.

Remark 3

In³⁵, the authors assume that migrations from one patch to other patches are random with a probability of migration (or migration constant) equal to the reciprocal of the number of dispersal links from a patch to other patches. They thus define a dispersal graph to be homogeneous if all nodes have the same degree (number of links), otherwise the graph is heterogeneous. With this definition, homogeneity, in general, is equivalent to the existence of cycles in the dispersal graph, whereas heterogeneity is equivalent to their absence. However, with our new definition, it is clear that this is not necessary. An example of such a case is shown in Fig. 2.

Coexistence equilibrium and its uniqueness

In this section, we present a theorem that gives an expression for a coexistence equilibrium of a balanced metapopulation model. Before we state our main theorem in this section, we need the following lemma.

Lemma 4

Let ({mathbf {B}}in {mathbb {R}}^{mtimes r}) denote the incidence matrix of a finite connected directed graph (G_2) and let ({mathcal {K}}in {mathbb {R}}^{rtimes r}) denote a diagonal matrix with positive diagonal entries. For any ({mathbf {w}}in {mathbb {R}}_+^m), it holds that (-{mathbf {w}}^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{mathbb {1}^m}{{mathbf {w}}}right) ge 0). Moreover (-{mathbf {w}}^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{mathbb {1}^m}{{mathbf {w}}}right) = 0) if and only if ({mathbf {w}}=qmathbb {1}^m), where (qin {mathbb {R}}_+).

Proof

Assume that the (p{ {text {th}}}) edge of the graph (G_2) is directed from vertex (i_p) to vertex (j_p). Hence, (B_{i_pp}=-1), (B_{j_pp}=1) and (B_{kp}=0) for (i_pne kne j_p). Thus,

Moreover, (-{mathbf {w}}^{top }{mathbf {B}}{mathcal {K}}{mathbf {B}}^{top }left( frac{mathbb {1}^m}{{mathbf {w}}}right) =0) if and only if (w_{j_p}=w_{i_p}) for (p=1,ldots ,m), which is equivalent with ({mathbf {B}}^{top }{mathbf {w}}={mathbf {0}}^r).

Since the graph (G_2) is connected, we recall from⁴⁸ that (text {rank}({mathbf {B}})=m-1) and furthermore (text {ker}({mathbf {B}}^{top })=mathbb {1}^m). Therefore ({mathbf {B}}^{top }{mathbf {w}}={mathbf {0}}^r) if and only if ({mathbf {w}}=qmathbb {1}^m), where (qin {mathbb {R}}_+). This completes the proof. (square )

We now state the main theorem of this section.

Theorem 5

A balanced metapopulation model described by System (14) admits a unique coexistence equilibrium ({mathbf {x}}^*in S^{mn}_{+}). The proportion (x_{i,j}^{*}) of species i in patch j at the unique coexistence equilibrium is given by

for (i=1,ldots ,n) and (j=1,ldots ,m).

Proof

We divide the proof into two parts. In the first part we prove that System (15) indeed yields an equilibrium for the model. In the second part, we prove that this coexistence equilibrium is unique.

Let us define

For ({mathbf {x}}^*) to be an equilibrium of System (14), it should render the right-hand side equal to zero. Note that

and

In addition,

Thus, ({mathbf {x}}^*) is a coexistence equilibrium of System (14).

Assume that there exists another coexistence equilibrium ({mathbf {x}}^{**}in , S^{mn}_{+}). Let (x_{i,j}^{**}) denote the corresponding proportion of species i in patch j and define

It follows that for any i, j it holds that

Multiplying both sides of this equality with (frac{x_{i,j}^*}{x_{i,j}^{**}}), we get

Summing the left-hand side of the above expression over the different species and patches, we get

Now consider the two terms in the left-hand side of the above equality separately. For the first term, note that for any j it holds that

Hence,

For the second term, we find

Thus, Eq. (17) can be simplified as

Since (y_i^*>0) for (i=1,ldots ,n), it holds for any (i=1,ldots ,n) that

From Eq. (18) and Lemma 4, it follows that ({mathbf {r}}_{i}^{**}=q_i{mathbf {z}}^*) with (q_iin {mathbb {R}}_+) for (i=1,ldots ,n). Thus, (x_{i,j}^{**}=q_iz_{j}^*) and ({mathbf {p}}_{j}^{**}=z_j^*{mathbf {q}}) for (i=1,ldots ,n) and (j=1,ldots ,m). Substituting the latter in the left-hand side of Eq. (16), we get

Since (q_i>0) for (i=1,ldots ,n), for Eq. (16) to hold, we should have (mathbf {Tq}={mathbf {0}}^n). Also note that

Since the metapopulation model is balanced, it follows that ({mathbf {q}}={mathbf {y}}^*). Thus, (x_{i,j}^{**}=y_i^*z_j^*=x_{i,j}^*) for (i=1,ldots ,n) and (j=1,ldots ,m). This proves the uniqueness of the coexistence equilibrium ({mathbf {x}}^*). (square )

We now give examples of two balanced metapopulation models.

Example 1

It is easy to verify that the network shown in Fig. 3 corresponds to a balanced metapopulation model governed by System (14) with

({mathbf {y}}^*=left( frac{1}{3}, frac{1}{3}, frac{1}{3} right) ^{top }), ({mathbf {z}}^*=left( frac{1}{5}, frac{2}{5}, frac{2}{5} right) ^{top }) and ({mathcal {K}}=text { diag}left( frac{1}{10},frac{3}{10},frac{1}{10}right) ). Note that this metapopulation model is heterogeneous. From Theorem 5, it follows that the species proportions at the unique coexistence equilibrium for this model are given by (x_{i,1}^*=frac{1}{15}) and (x_{i,2}^*=x_{i,3}^*= frac{2}{15}).

Example 2

It is easy to verify that the network shown in Fig. 4 corresponds to a balanced metapopulation model governed by System (14) with

({mathbf {y}}^{*}={mathbf {z}}^*=left( frac{1}{3}, frac{1}{3}, frac{1}{3} right) ^{top }) and ({mathcal {K}}=frac{1}{3}text { diag}(mathbb {1}_2)). Note that this metapopulation model is homogeneous. From Theorem 5, it follows that the species proportions at the unique coexistence equilibrium in this case are all given by (x_{i,j}^*=frac{1}{9}) for (i,j= 1,2,3).

Stability

We now prove the local stability of the unique coexistence equilibrium corresponding to the balanced metapopulation model (14). For the proof, we make use of the same Lyapunov function as in “Neutral stability” section, coupled with LaSalle’s invariance principle⁵¹, (⁵², Section 4.2), (⁵³, pp. 188–189).

Theorem 6

Consider the balanced metapopulation model (14) with coexistence equilibrium ({mathbf {x}}^*in , S^{mn}_{+}).

1. 1.

   If the model is heterogeneous, then ({mathbf {x}}^*) is locally asymptotically stable w.r.t. all initial conditions in (S^{mn}_{+}) in the neighbourhood of ({mathbf {x}}^*). Furthermore, if the model is persistent, then ({mathbf {x}}^*) is globally asymptotically stable w.r.t. all initial conditions in (S^{mn}_{+}).

2. 2.

   If the model is homogeneous and persistent, then as (trightarrow infty ), the solution trajectories converge to a limit cycle satisfying the equation ({dot{x}}_{i,j}=x_{i,j}({mathbf {T}}{mathbf {p}}_{j})_i) with (x_{i,j}=x_{i,k}), for (i=1,ldots ,n) and (j,k=1,ldots ,m).

Proof

Let (x_{i,j}) denote the proportion of species i in patch j. Assuming that ({mathbf {x}}in S^{mn}_{+}), consider the Lyapunov function

By Gibbs inequality, V(x) is positive on (S^{mn}_{+}) and is equal to zero only if ({mathbf {x}}={mathbf {x}}^*). Taking the time derivative of V, we have

From Eq. (14), it follows that

As in the proof of Theorem 5, it can be verified that

and

Thus,

Since (y_i^*>0) for (i=1,ldots ,n), it follows from Lemma 4 that ({dot{V}}({mathbf {x}})le 0) and ({dot{V}}({mathbf {x}})=0) if and only if ({mathbf {r}}_i=q_i{mathbf {z}}^*) with (q_iin {mathbb {R}}_+), for (i=1,ldots ,n). Thus,

for (i= 1,ldots ,n) and (j=1,ldots ,m). Since ((mathbb {1}^{mn})^{top }{mathbf {x}}=1), we obtain

Let ({mathcal {E}}subset S^{mn}_{+}) be the set of all vectors ({mathbf {x}}) for which condition (20) is satisfied with ((mathbb {1}^n)^{top }{mathbf {q}}=1). We now determine the largest subset of ({mathcal {E}}) that is positively invariant w.r.t. System (14). Assume that ({mathbf {x}}) continuously takes values from ({mathcal {E}}) and satisfies System (14). Since ({mathbf {x}}) takes values from ({mathcal {E}}), we have ({dot{x}}_{i,j}=z_j^*{dot{q}}_i). Since ({mathbf {x}}) also satisfies System (14), we have

Thus, (z_j^*{dot{q}}_i=q_i{z_j^*}^2(mathbf {Tq})_i) which implies that

for (i=1,ldots ,n) and (j=1,ldots ,m). We now consider two cases.

Case 1: The model is heterogeneous, i.e., the vector ({mathbf {z}}^*) is not parallel to (mathbb {1}^m).

In this case, Eq. (21) will be satisfied only if (q_i(mathbf {Tq})_i=0) for (i=1,ldots ,n). Since (q_iin {mathbb {R}}_+) for (i=1,ldots ,n), it follows that (mathbf {Tq}={mathbf {0}}^n). Since ((mathbb {1}^n)^{top }{mathbf {q}}=1), we have ({mathbf {q}}={mathbf {y}}^*). This implies that (x_{i,j}=y_i^*z_j^*=x_{i,j}^*) for (i=1,ldots ,n) and (j= 1,ldots ,m). Thus, the largest subset of ({mathcal {E}}) that is positively invariant w.r.t. System (14) consists of just the unique equilibrium ({mathbf {x}}^*in S^{mn}_{+}). By LaSalle’s invariance principle, it follows that the equilibrium ({mathbf {x}}^*) is locally asymptotically stable w.r.t. all initial conditions in (S^{mn}_{+}) in the neighbourhood of ({mathbf {x}}^*), and globally asymptotically stable w.r.t. all initial conditions in (S^{mn}_{+}) provided that System (14) is persistent.

Case 2: The model is homogeneous, i.e.
({mathbf {z}}^*=frac{1}{m}mathbb {1}^m)

In this case, Eq. (21) takes the form ({dot{q}}_i=frac{q_i}{m}(mathbf {Tq})_i). We have (x_{i,j}=q_iz_j^*=frac{q_i}{m}) and

Consequently, the largest subset of ({mathcal {E}}) that is positively invariant w.r.t. System (14) consists of all vectors ({mathbf {x}}(t)in , S^{mn}_{+}) satisfying ({dot{x}}_{i,j}=x_{i,j}({mathbf {T}}{mathbf {p}}_{j})_i) with (x_{i,j}=x_{i,k}) for (i=1,ldots ,n) and (j,k=1,ldots ,m). The proof for Case 2 again follows from LaSalle’s invariance principle. (square )

The above results can be illustrated by simulating System (14) for the metapopulation models shown in Fig. 3 and 4 in Examples 1 and 2, respectively. The results of the simulations are shown in Figs. 5 and 6, respectively.

Source: Ecology - nature.com