Stability and steady state of complex cooperative systems: a diakoptic approach

TL;DR

This paper introduces a graph-theoretical diakoptic approach to analyze the stability of complex cooperative systems by decomposing their dependence graphs into strongly connected components and assessing eigenvalues.

Contribution

It provides a novel divide-and-conquer method for stability analysis of large cooperative systems using graph decomposition and eigenvalue criteria.

Findings

A linear cooperative system is Lyapunov stable if all SCCs have non-positive dominant eigenvalues.

Stability requires no path connecting SCCs with zero eigenvalues.

The approach simplifies stability analysis of complex ecological and population systems.

Abstract

Cooperative dynamics are common in ecology and population dynamics. However, their commonly high degree of complexity with a large number of coupled degrees of freedom renders them difficult to analyse. Here we present a graph-theoretical criterion, via a diakoptic approach (`divide-and-conquer') to determine a cooperative system's stability by decomposing the system's dependence graph into its strongly connected components (SCCs). In particular, we show that a linear cooperative system is Lyapunov stable if the SCCs of the associated dependence graph all have non-positive dominant eigenvalues, and if no SCCs which have dominant eigenvalue zero are connected by a path.

Equations

$$\frac{d}{dt}\mathbf{m}(t)=A\mathbf{m}(t)$$

$$A=[a_{ij}],\text{mbox}with{a}_{ij}\ge 0\text{mbox}fori\ne j.$$

$$\displaystyle A=\left(\begin{array}[]{ccccc}B_{1}&0&0&0&...\\ C_{21}&B_{2}&0&0&...\\ C_{31}&C_{32}&B_{3}&0&...\\ \vdots&\vdots&\vdots&\ddots&0\\ ...&...&...&...&B_{h}\end{array}\right),$$

$$p_A(\lambda )=p_{B_1}(\lambda )\cdot \dots \cdot p_{B_h}(\lambda ).$$

$$\displaystyle\left(\begin{array}[]{ccccc}B_{1}&0&0&0&...\\ C_{21}&B_{2}&0&0&...\\ C_{31}&C_{32}&B_{3}&0&...\\ \vdots&\vdots&\vdots&\ddots&0\\ ...&...&...&...&B_{h}\end{array}\right)\left(\begin{array}[]{c}\mathbf{m}^{*}_{1}\\ \mathbf{m}^{*}_{2}\\ \mathbf{m}^{*}_{3}\\ \vdots\\ \mathbf{m}^{*}_{h}\end{array}\right)=0,$$

$$\sum _{\begin{array}{c}l<k\end{array}}{C}_{kl}{\mathbf{m}}_l^{\ast }+B_k{\mathbf{m}}_k^{\ast }=0$$

$${\mathbf{m}}_k^{\ast }=-B_k^{-1}\left[\sum _{l<k}{C}_{kl}{\mathbf{m}}_l^{\ast }\right].$$

$${\mathbf{m}}_k^{\ast }=-B_k^{-1}\left[\sum _{l\in I_f}{P}_{kl}{\mathbf{m}}_l^{\ast }\right],$$

$$P_{kl}=\sum _{(l_1,l_2,...,l_n)\in {\mathcal{P}}_{kl}}(-1{)}^{n-1}{C}_{k{l}_1}{B}_{l_1}^{-1}{C}_{l_1{l}_2}{B}_{l_2}^{-1}\cdots C_{l_{n}l}$$

$$e^{B_{k}t}{B}_k{\mathbf{m}}_k^{\ast }=B_k\left[e^{B_{k}t}{\mathbf{m}}_k^{\ast }\right]=-e^{B_{k}t}\left[\sum _{l<k}{C}_{kl}{\mathbf{m}}_l^{\ast }\right],$$

$$J_i=\left(\begin{array}{cccc}{\lambda }_i& 1& & \\ & {\lambda }_i& \ddots & \\ & & \ddots & 1\\ & & & {\lambda }_i\end{array}\right)\text{ then }{f}_t(J_i)=\left(\begin{array}{cccc}f({\lambda }_i)& f_t^{\prime }({\lambda }_i)& \dots & \frac{f_t^{(m_i-1)}({\lambda }_i)}{(m_i-1)!}\\ & f_t({\lambda }_i)& \ddots & ⋮\\ & & \ddots & f_t^{\prime }({\lambda }_i)\\ & & & f_t({\lambda }_i)\end{array}\right).$$

$$\displaystyle\mathbf{m}^{*}_{k}=\left\{\begin{array}[]{ccc}0&\mbox{ if }&B_{k}\mbox{ is upstream of any critical block,}\\ \alpha_{k}\bm{\phi}_{k}&\mbox{ if }&k\in I_{f},\\ B^{-1}_{k}\left[\sum_{l\in I_{f}}P_{kl}\alpha_{l}\bm{\phi}_{l}\right]&\mbox{ if }&B_{k}\mbox{ is sub-critical. }\end{array}\right.$$

$$\displaystyle\mathbf{m}^{*}=\sum_{k\in I_{f}}\alpha_{k}\mathbf{m}^{*(k)}_{l},\mbox{ with }\mathbf{m}^{*(k)}_{l}=\left\{\begin{array}[]{ccc}0&\mbox{ if }&B_{l}\mbox{ is upstream of }B_{k},\\ 0&\mbox{ if }&l\in I_{f}\mbox{ and }l\neq k,\\ \bm{\phi}_{k}&\mbox{ for }&l\in I_{f},\\ B^{-1}_{l}\left[P_{lk}\bm{\phi}_{k}\right]&\mbox{ if }&B_{l}\mbox{ is sub-critical. }\end{array}\right.$$

Full text

1 Introduction

Cooperative systems are a wide class of dynamical systems characterised by a non-negative dependence between components [1]. Common examples are (bio-)chemical reaction networks with mutually activating interactions and compartmental dynamics, where a conserved quantity transits between different compartments or states [2, 3]. However, cooperative systems also include non-conserved replicator dynamics, such as (multi-species) population dynamics, where a population of replicators transits between different states/compartments. Examples for the latter are organisms which transit through life cycles or tissue cells (e.g. stem cells) which proliferate, switch between different phenotypes [4], and differentiate during biological development and in renewing tissues. If one considers the dynamics of sub-populations embedded in a larger population, then the equations describing the system are linear: while a population as a whole may be subject to a non-linear feedback (for example by a finite carrying capacity), smaller embedded sub-populations compete neutrally with each other without affecting the population as a whole. This renders the dynamics linear.

In this article we find conditions for the stability of linear cooperative systems, based on graphical criteria of the underlying dependence graph. In an ecological or biological context, stability of populations is required to maintain ecological equilibrium (population of individuals) or a functional biological tissue (population of tissue cells). In particular, instability of a tissue cell population may lead to cancer, thus the study of a cell population’s stability is of high biomedical importance. However, the commonly applied property of asymptotic stability is not viable for linear systems in a biological context, since the only asymptotically stable state is extinction. In these contexts it is therefore more appropriate to study marginally stable steady states, a form of Lyapunov stability.

While general criteria for a cooperative system’s stability are well established [5], real-world systems can be very complex, with a large number of variables and complex interactions, in which case their analysis is a highly challenging endeavour. Topological features of trajectories, such as compactness, can theoretically be used to determine stability [6, 7], but are in practice difficult to apply without explicitly solving the underlying differential equation. To simplify the analysis of a system, it is useful to represent it as a directed graph in which dependent variables $x_{i}(t),i\in\mathbb{N}$, are nodes and links denote dependence relations between those variables. The Jacobian matrix $J=[\frac{\partial\dot{x}_{i}}{\partial x_{j}}]$ of such a system can be interpreted as an adjacency matrix of an underlying graph representing the mutual dependence of components. Cooperative systems are then defined by non-negativity of the Jacobian’s off-diagonal entries, which corresponds to positive-only weights of links in the network. The corresponding Jacobian matrix is a Metzler matrix and thus methods based on non-negative matrices (and the Perron-Fobenius theorem) can be applied to study them [8, 9, 10]. A paradigm to study complex systems is the diakoptic view (“divide-and-conquer”) [11]: a large interacting system is decomposed into suitable small subsystems, which are studied in isolation, a task which is usually easier to perform. Then a synthesis of subsystems yields the features of the whole system. The analogy between dynamical systems and graphs may be used to apply graph-theoretical tools to perform such a diakoptic decomposition of the system (i.e. graph) into smaller sub-systems which can significantly simplify the analysis of stability features of the respective system.

A diakoptic approach, based on the decomposition of the underlying graph into its strongly connected components (SCC) has been used to determine asymptotic stability of cooperative systems [12]. An SCC is a subset of nodes which are all mutually reachable by directed paths. Simply speaking, a system is asymptotically stable if, and only if, all its SCCs, when decoupled from each other, are asymptotically stable. This holds, since the eigenvalues of the system’s adjacency matrix are the union of the SCCs’ eigenvalues [12, 9], which can also easily be checked by evaluating row sums of the dynamical matrix [13]. However, this criterion cannot be straight-forwardly generalised to determine marginal stability; conclusions about marginal stability vs. instability can in general not be drawn just by considering SCCs in isolation and a system may be unstable even if there is no unstable individual SCC. Only for linear compartmental systems – i.e. cooperative systems that feature a conserved quantity – the existence of a marginally state can be found through analysis of SCCs in isolation: if there is at least one singular SCC, a so-called trap, then a non-trivial marginally stable steady state exists [14, 15]. However, this criterion cannot be applied to linear cooperative systems in general, when dynamics are not conserved.

Here, we introduce a diakoptic approach to determine the stability of general linear cooperative systems, which is also applicable for non-conserved systems, and allows us to identify conditions for marginal stability. This approach is based on graphical criteria of the underlying dependence graph when decomposed for its SCCs. The stability can then be inferred from (i) the spectrum of the Jacobian matrices of isolated SCCs, and (ii) the hierarchical arrangement of the SCCs. Our main result is Theorem 3 (illustrated in Fig. 2) which states that for (marginal) stability to prevail, no SCC may have positive eigenvalues, and any otherwise singular SCCs may not stand in any hierarchical relation to each other, i.e. there may be no (directed) path connecting them. This reflects the principle that the larger and more connected complex systems are, the more likely they are to become unstable [16].

2 Results

We consider a generic cooperative linear dynamical system of a positive quantity (‘mass’) $m$ on a directed weighted graph with $n$ nodes, whereby we denote $m_{i}=m_{i}(t)$ as the mass on node $i=1,...,n$ at time $t$. The state vector of the system is $\mathbf{m}=(m_{1},m_{2},...,m_{n})^{T}$ and the system is written as

[TABLE]

for a $n\times n$ real square matrix

[TABLE]

The condition $a_{ij}\geq 0$ for $i\neq j$, defines the system as cooperative, since $A$ is the Jacobian matrix of the system (1). We note that the system is not necessarily conserved, i.e. the ‘mass’ could replicate, such as a population of biological individuals, cells, or viruses.

We consider the underlying directed weighted graph $G(A)$ with transposed adjacency matrix $A$, that is, the graph with $n$ nodes and a link from $j$ to $i$, weighted by $a_{ij}$, only if $a_{ij}\neq 0$. This is a finite simple graph with positively weighted edges and arbitrarily weighted self-loops. We wish to relate the stability of the fixed points of (1) to the network structure of $G(A)$.

Since $A$ is the Jacobian of (1), the stability of a fixed point $\mathbf{m}^{*}$, defined by $A\,\mathbf{m}^{*}=\mathbf{0}$ (that is, a 0-eigenvector of $A$), is determined by the spectral properties of $A$. For the system to be asymptotically stable, all the real parts of the eigenvalues of $A$ must be negative. In this case, however, $\det(A)\neq 0$ and the only fixed point $A\,\mathbf{m}^{\ast}=\mathbf{0}$ is trivial, $\mathbf{m}^{\ast}=\mathbf{0}$. As we are interested in non-trivial solutions, we focus instead on Lyapunov stable fixed points which are at least marginally stable (also called semi-stable [5]). This is the case if the eigenvalue of $A$ with largest real part is zero and its geometric multiplicity is equal to its algebraic multiplicity [17]. Our main result is a necessary and sufficient condition on the structure of the graph $G(A)$, for the dynamical system to have non-trivial, marginally stable, non-negative solutions. Note that we call a vector $\mathbf{m}$ (or, similarly, a matrix) non-negative, written $\mathbf{m}\geq 0$, if all entries are real and non-negative, and positive, written $\mathbf{m}>0$, if all entries are real and positive.

First, we decompose $G(A)$ into its strongly connected components, as follows. A (sub-)graph is strongly connected if for any pair of nodes $i$ and $j$ in the graph there is a directed path from $i$ to $j$ and a directed path from $j$ to $i$, that is, every pair of nodes is mutually reachable. Every directed graph can be partitioned into maximal strongly connected subgraphs, the graph’s strongly connected components (SCCs). The SCCs of a directed graph $G$ form another graph called the condensation of $G$: in it, each node represents an SCC, and if two SCCs in $G$ are connected by at least one link, then the condensation possesses a link between them, in the same direction as in $G$ (see Fig. 1). The condensation of a directed graph is always a directed acyclic graph and, hence, its nodes (the SCCs of $G$) admit a topological ordering [18]: an ordering $B_{1},B_{2},\ldots,B_{h}$ (from now on, we will identify the $k$th connected component of $G$ with its adjacency matrix $B_{k}$) such that if there is a link from $B_{i}$ to $B_{j}$ then $i\leq j$ (see Fig. 1 for an example). We can extend the ordering to the nodes of $G$ so that node $u\in B_{i}$ appears before node $v\in B_{j}$ whenever $i\leq j$. With respect to this this re-ordering and re-labelling of the nodes of $G$, the adjacency matrix $A$ of $G(A)$ becomes a lower triangular matrix

[TABLE]

where $h$ is the number of SCCs of $G(A)$, $B_{k}$ is the adjacency matrix of the $k$-th SCC ($1\leq k\leq h$), and $C_{kl}$ encodes the connectivity from $B_{l}$ to $B_{k}$. This is sometimes called the normal form of a reducible matrix [19]. If there exist a path from $k$ to $l$ (thus $k\leq l$), we call $B_{k}$ upstream of $B_{l}$, and $B_{l}$ is downstream of $B_{k}$. If $B_{k}$ is connected by a single (directed) link to $B_{l}$ then we also call $B_{k}$ immediately upstream of $B_{l}$, and $B_{l}$ immediately downstream of $B_{k}$. From now on, we will implicitly assume a topological ordering and notation as above.

Since $A$, written in the form (8), is a lower triangular block matrix, the characteristic polynomial of $A$, $p_{A}(\lambda)=\det(\lambda I-A)$, is the product of the characteristic polynomials of the $B_{k}$’s

[TABLE]

Thus the spectrum of $A$ – seen as a multiset – is the union of the spectra of the $B_{k}$’s, and the algebraic multiplicity of the eigenvalues is preserved.

Since all off-diagonal elements of $A$ (and hence of each $B_{k}$) are non-negative, and each $B_{k}$ is the adjacency matrix of a strongly connected graph, the matrices $B_{k}$ are irreducible Metzler matrices, for which the Perron-Frobenius Theorem applies, to the shifted eigenvalues [20]. Therefore, each matrix $B_{k}$ has a real eigenvalue $\mu_{k}$ with (strictly) largest real part, which is simple and has a positive eigenvector $\mathbf{m}>0$. We call $\mu_{k}$ the dominant eigenvalue of the matrix $B_{k}$.

We now introduce some further terminology. We call each SCC, and equivalently its adjacency matrix $B_{k}$, a block of the system (we use the term ‘block’ and the notation $B_{k}$ for both the matrix and its graph). We call a block critical if its dominant eigenvalue $\mu_{k}=0$, sub-critical if $\mu_{k}<0$ and super-critical if $\mu_{k}>0$. Correspondingly, we define the index subsets $I_{c}=\{k\in 1,...,h\mid B_{k}\mbox{ critical}\}$, $I_{s}=\{k\in 1,...,h\mid B_{k}\mbox{ sub-critical }\}$, $I_{sp}=\{k\in 1,...,h\mid B_{k}\mbox{ super-critical }\}$. The first things we note are (see for example [12]) are

Lemma 1**.**

If at least one block $B_{k}$ of $A$ is super-critical, then the system (1) is unstable.

Lemma 2**.**

The system (1) is asymptotically stable if and only if all blocks $B_{k}$ of $A$ are sub-critical.

These Lemmas follow immediately from the the fact that a system is unstable if at least one real part of an eigenvalue of $A$ is positive and it is asymptotically stable if and only if all real parts of eigenvalues are negative, together with the property that the spectrum of $A$ is the multi-set union of spectra of the $B_{k}$ (note, however, that the ‘if and only if’ statement only holds for Lemma 2) [12]. In the situation of Lemma 2, observe that $\det(A)\neq 0$ and hence $\mathbf{m}^{\ast}=\mathbf{0}$ is the only fixed point of the system.

Lemmas 1 and 2 cover all cases where any super-critical blocks exist, or only sub-critical ones. In these cases, the system is either unstable, or has only a trivial (zero) fixed point. In the following, we will consider only the remaining cases when no super-critical blocks exist, but there is at least one critical block, and investigate the existence of non-trivial, non-negative (so that each node supports a non-negative fraction of the ‘mass’) marginally stable fixed points.

If no super-critical, and at least one critical, block exists, the dominant eigenvalue of $A$ is zero and, according to the Perron-Frobenius theorem, there exist non-trivial eigenvectors $\mathbf{m}^{*}$ for the eigenvalue zero. It is assured that all such $\mathbf{m}^{*}$ are equilibrium points of the system (1), however, to be a (Lyapunov) stable equilibrium it is required that the algebraic multiplicity of eigenvalue zero is equal to its geometric one, or equivalently, equal to the dimension to the nullspace of $A$. We will approach the latter question by explicitly constructing such equilibrium sets.

Let us first write the equilibrium condition of the dynamical system (1), using (8), as

[TABLE]

i.e. the equilibrium vector $\mathbf{m}^{*}$ is decomposed in the projections $\mathbf{m}^{*}_{k}$ on the sub-space of $B_{k}$, in the form $\mathbf{m}^{*}=(\mathbf{m}_{1},\mathbf{m}_{2},...,\mathbf{m}_{h})^{T}$. For simplicity, we call $\mathbf{m}^{*}_{k}$ the steady state on $B_{k}$. We further call a block $B_{k}$ trivial if $\mathbf{m}^{\ast}_{k}=\mathbf{0}$ for all non-negative marginally stable fixed points $\mathbf{m}^{*}$ of the system (20), and non-trivial otherwise111Note that $\mathbf{m}^{*}_{k}$ is the $k$-th subspace component of the global steady state $\mathbf{m}^{*}$ of $A$, but not necessarily the steady state of the isolated subsystem of $B_{k}$.. In other words, a trivial block is one that does not support any positive fraction of the ‘mass’ for any non-negative fixed point.

Our first result is a formula for the steady states $\mathbf{m}^{*}_{k}$ on sub-critical blocks $B_{k}$. Let us consider the $k$-th row of (20),

[TABLE]

where $B_{k}$ is a sub-critical block. Since all eigenvalue real parts of $B_{k}$ are negative, $\det(B_{k})\neq 0$ and thus $B_{k}$ is invertible, so that we obtain a recursive formula for the steady state:

[TABLE]

Let $I_{f}\subseteq I_{c}$ denote the indices of the critical SCCs for which there are no other critical SCCs downstream. Let us call them final critical blocks. With this terminology, we have, from the recursion relation above, the following:

Theorem 1**.**

If $B_{k}$ is a sub-critical block of $A$ and Eq. (20) holds, then $\mathbf{m}^{*}_{k}$, the steady state on $B_{k}$, is uniquely determined by the final critical blocks upstream of $B_{k}$, namely

[TABLE]

where

[TABLE]

and $\mathcal{P}_{kl}$ is the set of all paths from $B_{l}$ ($l\in I_{f}$) to $B_{k}$, written as a sequence of nodes $(l_{1},l_{2},...,l_{n})$, where $n$ is the length of the path.

This follows directly if we apply the relation Eq. (22) recursively to all steady states $\mathbf{m}^{*}_{k}$ of sub-critical blocks $B_{k}$ on the right hand side of Eq. (22), using that when propagating upstream, no critical block can be encountered before a final critical block is encountered.

Theorem 1 assures that the steady state on any sub-critical block is uniquely defined by the steady states on all critical blocks upstream of the former. Furthermore, we can conclude:

Corollary 1**.**

If $B_{k}$ is a sub-critical block of $A$, then $B_{k}$ is trivial if and only if all $B_{l}$ immediately upstream of $B_{k}$ are trivial.

Proof.

From Eq. (22) it directly follows that if all $B_{l}$ immediately upstream are trivial ($\mathbf{m}^{*}_{l}=0$), then $B_{k}$ is trivial ($\mathbf{m}^{*}_{k}=0$). Now let us consider the case that at least one $B_{l}$ immediately upstream has $\mathbf{m}^{*}_{l}\neq 0$. We first note that since $B_{k}$ is a Metzler matrix with $\det(B_{k})\neq 0$, $-B_{k}$ is a non-singular M-matrix, and its inverse is a positive matrix (shown in [21]). Thus $\mathbf{m}^{*}_{k}$ is positive if at least one $\mathbf{m}^{*}_{l}\neq 0$ (recall that $\mathbf{m}^{\ast}$ and the $C_{kl}$’s are non-negative). Therefore it follows: if $B_{k}$ is trivial, i.e. $\mathbf{m}^{*}_{k}=0$, then for all immediately upstream $B_{l}$, $\mathbf{m}_{l}=0$, and hence $B_{l}$ is trivial. ∎

Now we make a topological characterisation of the trivial blocks.

Theorem 2**.**

A block is trivial if and only if

• (i)

it is upstream of a critical block, or

• (ii)

it is a sub-critical block which is not downstream of a critical block.

Thereby all trivial blocks can be easily identified by inspecting the condensed graph and its critical blocks (Fig. 1).

Proof.

To prove Theorem 2, consider an equilibrium point $\mathbf{m}^{*}$. Then Eq. (20) holds, and in particular its $k$-th row Eq. (21). Now $B_{k}$ is critical and hence has an eigenvalue zero ($\mu_{k}=0$ by definition), $\det(B_{k})=0$ and thus $B_{k}$ is not invertible. Let us multiply both sides of Eq. (21) with the matrix exponential $e^{B_{k}t}=\sum_{n=1}^{\infty}\frac{(B_{k}t)^{n}}{n!}$ to yield,

[TABLE]

where we used that a square matrix $M$ commutes with its exponential, $e^{M}M=Me^{M}$. In general, $e^{Mt}\mathbf{x}$ is a solution of the linear ODE $\dot{\mathbf{x}}=M\mathbf{x}$ and thus converges to a linear combination of dominant eigenvectors (eigenvectors of the dominant eigenvalues) of $M$. Since $B_{k}$ is critical, the corresponding dominant eigenvalue is zero and thus $B_{k}[e^{B_{k}t}\mathbf{m}^{*}_{k}]\to\mathbf{0}$ for $t\to\infty$. This means that $-e^{B_{k}t}[\sum_{l<k}C_{kl}\mathbf{m}^{*}_{l}]=\mathbf{0}$ for $t\to\infty$. Let us call $L=\lim_{t\to\infty}e^{B_{k}t}$ and $\mathbf{v}=\sum_{l<k}C_{kl}\mathbf{m}^{*}_{l}$. We then have $L\mathbf{v}=\mathbf{0}$ and we want to show that $\mathbf{v}=\mathbf{0}$. This is not true for a general vector $\mathbf{v}$ unless $L$ is invertible, but will hold for non-negative eigenvectors such as $\mathbf{v}$. In fact, the matrix $L$ is not invertible in general as all its eigenvalues are zero, except a simple eigenvalue 1 with a positive (left) eigenvector $\mathbf{u}$ (see Lemma 4 below). In this case, $\mathbf{u}L=\mathbf{u}$ and $\mathbf{u}\mathbf{v}=\mathbf{u}L\mathbf{v}=\mathbf{0}$, a contradiction, since $\mathbf{u}$ is positive and $\mathbf{v}$ is non-negative, unless $\mathbf{v}=\mathbf{0}$.

All in all, we conclude that $\mathbf{v}=\sum_{l<k}C_{kl}\mathbf{m}^{*}_{l}=0$. Note that all entries of the matrices $C_{kl}$ and of the vector $\mathbf{m}^{*}_{l}$ are non-negative, so this can only be the case if, for all $l<k$, $C_{kl}=0$ or $\mathbf{m}^{*}_{l}=\mathbf{0}$. Since, for all immediately upstream blocks, we have $C_{kl}\neq 0$, it follows that

Lemma 3**.**

All blocks immediately upstream of a critical block are trivial.

Crucially, from Theorem 1 and Lemma 1, it follows that all blocks $B_{m}$ immediately upstream of any trivial block $B_{l}$ are trivial (either $B_{l}$ is critical, or sub-critical and trivial). By applying this argument recursively to Lemma 1, the first part of Theorem 2 follows. The second part is an immediate consequence of Corollary 1. ∎

To complete the proof above, we state and prove the following.

Lemma 4**.**

Let $B$ be an irreducible Metzler matrix with shifted Perron-Frobenius eigenvalue 0 and positive left eigenvector $\mathbf{u}$. Then the limit matrix $L=\lim_{t\to\infty}e^{tB}$ exists and satisfies $\mathbf{u}L=\mathbf{u}$.

Proof.

Define $f_{t}(x)=e^{tx}$ for $t>0$ and write $B$ in Jordan normal form as $B=PJP^{-1}$. By definition [22], the matrix function $f_{t}(B)=e^{tB}$ equals $Pf(J)P^{-1}$, where $f(J)$ is the block diagonal matrix obtained by applying $f$ to each diagonal Jordan block, $J_{i}$, of $J$ as

[TABLE]

The eigenvalue of $B$ with largest real part is [math] (dominant eigenvalue), hence $\lim_{t\to\infty}f_{t}^{(k)}(\lambda_{i})=\lim_{t\to\infty}t^{k}e^{t\lambda_{i}}=0$, if $\lambda_{i}\neq 0$, and $1$ otherwise, for all $k\geq 0$. All in all, $L=PMP^{-1}$ where $M$ is the zero matrix except a single 1 in the diagonal. Its eigenvectors (the columns of $P$) are the same as the (generalised) eigenvectors of $B=PJP^{-1}$. Hence the left 0-eigenvector $\mathbf{u}$ of $B$ becomes a left 1-eigenvector of $L$, $\mathbf{u}L=\mathbf{u}$. ∎

We can also easily follow from Theorem 2 and Eq. (25):

Corollary 2**.**

The steady states $\mathbf{m}^{*}_{k}$ on a non-trivial critical block $B_{k}$ (called a free block) is the one-dimensional family of dominant eigenvectors (of eigenvalue zero) of $B_{k}$. We can write these as $\alpha_{k}\bm{\phi}_{k}$ where $\alpha_{k}\in\mathbb{R}$ is a free parameter, and $\bm{\phi}_{k}$ is a (normalised) dominant eigenvector of $B_{k}$.

Theorems 1 and 2, and Corollary 2, allow us to construct the most generic steady state of the system (1), that is, the nullspace of $A$. From Theorem 2, it follows that the set of non-trivial SCCs is exactly the set of final critical blocks, as defined before Theorem 1. Hence $I_{f}\subseteq I_{c}$ is also the index set of non-trivial critical blocks, that is, $I_{f}=\{1\leq k\leq h\mid B_{k}\mbox{ critical and non-trivial}\}$. All in all, a steady state vector $\mathbf{m}^{*}=(\mathbf{m}_{1},\mathbf{m}_{2},...,\mathbf{m}_{h})^{T}$ has the form

[TABLE]

This can also be written as

[TABLE]

Hence, the dimension of the nullspace of $A$ (equivalently, the geometric multiplicity of the eigenvalue zero) is equal to the number of non-trivial critical blocks. The algebraic multiplicity, on the other hand, is the number of all critical blocks, since according to the Perron-Frobenius Theorem for Metzler matrices, each critical block has a simple (algebraic multiplicity 1) eigenvalue zero and thus contributes once to the multiset of eigenvalues of $A$, by Eq. (9). Recall that the system (1) is marginally stable if and only if the dominant eigenvalue of $A$ is zero and its geometric multiplicity is equal to the algebraic multiplicity; this is thus the case only if there are no super-critical blocks and all critical blocks are non-trivial. According to Theorem 2, this is the case if no critical block is upstream of another critical block, or alternatively, if there are no paths between any two critical blocks. We thereby arrive at

Theorem 3**.**

A dynamical system in the form of Eq. (1) with Jacobian matrix $A$ is marginally stable if these conditions hold:

• (a)

there are no super-critical blocks;

• (b)

there is at least one critical block;

• (c)

there are no (directed) paths in $G(A)$ which connect two critical blocks.

Part (a) and (b) follow from Lemmas 1 and 2, while part (c) follows from Eq. (34) and Theorem 2, that if there is a path between two critical blocks, one of them must be trivial. This theorem is illustrated in Fig. 2. To summarise our findings, including the necessary definitions, we can express the stability criteria of cooperative dynamical systems as

Theorem 4**.**

Let $A=[a_{ij}]$, with $a_{ij}\geq 0$ if $i\neq j$, be the Jacobian matrix of the linear cooperative system in Eq. (1), and $G(A)$ its weighted graph, defined by the edge weights $a_{ij}$ for all $i,j$. Let $B_{1},\ldots,B_{h}$ be the adjacency matrices of the strongly connected components (SCCs) of $G(A)$. We define an SCC as critical if its dominant eigenvalue is zero, sub-critical, if its dominant eigenvalue is negative, and super-critical if its dominant eigenvalue is positive. Then:

The system is asymptotically stable if and only if all SCCs are sub-critical. In that case the steady state vanishes (is the zero-vector). 2. 2.

Otherwise, the system is marginally stable if

• (a)

there are no super-critical SCCs, and

• (b)

there are no paths in $G(A)$ which connect two critical SCCs. 3. 3.

Otherwise, the system is unstable.

The corresponding equilibrium set of the system is given by Eq. (29).

3 Conclusions

The conditions stated in Theorem 4 prescribe a way to simplify the analysis of a high-dimensional linear cooperative system by decomposition into lower dimensional subsystems, the strongly connected components (SCCs) of the dynamical system’s dependence graph. By spectral analysis of these SCCs and checking whether the topological conditions of Theorem 4 are fulfilled, the system’s stability can be determined. In particular, marginal stability is of importance for linear systems, since marginally stable states represent the only possible non-vanishing stable states – i.e. not identical to the zero-vector – called steady states. Such a steady state features conservation of the quantity of interest “on average”, i.e. the mean value stays constant even if the quantity itself is not strictly conserved222We note that since a marginal steady state $\bm{m}^{*}$ is a right 0-eigenvector, by fulfilling $A\bm{m}^{*}=0$, there must also exist a left 0-eigenvector $\bm{c}$, fulfilling $\bm{c}A=0$. The latter equation defines a generalised conservation law with coefficients $\bm{c}$; however, this conservation law may be non-trivial and cannot be directly derived from $\bm{m}^{*}$.. In contrast, asymptotically stable states are trivially vanishing for linear systems.

Moreover, our analysis revealed that a critical SCC, i.e. one with dominant eigenvalue zero, uniquely determines the steady state of the (necessarily sub-critical) SCCs upstream and downstream of it. In particular, the steady state configurations of all SCCs downstream of a critical SCC do generally not vanish (Theorem 1), and are uniquely determined by Eqs. (23, 24), while the steady state configurations of all SCCs upstream of it must vanish (Theorem 2). This leads to an explicit formula (Eqs. (29,34)) to construct the steady state of the whole system by the knowledge of the steady states on the critical SCCs only.

The results, Theorems 1-4, can be seen as a generalisation of a similar condition found for linear compartmental systems, i.e linear cooperative systems where the quantity of interest is strictly conserved (apart from external sources and sinks) [2]. For those systems, it has been found that the existence of at least one singular SCC (having an eigenvalue zero) is sufficient to ensure a non-trivial steady state [14, 15]. Notably, due to the conservation law in compartmental systems, no SCCs with positive eigenvalue may exist, meaning that all singular SCCs are critical, and furthermore, no critical SCCs can have any outgoing links. Hence, the existence of a singular SCC in a compartmental system automatically implies the conditions of our Theorem 3. The stability conditions of Theorems 3 and 4 therefore represent a generalisation to cooperative systems where the quantity of interest is not necessarily conserved, and can therefore also be applied to population dynamics where individuals can replicate and transit between different states. An example are populations of stem cells in animal tissues which differentiate, thereby changing their cell type. Note that while (cell) populations as a whole are often subject to feedback and thus follow non-linear dynamics, when considering subpopulations therein, which compete neutrally, the corresponding subsystem is linear.

In general, cooperative systems can be highly complex, with a large number of variables and very complex interactions, hence represented by large and often irregular graphs. The method presented here is a way to significantly simplify the analysis of a wide range of systems, ranging from cooperative (bio-)chemical reactions to complex population dynamics, by decomposing the systems into its strongly connected components. We have shown that a spectral analysis of each SCC, and a simple graphical criterion of the connectivity between SCCs (Theorem 4, Figure 2) completely determine the stability of any linear cooperative system. This provides a unique insight into the possible configurations of cooperative systems and demonstrates the power of graph theoretic techniques in the analysis of complex dynamical systems.

Funding Statement

PG was supported by a Medical Research Council New Investigator Research Grant MR/R026610/1.

Author contributions

PG and RJSG carried out the mathematical analysis, PG, BDM, and CP conceived and designed the project. All authors helped draft the manuscript. All authors gave final approval for publication.

Acknowledgements

We thank David Chillingworth for hinting us towards some literature on cooperative systems.