Dynamic Ecological System Analysis

TL;DR

This paper introduces a comprehensive mathematical framework for analyzing nonlinear ecological systems dynamically, enabling detailed tracking of flows, storages, and interactions within complex ecological networks.

Contribution

It develops a novel system decomposition and flow characterization method that extends static network analysis to nonlinear dynamic ecological systems.

Findings

Analytical characterization of flow paths and interactions.

Quantitative classification of interspecific interactions.

Extension of static network concepts to dynamic settings.

Abstract

This article develops a new mathematical method for holistic analysis of nonlinear dynamic compartmental systems through the system decomposition theory. The method is based on the novel dynamic system and subsystem partitioning methodologies through which compartmental systems are decomposed to the utmost level. The dynamic system and subsystem partitioning enable tracking the evolution of the initial stocks, environmental inputs, and intercompartmental system flows, as well as the associated storages derived from these stocks, inputs, and flows individually and separately within the system. Moreover, the transient and the dynamic direct, indirect, acyclic, cycling, and transfer (diact) flows and associated storages transmitted along a given flow path or from one compartment, directly or indirectly, to any other are analytically characterized, systematically classified, and mathematically formulated. Further, the article develops a dynamic technique based on the diact transactions for the quantitative classification of interspecific interactions and the determination of their strength within food webs. Major concepts and quantities of the current static network analyses are also extended to nonlinear dynamic settings and integrated with the proposed dynamic measures and indices within the proposed unifying mathematical framework. Therefore, the proposed methodology enables a holistic view and analysis of ecological systems. We consider that this methodology brings a novel complex system theory to the service of urgent and challenging environmental problems of the day and has the potential to lead the way to a more formalistic ecological science.

Tables

Table 1. Table 1 : The dynamic diact flow distribution and the simple and composite diact (sub)flow matrices. The superscript ( * * {}^{\texttt{*}} ) in each equation represents any of the diact symbols. For the sake of readability, the function arguments are dropped.

diact	flow distribution matrix	flows
d	$N^{\text{d}}=F{𝒯}^{-1}$
i	$N^{\text{i}}=\tilde{T}{\widehat{𝖳}}^{-1}-F{𝒯}^{-1}$
a	$N^{\text{a}}=\tilde{T}{\widehat{𝖳}}^{-1}-\widetilde{𝖳}{\widehat{𝖳}}^{-1}\hat{T}{\widehat{𝖳}}^{-1}$
c	$N^{\text{c}}=\widetilde{𝖳}{\widehat{𝖳}}^{-1}\hat{T}{\widehat{𝖳}}^{-1}$
t	$N^{\text{t}}=\tilde{T}{\widehat{𝖳}}^{-1}$

Table 2. Table 2 : The input-oriented, flow-based diact flow and storage distribution and the simple and composite diact (sub)flow and (sub)storage matrices. The superscript ( * * {}^{\texttt{*}} ) in each equation represents any of the diact symbols.

diact	flow and storage distribution matrices		flows	storages
d	$N^{\text{d}}=F{𝒯}^{-1}$	$S^{\text{*}}=ℛN^{\text{*}}$	$T^{\text{*}}=N^{\text{*}}𝒯$ $T_{\ell }^{\text{*}}=N^{\text{*}}{𝒯}_{\ell }$ ${\tilde{T}}^{\text{*}}=N^{\text{*}}𝖳$	$X^{\text{*}}=S^{\text{*}}𝒯$ $X_{\ell }^{\text{*}}=S^{\text{*}}{𝒯}_{\ell }$ ${\tilde{X}}^{\text{*}}=S^{\text{*}}𝖳$
i	$N^{\text{i}}=\left(N-I\right){𝒩}^{-1}-F{𝒯}^{-1}$
a	$N^{\text{a}}=\left({𝒩}^{-1}N-I\right){𝒩}^{-1}$
c	$N^{\text{c}}=\left(N-{𝒩}^{-1}N\right){𝒩}^{-1}$
t	$N^{\text{t}}=\left(N-I\right){𝒩}^{-1}$

Equations

$${\dot{x}}_i(t)={\check{\tau }}_i(t,x)-{\hat{\tau }}_i(t,x)$$

$${\check{\tau }}_i(t,x)=\sum _{\begin{array}{c}j=0\end{array}}^n{f}_{ij}(t,x)\text{mbox}and{\hat{\tau }}_i(t,x)=\sum _{\begin{array}{c}j=0\end{array}}^n{f}_{ji}(t,x)$$

$$f_{ij}(t,x)=q_{ij}^x(t,x)x_j(t)$$

$${\dot{x}}_i(t)=\left(z_i(t,x)+\sum _{\begin{array}{c}j=1\end{array}}^n{f}_{ij}(t,x)\right)-\left(y_i(t,x)+\sum _{\begin{array}{c}j=1\end{array}}^n{f}_{ji}(t,x)\right)$$

$$\sum _{i=1}^n{\dot{x}}_i(t)=0\text{mbox}whenz(t,x)=y(t,x)=0\text{mbox}forallt$$

$$F(t,x)=\left(f_{ij}(t,x)\right)$$

$$\check{\tau }(t,x)$$

$$\hat{\tau }(t,x)$$

$$z(t,x)=[z_1(t,x),\dots ,z_n(t,x){]}^T\text{mbox}andy(t,x)=[y_1(t,x),\dots ,y_n(t,x){]}^T$$

$$x_i(t)=\sum _{k=0}^n{x}_{i_k}(t),i=1,\dots ,n.$$

$$\mathrm{x}(t)={\left[x_{1_0}(t),\dots ,x_{n_0}(t),x_{1_1}(t),\dots ,x_{n_1}(t),\dots ,x_{1_n}(t),\dots ,x_{n_n}(t)\right]}^T.$$

$$z_{i_k}(t,\mathrm{x})={\delta }_{ik}{z}_k(t,x)=\{\begin{array}{rl}& z_{k_k}(t,\mathrm{x})=z_k(t,x),i=k\\ & 0.i\ne k\end{array}$$

$$f_{i_k{j}_k}(t,\mathrm{x})=x_{j_k}(t)\frac{f_{ij}(t,x)}{x_j(t)}=x_{j_k}(t)q_{ij}^x(t,x)=d_{j_k}(\mathrm{x})f_{ij}(t,x)$$

$$f_{ij}(t,x)=\sum _{k=0}^n{f}_{i_k{j}_k}(t,\mathrm{x}),i,j=1,\dots ,n.$$

$${\dot{x}}_{i_k}(t)=\left(z_{i_k}(t,\mathrm{x})+\sum _{\begin{array}{c}j=1\end{array}}^n{f}_{i_k{j}_k}(t,\mathrm{x})\right)-\left(y_{i_k}(t,\mathrm{x})+\sum _{\begin{array}{c}j=1\end{array}}^n{f}_{j_k{i}_k}(t,\mathrm{x})\right)$$

$$x_{i_k}(t_0)=\{\begin{array}{r}{x}_{i,0},k=0\\ 0.k\ne 0\end{array}$$

$${\check{\tau }}_{i_k}(t,\mathrm{x})=z_{i_k}(t,\mathrm{x})+\sum _{\begin{array}{c}j=1\end{array}}^n{f}_{i_k{j}_k}(t,\mathrm{x})\text{mbox}and{\hat{\tau }}_{i_k}(t,\mathrm{x})=y_{i_k}(t,\mathrm{x})+\sum _{\begin{array}{c}j=1\end{array}}^n{f}_{j_k{i}_k}(t,\mathrm{x})$$

$$X_{ik}(t)=x_{i_k}(t),{\check{T}}_{ik}(t,\mathrm{x})={\check{\tau }}_{i_k}(t,\mathrm{x}),\text{mbox}and{\hat{T}}_{ik}(t,\mathrm{x})={\hat{\tau }}_{i_k}(t,\mathrm{x}),$$

$$x_0(t)$$

$${\check{\tau }}_0(t,\mathrm{x})$$

$$\mathcal{X}(t)=\mathrm{diag}(x(t)),\mathcal{Y}(t,x)=\mathrm{diag}(y(t,x)),\text{mbox}and\mathcal{Z}(t,x)=\mathrm{diag}(z(t,x)).$$

$$\check{T}(t,\mathrm{x})$$

$$\hat{T}(t,\mathrm{x})$$

$$=\mathcal{T}(t,x){\mathcal{X}}^{-1}(t)X(t)$$

$$x(t)$$

$$\check{\tau }(t,x)$$

$$\dot{X}(t)$$

$${\dot{x}}_0(t)$$

$$A(t,x)$$

$$=\left(F(t,x)-\mathcal{T}(t,x)\right){\mathcal{X}}^{-1}(t)$$

$$=Q^x(t,x)-{\mathcal{R}}^{-1}(t,x)$$

$$\dot{X}(t)$$

$${\dot{x}}_0(t)$$

$$T(t,\mathrm{x})$$

$$\tau (t,x)$$

Full text

{AMS}

34A34, 35A24, 37C60, 37N25, 37N40, 70G60, 91B74, 92B20, 92C42, 92D30, 92D40, 93C15, 94A15

1 Introduction

Compartmental systems are mathematical abstractions of networks that model behaviors of continuous physical systems composed of discrete living and nonliving homogeneous components. Based on conservation principles, system compartments are interconnected through the flow of energy, matter, or currency between them and their environment. Therefore, formulating flows and associated storages accurately and explicitly is critically important in quantifying compartmental system function. Various mathematical aspects of compartmental systems are studied in the literature [27, 2]. While many fields utilize compartmental modeling, this approach proves particularly well-suited for analysis of ecological systems to address environmental phenomena.

Due to the current technological advancements as well as scientific understandings of population and industrial growth and resource demands, environmental issues have assumed center stage in human communities. On the other hand, in spite of this increased attention to the environment, traditional ecology has an applied nature and is still in the empirical stage of development. In the mainstream framework of traditional ecology, a first principles-based formal theory has yet to emerge. This disconnect narrows the scope of applicability of the field and reduces its ability to deal with complex organism-environment relationships. To that extent, ecology and environmental science are limited in their capacity to realistically model and analyze complex systems. Mathematical theories and modeling have significant potential to lead the way to a more formalistic and theoretical ecoscience devoted to the discovery of basic scientific laws. More exact, precise, and incisive environmental applications can then be materialized based on this understanding.

Sound rationales have been offered in the literature for ecological network analysis, but these are for special cases, such as linear and static models. One such static approach called the environ theory has been developed by [38, 34] based on economic input-output analysis of [31, 32] introduced into ecology by [20]. Ecological networks and complexity in living systems are analyzed also in the context of information theory, thermodynamics [44, 24, 45, 46], and hierarchy theory [1], yet only for static systems. Several software have been developped to computerize these static methods [47, 8, 15, 28, 40, 5].

Although the steady-state analysis is well-established, dynamic analysis of nonlinear compartmental systems has remained a long-standing, open problem. For example, Finn’s cycling index—a celebrated ecosystem measure that quantifies cycling system flows defined in static ecological network analyses over four decades ago—has still not been made applicable to ecosystem models that change over time [17]. The indirect effects in ecosystems have also long been a well-established empirical fact [37, 42, 49, 36, 35, 48]. Theoretical explorations of the concept began as early as the 1970s, and it has been a topic of scholarly conversation for the past five decades [26, 39, 16, 33, 13]. Despite the urgent need, the indirect flow and storage transfers have never been formulated before. There are earlier approaches in the literature for the analysis of dynamic ecosystems, but these are either essentially closed-form abstract formulations [19], or designed for special cases, such as linear systems with time-dependent inputs [23]. In addition, there are also agent-based techniques for dynamic compartmental system analysis [41, 30, 29]. These are, however, computational methods that rely on network particle tracking simulations.

In ecosystem ecology, food webs provide a framework to link community structure with flows of energy and material through trophic interactions and, therefore, relate biodiversity with ecosystem function. Temporal variation in web architecture and nonlinearity are discussed in the literature [14, 51]. It is suggested that the dynamic nature of food webs is affecting ecosystem attributes. Nonlinearity and dynamic behavior, such as extinction in food webs, however, has yet to be addressed methodologically. Not only food webs, but today’s major environmental and ecological phenomena and problems–human impact, climate change, biodiversity loss, etc.– all involve change, which demonstrates that the need for dynamic methods for nonlinear system analysis is not only appropriate, but also urgent [7, 22].

This is the first manuscript in the literature that potentially addresses the disconnect between the current static and computational methods and applied ecological needs. We consider that the methodology proposed herein, in effect, brings a novel complex system theory to the service of pressing and challenging environmental problems of the day. Due to its theoretical and mathematical nature, it has the potential to lead the way to a more formalistic ecological science. The proposed methodology is a comprehensive approach in the sense that the major concepts and quantities in the current static ecological network analyses are extended from static to nonlinear dynamic settings, as well as integrated effectively with the proposed dynamic measures in this unifying mathematical framework. This novel and unifying approach leads to a holistic analysis of ecosystems.

Aligned with the mathematical theory, called the system decomposition theory, introduced recently by [10], the proposed comprehensive method is composed of the novel dynamic system and subsystem partitioning methodologies. The system partitioning methodology yields the subthroughflow and substorage vectors and matrices that represent the flows and storages generated by the initial stocks and individual environmental inputs in each compartment separately. Therefore, the system partitioning enables dynamically decomposing composite compartmental flows and storages into subcompartmental segments based on their constituent sources from the initial stocks and environmental inputs. In other words, this methodology enables dynamically tracking the evolution of the initial stocks and environmental inputs, as well as the associated storages derived from these stocks and inputs individually and separately within the system.

The transient flows transmitted along a given flow path and the associated storages generated by these flows in each compartment on the path are then formulated through the subsystem partitioning methodology. Therefore, this methodology allows for the dynamic decomposition of arbitrary composite intercompartmental flows and associated storages into the transient subflow and substorage segments along a given set of subflow paths. Consequently, the subsystem partitioning enables dynamically tracking the fate of arbitrary intercompartmental flows and associated storages within the subsystems. Moreover, the spread of an arbitrary flow or storage segment from one compartment to the entire system can be determined and monitored. For the quantification of intercompartmental flow and storage transfer dynamics, the direct, indirect, acyclic, cycling, and transfer (diact) flows and associated storages transmitted from one compartment, directly or indirectly, to any other are also analytically characterized, systematically classified, and mathematically formulated.

In a nutshell, the system and subsystem partitioning methodologies dynamically determine the distribution of the initial stocks, environmental inputs, and arbitrary intercompartmental flows, as well as the organization of the associated storages derived from these stocks, inputs, and flows individually and separately within the system. In other words, the proposed method as a whole enables tracking the evolution of the initial stocks, environment inputs, and arbitrary intercmopartmental system flows, as well as associated storages individually and separately. The dynamic quantities such as the subthroughflows, substorages, and transient and diact flows and storages are systematically introduced through the proposed method for the first time in the literature. Equipped with these measures, the proposed methodology serves as a quantitative platform for testing empirical hypotheses, ecological inferences, and, potentially, theoretical developments. The method also constructs a foundation for the development of new mathematical system analysis tools as quantitative ecological indicators. Multiple such dynamic diact measures and indices of matrix, vector, and scalar types which may prove useful for environmental assessment and management were systematically introduced by [9].

The temporal variations of trophic interactions in food webs is an important topic in ecology as outlined above [14, 51]. The conditions or states of communities in food webs, such as extinction, can be dynamically regulated by the temporal variations and seasonal shifts. The present manuscript develops also a novel mathematical technique based on the diact transactions for the dynamic classification of interspecific interactions, and notably, for the determination of their strength within food webs. This technique effectively addresses the nonlinearity in and dynamic architecture of the food chains and webs.

The proposed methodology is applicable to any conservative compartmental system of naturogenic or anthropogenic nature. The method can be used, for example, to analyze models designed for material flows in industry [3]. It can also be used to analyze mass or energy transfers between species of different trophic levels in a complex network or along a given food chain of a food web in nonlinear dynamic settings [21, 4, 18]. Although the motivating applications are ecological and environmental for this paper, the applicability of the proposed method extends to other realms, such as economics, pharmacokinetics, chemical reaction kinetics, epidemiology, biomedical systems, neural networks, social networks, and information science—in fact, wherever dynamical compartmental models of conserved quantities can be constructed. An input-output analysis in economics was developed several decades ago, but only for static systems [31, 32]. The proposed methodology in the context of economics, in particular, can be considered as the mathematical foundation of the dynamic input-output economics.

The proposed method is applied to two models in Section 3 to illustrate its efficiency and wide applicability. In the first case study, a linear ecosystem model introduced by [23] is analyzed. The second case study concerns nutrient transfer within a nutrient-producer-consumer ecosystem [19]. Analytical and numerical solutions for the substorages, subthroughflows, transient and diact transactions, and residence times are presented for both models. The interspecific interactions in the nonlinear model and their strength are also analyzed through the proposed mathematical classification technique.

This paper is organized as follows: the mathematical method is introduced in Section 2.1, the transient and diact flows and storages are formulated in Section 2.5, system analysis and measures are discussed in Section 2.7, results and case studies are presented in Section 3, and discussion and conclusions follow in Section 4 and 5.

2 Methods

A new mathematical theory for the dynamic decomposition of nonlinear compartmental systems has recently been introduced by [10]. In line with this theory, a mathematical method for the dynamic analysis of nonlinear ecological systems is developed in the present paper.

The proposed theory is based on the novel system and subsystem partitioning methodologies. The system and subsystem partitioning determine the distribution of the initial stocks, environmental inputs, and intercompartmental flows, as well as the organization of the associated storages derived from these stocks, inputs, and flows individually and separately within compartmental systems. The proposed method, therefore, as a whole, yields the decomposition of all system flows and storages to the utmost level. The method together with the corresponding concepts and quantities will be introduced in this section.

The terminology and notations used in this paper are adopted from [10] as follows:

[TABLE]

The governing equations for the compartmental dynamics are

[TABLE]

for $i=1,\ldots,n$. The state vector $x(t)=[x_{1}(t),\ldots,x_{n}(t)]^{T}$ is a differentiable function of compartmental storages with the initial conditions of $x(t_{0})=x_{0}=[x_{1,0},\ldots,x_{n,0}]^{T}$ where the superscript $T$ represents the matrix transpose. The total inflow, $\check{\tau}_{i}(t,x)$, and outflow, $\hat{\tau}_{i}(t,x)$, are called the inward and outward throughflows at compartment $i$, respectively, and formulated as

[TABLE]

for $i=1,\ldots,n$. The nonlinear differentiable function $f_{ij}(t,x)\geq 0$ represents nonnegative flow rate from compartment $j$ to $i$ at time $t$. In general, it is assumed that $f_{ii}(t,x)=0$, but the following analysis is also valid for nonnegative flow from a compartment into itself. Index $j=0$ stands for the environment. We further assume that $f_{ij}(t,x)$ has the following special form:

[TABLE]

where $q^{x}_{ij}(t,x)$ is a nonlinear function of ${x}$ and $t$, and has the same properties as $f_{ij}(t,x)$. We will call $q^{x}_{ij}(t,x)={f_{ij}(t,x)}/{x_{j}(t)}$ the flow intensity directed from compartment $j$ to $i$ per unit storage or the flow distribution factor for system storages in the context of the proposed methodology [11].

Combining Eqs. 1 and 2 and separating environmental inputs and outputs, the system of governing equations takes the following standard form:

[TABLE]

with the initial conditions $x_{i}(t_{0})=x_{i,0}$, for $i=1,\ldots,n$. There are $n$ equations; one for each compartment. The condition, Eq. 3, guaranties non-negativity of the compartmental storages, that is $x_{i}(t)\geq 0$ for all $i$. If an environmental input or initial condition is positive, that is, $z_{i}(t,x)>0$ or $x_{i,0}>0$, the corresponding storage value is always strictly positive, $x_{i}(t)>0$.

The proposed methodology is designed for conservative compartmental systems. A dynamical system is called compartmental if it can be expressed in the form of Eq. 4. The compartmental systems will be called conservative if all internal flow rates add up to zero when the system is closed, that is, when there is neither environmental input nor output:

[TABLE]

where $\mathbf{0}$ is used for both the $n\times n$ zero matrix and zero vector of size $n$ [10].

For notational convenience, we define a direct flow matrix function $F$ of size $n\times n$ as

[TABLE]

and the inward and outward throughflow vector functions as

[TABLE]

respectively, where

[TABLE]

are the input and output vector functions, and ${\mathbf{1}}$ denotes the column vector of size $n$ whose entries are all one.

2.1 System partitioning methodology

In this section, we introduce the dynamic system partitioning methodology for analytically partitioning the governing system into mutually exclusive and exhaustive subsystems, as a simplified version of the system decomposition methodology recently proposed by [10]. By mutual exclusiveness, we mean that transactions are possible only among corresponding subcompartments belong to the same subsystem. By exhaustiveness, we mean that all generated subsystems sum to the entire system so partitioned. The system partitioning enables dynamically partitioning composite compartmental flows and storages into subcompartmental segments based on their constituent sources from the initial stocks and environmental inputs of the same conserved quantity. The system partitioning methodology, consequently, yields the subthroughflow and substorage matrices representing the distribution of the initial stocks and environmental inputs, as well as the organization of the associated storages derived from these stocks and inputs individually and separately within the system. In other words, this methodology enables tracking the evolution of the initial stocks and environmental inputs, as well as associated storages individually and separately within the system.

The system partitioning involves the dynamic subcompartmentalization and flow partitioning components, whose mechanisms are explained in this section (Figs. 1 and 2). The related concepts and notations are summarized below:

[TABLE]

The system is partitioned explicitly and analytically into mutually exclusive and exhaustive subsystems as follows: Each compartment is partitioned into $n+1$ subcompartments; $n$ initially empty subcompartments for $n$ environmental inputs and $1$ subcompartment for the initial stock of the compartment. The notation $i_{k}$ is used to represent the $k^{th}$ subcompartment of the $i^{th}$ compartment for $i=1,\ldots,n$ and $k=0,\ldots,n$. The subscript index $k=0$ represents the initial subcompartment of compartment $i$ (see Fig. 1).

The storage in subcompartment $i_{k}$ will be called the substorage in $i_{k}$ and denoted by $x_{i_{k}}(t)$. More specifically, the substorage $x_{i_{k}}(t)$ is defined as the storage in compartment $i$ at time $t$ that is generated by the environmental input into compartment $k\neq 0$, $z_{k}(t)$, during the time interval $[t_{0},t]$ (see Fig. 2). Consequently, due to the exhaustiveness of the system partitioning, we have

[TABLE]

We define a new vector variable for the substorages as

[TABLE]

We assume that environmental input $z_{k}(t,\rm{x})$ enters the system at subcompartment $k_{k}$, for all $k$. Moreover, no other $k^{th}$ subcompartment of any other compartment $i$, that is, subcompartment $i_{k}$, receives environmental input. This input partitioning can be expressed as

[TABLE]

The intercompartmental flows are also partitioned in line with the subcompartmentalization (see Fig. 2). The composite intercompartmental direct flow, $f_{ij}(t,x)$, is partitioned based on the assumption that the subcompartmental flow segments, $f_{{i_{k}}{j_{k}}}(t,{\rm x})$, $k=0,\ldots,n$, are proportional to the corresponding substorages, $x_{j_{k}}(t)$, with the proportionality factor of the flow intensity in the flow direction, $q^{x}_{ij}(t,{x})$. The subcompartmental flow $f_{{i_{k}}{j_{k}}}(t,{\rm x})$ will be called the subflow from subcompartment $j_{k}$ to $i_{k}$ at time $t$. It can be formulated as follows:

[TABLE]

where the coefficients $d_{j_{k}}({\rm x})={x_{j_{k}}(t)}/{x_{j}(t)}$ will be called the decomposition factors. Consequently, due to the exhaustiveness of the system partitioning, we have

[TABLE]

In summary, the dynamic system partitioning methodology explicitly generates mutually exclusive and exhaustive subsystems running within the original system. The $k^{th}$ subsystem is composed of all $k^{th}$ subcompartments of each compartment together with the corresponding subflows and substorages. These subsystems have the same structure and dynamics as the original system itself, except for their environmental inputs and initial conditions (see Figs. 1 and 2). Each subsystem, except the initial one—which is driven by the initial stocks—is generated by a single environmental input. Therefore, the number of non-intersecting subcompartments in each compartment is equal to the number of inputs or compartments, plus one for the initial stocks. If an input or all initial conditions are zero, the corresponding subsystem becomes null. Consequently, for a system with $n$ compartments, each compartment has $n+1$ non-intersecting subcompartments, and therefore the system has $n+1$ mutually exclusive subsystems indexed by $k=0,\ldots,n$. The initial subsystem ($k=0$) represents the evolution of the initial stocks, receives no environmental input, and has the same initial conditions as the original system. The initial conditions for all the other subcompartments ($k\neq 0$) are zero, since they are initially assumed to be empty.

The governing equation for each subcompartment $i_{k}$ then becomes

[TABLE]

for $i=1,\ldots,n$, $k=0,\ldots,n$. There are $n\times(n+1)$ of such governing equations, one for each subcompartment. In order to track the evolution of environmental inputs within the system individually and separately, all except the initial subcompartments are assumed to be initially empty, as mentioned above. Therefore, the initial conditions become

[TABLE]

The governing system of equations, Eq. 10, is solved numerically with the initial conditions, Eq. 11. The result yields the substorages at any time $t$, that is, $x_{i_{k}}(t)$.

The total subcompartmental inflows and outflows at compartment $i$ at time $t$ generated by the environmental input into compartment $k$, $z_{k}(t)$, during $[t_{0},t]$ can then be defined, respectively, as

[TABLE]

for $k=0,1,\ldots,n$. The functions $\check{\tau}_{i_{k}}(t,{\rm x})$ and $\hat{\tau}_{i_{k}}(t,{\rm x})$ will respectively be called inward and outward subthroughflow at subcompartment $i_{k}$ at time $t$ (see Fig. 4). Therefore, the system partitioning enables dynamically partitioning composite compartmental flows and storages into subcompartmental segments based on their constituent sources from the initial stocks and environmental inputs.

We define the $n\times n$ substorage and associated inward and outward subthroughflow matrix functions, $X(t)$, $\check{T}(t,{\rm x})$, and $\hat{T}(t,{\rm x})$, respectively, as follows:

[TABLE]

for $i,k=1,\ldots,n$. The substorage and associated inward and outward subthroughflow vector functions of size $n$ for the initial subsystem, ${x}_{0}(t)$, $\check{\tau}_{0}(t,{\rm x})$, and $\hat{\tau}_{0}(t,{\rm x})$, can also be defined, respectively, as

[TABLE]

We use the constant vector notation $x_{0}$ for the initial conditions and the function notation $x_{0}(t)$ for the evolution of these initial stoks for $t>t_{0}$ with $x_{0}(t_{0})=x_{0}$.

The notation $\operatorname{diag}({x(t)})$ will be used to represent the diagonal matrix whose diagonal elements are the elements of vector $x(t)$, and $\operatorname{diag}({X(t)})$ to represent the diagonal matrix whose diagonal elements are the same as the diagonal elements of matrix $X(t)$. The $n\times n$ diagonal storage, output, and input matrix functions, $\mathcal{X}(t)$, $\mathcal{Y}(t,x)$, and $\mathcal{Z}(t,x)$ will be defined, respectively, as

[TABLE]

Using Eq. 8, the subthroughflow matrices can then be formulated as follows:

[TABLE]

where $\mathcal{T}(t,x)=\operatorname{diag}{(\hat{\tau}(t,x))}=\mathcal{Y}(t,x)+\operatorname{diag}\left(F^{T}(t,x)\,\mathbf{1}\right)$. Note that,

[TABLE]

The governing equations for the decomposed system, Eq. 10, can be expressed in terms of the vector and matrix functions introduced above as follows:

[TABLE]

We define an $n\times n$ matrix function $A(t,x)$ as

[TABLE]

where $Q^{x}(t,{x})=F(t,{x})\,\mathcal{X}^{-1}(t)$ and $\mathcal{R}^{-1}(t,{x})=\mathcal{T}(t,{x})\,\mathcal{X}^{-1}(t)$, assuming $\mathcal{R}(t,{x})$ is invertible. Note that the first term in the definition of $A(t,x)$, $Q^{x}(t,{x})$, represents the intercompartmental flow intensity defined in Eq. 8, and the second term, $\mathcal{R}^{-1}(t,x)$, represents the outward throughflow intensity. Consequently, we will call $A(t,x)$ the flow intensity matrix per unit storage. It is sometimes called the compartmental matrix [27]. The matrix measure $\mathcal{R}(t,{x})$ will be called the residence time matrix, and the matrix measure $Q^{x}(t,{x})$ will be called the flow intensity matrix per unit storage or the flow distribution matrix for system storages [10, 11]. The governing equations, Eq. 16, can be expressed using the flow intensity matrix in the following from:

[TABLE]

The dynamic system partitioning methodology that yields a decomposed system of $n^{2}+n$ governing equations for all subcompartments, Eq. 18, from the original system of $n$ governing equations for all compartments, Eq. 1, can algebraically be schematized as follows (see Figs. 1 and 2 for graphical illustrations):

{\dot{x}(t)={\tau}(t,x)}$${\mathrel{\rotatebox{90.0}{=}}}$${{\begin{bmatrix}\dot{x}_{1}\\ \dot{x}_{2}\\ \vdots\\ \dot{x}_{n}\end{bmatrix}=\begin{bmatrix}{\tau}_{1}\\ {\tau}_{2}\\ \vdots\\ {\tau}_{n}\end{bmatrix}}}$${{\begin{bmatrix}\dot{x}_{1_{0}}\\ \dot{x}_{2_{0}}\\ \vdots\\ \dot{x}_{n_{0}}\end{bmatrix}=\begin{bmatrix}{\tau}_{1_{0}}\\ {\tau}_{2_{0}}\\ \vdots\\ {\tau}_{n_{0}}\end{bmatrix}}}$${{\quad\mbox{and}\quad\begin{bmatrix}\dot{x}_{1_{1}}&&\cdots&&\dot{x}_{1_{n}}\\ \dot{x}_{2_{1}}&&\cdots&&\dot{x}_{2_{n}}\\ \vdots&&\ddots&&\vdots\\ \dot{x}_{n_{1}}&&\cdots&&\dot{x}_{n_{n}}\end{bmatrix}=\begin{bmatrix}{\tau}_{1_{1}}&&\cdots&&{\tau}_{1_{n}}\\ {\tau}_{2_{1}}&&\cdots&&{\tau}_{2_{n}}\\ \vdots&&\ddots&&\vdots\\ {\tau}_{n_{1}}&&\cdots&&{\tau}_{n_{n}}\end{bmatrix}}}$${\mathrel{\rotatebox{90.0}{=}}}$${\mathrel{\rotatebox{90.0}{=}}}$${{\dot{x}_{0}(t)={\tau}_{0}(t,{\rm x})}}$${\dot{X}(t)=T(t,{\rm x})}dynamic system partitioning

In the diagram above, the net subthroughflow matrix, $T(t,{\rm x})$, as well as the net throughflow and initial throughflow vectors, ${\tau}(t,{\rm x})$ and ${\tau}_{0}(t,{\rm x})$, are defined as the difference between the corresponding inward and outward throughflows. That is,

[TABLE]

The system partitioning introduced in this section is input-oriented. The governing system can be partitioned based on environmental outputs instead of inputs, by conceptually reversing all system flows. The following condition on the direct flows, instead of Eq. 3, ensures the possibility of the system partitioning and analysis in both the input- and output-orientations:

[TABLE]

This form of flow rates makes both the original and reversed decomposed systems well-defined.

2.2 Subsystem flows and storages

The system partitioning methodology dynamically decomposes a nonlinear system into mutually exclusive and exhaustive subsystems through the formulation of a set of governing equations derived from subcompartmentalization and flow partitioning components, as introduced above. This methodology enables the dynamic analysis of subsystems generated by the initial stocks and environmental inputs individually and separately. The subsystem flows and storages are formulated in matrix form in this section.

We define the $k^{th}$ direct subflow matrix function for the $k^{th}$ subsystem as $F_{k}(t,{\rm x})=\left(f_{i_{k}j_{k}}(t,{\rm x})\right)$, $k=0,\ldots,n$. Using the relationships formulated in Eq. 8, this matrix can be expressed as follows:

[TABLE]

where $\mathcal{X}_{k}(t)=\operatorname{diag}\left([x_{1_{k}}(t),\ldots,x_{n_{k}}(t)]\right)$ is the diagonal matrix of the substorage functions in the $k^{th}$ subsystem. The matrix $\mathcal{X}_{k}(t)$ will accordingly be called the $k^{th}$ substorage matrix function. The $k^{th}$ output and input matrix functions then become

[TABLE]

where $\bm{e}_{k}$ is the elementary unit vector whose components are all zero except the $k^{th}$ element, which is $1$, and we set $\bm{e}_{0}={\bm{0}}$. The $k^{th}$ direct subflow matrix, $F_{k}(t,{\rm x})$, and the $k^{th}$ input and output vector functions, defined as $\check{z}_{k}(t,{\rm x})=\mathcal{Z}_{k}(t,{\rm x})\,\bm{1}$ and $\hat{y}_{k}(t,{\rm x})=\mathcal{Y}_{k}(t,{\rm x})\,\bm{1}$, are the counterparts for the $k^{th}$ subsystem of the direct flow matrix, $F(t,x)$, and the input and output vectors, $z(t,x)$ and $y(t,x)$, for the original system. Altogether, they represent the subflow regime of the $k^{th}$ subsystem.

Using the notations and definitions of Eqs. 15 and 21, the $k^{th}$ inward and outward subthroughflow matrices, $\check{\mathcal{T}}_{k}(t,{\rm x})=\operatorname{diag}{\left([\check{\tau}_{1_{k}}(t,{\rm x}),\ldots,\check{\tau}_{n_{k}}(t,{\rm x})]\right)}$ and $\hat{\mathcal{T}}_{k}(t,{\rm x})=\operatorname{diag}{\left([\hat{\tau}_{1_{k}}(t,{\rm x}),\ldots,\hat{\tau}_{n_{k}}(t,{\rm x})]\right)}$, for the $k^{th}$ subsystem can be expressed as follows:

[TABLE]

We also define the decomposition and $k^{th}$ decomposition matrices, ${D}({\rm x})=(d_{i_{k}}({\rm x}))$ and $\mathcal{D}_{k}({\rm x})=\operatorname{diag}{([d_{1_{k}}({\rm x}),\ldots,d_{n_{k}}({\rm x})])}$, as

[TABLE]

The second equalities in the definitions of ${D}({\rm x})$ and $\mathcal{D}_{k}({\rm x})$ are due to Eq. 15 and 22, respectively. Note that ${D}({\rm x})$ and $\mathcal{D}_{k}({\rm x})$ decompose the compartmental throughflow matrix, $\mathcal{T}(t,{x})$, into the subthroughflow and $k^{th}$ subthroughflow matrices as indicated in Eqs. 15 and 22, similar to the decomposition of ${F}_{k}(t,{\rm x})$ as formulated below in Eq. 25. That is,

[TABLE]

The $k^{th}$ direct subflow and substorage matrices, $F_{k}(t,{\rm x})$ and $\mathcal{X}_{k}(t)$, can then be written in the following various forms:

[TABLE]

where $Q^{\tau}(t,{x})=F(t,{x})\,\mathcal{T}^{-1}(t,{x})$ will be called the flow intensity matrix per unit throughflow or the flow distribution matrix for system flows in the context of the proposed methodology [11]. Note that the elements of $Q^{\tau}(t,{x})$, $q^{\tau}_{ij}(t,{x})$, are sometimes called transfer coefficients, technical coefficients in economics, or stoichiometric coefficients in chemistry. The system level counterpart of Eq. 25 can also be formulated as follows:

[TABLE]

using Eq. 15. The matrix function $\tilde{T}(t,{\rm x})=\check{T}(t,{\rm x})-\mathcal{Z}(t,{x})$ will be called the intercompartmental subthroughflow matrix. Componentwise, $\tilde{T}(t,{\rm x})=(\tilde{\tau}_{i_{k}}(t,{\rm x}))$ can be expressed as $\tilde{\tau}_{i_{k}}(t,{\rm x})=\check{\tau}_{i_{k}}(t,{\rm x})-z_{i_{k}}(t)$. Consequently, for any given storage organization within the system, $X(t)$, the corresponding intercompartmental flow distributions—that is, intercompartmental inward and outward subthroughflow matrices, $\tilde{T}(t,{\rm x})$ and $\hat{T}(t,{\rm x})$—can be determined at any time $t$ as follows:

[TABLE]

It is worth noting that the residence time matrix can be written in the following various forms:

[TABLE]

as formulated in Eqs. 17, 25, and 26. For later use, we also define three diagonal matrices as follows:

[TABLE]

2.3 Analytic solution to linear systems

This section formulates analytic solutions to linear systems with time-dependent inputs. The system partitioning methodology yields a linear system, if the original system is linear. That is, if Eq. 4 is linear, Eq. 18 is also linear. Due to the cancellations in matrix $A(t,x)$ defined in Eq. 17, the decomposed linear system, Eq. 18, can be expressed in the following matrix form:

[TABLE]

Let $V(t)$ be the fundamental matrix solution to the system Eq. 30, as defined by [10]. That is, let $V(t)$ be the unique solution of the system

[TABLE]

The solutions to Eq. 30 for substorage matrix, $X(t)$, and initial substorage vector, $x_{0}(t)$, in terms of $V(t)$, become

[TABLE]

as formulated by [10]. Therefore, the solution to the original system, Eq. 1 can be written as

[TABLE]

For the special case of constant diagonalizable flow intensity matrix $A$, we have

[TABLE]

where $\Omega$ is the matrix whose columns are the eigenvectors of A, and $\Lambda$ is the diagonal matrix whose diagonal elements are the eigenvalues of $A$. For this particular case, Eq. 31 takes the following form:

[TABLE]

Consequently,

[TABLE]

A subsystem scaling argument is proposed to analyze static system behavior per unit input by [11]. The scaled substorage matrix, $S(t)=X(t)\,\mathcal{Z}^{-1}$, can be expressed for constant invertible input matrix, $\mathcal{Z}(t)=\mathcal{Z}>{\bm{0}}$, as follows:

[TABLE]

using Eq. 31 and 33. The static version of this measure $S(t)$ is widely used in static ecological network analyses as outlined in the next section [11].

An example of the analytic solution to a linear ecosystem model with time dependent environmental input is presented in Section 3.1.

2.4 Static ecological system analysis

At steady state, the time derivatives of the state variables are zero, and all system flows and storages are constant. That is,

[TABLE]

The constant static quantities will be denoted by the same symbols without the time argument. The constant substorage matrix, for example, will be denoted by ${X}(t)=X$.

Summing up the equations in Eq. 10 side by side over index $k$ yields Eq. 4 because of the relationship

[TABLE]

deduced from Eq. 7 and the definition of the decomposition factors, $d_{i_{k}}({\rm x})$, given in Eq. 8. Therefore, if the partitioned system, Eq. 10, is at steady state, the original system, Eq. 4, is also at steady state. The static version of the proposed dynamic methodology is introduced by [11, 12], as summarized below in this section.

Since $A$ is a strictly diagonally dominant constant matrix, it is invertible. It can be expressed as

[TABLE]

We then have the following solutions to Eq. 18 for the substorage matrix, $X(t)$, and initial substorage vector, $x_{0}(t)$, at steady state:

[TABLE]

From Eq. 15 and the fact that ${\tau}=\hat{\tau}=\check{\tau}$ and $T=\check{T}=\hat{T}$ at steady state, the throughflow matrix can be written in terms of system flows only:

[TABLE]

The residence time matrix $\mathcal{R}$ can also be expressed as

[TABLE]

similar to Eq. 28.

The scaled substorage and subthroughflow matrices are defined for system analysis per unit input by [11]. They are formulated as $S=X\,\mathcal{Z}^{-1}$ and $N=T\,\mathcal{Z}^{-1}$, where $\mathcal{Z}$ is invertible. Using Eqs. 36 and 37, these matrix measures can be expressed as follows:

[TABLE]

Note that $S(t)$ formulated in Eq. 34 is equivalent to $S$ at steady state, that is

[TABLE]

The matrices $N$ and $S$ are called the cumulative flow and storage distribution matrices in the context of the proposed methodology [11].

Although the derivation rationales are different, the proposed matrix measures $S$ and $N$ are equivalent to the ones formulated in the current static ecological network analyses, as shown by [11]. These matrices are treated separately in the current static methodologies, although they are naturally related by a factor of the residence time matrix. Equation 38 implies that

[TABLE]

This relationship enables the holistic view of static ecological networks:

[TABLE]

as introduced by [12].

The substorage and subthroughflow matrices can be scaled by output matrix instead, for the output-oriented system analysis. We use a bar notation over the output-oriented counterparts of the input oriented quantities. In the output-oriented analysis, we assume that all system flows are conceptually reversed. That is,

[TABLE]

The counterparts of the Eq. 40 for the output-oriented analysis then become

[TABLE]

where $\bar{S}=\bar{X}\,\mathcal{Y}^{-1}$, $\bar{N}=\bar{T}\,\mathcal{Y}^{-1}$, and the diagonal output matrix, $\mathcal{Y}$, is assumed to be invertible. Since $x=\bar{x}$ and $\tau=\bar{\tau}$ at steady state, $\mathcal{\bar{R}}=\mathcal{R}$.

It is worth noting that the input- and output-oriented cumulative flow and storage distribution matrices are similar. This duality can be expressed as

[TABLE]

The holistic input- and output-oriented static ecological system analyses and their duality have recently been introduced by [11, 12].

2.5 Subsystem partitioning methodology

In this section, we introduce the dynamic subsystem partitioning methodology for further partitioning or segmentation of subsystems along a given set of mutually exclusive and exhaustive subflow paths, as a simplified version of the subsystem decomposition methodology recently proposed by [10]. The subsystem partitioning methodology dynamically apportions arbitrary composite intercompartmental flows and the associated storages generated by these flows into transient subflow and substorage segments along given subflow paths. The subsystem partitioning, therefore, determines the distribution of arbitrary intercompartmental flows and the organization of the associated storages generated by these flows within the subsystems. In other words, this methodology enables tracking the evolution of arbitrary intercompartmental flows and associated storages within and monitoring their spread throughout the system.

The natural subsystem decomposition is defined as the set of mutually exclusive and exhaustive subflow paths whose local inputs and outputs, except for the closed paths, are environmental inputs and outputs, respectively [10, 11]. By mutually exclusive subflow paths, we mean that no given subflow path is a subpath, that is, completely inside of another path in the same subsystem. Exhaustiveness in this context means that such mutually exclusive subflow paths—together with the corresponding transient subflows and substorages along the paths—all together sum to the entire subsystem so partitioned. The natural subsystem decomposition of each subsystem then results in a mutually exclusive and exhaustive decomposition of the entire system.

We will first introduce the transient flows and storages below. Thereafter, they will then be used for the formulation of the diact flows and storages in the subsequent section.

*No man ever steps in the same river twice. – Heraclitus (535-475 BC) *

2.5.1 Transient flows and storages

As indicated in the famous dictum by Heraclitus that “everything flows,” flows are one of the most important physical phenomena of existence. In this section, we formulate the transient flows and the associated storages generated by these flows.

The transient and cumulative transient subflows along a subflow path are defined as follows: Along a given subflow path $p^{w}_{n_{k}j_{k}}=i_{k}\mapsto j_{k}\to\ell_{k}\to n_{k}$, the transient inflow at subcompartment ${\ell_{k}}$, $f^{w}_{\ell_{k}j_{k}i_{k}}(t)$, generated by the local input from $i_{k}$ to ${j_{k}}$ during $[t_{1},t]$, $t_{1}\geq t_{0}$, is the input segment that is transmitted from $j_{k}$ to ${\ell_{k}}$ at time $t$. Similarly, the transient outflow generated by the transient inflow at ${\ell_{k}}$ during $[t_{1},t]$, $f^{w}_{n_{k}\ell_{k}j_{k}}(t)$, is the inflow segment that is transmitted from ${\ell_{k}}$ to the next subcompartment, ${n_{k}}$, along the path at time $t$. The associated transient substorage in subcompartment ${\ell_{k}}$ at time $t$, $x^{w}_{n_{k}\ell_{k}j_{k}}(t)$, is the substorage segment governed by the transient inflow and outflow balance during $[t_{1},t]$ (see Fig. 3).

The transient outflow at subcompartment ${\ell_{k}}$ at time $t$ along subflow path $p^{w}_{n_{k}j_{k}}$ from ${j_{k}}$ to ${n_{k}}$, $f^{w}_{n_{k}\ell_{k}j_{k}}(t)$, can be formulated as follows:

[TABLE]

similar to Eq. 8, due to the equivalence of flow and subflow intensities, where the transient substorage, $x^{w}_{n_{k}\ell_{k}j_{k}}(t)$, is determined by the governing mass balance equation

[TABLE]

The equivalence of the throughflow and subthroughflow intensities, as well as the flow and subflow intensities in the same direction, that is

[TABLE]

are given by Eqs. 8 and 15, for $\ell,n=1,\ldots,n$, and $k=0,1,\ldots,n$ [10]. Therefore, since the intensities in Eqs. 43 and 44 can be expressed at both the compartmental and subcompartmental levels, the subsystem partitioning is actually independent from the system partitioning. That is, the same analysis can be done along flow paths within the system, instead of subflow paths within the subsystems. This allows the flexibility of tracking arbitrary intercompartmental flows and storages generated by all or individual environmental inputs within the system. The governing equations, Eqs. 43 and 44, establish the foundation of the dynamic subsystem partitioning. These equations for each subcompartment along a given flow path of interest will then be coupled with the partitioned system, Eq. 10, or the original system, Eq. 4, and be solved simultaneously. The equations can alternatively be solved individually and separately once the original or partitioned system is solved.

The transient subflows and substorages are defined for linear subflow paths above. The sum of the transient inflows from subcompartment $j_{k}$ to $\ell_{k}$ and the outflows from $\ell_{k}$ to $n_{k}$ generated at subcompartment $\ell_{k}$ at time $t$ by the local input into the connection of a given non-self-intersecting closed subflow path $p^{w}_{n_{k}j_{k}}$ during $[t_{1},t]$, $t_{1}\geq t_{0}$, will respectively be called the inward and outward cumulative transient subflow at subcompartment ${\ell_{k}}$ at time $t$. The associated storage generated by the inward cumulative transient subflow will be called cumulative transient substorage. These inward and outward cumulative transient subflows will be denoted by $\check{\tau}^{w}_{\ell_{k}}(t)$ and $\hat{\tau}^{w}_{\ell_{k}}(t)$, respectively, and associated cumulative transient substorage by ${x}^{w}_{\ell_{k}}(t)$. They can be formulated as

[TABLE]

for $k=0,1,\ldots,n$, where the superscript $m$ represents the cycle number, and ${m_{w}}$ is the number of cycles, that is, the number of times the path $p^{w}_{n_{k}i_{k}}$ pass through subcompartment $\ell_{k}$. Large number of terms, ${m_{w}}$, in computation of these summations reduce truncation errors and, thus, improve the approximations.

Using the equivalence of flow intensities as formulated in Eq. 45, it has been recently shown for compartmental systems that the parallel subflows and the corresponding subthroughflows and substorages are proportional [10]. By parallel subflows, we mean the intercompartmental flows that transit through different subcompartments of the same compartment along the same flow path at the same time. This proportionality can be formulated as

[TABLE]

for $k=1,\ldots,n$ and $\ell=0,\ldots,n$, where the denominators are nonzero.

2.5.2 The diact flows and storages

In this section, we formulate five main transaction types for nonlinear systems at both the subcompartmental and compartmental levels using two approaches: the direct (d), indirect (i), cycling (c), acyclic (a), and transfer (t) flows and the associated storages generated by these diact flows. The first approach based on the subsystem partitioning methodology will be called the path-based approach, and the second approach based on the system partitioning methodology will be called the dynamic approach.

The composite transfer flow will be defined as the total intercompartmental transient flow that is generated by all environmental inputs from one compartment, directly or indirectly through other compartments, to another. The composite direct, indirect, acyclic, and cycling flows from the initial compartment to the terminal compartment are then defined as the direct, indirect, non-cycling, and cycling segments at the terminal compartment of the composite transfer flow (see Fig. 4). The cycling and acyclic flows can, therefore, be interpreted as the flows that visit the terminal compartment multiple times and only once, respectively, after being transmitted from the initial compartment.

The composite transfer subflow within the initial subsystem can also be defined as the total intercompartmental transient subflow that is derived from all initial stocks from one initial subcompartment, directly or indirectly through other initial subcompartments, to another. The composite direct, indirect, acyclic, and cycling subflows within the initial subsystem from the initial subcompartment to the terminal subcompartment are then defined as the direct, indirect, non-cycling, and cycling segments at the terminal subcompartment of the composite transfer subflow.

The simple transfer flow will be defined as the total intercompartmental transient subflow that is generated by the single environmental input from an input-receiving subcompartment, directly or indirectly through other compartments, to another subcompartment. The simple direct, indirect, acyclic, and cycling flows from the initial input-receiving subcompartment to the terminal subcompartment are then defined as the direct, indirect, non-cycling, and cycling segments at the terminal subcompartment of the simple transfer flow (see Fig. 4). The associated simple and composite diact storages are defined as the storages generated by the corresponding diact flows.

Let $P^{\texttt{t}}_{i_{k}j_{k}}$ be the set of mutually exclusive subflow paths $p^{w}_{i_{k}j_{k}}$ from subcompartment $j_{k}$ directly or indirectly to $i_{k}$ in subsystem $k$. The sets $P^{\texttt{d}}_{i_{k}j_{k}}$ and $P^{\texttt{i}}_{i_{k}j_{k}}$ are also defined as the sets of mutually exclusive direct and indirect subflow paths $p^{w}_{i_{k}j_{k}}$ from subcompartment $j_{k}$ directly and indirectly to $i_{k}$, respectively. Similarly, the sets $P^{\texttt{c}}_{i_{k}j_{k}}$ and $P^{\texttt{a}}_{i_{k}j_{k}}$ are defined as the sets of mutually exclusive cyclic and acyclic subflow paths $p^{w}_{i_{k}j_{k}}$ from $j_{k}$ to $i_{k}$ with a closed and linear subpath at terminal subcompartment $i_{k}$, respectively (see Fig. 4). The cyclic subflow set, ${P}^{\texttt{c}}_{i_{k}}$, can alternatively be defined as the set of mutually exclusive subflow paths ${p}^{w}_{i_{k}}$ from subcompartment $i_{k}$ indirectly back to itself. The number of subflow paths in $P^{\texttt{*}}_{i_{k}j_{k}}$ will be denoted by $w_{k}$, where the superscript $(^{\texttt{*}})$ represent any of the diact symbols.

The composite diact subflow from subcompartment $j_{k}$ to $i_{k}$, $\tau^{\texttt{*}}_{i_{k}j_{k}}(t)$, is defined as the sum of the cumulative transient subflows, $\check{\tau}_{i_{k}}^{w}(t)$, generated by the outward subthroughflow at subcompartment $j_{k}$, $\hat{\tau}_{j_{k}}(t,{\rm x})$, during $[t_{1},t]$, $t_{1}\geq t_{0}$, and transmitted into $i_{k}$ at time $t$ along all subflow paths $p^{w}_{i_{k}j_{k}}\in P^{\texttt{*}}_{i_{k}j_{k}}$. The associated composite diact substorage, $x^{\texttt{*}}_{i_{k}j_{k}}(t)$, at subcompartment $i_{k}$ at time $t$ is the sum of the cumulative transient substorages, $x^{w}_{i_{k}}(t)$, generated by the cumulative transient inflows, $\check{\tau}_{i_{k}}^{w}(t)$, during $[t_{1},t]$. Alternatively, $x^{\texttt{*}}_{i_{k}j_{k}}(t)$ can be defined as the storage segment generated by the composite diact inflow $\tau^{\texttt{*}}_{i_{k}j_{k}}(t)$ in subcomaprtment $i_{k}$ during $[t_{1},t]$. Note that, for the cycling case, the first entrance of the transient subflows and substorages into $i_{k}$ are not considered as cycling subflows and substorages.

The composite diact subflows and substorages can then be formulated as follows:

[TABLE]

The sum of all composite diact subflows and substorages from subcompartment $j_{k}$ to $i_{k}$ within each subsystem $k\neq 0$ will be called the composite diact flow and storage from compartment $j$ to $i$ at time $t$, ${\tau}^{\texttt{*}}_{ij}(t)$ and $x^{\texttt{*}}_{ij}(t)$, generated by all environmental inputs during $[t_{1},t]$. They can be formulated as

[TABLE]

For notational convenience, we define $n\times n$ matrix functions ${T}^{\texttt{*}}_{k}(t)$ and ${X}^{\texttt{*}}_{k}(t)$ whose $(i,j)-$elements are $\tau^{\texttt{*}}_{i_{k}j_{k}}(t)$ and $x^{\texttt{*}}_{i_{k}j_{k}}(t)$, respectively. That is,

[TABLE]

for $k=0,\ldots,n$. These matrix measures ${T}^{\texttt{*}}_{k}(t)$ and ${X}^{\texttt{*}}_{k}(t)$ are called the $k^{th}$ composite diact subflow and associated substorage matrix functions. The corresponding composite diact flow and associated storage matrix functions generated by environmental inputs are ${T}^{\texttt{*}}(t)=\left(\tau^{\texttt{*}}_{ij}(t)\right)$ and ${X}^{\texttt{*}}(t)=\left(x^{\texttt{*}}_{ij}(t)\right)$, respectively [10].

The simple dicat flows and storages generated by single environmental inputs can be formulated in terms of their composite counterparts as follows:

[TABLE]

To distinguish the composite and simple diact flow and storage matrices, we use a tilde notation over the simple versions. That is, the simple diact flow and storage matrices, for example, will be denoted by $\tilde{T}^{\texttt{*}}(t)=\left(\tau^{\texttt{*}}_{i_{k}}(t)\right)$ and $\tilde{X}^{\texttt{*}}(t)=\left(x^{\texttt{*}}_{i_{k}}(t)\right)$.

The simple diact throughflow and compartmental storage matrices and vectors can then be formulated as

[TABLE]

The composite counterparts of these quantities can similarly be formulated in parallel.

The difference between the composite and simple diact flows, $\tau^{\texttt{*}}_{ik}(t)$ and $\tau^{\texttt{*}}_{i_{k}}(t)$, and storages, $x^{\texttt{*}}_{ik}(t)$ and $x^{\texttt{*}}_{i_{k}}(t)$, is that the composite flow and storage from compartment $k$ to $i$ are generated by outward throughflow $\hat{\tau}_{k}(t,{x})-\hat{\tau}_{k_{0}}(t,{\rm x})$ derived from all environmental inputs and their simple counterparts from input-receiving subcompartment $k_{k}$ to $i_{k}$ are generated by outward subthroughflow $\hat{\tau}_{k_{k}}(t,{\rm x})$ derived from single environmental input $z_{k}(t,{x})$ (see Fig. 4). In that sense, the composite and simple diact flows and storages measure the influence of one compartment on another induced by all and a single environmental input, respectively. The composite diact subflows and substorages within the initial subsystem, $\tau^{\texttt{*}}_{i_{0}k_{0}}(t)$ and $x^{\texttt{*}}_{i_{0}k_{0}}(t)$, from compartment $k$ to $i$ are then generated by outward throughflow $\hat{\tau}_{k_{0}}(t,{\rm x})$ derived from all initial stocks.

The simple and composite diact flows have been explicitly formulated for static systems through the system partitioning methodology in a recent study by [11]. This static approach will be extended to dynamic systems pointwise in time, that is, at each time step, in what follows. In addition to the path-based approach introduced above, this alternative approach to formulate the dynamic diact flows and storages will be called the dynamic approach.

The simple transfer flow from an input-receiving subcompartment $k_{k}$ to $i_{k}$ can be expressed as follows:

[TABLE]

for $i,k=1,\ldots,n$. Note that the simple transfer flow at subcompartment $i_{k}$ is equal to the intercompartmental subthroughflow at that subcompartment. The simple direct flow from $k_{k}$ to $i_{k}$ is

[TABLE]

The simple indirect flow from $k_{k}$ to $i_{k}$ at time $t$ can then be formulated as the simple transfer flow, $\tau^{\texttt{t}}_{i_{k}}(t)$, diminished by the simple direct flow, $\tau^{\texttt{d}}_{i_{k}}(t)$. That is,

[TABLE]

Due to the reflexivity of the cycling flow, the simple cycling flow from an input-receiving subcompartment $k_{k}$ back into itself can be formulated in terms of the simple indirect or transfer flows as follows:

[TABLE]

That is, the simple cycling flow is the simple transfer or indirect flow from a subcompartment back into itself. The complementary nature of the cycling and indirect flows are schematized in Fig. 5. The proportionality given in Eq. 60 implies that the simple (composite) cycling (sub)flow from an input-receiving (arbitrary) subcompartment $k_{k}$ ($i_{k}$) into $i_{k}$ is

[TABLE]

The simple acyclic flow can, therefore, be formulated as

[TABLE]

Note that the simple acyclic flow from subcompartment $k_{k}$ back into itself is zero:

[TABLE]

This means that there is, obviously, no acyclic flow from a subcompartment back into itself.

The proportionality of the parallel subflows and the corresponding subthroughflows and substorages formulated in Eq. 47 can be expressed for the diact subflows as follows:

[TABLE]

for $i,k=1,\ldots,n$ and $\ell=0,\ldots,n$, where the denominator is nonzero, ${\hat{\tau}}_{k_{k}}(t,{\rm x})\neq 0$. Using this proportionality and the simple diact flows formulated above, the composite diact subflows from subcompartment $k_{\ell}$ to $i_{\ell}$ can also be formulated as follows:

[TABLE]

for $t>t_{0}$. Note that $\hat{\tau}_{k_{k}}(t_{0},{\rm x})=0$ and we assume that $\hat{\tau}_{k_{k}}(t,{\rm x})$ is nonzero for all $t>t_{0}$. The second equality of the first equation in Eq. 61 for the composite direct subflow is due to the equivalence of flow and subflow intensities in the same direction, as formulated in Eq. 45 [10, 11].

The composite diact flows generated by environmental inputs then become

[TABLE]

where the diact flow distribution factor, $n^{\texttt{*}}_{ik}(t)$, is $n^{\texttt{*}}_{ik}(t)={\tau^{\texttt{*}}_{i_{k}}(t)}/{{\hat{\tau}}_{k_{k}}(t,{\rm x})}$. The dynamic diact flow distribution matrix function can be defined as $N^{\texttt{*}}(t)=(n^{\texttt{*}}_{ik}(t))$. The dynamic diact flow distribution matrices, as well as the dynamic simple and composite diact flow matrices are explicitly formulated in Table 1 based on their componentwise definitions in Eq. 61, similar to their static counterparts introduced by [11]. The inverted matrices in the table are assumed to be invertible. For the sake of readability, the function arguments are dropped in the table.

The composite diact substorages can also be formulated using the corresponding diact subflows as transient inflows in Eq. 44 as follows:

[TABLE]

for $t_{1}>t_{0}$, $i,k=1,\ldots,n$ and $\ell=0,\ldots,n$. The solution to this governing equation, ${x}^{\texttt{*}}_{i_{\ell}k_{\ell}}(t)$, represents the diact substorage in subcompartment $i_{\ell}$ at time $t\geq t_{1}$ generated by the corresponding diact subflow, $\tau^{\texttt{*}}_{i_{\ell}k_{\ell}}(t)$, during $[t_{1},t]$ (see Fig. 3).

The analytic solutions for linear systems are introduced in Section 2.3. In this case, the governing equation for the diact substorages, Eq. 63, can be solved explicitly for $x^{\texttt{*}}_{i_{\ell}k_{\ell}}(t)$ as well. The solution becomes

[TABLE]

where $r_{i}^{-1}(t,{x})={{\hat{\tau}}_{i}(t,{x})}/{x_{i}(t)}$ is the outward throughflow intensity function.

The supplementary relationship between the simple direct and indirect subflows given in Eq. 55 can be expressed in matrix form, in terms of both flows and storages, as follows:

[TABLE]

A similar supplementary relationship can be formulated between the simple cycling and acyclic flows and storages:

[TABLE]

due to Eq. 58. The reflexivity of the simple cycling flows and storages, formulated in Eq. 56, can also be written in matrix form, as follows:

[TABLE]

The matrix relationships formulated in Eqs. 65, 66, and 67 can similarly be expressed for the composite diact flows and storages.

2.6 Static subsystem partitioning and diact transactions

The static version of the dynamic subsystem partitioning methodology given in Eqs. 43 and 44 has recently been formulated by [11]. This static partitioning is summarized in this section.

Since time derivatives are zero at steady state, we set $\dot{x}^{w}_{n_{k}\ell_{k}j_{k}}(t)=0$ in Eq. 44. The static transient outflow at subcompartment $\ell_{k}$ along subflow path $p^{w}_{n_{k}j_{k}}=i_{k}\mapsto j_{k}\to\ell_{k}\to n_{k}$ from $j_{k}$ to $n_{k}$, $f^{w}_{n_{k}\ell_{k}j_{k}}$, and the associated transient substorage generated in $\ell_{k}$, $x^{w}_{n_{k}\ell_{k}j_{k}}$, by the transient inflow, $f^{w}_{\ell_{k}j_{k}i_{k}}$, are then formulated as follows:

[TABLE]

The second equality for $f^{w}_{n_{k}\ell_{k}j_{k}}$ is obtained by using the first equation for ${x}^{w}_{n_{k}\ell_{k}j_{k}}$ in the first equality. The relationships in Eq. 68 establish the foundation of the static subsystem partitioning.

Through the system partitioning methodology, the static diact flows and storages are formulated in matrix form by [11], as presented in Table 2. Note that the matrix $\mathcal{N}$ used in the table is defined as $\mathcal{N}=\operatorname{diag}(N)$.

2.7 System measures and indices

The dynamic system partitioning methodology yields the subthroughflow and substorage vectors and matrices that measure the influence of the initial stocks and environmental inputs on system compartments in terms of the flow and storage generation. These vector and matrix measures enable tracking the evolution of the initial stocks, environmental inputs, as well as the associated storages sourced from these stocks and inputs individually and separately within the system. For the analysis of intercompartmental flow and storage dynamics, the dynamic subsystem partitioning methodology then formulates the transient and dynamic diact flows and associated storages. The transient and diact transactions enable tracking of arbitrary intercompartmental flows and storages along a particular and all possible flow paths within the system and determining the influence of system compartments directly and indirectly on one another.

The proposed methodology constructs a foundation for the development of new mathematical system analysis tools as quantitative ecosystem indicators. In addition to the measures summarized above, multiple novel dynamic and static system analysis tools of matrix, vector, and scalar types have recently been introduced in separate works, based on the proposed methodology [9, 12]. More specifically, these manuscripts introduce measures and indices for the diact effects, utilities, exposures, and residence times, as well as the system efficiencies, stress, and resilience.

2.8 Quantitative classification of interspecific interactions

An immediate ecological application of the proposed methodology is the quantitative analysis of food webs and chains. The system compartments of a food web ecosystem represent species, the conserved quantity in question becomes nutrient or energy, and flow paths correspond to food chains in the web. In this setting, direct flow rate $f_{ik}(t,x)$ represents an interspecific interaction, such as predation, between species $i$ and $k$ and measures the rate of nutrient or energy flow from species $k$ in a lower trophic level to $i$ in the next level at time $t$. Nutrient or energy stored in species $i$ through all trophic interactions is represented by $x_{i}(t)$.

Community ecology classifies interspecific interactions qualitatively using network topology without regard for system flows [36, 35]. This structural determination, however, gets more complicated, if at all possible, with the increasing complexity of intricate food webs [25, 50]. Multiple food chains of potentially different lengths between two species, for example, disallow the classification based on the length of the chains [6]. A mathematical characterization and classification technique for the analysis of the nature and strength of food chains has recently been proposed by [11] for static systems. This section introduces a dynamic version of this technique with slight modifications.

The proposed methodology can quantitatively determine the net benefit in terms of flow and storage transfers received by the involved species from each other. The sign analysis of the diact interspecific interactions determines the neutral and antagonistic nature of the interactions—whether the interaction is beneficial or harmful to the species involved. The strength analysis then quantifies the strength of these diact interactions. The sign and strength of the diact interactions induced by environmental inputs between species $i$ and $j$ will be defined respectively as follows:

[TABLE]

where $\operatorname{sgn}(\cdot)$ is the sign function, and the superscript (${}^{\texttt{*}}$) represents any of the diact symbols. Following the convention of community ecology, instead of $(+1)$ and $(-1)$, $(+)$ and $(-)$ notations will be used for the sign of the diact interactions. The strength, $0\leq\mu^{\texttt{*}}_{ij}(t)\leq 1$, is defined to be zero, if both terms in its denominator are zero.

For the analysis of diact interactions ranging from the individual and local to the system-wide and global scale, the strength of the interactions can be formulated with the normalization by $\tau^{\texttt{*}}_{ij}(t)+\tau^{\texttt{*}}_{ji}(t)$, $\tau^{\texttt{t}}_{ij}(t)+\tau^{\texttt{t}}_{ji}(t)$, $\check{\tau}_{i}(t,x)+\check{\tau}_{j}(t,x)$, as in Eq. 69, or $\check{\sigma}^{\tau}(t)=\bm{1}^{T}\,\check{\tau}(t,x)$ in the given order, where $\check{\sigma}^{\tau}(t)$ will be called the total inward system throughflow. As an example, for the analysis of local diact interactions at the global scale, $\mu^{\texttt{*}}_{ij}(t)$ can be defined as $\mu^{\texttt{*}}_{ij}(t)={|\tau^{\texttt{*}}_{ij}(t)-\tau^{\texttt{*}}_{ji}(t)|}/{\check{\sigma}^{\tau}(t)}$.

The diact neutral relationship between species $i$ and $j$ and “predation” of species $i$ on $j$ can quantitatively be characterized, respectively, as follows:

[TABLE]

This classification can also be extended to the transient interactions between two species along a given food chain using the transient inflows. We will use the notations $\delta^{w}_{ij}(t)$ and $\mu^{w}_{ij}(t)$ for the sign and strength of the net flow from compartment $j$ indirectly to $i$ through trophic interactions along a given food chain, $p^{w}_{ij}$.

The classification of the diact interactions induced by the initial stocks can be determined by using the composite diact subflows for the initial subsystem, $\tau^{\texttt{*}}_{i_{0}j_{0}}(t)$, in Eqs. 69 and 70. The storage-based quantitative definition of the diact interactions can be formulated in parallel by substituting the diact storages for the corresponding diact flows in the definitions above as well. The storage-based formulations represent the history of the diact interactions during $[t_{1},t]$ while the flow-based formulations represent simultaneous interactions at time $t$. We will use the superscript $x$ to distinguish the storage-based measures, $\delta^{\texttt{*},x}_{ij}(t)$ and $\mu^{\texttt{*},x}_{ij}(t)$. For the classification of the diact interactions induced by individual environmental inputs, the simple diact flows and storages can be used instead of their composite counterparts in Eqs. 69 and 70. A tilde notation will be used over the measures for simple diact interspecific interactions, $\tilde{\delta}^{\texttt{*}}_{ij}(t)$ and $\tilde{\mu}^{\texttt{*}}_{ij}(t)$.

A mathematical technique for the dynamic characterization and classification of the main interspecific interaction types, such as neutralism, mutualism, commensalism, competition, and exploitation, has also been developed recently in a separate paper [9].

3 Results

The proposed dynamic methodology is applied to a linear and nonlinear dynamic ecosystem models. Numerical results for the system analysis tools developed in this manuscript, such as the substorage and subthroughflow matrix measures, as well as the transient and dynamic diact flows and storages, are presented in this section.

The results indicate that the proposed methodology precisely quantifies system functions, properties, and behaviors, enables tracking the evolution of the initial stocks, environmental inputs, and intercompartmental flows, as well as associated storages individually and separately within the system, is sensitive to perturbations due to even a brief unit impulse, and, thus, can be used for rigorous dynamic analysis of nonlinear ecological systems. It is worth noting, however, that this present work proposes a mathematical method—a systematic technique designed for solving and analyzing any nonlinear dynamic compartmental model—and it is not itself a model. Therefore, we focus more on demonstrating the efficiency and wide applicability of the method. It is expected that once the method is accessible to a broader community of environmental ecologists, it can be used for ecological inferences and the holistic analysis of specific models of interest.

3.1 Case study

A linear dynamic ecosystem model introduced by [23] is analyzed through the proposed methodology in this case study to demonstrate the capability of the method to analytically solve linear systems with time-dependent inputs. The graphical representations of the results are presented.

The model has two compartments, $x_{1}(t)$ and $x_{2}(t)$ (see Fig. 6). The system flows are described as

[TABLE]

The governing equations take the following form:

[TABLE]

with the initial conditions $[x_{1,0},x_{2,0}]^{T}=[3,3]^{T}$.

The subcompartmentalization step yields the substate variables that represent the substorage values as follows:

[TABLE]

The flow partitioning then yields the subflows for the subsystems:

[TABLE]

where the decomposition factors $d_{i_{k}}({\rm x})$ are defined by Eq. 8. The dynamic system partitioning methodology then yields the following governing equations for the decomposed system:

[TABLE]

with the initial conditions

[TABLE]

for $i=1,2$. There are $n\times(n+1)=2\times 3=6$ equations in the system.

The system can be written in matrix form as formulated in Eq. 30:

[TABLE]

where the constant flow intensity matrix $A$, as defined in Eq. 35, becomes

[TABLE]

The governing decomposed system, Eq. 72, is linear. We can, therefore, solve it analytically, as formulated in Section 2.3. Since the flow intensity matrix $A$ is constant, we have the following fundamental matrix solution as given in Eq. 32:

[TABLE]

For $z=[1,1]^{T}$, the solutions for the matrix equation, Eq. 72, then become

[TABLE]

as given in Eq. 33. Therefore, the solution to the original system, Eq. 71, in vector form, is

[TABLE]

The subthroughflow matrices, $\check{T}(t)$ and $\hat{T}(t)$, can also be expressed as formulated in Eq. 15, using the solution for the substorage matrix, $X(t)$.

The steady state solutions can also be computed as formulated in Eq. 36:

[TABLE]

It can easily be seen that, this steady-state solution is the same as the limit of the dynamic solution, Eq. 73, as $t$ tends to infinity. That is, $\lim_{t\rightarrow\infty}X(t)=X$.

We also analyze the system with a time dependent input $z(t)=[3+\sin(t),3+\sin(2t)]^{T}$. Similar computations with the same fundamental matrix, $V(t)$, lead us to the following initial substorage vector, $x_{i_{0}}(t)$, and substorage matrix components, $x_{i_{k}}(t)$:

[TABLE]

Using these solutions, we can express the solutions to the original system, Eq. 71, as

[TABLE]

The proposed dynamic method solves linear systems analytically. Explicit solutions can be used to compute the quantities in question at any time $t$.

The elements of the inward initial throughflow vector, $\check{\tau}_{0}(t)$, and subthroughflow matrix, $\check{T}(t)$, can be computed using Eq. 15 as follows:

[TABLE]

The outward initial throughflows and throughflows can also be obtained similarly, using Eq. 15. The substorage and inward subthroughflow matrix functions, $X(t)$ and $\check{T}(t)$, given in Eqs. 75 and 76, determine the dynamic distribution of the environmental inputs and the organization of the associated storages generated by these inputs individually and separately within the system. In other words, using these functions, the evolution of the environmental inputs and associated storages can be tracked individually and separately throughout the system. The graphical representation of $X(t)$ and $\check{T}(t)$ for time dependent input $z(t)=\left[3+\operatorname{sin}(t),3+\operatorname{sin}(2\,t)\right]^{T}$ are depicted in Fig. 7. The evolution of initial stocks and inward throughflows, $x_{i_{0}}(t)$ and $\check{\tau}_{i_{0}}(t)$, are also presented in Fig. 7.

The residence time matrix for this model, defined in Eq. 28, becomes

[TABLE]

The residence time of compartment $2$ is constantly smaller than that of compartment 1. That is, $r_{2}(t,x)=0.43<0.6=r_{1}(t,x)$. This result ecologically indicates that compartment $2$ is more active than $1$.

The subsystem partitioning methodology allows for the further analysis of the system and brings out additional insights that are not available through the state-of-the-art techniques. The composite transfer flow and storage from compartment $2$ to $1$, $\tau^{\texttt{t}}_{12}(t)$ and $x^{\texttt{t}}_{12}(t)$, and composite transfer subflow and substorage from initial subcompartment $2_{0}$ to $1_{0}$, $\tau^{\texttt{t}}_{1_{0}2_{0}}(t)$ and $x^{\texttt{t}}_{1_{0}2_{0}}(t)$, are computed below as an application of the proposed subsystem partitioning methodology. They can be expressed using Eq. 49 as

[TABLE]

The sets of mutually exclusive subflow paths from $2_{k}$ to $1_{k}$, $P_{1_{k}2_{k}}$, for $k=0,1,2$, can be formulated as follows: $P_{1_{0}2_{0}}=\{p^{1}_{1_{0}2_{0}},p^{2}_{1_{0}2_{0}}\}$, $P_{1_{1}2_{1}}=\{p^{1}_{1_{1}2_{1}}\}$, $P_{1_{2}2_{2}}=\{p^{1}_{1_{2}2_{2}}\}$ where $p_{1_{0}2_{0}}^{1}=0_{0}\mapsto 1_{0}\rightsquigarrow 2_{0}\to 1_{0}$, $p_{1_{0}2_{0}}^{2}=0_{0}\mapsto 2_{0}\to 1_{0}\rightsquigarrow 2_{0}$, $p_{1_{1}2_{1}}^{1}=0_{1}\mapsto 1_{1}\rightsquigarrow 2_{1}\to 1_{1}$, and $p_{1_{2}2_{2}}^{1}=0_{2}\mapsto 2_{2}\to 1_{2}\rightsquigarrow 2_{2}$.

There are two subflow paths in the initial subsystem, $p_{1_{0}2_{0}}^{1}$ and $p_{1_{0}2_{0}}^{2}$, and therefore, $w_{0}=2$. The corresponding transfer subflow and associated substorage functions, as formulated in Eq. 48, become

[TABLE]

Similarly, we have

[TABLE]

for $k=1,2$, since there is only one subflow path in these subsystems ($w_{k}=1$).

The links that directly contribute to the cumulative transient inflow and substorage, $\check{\tau}^{1}_{1_{1}}(t)$ and $x^{1}_{1_{1}}(t)$, at subcompartment $1_{1}$ along $p^{1}_{1_{1}2_{1}}$ are marked with cycle numbers, $m$, in the extended subflow path diagram below:

[TABLE]

The cumulative transient inflow and substorage will be approximated by two terms ($m_{1}=2$) using Eq. 46:

[TABLE]

The governing equations, Eqs. 43 and 44, for the transient subflows and associated substorages, $f^{1,m}_{1_{1}2_{1}1_{1}}(t)$ and $x^{1,m}_{2_{1}1_{1}2_{1}}(t)$, and the other transient subflows and substorages involved in Eqs. 78 and 79, are solved simultaneously, together with the decomposed system, Eq. 10. Numerical results for the transfer subflows and associated substorages are presented in Fig. 8.

The subflow paths in $P_{1_{k}2_{k}}$ for each subsystem $k$ are mutually exclusive and exhaustive. Therefore, $x_{1}(t)$ and $x_{1_{0}}(t)+x_{2_{1}1_{1}0_{1}}(t)+x^{\texttt{t}}_{1_{0}2_{0}}(t)+x^{\texttt{t}}_{12}(t)$ must be the same, as well as $f_{12}(t)$ and $\tau^{\texttt{t}}_{1_{0}2_{0}}(t)+\tau^{\texttt{t}}_{12}(t)$. The terms added to $x^{\texttt{t}}_{1_{0}2_{0}}(t)+x^{\texttt{t}}_{12}(t)$ for a comparison, $x_{2_{1}1_{1}0_{1}}(t)$ and $x_{1_{0}}(t)$, are the transient substorage generated by environmental input in $1_{1}$ ($0_{1}\mapsto 1_{1}$) and the initial substorage in compartment 1 (see Fig. 7(a)). Therefore, they are not included in the transfer storage and initial substorage, $x^{\texttt{t}}_{12}(t)$ and $x^{\texttt{t}}_{1_{0}2_{0}}(t)$. These quantities, however, are approximately equal as presented in Fig. 8:

[TABLE]

The small differences are caused by the truncation errors in the computation of cumulative transient subflows, and larger ${m_{w}}$ values further improve the approximations. These close approximations demonstrate the accuracy and consistency of both the system and subsystem partitioning methodologies.

Instead of the path-based approach used in the numerical computations above, the diact flows can also be obtained analytically and explicitly using the dynamic approach as introduced in Section 2.5.2. The composite transfer subflows $\tau^{\texttt{t}}_{1_{k}2_{k}}(t)$ for $k=0,1,2$, and transfer flow $\tau^{\texttt{t}}_{12}(t)$ become:

[TABLE]

as formulated in Eq. 61. The composite transfer substorages can then be obtained by coupling Eq. 63 for the transfer subflows with the decomposed system, Eq. 10, and solving them simultaneously. Alternatively, the corresponding transfer substorages can be obtained analytically as formulated in Eq. 64:

[TABLE]

The graphs of these explicit transfer subflow and substorage functions in Eqs. 81 and  82 obtained through the dynamic approach are exactly the same as the ones obtained by numerical computation through the path-based approach, Eq. 77, as depicted in Fig. 8.

The cycling flows and the associated storages generated by these flows are also calculated below for both compartments. The sets of mutually exclusive subflow paths from subcompartment $k_{k}$ to $1_{k}$ with a closed subpath at $1_{k}$, $P^{\texttt{c}}_{1_{k}k_{k}}$, are given as $P^{\texttt{c}}_{1_{0}0_{0}}=\{p^{1}_{1_{0}1_{0}},p^{2}_{1_{0}2_{0}}\}$, $P^{\texttt{c}}_{1_{1}1_{1}}=\{p^{1}_{1_{1}1_{1}}\}$, $P^{\texttt{c}}_{1_{2}2_{2}}=\{p^{1}_{1_{2}2_{2}}\}$, where $p_{1_{0}1_{0}}^{1}=0_{0}\mapsto 1_{0}\rightsquigarrow 2_{0}\to 1_{0}$, $p_{1_{0}2_{0}}^{2}=0_{0}\mapsto 2_{0}\rightsquigarrow 1_{0}\rightsquigarrow 2_{0}\to 1_{0}$, $p_{1_{1}1_{1}}^{1}=0_{1}\mapsto 1_{1}\rightsquigarrow 2_{1}\to 1_{1}$, and $p_{1_{2}2_{2}}^{1}=0_{2}\mapsto 2_{2}\rightsquigarrow 1_{2}\rightsquigarrow 2_{2}\to 1_{2}$. For the subflow paths in $P^{\texttt{c}}_{1_{0}0_{0}}$, the composite cycling subflows are derived from the initial stocks, and for the ones in $P^{\texttt{c}}_{1_{1}1_{1}}$ and $P^{\texttt{c}}_{1_{2}2_{2}}$, the simple cycling flows are generated by the respective environmental inputs of $z_{1}(t)$ and $z_{2}(t)$. The sets of subflow paths $P^{\texttt{c}}_{2_{k}k_{k}}$ for $k=0,1,2$ can similarly be defined.

The simple cycling subflow at subcompartment $1_{2}$ along the only subflow path ($w_{2}=1$) in subsystem $2$, $p^{1}_{1_{2}2_{2}}\in P^{\texttt{c}}_{1_{2}2_{2}}$, and associated substorage are

[TABLE]

as formulated in Eq. 48. The links contributing to the cycling subflow along the path are marked with cycle numbers in the extended subflow diagram below:

[TABLE]

Note that the first flow entrance into $1_{2}$ is not considered as cycling flow. The cumulative transient inflow $\check{\tau}^{1}_{1_{2}}(t)$ and substorage $x^{1}_{1_{2}}(t)$ can be approximated by two terms ($m_{1}=2$) along the closed subpath as formulated in Eq. 46:

[TABLE]

The governing equations for the transient subflows and associated substorage functions, $f^{w,m}_{1_{2}2_{2}1_{2}}(t)$ and $x^{w,m}_{2_{2}1_{2}2_{2}}(t)$, as well as the other transient subflows and substorages involved in Eq. 83, as formulated in Eqs. 43 and 44, are coupled and solved simultaneously together with the decomposed system, Eq. 10. The numerical results for the composite cycling flows and associated storages induced both by the environmental inputs and initial stocks,

[TABLE]

for $i=1,2$, are presented in Fig. 8(c).

Note that, due to the reflexivity of cycling flows, the same computations can be done more practically in only two steps using the sets of closed subflow paths, $P^{\texttt{c}}_{i_{k}}$, instead, with the local inputs being the corresponding outwards subthroughflows. The subflow path $p_{1_{1}}^{1}=1_{1}\mapsto 1_{1}\rightsquigarrow 2_{1}\to 1_{1}$ in subcompartment $1_{1}$ with local input $\hat{\tau}_{1_{1}}(t,{\mathbf{x}})$, for example, is depicted in Fig. 6. The cycling flows can also be computed along closed paths at the compartmental level, where the local inputs are the outwards throughflows.

The composite cycling subflows can also be computed analytically through the dynamic approach as formulated in Eq. 57. As examples, $\tau^{\texttt{c}}_{1_{0}1_{0}}(t)$ and $\tau^{\texttt{c}}_{2_{2}}(t)=\tau^{\texttt{c}}_{2_{2}2_{2}}(t)$ become

[TABLE]

The composite cycling storages can then be obtained by coupling Eq. 63 for the cycling flows and storages with the decomposed system, Eq. 10, and solving them simultaneously. Alternatively, they can be obtained analytically by using Eq. 64, similar to the transfer storages presented above in this example. Because of the lengthy analytical formulations of the other cycling subflows and substorages, only $\tau^{\texttt{c}}_{1_{0}1_{0}}(t)$ and $\tau^{\texttt{c}}_{2_{2}}(t)$ are presented in Eq. 84 as examples.

3.2 Case study

In this section, a nonlinear resource-producer-consumer ecosystem model introduced by [19] is analyzed through the proposed methodology. A comparison of the results is not possible, as the authors did not provide any computational or explicit results in the article. They only provided some results at steady state. Besides a constant environmental input, the system is also examined for a time dependent, symmetric Gaussian impulse to illustrate the efficiency of the proposed method in capturing the system response to disturbances. Such analysis can be used to quantify the system resistance and resilience in the face of disturbances and perturbations.

The resource-producer-consumer model by [19] consists of the dynamics for three components: $x_{1}(t)=r(t)$ is the nutrient storage (such as phosphorus or nitrogen) present at time t; $x_{2}(t)=s(t)$ represents the nutrient storage in the producer (such as phytoplankton) population; and $x_{3}(t)=c(t)$ denotes the nutrient storage in the consumer (such as zooplankton) population (see Fig. 9). The conservation of nutrient is the basic model assumption. The system flows are described as follows:

[TABLE]

where the constant input is $z(t)=[1,1,1]^{T}$, and the parameters are given as

[TABLE]

The value for $\alpha_{1}$ was not provided in [19] and was chosen arbitrarily for this example. The governing equations take the following form:

[TABLE]

with the initial conditions of $[r_{0},s_{0},c_{0}]=[1,1,1]$.

The system partitioning methodology is composed of the subcompartmentalization and flow partitioning components. The subcompartmantalization yields

[TABLE]

The flow partitioning then gives the flow regime for each subsystem as follows:

[TABLE]

where $F_{k}$, $\check{z}_{k}$, and $\hat{y}_{k}$ describe the $k^{th}$ direct flow matrix, input, and output vectors for the $k^{th}$ subsystem, and the decomposition factors $d_{i_{k}}({\rm x})$ are defined by Eq. 8. Therefore, the dynamic system partitioning methodology yields the following governing equations for the decomposed system:

[TABLE]

with the initial conditions

[TABLE]

for $i=1,\ldots,3$. There are $n\times(n+1)=3\times 4=12$ equations in this system.

The system is solved numerically and the graphs for selected elements of the substorage and subthroughflow matrices are depicted in Fig. 10. As seen from the graphs, the system converges to a steady-state at about $t\approx 6$. The results show, for example, that the nutrient storage in the resource compartment ($i=1$) derived from nutrient input into the consumer compartment ($i=3$), $x_{1_{3}}(t)$, increases from [math] to $0.64$ units until the system reaches the steady state, while the initial nutrient storage, $x_{1_{0}}$, first increases from $1$ to $1.38$ units and then vanishes. The throughflow into the resource compartment generated by nutrient input into the producer compartment ($i=2$), $\check{\tau}_{1_{2}}(t,{\rm x})$, increases until about $t\approx 2$. The outward throughflow at the same subcompartment, $\hat{\tau}_{1_{2}}(t,{\rm x})$, is slightly smaller than inward throughflow, $\check{\tau}_{1_{2}}(t,{\rm x})$, but has a similar behavior. As seen from these results, the distribution of environmental nutrient inputs and the organization of the associated nutrient storages generated by the inputs can be analyzed individually and separately within the system.

In general terms, the state variable $x_{i}(t)$ of the original system for the resource-producer-consumer dynamics, Eq. 85, gives the nutrient storage in compartment $i$ at time $t$ based on its initial stock, $x_{i}(t_{0})$. It cannot be used to distinguish the nutrient storage derived from individual environmental nutrient inputs. On the other hand, the state variable $x_{i_{k}}(t)$ of the decomposed system, Eq. 86, represents the nutrient storage in compartment $i$ that is derived from the specific environmental nutrient input into compartment $k$, $z_{k}(t)$. Similarly, the state variable ${x}_{i_{0}}(t)$ of the decomposed system represents the dynamics of the initial nutrient stocks in compartment $i$. Parallel interpretations are possible for the inward and outward throughflow functions of the original system, $\check{\tau}_{i}(t,x)$ and $\hat{\tau}_{i}(t,x)$, and the inward and outward subthroughflow functions of the decomposed system, $\check{\tau}_{i_{k}}(t,{\rm x})$ and $\hat{\tau}_{i_{k}}(t,{\rm x})$, as well.

The proposed dynamic system partitioning methodology, consequently, enables partitioning the compartmental composite nutrient flows and storages into subcompartmental segments based on their constituent sources from the initial stocks and environmental inputs. In other words, the system partitioning enables tracking the evolution of the initial nutrient stocks and environmental nutrient inputs a well as the associated storages generated by the stocks and inputs individually and separately within the system. This partitioning also allows for compiling a history of compartments visited by individual nutrient inputs separately.

The system is also perturbed with a Gaussian input $z_{2}(t)={\mathrm{e}}^{\frac{-(t-15)^{2}}{2}}+0.1$, which represents a brief, unit local impulse at about $t=15$ to demonstrate the capability of the proposed method to analyze the influence of time dependent inputs on the system. The other two environmental nutrient inputs are kept constant as before for a comparison, that is, $z_{1}(t)=z_{3}(t)=1$. The graphical representations for the selected elements of the substorage and subthroughflow matrices are given in Fig. 11. It is clear from the graphs that the dynamic substorage and subthroughflow matrix measures reflect the impact of the unit impulse at about $t=15$. Note that, the system completely recovers after the disturbance in about $10$ time units. This time interval can be taken as a quantitative measure for the restoration time or system resilience. Therefore, the proposed measures can be used as quantitative ecological indicators for various ecosystem characteristics and behaviors.

The subsystem partitioning methodology is also applied to this model to track the fate of arbitrary intercompartmental flows and the associated storages generated by these flows within the subsystems. Along the subflow path $p_{3_{1}1_{1}}^{1}=0_{1}\mapsto 1_{1}\to 2_{1}\rightsquigarrow 3_{1}$ from subcompartment $1_{1}$ to $3_{1}$ in subsystem $1$, the transient subflows and associated substorages are computed as formulated in Eqs. 43 and 44. The numerical results for the transient subflows, $f^{1}_{2_{1}1_{1}0_{1}}(t)$ and $f^{1}_{3_{1}2_{1}1_{1}}(t)$, and associated substorage functions, $x^{1}_{2_{1}1_{1}0_{1}}(t)$ and $x^{1}_{3_{1}2_{1}1_{1}}(t)$, are presented in Fig. 12(b) and 12(a).

The subflow path $p_{3_{1}1_{1}}^{1}$ is extended to path

[TABLE]

to compute the local output $f^{1}_{0_{1}1_{1}2_{1}}(t)$ (a segment of environmental output $y_{1}(t)$) derived from the local (and environmental) input $z_{1}(t)=1$ along that particular path (see Fig. 9). That is, the fate of $z_{1}(t)$ along path $p_{0_{1}1_{1}}^{1}$ within the system is determined. The corresponding transient subflow and associated substorage functions at each step (subcompartment) along the path are also presented in Fig. 12(d) and 12(c). Since $f^{1}_{0_{1}1_{1}2_{1}}(t)\leq 6.28\times 10^{-5}$, at most about $\%0.006$ of $z_{1}(t)$ exits the system through the given subflow path $p^{1}_{0_{1}1_{1}}$ at any time $t$.

These results indicate that the proposed dynamic subsystem partitioning methodology enables dynamically tracking the fate of an arbitrary amount of nutrient flow and associated nutrient storage along a given flow path. Consequently, the spread of an arbitrary amount of nutrient from one compartment to the entire system can be monitored. Moreover, the effect of one compartment on any other in terms of the nutrient transfer, through not only direct but also indirect interactions, can be determined.

The diact flows and storages are introduced in Section 2.5.2. The indirect flow and storage from compartment 1 to 3, $\tau^{\texttt{i}}_{31}(t)$ and $x^{\texttt{i}}_{31}(t)$, the acyclic flow and storage from compartment 2 to 1, $\tau^{\texttt{a}}_{12}(t)$ and $x^{\texttt{a}}_{12}(t)$, and the cycling flow and storage at compartment 1, $\tau^{\texttt{c}}_{11}(t)$ and $x^{\texttt{c}}_{11}(t)$, generated by the environmental inputs, as well as the corresponding initial subflows and substorages derived from the initial stocks transmitted in the same directions are depicted in Fig. 13. As seen from the graphs, all initial diact subflows and substorages vanish as the system converges to a steady-state and, then, the system behavior is eventually dominated by the environmental inputs. Ecologically, the acyclic flow and storage, $\tau^{\texttt{a}}_{12}(t)$ and $x^{\texttt{a}}_{12}(t)$, represent the nutrient flow at time $t$ and the associated nutrient storage generated by this flow during $[t_{1},t]$ that visit the resource compartment only once—do not return to this compartment for a second time later—after being directly or indirectly transmitted from the producer compartment. The initial acyclic subflow and substorage, $\tau^{\texttt{a}}_{1_{0}2_{0}}(t)$ and $x^{\texttt{a}}_{1_{0}2_{0}}(t)$, represent the same phenomena within the initial subsystem. Similarly, the indirect flow and storage, $\tau^{\texttt{i}}_{31}(t)$ and $x^{\texttt{i}}_{31}(t)$, represent the nutrient flow and storage transmitted indirectly from the resource compartment through the producer to the consumer compartment. The cycling flow and storage, $\tau^{\texttt{c}}_{11}(t)$ and $x^{\texttt{c}}_{11}(t)$, represent the nutrient flow and storage transmitted indirectly from the resource compartment through other compartments back into itself. The other diact (sub)flows and (sub)storages can be interpreted similarly, for both the subsystems and initial subsystem to analyze the intercompartmental dynamics generated respectively by the environmental inputs and initial stocks, individually and separately.

The residence time matrix is another novel mathematical measure proposed for quantitative system analysis [10, 9]. The $i^{th}$ diagonal element of $\mathcal{R}(t,{x})$ at time $t_{1}$, ${r}_{i}(t_{1},{x})$, can be interpreted as the time required for the outward throughflow, at the constant rate of $\hat{\tau}_{i}(t_{1},{x})$, to completely empty compartment $i$ with the storage of $x_{i}(t_{1})$. The diagonal structure of the residence time matrix indicates that all subcompartments of compartment $i$ vanish simultaneously. The residence times measure compartmental activity levels [11]. The smaller the residence time the more active the corresponding compartment. The derivative of the residence time matrix is called the reverse activity rate matrix [9].

The residence time functions for this model with both the constant and time-dependent environmental inputs are depicted in Fig. 14, for a comparison. The residence times of both the consumer and producer compartments are almost constant and the same in both cases. Interestingly, the Gaussian impulse at the producer compartment, $z_{2}(t)$, has no significant impact on the activity level of the consumer and even that of the producer compartment itself. However, the decrease in the input into the producer compartment from constant $z_{2}(t)=1$ to $z_{2}(t)={\mathrm{e}}^{\frac{-(t-15)^{2}}{2}}+0.1$ results in an overall increase in the residence time of the resource compartment (and all of its subcompartments) from the steady state value of ${r}_{1}=0.87$ days to ${r}_{1}=0.98$ days. Moreover, the maximum impulse at $t=15$ decreases this residence time, $r_{1}(t)$, locally in time. That is,

[TABLE]

Consequently, the residence time of the resource compartment adversely impacted by the environmental input into the producer compartment.

The mathematical classification of the diact interspecific interactions is also introduced in Section 2.8. The sign and strength of the diact interactions, induced by environmental inputs, between the producer and consumer compartments become

[TABLE]

for the storage-based analysis. The numerical results for the net diact storage transactions and the strengths of the diact interactions are presented in Fig. 15.

As seen from the sign analysis in Fig. 15(a),

[TABLE]

These results indicate that the diact interactions induced by the environmental inputs between the producer and consumer compartments are all antagonistic. Although the consumer compartment directly benefits from the producer compartment as expected, interestingly, their indirect, cycling, acyclic, and total interactions are detrimental to the consumer compartment. The strengths of the diact iterations are ordered as follows:

[TABLE]

for $t\geq 0.39$. Therefore, the indirect interaction between the producer and consumer compartments is the strongest of all diact interactions. Since the indirect interaction dominates the direct interaction, their overall interactions is counterintuitively detrimental to the consumer compartment after $t=0.39$.

The detailed information and inferences enabled by the proposed methodology cannot be obtained through the analysis of the original system by the state-of-the-art techniques, as demonstrated in these case studies.

4 Discussion

Environment is not an easy concept to define in general and, in particular, to analyze mathematically. One reason for this is that nature is always on the move and ecological systems struggle to adapt to constantly changing circumstances. Although sound rationales are offered in the literature for the analysis of natural system dynamics, they are only for special cases, such as linear models and static systems. In recent decades, there has been several attempts to analyze dynamic ecological networks, but each of them bears disadvantages. The need for dynamic and nonlinear methodologies has always existed. This manuscript proposes a novel mathematical methodology for the analysis of nonlinear dynamic compartmental systems to comprehensively address these shortcomings.

Considering a hypothetical ecosystem with several interacting species for which the effect of a specific pollutant needs to be investigated, monitoring the evolution of that pollutant within this food web would be critical to addressing the potential harm. Given the initial pollutant stocks in each of the species, current deterministic mathematical methods can analyze the composite throughflow and storage of the toxin in the species. The evolution of each environmental pollutant input separately within the web, however, cannot be determined through the current methodologies. In the case that multiple species exposed to the same pollutant from the environment, the proposed system partitioning methodology enables dynamically partitioning the composite pollutant flow at and storage in any species into subcompartmental segments based on their constituent sources from the initial pollutant stocks and environmental pollutant inputs. In other words, the system partitioning enables tracking the evolution of the initial pollutant stocks and environmental pollutant inputs in each species, as well as the associated pollutant storages derived from these inputs and stocks individually and separately within the food web.

The proposed subsystem partitioning methodology can then dynamically track the fate of arbitrary intercompartmental pollutant flows and associated storages in each species along a given food chain in the web as well. Therefore, the spread of an arbitrary amount of toxin from one species to the entire web or along a specific food chain can be monitored. Such information can help, for example, to identify critical pathways in which intervention helps to alleviate the effects of the pollutant. The direct, indirect, acyclic, cycling, and transfer (diact) flows and associated storages of the pollutant from one species directly or indirectly to any other—including itself—can also be determined.

More technically, the proposed system partitioning methodology explicitly generates mutually exclusive and exhaustive subsystems. Except the initial subsystem—which is driven by the initial stocks—each subsystem is driven by a single environmental input. The subsystems are running within the original system and have the same structure and dynamics as the system itself, except for their initial stocks and environmental inputs. The system partitioning yields the subthroughflows and substorages that represent flows and storages derived from the initial stocks and generated by environmental inputs in each compartment. That is, the composite compartmental storages and throughflows, $x_{i}(t)$ and $\tau_{i}(t)$, are dynamically partitioned into the subcompartmental substorage and subthroughflow segments, $x_{i_{k}}(t)$ and $\tau_{i_{k}}(t)$, based on their constituent sources from the initial stocks and environmental inputs, $x_{i_{0}}$ and $z_{k}(t)$. Equipped with these measures, the system partitioning ascertains the dynamic distribution of the initial stocks and environmental inputs and the organization of the associated storages derived from these stocks and inputs within the system. In other words, the system partitioning enables dynamically tracking the evolution of the initial stocks and environmental inputs individually and separately within the system. The system partitioning methodology, therefore, refines system analysis from the current static, linear, compartmental to the dynamic, nonlinear, subcompartmental level to explore the full complexity of the ecological systems.

The subsystems are then further decomposed into subflows and substorages along a set of mutually exclusive and exhaustive directed subflow paths. The subsystem partitioning methodology yields the transient subflows and substorages in each subcompartment along a given subflow path within a subsystem, generated by or derived from an arbitrary intercompartmental flow or storage. Therefore, arbitrary composite intercompartmental flows and associated storages can be apportioned dynamically into transient subflow and substorage segments along a given set of subflow paths. That is, the transient subflows and substorages determine the dynamic distribution of arbitrary intercompartmental flows and the organization of the associated storages generated by these flows within the subsystems. As a result, the spread of an arbitrary flow or storage segment from one compartment to the entire system can be monitored. Moreover, an archive of compartments visited by arbitrary system flows and storages can be also compiled. In brief, the subsystem partitioning methodology enables tracking the fate of arbitrary intercompartmental flows and associated storages within the subsystems.

The proposed mathematical method, as a whole, enables the decomposition to the utmost level, or “atomization,” of the system flows and storages. In addition to the subthroughflows, substorages, transient flows and storages, the dynamic direct, indirect, acyclic, cycling, and transfer (diact) flows and storages from one compartment directly or indirectly to any other—including itself—are also systematically formulated for the quantification of intercompartmental flow and storage dynamics. The diact flows and storages are derived explicitly through the dynamic and path-based approaches, which are based on the system and subsystem partitioning methodologies, respectively. As an immediate application, a mathematical technique that characterizes and classifies the neutral and antagonistic nature of diact interspecific interactions in food webs and determines the strength of these interactions is also developed based on the diact flows and storages. The illustrative case studies in Section 3 demonstrate the rigor and efficiency of these mathematical system analysis tools as ecological system indicators.

For a comparison of the proposed methodology with the state-of-the-art techniques, we first note that, at a steady state, the proposed dynamic methodology agrees with the current techniques for static ecological network analysis, as shown by [11]. In recent decades, there have been several attempts to analyze dynamic ecological networks. The first actual dynamic analysis was limited to linear systems with time-dependent input [23]. The proposed method is applied to the linear ecosystem model introduced by [23] as an illustrative case study in Section 3. It is shown that the analytic solutions obtained by the proposed methodology agree with those obtained by Hippe’s approach. Further results that are not available through Hippe’s approach, such as the diact flows and storages, are also presented for this linear system.

The dynamic approach is extended from linear to nonlinear systems by [19]. The authors provided, however, only closed-form, abstract formulations that are difficult to apply to real cases. The proposed methodology is also applied to the nonlinear ecosystem model analyzed by [19], and the results, together with their ecological interpretations, are presented in Section 3. A comparison of our results with the ones provided by the authors was not possible because, unlike the comprehensive dynamic analysis enabled by the proposed methodology for nonlinear systems, the authors could only provide asymptotic solutions to the model at steady state through their methodology.

Individual-based algorithms that rely on particle tracking simulations are also proposed for dynamic nonlinear ecological models in the literature. A truncated infinite series formulation, for example, was proposed by [41]. However, the authors’ approach was approximate, computationally resource-intensive, and offered no guarantees of series convergence. While a guarantee of convergence is added by [29], the computation remained resource-intensive due to the individual-based simulation technique [30, 43].

This is the first manuscript in the literature that comprehensively addresses all the previously identified problems and shortcomings. The method’s primary limitation is that it is designed for the analysis of conservative models defined in Eq. 5. Since the conservation principles are fundamental laws of nature, a large class of real-world problems are formulated based on conservation principles in many fields. On the other hand, there are still various non-conservative systems that cannot be analyzed by the proposed methodology in its current form.

The proposed methodology can easily be extended for similar analyses to systems of higher order ordinary and partial differential equations whose source terms governing the intercompartmental interactions are in the form of the conservative compartmental systems, as defined in Eq. 5. Such an extension of the system partitioning methodology to partial differential equations enables spatiotemporal analysis of ecological systems.

5 Conclusions

In the present paper, we developed a comprehensive mathematical method for the analysis of nonlinear dynamic compartmental systems through the system decomposition theory. The proposed method is based on the novel analytical and explicit, mutually exclusive and exhaustive system and subsystem partitioning methodologies. While the proposed dynamic system partitioning provides the subthroughflow and substorage matrices to determine the distribution of the initial stocks and environmental inputs, as well as the organization of associated storages individually and separately within the system, the subsystem partitioning yields the transient flows and storages to determine the distribution of arbitrary intercompartmental flows and the organization of associated storages within the subsystems. Consequently, the evolution of the initial stocks, environmental inputs, and arbitrary intercompartmental system flows, as well as the associated storages derived from these stocks, inputs, and flows can be tracked individually and separately within the system. Moreover, the transient and the dynamic direct, indirect, acyclic, cycling, and transfer (diact) flows and associated storages transmitted along a given flow path or from one compartment directly or indirectly to any other within the system are systematically formulated to ascertain the intercompartmental dynamics.

Traditional ecology is still largely a descriptive empirical science. This narrows the field’s scope of applicability and compromises its capacity to deal with complex ecological networks. The proposed dynamic method enhances the strength and extends the applicability of the state-of-the-art techniques and provides significant advancements in theory, methodology, and practicality. It serves as a quantitative platform for testing empirical hypotheses, ecological inferences, and, potentially, theoretical developments. Therefore, this method has the potential to lead the way to a more formalistic ecological science. We consider that the proposed methodology brings a novel complex system theory to the service of urgent and challenging environmental problems of the day. Several case studies from ecosystem ecology are presented to demonstrate the accuracy and efficiency of the method.

The proposed methodology also lays groundwork for the development of new mathematical system analysis tools as quantitative ecological indicators. The time dependent nature of these quantities enables also their time derivatives and integrals to be formulated as novel system measures. Multiple such dynamic diact measures and indices of matrix, vector, and scalar types which may prove useful for environmental assessment and management are systematically introduced in a separate paper by [9].

Acknowledgments

The author would like to thank the members of the ecosystem ecology group in the Department of Mathematics at UGA, for useful discussions and valuable feedback on a prior draft of this paper during his visit in 2016. In particular, the author is indebted to Caner Kazanci for introducing the static environ theory to the author, identifying some open problems in the area, and reviewing a prior draft of the paper, and other group members, Bernard C. Patten for his detailed review of a prior draft of the paper and Malcolm Adams for useful discussions and comments. The author also thanks Hasan Coskun, Sangwon Suh, Stuart Borett, and Brian Fath for their review of a prior draft of the paper and helpful comments.