Exact Reduction of the Generalized Lotka–Volterra Equations via Integral and Algebraic Substitutions

Abstract

Systems of interacting species, such as biological environments or chemical reactions, are often described mathematically by sets of coupled ordinary differential equations. While a large number $\beta$ of species may be involved in the coupled dynamics, often only $\alpha <\beta$ species are of interest or of consequence. In this paper, we explored how to construct models that include only those given $\alpha$ species, but still recreate the dynamics of the original $\beta$-species model. Under some conditions detailed here, this reduction can be completed exactly, such that the information in the reduced model is exactly the same as the original one, but over fewer equations. Moreover, this reduction process suggests a promising type of approximate model—no longer exact, but computationally quite simple.

1. Introduction

Consider an environment in which a large number of species interact. This could be, for example, a chemical reaction (with reacting chemical species) or an ecological system (with interacting organisms). To model the behavior of interacting species in many such environments, we often use the generalized Lotka–Volterra equations—a set of coupled ordinary differential equations (ODEs) [1]. The generalized Lotka–Volterra (GLV) equations extend two-species predator-prey models common in undergraduate differential equations classes to any number of species, and all possible combinations of interspecific interactions. They are thus much more general than a predator-prey-type system: they provide a framework to describe the time dynamics of any number $\beta$ of interacting species, allowing for linear (growth rate) and quadratic (interaction) terms.

In many applications, a large number of species may play a role in the dynamical community behavior. However, from a modeling perspective, the species of interest or of consequence may be restricted to a much smaller subset of $\alpha$ species, where $\alpha <\beta$. This occurs in fields as diverse as ecology [2], epidemiology [3,4], chemical kinetics [5,6,7], and agriculture [8]. For example, crops and weeds coexist, but farmers are typically only concerned with the performance of the crops. There are good reasons to omit many species. First, a modeler may not have access to data about certain species’ concentrations, which would be needed to define initial conditions or calibrate additional model parameters. Second, the modeler may not actually know what species should be included. Third, for a real system, $\mathcal{O}\left({\beta }^2\right)$ interspecific interaction coefficients must be estimated, requiring as many experimental treatments. If $\beta$ is large, then this estimation problem may quickly become experimentally infeasible [9]. Fourth, the more species included in the model, then in general, the more computationally expensive the model becomes. Except in very special (and usually unrealistic) cases, a system of GLV equations will not admit a closed-form solution; instead, we solve these systems with computational models. Therefore, it is common to build reduced models that include only $\alpha <\beta$ given species.

One immediate question that arises in this context is the following: Given a system of $\beta$ species, suppose only $\alpha$ are known, or of interest. What is the best reduced deterministic model, in terms of only those given $\alpha$ species? Of course, to answer this question, we must define what is meant by “best.” To do so, let us first consider the landscape of computational models—one of science, mathematics, and computation. We will broadly classify this consideration into three major areas: model validation, model reduction, and computational implementation, and what the major questions are in each.

I

Model validation: Does the mathematical model adequately represent the scientific system in question?

II

Model reduction: How much error is incurred by the use of the reduced model compared to the original, or high-fidelity, model?

III

Computational implementation: Is the reduced model less computationally expensive than the original model and, if so, by how much?

(Point III is closely related to model verification, a process that checks that the computational model solves the mathematical model correctly [10].) Returning to the question of the best reduced model: the best reduced model would address all three areas, i.e., well represent the scientific system under study (I); recreate the dynamics of the high-fidelity model (II); and be computationally easier to solve than the high-fidelity model (III).

The first topic, model validation, is a recently growing and still very open field. For a general description, see, for example, [11]. For a few specific works, see [12,13]. In this paper, the original model of $\beta$ species represents the true (physical, chemical, biological, etc.) system under study.

In many types of model reduction of systems of differential equations, the goal is to reduce the computational cost (III) while controlling the incurred error (II). This certainly makes sense as an objective: one would expect a reduced model to offer computational savings over the original model. In order to gain computational efficiency, these types of reductions are not exact, but instead offer an approximation of the species behavior. For example, eigendecompositions may yield an approximation of the static equilibrium state and not the dynamical behavior. Other techniques reduce the computational complexity while maintaining time dynamics, such as volume averaging [14], perturbation theory [15], spectral analysis [16], and the separation of fine- and coarse-scale variables [17]. Related work identifies the exact dynamics on the relevant slow manifold, using, for example, fast-slow decompositions [18] or polynomial approximations up to a specified accuracy [19]. The identification of symmetries can reveal reduced exact dynamics in initial-value problems [20] and the KdV–Zakharov–Kuznetsov equation [21]. For a good overview of more methods, see [22].

In contrast, our goal was to reduce the number of coupled equations that make up the model, or equivalently, the number of species involved in the model. We investigated two possibilities for this type of dimension reduction in the context of the generalized LV equations, called integral and algebraic substitution. There are two defining characteristics of these methods. First, they preserve the correspondence between the set of $\alpha$ species of interest as they appear in the original model and the resulting set after the reduction occurs. This property is termed species correspondence, or simply correspondence. Second, they create a reduced model that contains the exact same information as the original one, but with fewer equations. Variables are eliminated without loss of information. This means that the resulting derivatives (after reduction) are equivalent to those of the original system and therefore that the solutions of the original and the reduced equations are the same dynamical system. In terms of the points above, (I) the model is an exact representation of the true system, and (II) the reduction incurs zero error. Although the resultant model is not necessarily better in a computational sense, this method reveals a path towards model reduction that preserves correspondence, closely approximates the dynamics, and is computationally quite simple.

The two methods discussed here do not automatically find the $\alpha$ species comprising the reduced set, but rather assume them given. Some techniques do choose this set as a step within the reduction process itself. Doing so may favor species with high relative concentrations, or those with slow dynamics, for example. However, consider a combustion model in which we are concerned about trace amounts of a contaminant or an ecological model built to track a species near extinction. The methods presented here allow the modeler to include these critical species a priori.

The rest of the paper is structured as follows. Section 2 outlines the generalized Lotka–Volterra equations. Section 3 presents our main results about two types of exact dimension reduction. Based on these results, Section 4 presents some possible approximate methods, and we conclude with Section 5.

2. Background: The Generalized Lotka–Volterra Equations

Let us briefly review the GLV equations. Let $\mathit{x}$ be the $\beta$-vector of species concentrations. Here, units refer to the number of specimens per unit area, but specific units were omitted for this paper. In the framework of generalized Lotka–Volterra equations, the system of ODEs is given as:

$$\frac{d\mathit{x}}{dt}=\mathrm{diag}\left(\mathit{x}\right)(\mathit{b}+\mathit{A}\mathit{x}),$$

(1)

where the $\beta$-vector $\mathit{b}$ is the intrinsic growth rate vector and $\mathit{A}$ is the $\beta \times \beta$ interaction matrix. Note that there is one differential equation for each species variable. There is coexistence, or a feasible equilibrium, if, at equilibrium, $x_i>0\forall i\in 1,\dots ,\beta$. In this case, then the equilibrium solution is given by:

$$\mathit{x}=-A^{-1}\mathit{b}.$$

(2)

If instead some species become extinct (i.e., $x_j=0$ for some $j\in \{1,\dots ,\beta \}$), then other equilibria are possible.

Here, it was assumed all $x_i\left(t\right)>0,t\ge 0$. Conditions on the entries of A and $\mathit{b}$ that guarantee coexistence were given in [23]; the authors studied the effects of intra- (within a species) and inter-specific (among different species) competition on the coexistence of large systems. In [24], the authors explored the conditions for stability under perturbations, asymptotic feasibility, and the equilibrium of large systems.

3. Exact Dimension Reduction

We begin with some definitions. Let $\alpha ,\beta \in {\mathbb{Z}}^{+}$.

Definition 1

($(\beta ,\alpha )$-reducible). A system of β differential equations that can be converted to a set of α differential equations, where $\alpha <\beta$ and without loss of information, is called ($\beta ,\alpha$)-reducible.

Definition 2

(Integral substitution). Integral substitution (IS) may allow for reductions from a system of β equations to α by eliminating the variables $x_{\alpha +1},\dots ,x_{\beta }$. During this process, we introduce the entire history, or memory, of a subset of the variables $x_1,\dots ,x_{\alpha }$.

This approach is similar to the Mori–Zwanzig method of model reduction [25].

Definition 3

(Algebraic substitution). Algebraic substitution (AS) may allow for reductions from a system of β equations to α by eliminating the variables $x_{\alpha +1},\dots ,x_{\beta }$. During this process, we introduced higher order derivatives of a subset of the variables $x_1,\dots ,x_{\alpha }$.

A classic example from physics in which AS allows for exact reduction is that of a harmonic oscillator, such as a mass on a spring or a simple pendulum [26]. Consider the two coupled first-order differential equations for the position and velocity of a mass on a spring in one dimension:

$$\begin{array}{cc}\hfill \frac{dx}{dt}\left(t\right)& =v\left(t\right)\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill m\frac{dv}{dt}\left(t\right)& =-kv\left(t\right)+f\left(t\right).\hfill \end{array}$$

(3b)

where m is the mass, x is the position, v the velocity, t time, k the spring constant, and $f\left(t\right)$ the forcing function. Notice that the variable $v=dx/dt$; this equation can be converted into a single second-order differential equation:

$$m\frac{d^{2}x}{d{t}^2}\left(t\right)=-k\frac{dx}{dt}\left(t\right)+f\left(t\right).$$

(4)

As seen here, a second-order system of one equation contains exactly the same information as two first-order equations.

This process is also similar to the algebraic reduction presented in [27].

Example 1.

The GLV equations are $(2,1)$-reducible via IS.

With $\beta =2$, the original model is:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =b_1{x}_1+(a_{11}{x}_1+a_{12}{x}_2)x_1\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill {\dot{x}}_2& =b_2{x}_2+(a_{21}{x}_1+a_{22}{x}_2)x_2.\hfill \end{array}$$

(5b)

The goal is to rewrite $x_2$ in terms of $x_1$ in Equation (5a), specifically in the term $a_{12}{x}_2{x}_1$.

First, rearranging Equation (5a) gives:

$$\begin{array}{cc}\hfill x_2& =\frac{1}{a_{12}}\left(\frac{1}{x_1}({\dot{x}}_1-b_1{x}_1)-a_{11}{x}_1\right)\equiv {\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}.\hfill \end{array}$$

(6)

We denote this quantity ${\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}$ because it is a representation of $x_2$ that only depends on $x_1$. With regards to the notation, bold type indicates any such introduced variables to more easily distinguish them from the $x_i$s. Furthermore, a subscript indicates which variable the new one is replacing, and the superscript shows on which variables this new one actually depends. Now, substituting ${\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}$ back into (5a) yields $0=0$. Instead, let us substitute ${\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}$ for $x_2$ in (5b):

$${\dot{x}}_2=b_2{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+(a_{21}{x}_1+a_{22}{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}){\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}.$$

(7)

Integrating, we have,

$$\begin{array}{cc}\hfill x_2& ={\int }_0^t{b}_2{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+(a_{21}{x}_1+a_{22}{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}){\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}\equiv {\mathit{\chi }}_{\mathbf{2}}^{\mathbf{1}}.\hfill \end{array}$$

(8)

Similarly, the symbol ${\mathit{\chi }}_{\mathbf{2}}^{\mathbf{1}}$ means this is a variable that is equivalent to $x_2$, but only in terms of $x_1$. Finally, (5a) becomes:

$${\dot{x}}_1=b_1{x}_1+(a_{11}{x}_1+a_{12}{\mathit{\chi }}_{\mathbf{2}}^{\mathbf{1}})x_1.$$

(9)

We now have a system of a single differential equation, in terms of $x_1$ and its memory. Note that this process has preserved species correspondence and that there is no loss of information: the variable $x_1$ and its derivative in Equation (9) are equivalent to those in (5a).

Example 2.

The GLV equations are also $(2,1)$-reducible via AS.

The first step here is the same as that of the previous subsection, yielding ${\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}$. Next, however, by differentiating ${\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}$, we have:

$${\dot{\mathit{y}}}_{\mathbf{2}}^{\mathbf{1}}=\frac{1}{a_{12}}\left(\frac{-{\dot{x}}_1}{x_1^2}\left({\dot{x}}_1-b_1{x}_1\right)+\frac{1}{x_1}\left({\ddot{x}}_1-b_1{\dot{x}}_1\right)-a_{11}{\dot{x}}_1\right).$$

(10)

Equations (23) and (10) express $x_2$ and ${\dot{x}}_2$, respectively, in terms of $x_1$, ${\dot{x}}_1$, and ${\ddot{x}}_1$. Finally, we can rewrite the second equation as:

$${\dot{\mathit{y}}}_{\mathbf{2}}^{\mathbf{1}}=b_2{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+(a_{21}{x}_1+a_{22}{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}){\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}.$$

(11)

Again, Equation (11) is a system of a single differential equation, but here in terms of $x_1$ and its derivatives, ${\dot{x}}_1$ and ${\ddot{x}}_1$.

Remark 1

(Resulting functional form). Both the integral and algebraic forms, manipulated in this way, yield a single differential equation in terms of $x_1$. Both methods respect species correspondence and are exact. A major difference between the two methods is revealed by inspection of the structure of the resulting sole equations. After reduction via the integral method, the structure of the final equation resembles that of the initial Equation (5a) for $x_1$. The functional form matches in the placement of the variables $x_1$ and ${\mathit{\chi }}_{\mathbf{2}}^{\mathbf{1}}$ (in place of $x_2$) and also the remaining constants $b_1$, $a_{11}$, and $a_{12}$ (which do not appear in the second equation of the original system). In contrast, after reduction via the algebraic method, the resultant equation has the structure of (5b). This suggests that in an applied setting, one type of reduction may be advantageous, depending on what information is known about the high-fidelity model, such as various model parameters.

For the remainder of the paper, we focused on IS, as opposed to AS. The techniques and results between the two methods are quite similar. An advantage of IS is its notational convenience (as described in Remark 1): the reduced equations preserve the structure of the original ones, which simplifies tracking the reductions, in theory.

Remark 2

($(\beta ,\beta -1)$-reducibility). In a similar way, we can reduce any generalized LV system of β species to one of $\beta -1$ species. Note that this level of reduction, from β to $\beta -1$ equations, can happen when the ODEs describe the dynamics of fractional concentrations. In that case, the extra constraint that ${\sum }_{i=1}^{\beta }{x}_i=1$ readily allows for a reduction to $\beta -1$ equations, since, for example, $x_{\beta }$ can be expressed as $x_{\beta }=1-{\sum }_{i=1}^{\beta -1}{x}_i$ (in other words, as shown by [28], the replicator equation in S variables is equivalent to the generalized LV equations in $S-1$ variables). In the current work, however, there was no such restriction on the concentrations, and yet, such a reduction is always possible.

Reductions from two to one equation and from $\beta$ to $\beta -1$ are similar.

Lemma 1.

The GLV equations are ($\beta ,\beta -1$)-reducible via IS.

Proof. 

The original model for $\beta$ species is:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =b_1{x}_1+(a_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1S}{x}_{\beta })x_1\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill {\dot{x}}_2& =b_2{x}_2+(a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2\beta }{x}_{\beta })x_2\hfill \\ \hfill & ⋮\hfill \end{array}$$

(12b)

$$\begin{array}{cc}\hfill {\dot{x}}_{\beta }& =b_{\beta }{x}_{\beta }+(a_{\beta 1}{x}_1+a_{\beta 2}{x}_2+\cdots +a_{\beta \beta }{x}_{\beta })x_{\beta }.\hfill \end{array}$$

(12c)

At this point, we want to rewrite $x_{\beta }$ in terms of the remaining $\beta -1$ variables, but we could use any of the first $\beta -1$ equations to do so. Without loss of generality, we choose the second to last one: -4.6cm0cm

$$\begin{array}{cc}\hfill x_{\beta }& =\frac{1}{a_{\beta -1,\beta }}\left(\frac{1}{x_{\beta -1}}({\dot{x}}_{\beta -1}-b_{\beta -1}{x}_{\beta -1})-a_{\beta -1,1}{x}_1-a_{\beta -1,2}{x}_2-\cdots -a_{\beta -1,\beta -1}{x}_{\beta -1}\right)\equiv {\mathit{y}}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}.\hfill \end{array}$$

(13)

The colon notation in ${\mathit{y}}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}$ signifies that this new variable is written in terms of variables $x_1$ through $x_{\beta -1}$.

Now, substituting ${\mathit{y}}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}$ into (12c) yields:

$$\begin{array}{cc}\hfill x_{\beta }& ={\int }_0^t\left(b_{\beta }{\mathit{y}}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}+\left(a_{\beta 1}{x}_1+a_{\beta 2}{x}_2+\cdots +a_{\beta \beta }{\mathit{y}}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right){\mathit{y}}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right)\equiv {\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}.\hfill \end{array}$$

Finally, substituting ${\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}$ back into the first $\beta -1$ equations:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =b_1{x}_1+\left(a_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right)x_1\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill {\dot{x}}_2& =b_2{x}_2+\left(a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right)x_2\hfill \\ \hfill & ⋮\hfill \end{array}$$

(14b)

$$\begin{array}{cc}\hfill {\dot{x}}_{\beta -1}& =b_{\beta -1}{x}_{\beta -1}+\left(a_{\beta -1,1}{x}_1+a_{\beta -1,2}{x}_2+\cdots +a_{\beta -1,\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right)x_{\beta -1}.\hfill \end{array}$$

(14c)

Thus, the set of $\beta$ ODEs is reduced to $\beta -1$ without loss of information. □

In general, a system of GLV equations is not ($\beta ,\alpha$)-reducible because of the coupled-ness, or “entanglement”, that occurs between the species in the reduced set via the introduced $\mathit{y}$ and/or $\mathit{\chi }$ variables. To see this, consider the case of $\beta =3$, $\alpha =1$. With both the integral and algebraic methods, the model cannot be exactly reduced, unless one of the coefficients is zero, thereby breaking the entanglement.

Lemma 2.

The GLV equations are $(3,1)$-reducible if only one $a_{ij}=0,i\ne j$.

Proof. 

We begin the process to reduce the system of three equations to one by eliminating the variables $x_2$ and $x_3$. When $\beta =3$, the original model is:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =b_1{x}_1+(a_{11}{x}_1+a_{12}{x}_2+a_{13}{x}_3)x_1\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill {\dot{x}}_2& =b_2{x}_2+(a_{21}{x}_1+a_{22}{x}_2+a_{23}{x}_3)x_2\hfill \end{array}$$

(15b)

$$\begin{array}{cc}\hfill {\dot{x}}_3& =b_3{x}_3+(a_{31}{x}_1+a_{32}{x}_2+a_{33}{x}_3)x_3.\hfill \end{array}$$

(15c)

Repeating the process described in the previous section, we can rewrite either Equation (15a) or (15b) for $x_3$. Using Equation (15b), then:

$$\begin{array}{cc}\hfill x_3& =\frac{1}{a_{23}}\left(\frac{1}{x_2}({\dot{x}}_2-b_2{x}_2)-a_{21}{x}_1-a_{22}{x}_2\right)\equiv {\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}.\hfill \end{array}$$

(16)

Next, substituting ${\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}$ into (15c),

$$\begin{array}{cc}\hfill {\dot{x}}_3& =b_3{\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}+(a_{31}{x}_1+a_{32}{x}_2+a_{33}{\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}){\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}.\hfill \end{array}$$

(17)

Integrating,

$$\begin{array}{cc}\hfill x_3& ={\int }_0^t{b}_3{\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}+(a_{31}{x}_1+a_{32}{x}_2+a_{33}{\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}){\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}\equiv {\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}.\hfill \end{array}$$

(18)

Now, we can substitute ${\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}$ into Equations (15a) and (15b):

$$\begin{array}{cc}\hfill {\dot{x}}_1& =b_1{x}_1+(a_{11}{x}_1+a_{12}{x}_2+a_{13}{\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}})x_1\hfill \end{array}$$

(19a)

$$\begin{array}{cc}\hfill {\dot{x}}_2& =b_2{x}_2+(a_{21}{x}_1+a_{22}{x}_2+a_{23}{\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}})x_2.\hfill \end{array}$$

(19b)

At this point, the model has been reduced from three equations to two. Consider if we now tried to remove another variable, say $x_2$. The first step would be to rewrite Equation (19a) so that $x_2$ is alone on the left-hand side (LHS). However, now that ${\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}$ has been introduced, we cannot cleanly separate $x_2$ out of the right-hand side (RHS) since it is embedded inside the integral term.

A similar problem occurs with the algebraic approach.

Now assume of the six off-diagonal interaction terms, it is set to zero. Without loss of generality, set $a_{13}=0$. Then, Equation (19a) becomes:

$${\dot{x}}_1=b_1{x}_1+(a_{11}{x}_1+a_{12}{x}_2)x_1$$

(20)

and so,

$$\begin{array}{cc}\hfill x_2& =\frac{1}{a_{12}}\left(\frac{1}{x_1}({\dot{x}}_1-b_1{x}_1)-a_{11}{x}_1\right)\equiv {\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}.\hfill \end{array}$$

(21)

Substituting ${\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}$ into Equation (19b),

$$\begin{array}{cc}\hfill {\dot{x}}_2& =b_2{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+(a_{21}{x}_1+a_{22}{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+a_{23}{\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}){\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & =b_2{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+(a_{21}{x}_1+a_{22}{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+a_{23}{\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}}){\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}.\hfill \end{array}$$

(23)

The second line above appears after replacing ${\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}$ with ${\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}}$. This variable ${\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}}$ is found by replacing any explicit dependence on $x_2$ with ${\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}$. Specifically,

$$\begin{array}{cc}\hfill {\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}}& ={\int }_0^t{b}_3{\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}+(a_{31}{x}_1+a_{32}{x}_2+a_{33}{\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}){\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & ={\int }_0^t{b}_3{\mathit{y}}_{\mathbf{3}}^{\mathbf{1}}+(a_{31}{x}_1+a_{32}{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+a_{33}{\mathit{y}}_{\mathbf{3}}^{\mathbf{1}}){\mathit{y}}_{\mathbf{3}}^{\mathbf{1}},\hfill \end{array}$$

(25)

where, to in turn find ${\mathit{y}}_{\mathbf{3}}^{\mathbf{1}}$, we modify Equation (16):

$$\begin{array}{cc}\hfill {\mathit{y}}_{\mathbf{3}}^{\mathbf{1}:\mathbf{2}}& =\frac{1}{a_{23}}\left(\frac{1}{x_2}({\dot{x}}_2-b_2{x}_2)-a_{21}{x}_1-a_{22}{x}_2\right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & =\frac{1}{a_{23}}\left(\frac{1}{{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}}(\dot{{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}}-b_2{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}})-a_{21}{x}_1-a_{22}{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}\right)\equiv {\mathit{y}}_{\mathbf{3}}^{\mathbf{1}}.\hfill \end{array}$$

(27)

Integrating the line (23),

$$\begin{array}{cc}\hfill x_2& ={\int }_0^t{b}_2{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+(a_{21}{x}_1+a_{22}{\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}+a_{23}{\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}}){\mathit{y}}_{\mathbf{2}}^{\mathbf{1}}\equiv {\mathit{\chi }}_{\mathbf{2}}^{\mathbf{1}},\hfill \end{array}$$

(28)

and finally,

$${\dot{x}}_1=b_1{x}_1+(a_{11}{x}_1+a_{12}{\mathit{\chi }}_{\mathbf{2}}^{\mathbf{1}})x_1.$$

(29)

□

Therefore, we can exactly reduce the system of three ODEs to one differential equation if we assume $a_{13}=0$. However, note that a nested set of integrals comprises the term ${\mathit{\chi }}_{\mathbf{3}}^{\mathbf{1}}$. Indeed, such a system of one differential equation with nested integral terms as Equation (29) would likely be more difficult to solve than the original set of three coupled ODEs shown in (15a)–(15c).

Remark 3

(Entanglement). Multiple variables can be eliminated if enough information (memory terms, higher derivatives) is known about the remaining variables. Lemma 3 shows the necessity of the assumption that some interaction term be zero in order to exactly reduce (in this demonstration, that $a_{13}=0$). This assumption begins to reveal the limitations of such reductions via substitutions—in particular, how, after a single substitution, we cannot escape the entanglement of the remaining species.

At the same time, the methods here give insight into how one might approximate the role of eliminated variables in terms of the reduced set. For example, in the case of $\beta =2,\alpha =1$, consider an approximation ${\widehat{\mathit{\chi }}}_{\mathbf{2}}^{\mathbf{1}}$ in terms of $x_1$ and this extra information:

$${\mathit{\chi }}_{\mathbf{2}}^{\mathbf{1}}\approx {\widehat{\mathit{\chi }}}_{\mathbf{2}}^{\mathbf{1}}\left(x_1,{\dot{x}}_1,{\ddot{x}}_1,K\left(x_1\right),\dots \right),$$

(30)

where $K\left(x_1\right)$ represents some memory kernel of $x_1$. We explore this type of approximation in Section 4.

However, first, we complete the reduction via IS for general $\beta$ and $\alpha$. Of course, in this setting, a larger collection of interaction terms must be assumed zero. Define the set of coefficients ${\mathcal{A}}_{\beta ,\alpha }$, so that:

$${\mathcal{A}}_{\beta ,\alpha }=\left\{a_{ij}\right\},i=\alpha ,\dots ,\beta -2;j>i+1.$$

(31)

Theorem 1.

The GLV equations are ($\beta ,\alpha$)-reducible via IS if $a=0\forall a\in {\mathcal{A}}_{\beta ,\alpha }$, and $a\ne 0\forall a\in {\mathcal{A}}_{\beta ,\alpha }^c$.

Proof. 

In Lemma 1, we completed the reduction from $\beta$ to $\beta -1$ equations. The next step reduces to $\beta -2$. Let us start with the $\beta -1$ equations:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =b_1{x}_1+(a_{11}{x}_1+\cdots +a_{1\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}})x_1\hfill \end{array}$$

(32a)

$$\begin{array}{cc}\hfill {\dot{x}}_2& =b_2{x}_2+(a_{21}{x}_1+\cdots +a_{2\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}})x_2\hfill \\ \hfill & ⋮\hfill \end{array}$$

(32b)

$$\begin{array}{cc}\hfill {\dot{x}}_{\beta -2}& =b_{\beta -2}{x}_{\beta -2}+\left(a_{\beta -2,1}{x}_1+\cdots +a_{\beta -2,\beta -2}{x}_{\beta -2}+a_{\beta -2,\beta -1}{x}_{\beta -1}+a_{\beta -2,\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right)x_{\beta -2}\hfill \end{array}$$

(32c)

$$\begin{array}{cc}\hfill {\dot{x}}_{\beta -1}& =b_{\beta -1}{x}_{\beta -1}+\left(a_{\beta -1,1}{x}_1+\cdots +a_{\beta -1,\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right)x_{\beta -1}.\hfill \end{array}$$

(32d)

Now, we rewrite $x_{\beta -1}$ in terms of the remaining $\beta -2$ variables; again, we could use any of the first $\beta -2$ equations to do so. Let us choose the ODE for $x_{\beta -2}$ with the assumption that $a_{\beta -2,\beta }=0$:

$$\begin{array}{cc}\hfill x_{\beta -1}& =\frac{1}{a_{\beta -2,\beta -1}}\left(\frac{1}{x_{\beta -2}}({\dot{x}}_{\beta -2}-b_{\beta -2}{x}_{\beta -2})-a_{\beta -2,1}{x}_1-\cdots -a_{\beta -2,\beta -3}{x}_{\beta -3}\right)\equiv {\mathit{y}}_{\mathit{\beta }-\mathbf{1}}^{\mathbf{1}:\mathit{\beta }-\mathbf{2}}.\hfill \end{array}$$

(33)

Substituting ${\mathit{y}}_{\mathit{\beta }-\mathbf{1}}^{\mathbf{1}:\mathit{\beta }-\mathbf{2}}$ into (32d),

$$\begin{array}{cc}\hfill x_{\beta -1}& ={\int }_0^t\left(b_{\beta -1}{\mathit{y}}_{\mathit{\beta }-\mathbf{1}}^{\mathbf{1}:\mathit{\beta }-\mathbf{2}}+(a_{\beta -1,1}{x}_1+\cdots +a_{\beta -1,\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}){\mathit{y}}_{\mathit{\beta }-\mathbf{1}}^{\mathbf{1}:\mathit{\beta }-\mathbf{2}}\right){\mathit{\chi }}_{\mathit{\beta }-\mathbf{1}}^{\mathbf{1}:\mathit{\beta }-\mathbf{2}}.\hfill \end{array}$$

(34)

Finally, substituting ${\mathit{\chi }}_{\mathit{\beta }-\mathbf{1}}^{\mathbf{1}:\mathit{\beta }-\mathbf{2}}$ back into the first $\beta -2$ equations,

$$\begin{array}{cc}\hfill {\dot{x}}_1& =b_1{x}_1+\left(a_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right)x_1\hfill \end{array}$$

(35a)

$$\begin{array}{cc}\hfill {\dot{x}}_2& =b_2{x}_2+\left(a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2\beta }{\mathit{\chi }}_{\mathit{\beta }}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right)x_2\hfill \\ \hfill & ⋮\hfill \end{array}$$

(35b)

$$\begin{array}{cc}\hfill {\dot{x}}_{\beta -2}& =b_{\beta -2}{x}_{\beta -2}+\left(a_{\beta -2,1}{x}_1+a_{\beta -2,2}{x}_2+\cdots +a_{\beta -2,\beta }{\mathit{\chi }}_{\mathit{\beta }-\mathbf{1}}^{\mathbf{1}:\mathit{\beta }-\mathbf{1}}\right)x_{\beta -2}.\hfill \end{array}$$

(35c)

This process can be repeated once more, yielding a system of $\beta -3$ equations, if we also assume $a_{\beta -3,\beta }=a_{\beta -3,\beta -1}=0$. We may find the equivalent system in $\beta -4$ equations if the set of zero interaction terms also includes $a_{\beta -4,\beta }$, $a_{\beta -4,\beta -1}$, and $a_{\beta -4,\beta -2}$. Continuing in this way, the model can be reduced from $\beta$ to s equations if $a=0\forall a\in {\mathcal{A}}_{\beta ,s}$. Note the assumption that $a\ne 0\forall a\in {\mathcal{A}}_{\beta ,\alpha }^c$ allows for division by these coefficients. □

Let us examine this set ${\mathcal{A}}_{\beta ,\alpha }$ in more detail. We compute the fraction $\rho$ of these “zeroed” terms out of all interaction coefficients.

Corollary 1.

Let $\lambda =\alpha /\beta$. Then, ${lim}_{\beta \to \infty }\rho =\frac{1}{2}{(1-\lambda )}^2.$

Proof. 

Given $\beta$ species in the original model, there are ${\beta }^2$ interaction terms, and $|{\mathcal{A}}_{\beta ,\alpha }|=\frac{(k-1)k}{2}$. In terms of $\alpha$, $\beta$, then $\rho$ is:

$$\begin{array}{cc}\hfill \rho & =\frac{(k-1)k}{2{\beta }^2}\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & =\frac{(\beta -s-1)(\beta -s)}{2{\beta }^2}\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & =\frac{{\beta }^2+s^2-2\beta s-\beta -s}{2{\beta }^2}.\hfill \end{array}$$

(38)

Then, the fraction $\rho$ may be written as:

$$\begin{array}{cc}\hfill \rho & =\frac{1}{2}(1+{\lambda }^2-2\lambda )-\frac{1}{2\beta }(1+\lambda )\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill & =\frac{1}{2}{(1-\lambda )}^2-\frac{1}{2\beta }(1+\lambda ).\hfill \end{array}$$

(40)

In the limit as $\beta \to \infty$ and for a fixed $\lambda$, we have:

$$\begin{array}{c}\hfill \underset{\beta \to \infty }{lim}\rho =\frac{1}{2}{(1-\lambda )}^2.\end{array}$$

(41)

□

This limiting value is plotted for $\lambda \in [0,1]$ in Figure 1. Note that in this limit, when $\lambda =0$, i.e., there is a complete reduction, then half of the interaction terms vanish. At the other extreme, when $\lambda =1$, i.e., $\alpha =\beta$ so there is no reduction, then of course, all interaction terms can be nonzero.

Note that the conditions on the set ${\mathcal{A}}_{\beta ,\alpha }$ and its complement are not necessary, but sufficient. They are not necessary at least in the sense of non-uniqueness: reordering the variables or making different substitutions could lead to a different set of zeroed interactions terms. However, the current results suggest at least some number $z^{\ast }$ of coefficients must be zero; otherwise, the entanglement of the species prevents exact reduction. We thus end this section with the following conjecture.

Conjecture 1

For $(\beta ,\alpha )$-reducibility, it is necessary that at least $z^{\ast }=\left|{\mathcal{A}}_{\beta ,\alpha }\right|$ terms are zero, where, for $\alpha =\beta -k$, $|{\mathcal{A}}_{\beta ,s}|=\frac{1}{2}k(k-1)$.

4. Approximate Model Reduction

Motivated by the previous results, the final topic of this paper explored approximate model reduction. In this section, we consider the situation that information—interaction coefficients, growth rates, and some observations—is only available about a subset of $\alpha$ out of $\beta$ species. We aimed to create an approximate model for those $\alpha$ species without including any information about the excluded $\beta -\alpha$ variables.

Interestingly, another recently developed method tackles the same situation: the focal species approximation (FSA) [9,29]. For an environment of S species, the FSA estimates interaction coefficients (from experiments) involving the species of interest and ignores the rest, thereby dropping the number of needed experiments from $\mathcal{O}\left(S^2\right)$ to $\mathcal{O}\left(S\right)$. While the accuracy of this approximation is lower than that of the complete GLV model, the approximation still provides useful predictions and can thus be employed in applied settings for which full parameter estimation is infeasible. In a sense, the approximate AS and IS and the FSA demonstrate that (severely) reduced models may still capture the dynamics of interest within reasonable error.

4.1. Algebraic Substitutions

In the algebraic method, eliminated variables are exchanged for higher derivatives about the remaining variables.

Example 3

(Approximate AS for $\beta =2$, $\alpha =1$). The goal is to find an approximate model for $x_1$, in terms of only $x_1$.

The original model is, as before:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =r_1{x}_1+(a_{11}{x}_1+a_{12}{x}_2)x_1\hfill \end{array}$$

(42a)

$$\begin{array}{cc}\hfill {\dot{x}}_2& =r_2{x}_2+(a_{21}{x}_1+a_{22}{x}_2)x_2.\hfill \end{array}$$

(42b)

Motivated by the AS above, we write $x_2$ as an expansion in $x_1$ and its higher derivatives:

$$x_2\to {\kappa }_0{x}_1+{\kappa }_1{\dot{x}}_1.$$

(43)

This leads to an approximate model for $x_1$:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =r_1{x}_1+(a_{11}{x}_1+a_{12}{x}_2)x_1\hfill \\ \hfill & \approx r_1{x}_1+\left(a_{11}{x}_1+a_{12}({\kappa }_0{x}_1+{\kappa }_1{\dot{x}}_1)\right)x_1\hfill \\ \hfill & =r_1{x}_1+a_{11}{x}_1^2+{\delta }_0{x}_1^2+{\delta }_1{\dot{x}}_1{x}_1,\hfill \end{array}$$

(44)

where ${\delta }_0=a_{12}{\kappa }_0$ and ${\delta }_1=a_{12}{\kappa }_1$.

It may also be of interest to see the behavior of $x_1$ given the remaining terms from the original model, i.e., $\dot{x}=r_1{x}_1+a_{11}{x}_1^2$. We call this the naive reduced model. Note that this is the standard LV model if one only has access to the interaction coefficients of the reduced set of species.

Let us specify the coefficients of the original model as:

$$\mathit{r}=\left[\begin{array}{c}5\\ 3\end{array}\right],\mathit{A}=\left[\begin{array}{cc}-3& -1\\ -1& -2\end{array}\right],$$

(45)

and set the initial conditions as $x_1\left(0\right)=0.289$ and $x_2\left(0\right)=0.665$. (These values are approximate. The precise initial conditions were chosen randomly from a standard lognormal distribution.)

The introduced coefficients ${\delta }_0$ and ${\delta }_1$ are unknown and must be calibrated. We assumed data about $x_1$ over time. Calibration was performed under a Bayesian framework; the posterior means were ${\overline{\delta }}_1\approx -0.280$ and ${\overline{\delta }}_2\approx -0.357$. The code to run Bayesian inverse and forward problems is available here: github.libqueso.com (accessed on 21 April 2021) [30]; specific implementations are provided in the supplementary material. The trajectory of $x_1$ is shown in Figure 2 from the three models: original, naive reduced, and approximate.

The approximate model is quite close to the original. Importantly, the approximate model output completely covers the original model output, which lies entirely within the approximate model 50% confidence interval.

Example 4

(Approximate AS for $\beta =3$, $\alpha =1$.). Here, we tried the same framework with a slightly larger reduction: from three to one. In this case, note that we still replaced $x_2$ and $x_3$ with our approximate expression, but still only two new coefficients need be introduced.

Recall the original model for $x_1$ is:

$${\dot{x}}_1=r_1{x}_1+(a_{11}{x}_1+a_{12}{x}_2+a_{13}{x}_3)x_1.$$

The approximate model is then:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =r_1{x}_1+\left(a_{11}{x}_1+a_{12}({\kappa }_1{x}_1+{\kappa }_2{\dot{x}}_1)+a_{13}({\gamma }_1{x}_1+{\gamma }_2{\dot{x}}_1)\right)x_1\hfill \\ \hfill & =r_1{x}_1+a_{11}{x}_1^2+(a_{12}{\kappa }_1+a_{13}{\gamma }_1)x_1^2+(a_{12}{\kappa }_2+a_{13}{\gamma }_2)x_1{\dot{x}}_1\hfill \\ \hfill & =r_1{x}_1+a_{11}{x}_1^2+{\delta }_1{x}_1^2+{\delta }_2{x}_1{\dot{x}}_1\hfill \end{array}$$

(46)

where ${\delta }_1=a_{12}{\kappa }_1+a_{13}{\gamma }_1$ and ${\delta }_2=a_{12}{\kappa }_2+a_{13}{\gamma }_2$. Thus, we can again calibrate simply ${\delta }_1$ and ${\delta }_2$.

Note the naive reduced model is again just:

$${\dot{x}}_1=r_1{x}_1+a_{11}{x}_1^2.$$

Specifically, we set:

$$\mathit{r}=\left[\begin{array}{c}2.2\\ 1.7\\ 3.1\end{array}\right],\mathit{A}=\left[\begin{array}{ccc}-1.2& -0.5& -0.1\\ 0.2& -1.5& -0.8\\ -0.3& -0.7& -1.1\end{array}\right],$$

(47)

and the initial conditions to $x_1\left(0\right)=0.289$, $x_2\left(0\right)=0.665$, and $x_3\left(0\right)=5.01$.

The posterior means were ${\overline{\delta }}_1\approx -0.311$ and ${\overline{\delta }}_2\approx -0.432$. The trajectories from the three models are plotted in Figure 3.

Again, the approximate model is a good match to the original for $x_1$. Note that ten coefficients from the original model were omitted (eight interaction and two growth rate coefficients), while only two new parameters were introduced, ${\delta }_1$ and ${\delta }_2$.

4.2. Integral Substitutions

Another approach is instead inspired by the IS above, so that information about eliminated variables is replaced by memory, or integral, information of the remaining species.

Example 5

(Approximate IS for $\beta =2$, $\alpha =1$). Starting, as usual, with the original formulation:

$${\dot{x}}_1=r_1{x}_1+a_{11}{x}_1^2+a_{12}{x}_2{x}_1$$

we aimed for an approximate model by eliminating $x_2$.

In this example, we also tried a slight variation from the algebraic one in that we replaced the entire term containing $x_2$:

$$a_{12}{x}_1{x}_2\to {\delta }_1{x}_1+{\delta }_2{\int }_0^t{x}_1\left(\tau \right)d\tau .$$

(48)

Moreover, since, in all of these examples, we assumed some data about the species of interest, we can approximate this integral directly from the data, which is generated according to the original model (and not the $x_1$ as given by the approximate model). In these numerical results, this integral was estimated numerically using a simple trapezoid rule. Call this numerical approximation $\mathcal{I}\left(t\right)$:

$$\mathcal{I}\left(t\right)\approx {\int }_0^t{x}_1^{\ast }\left(\tau \right)d\tau ,$$

(49)

where the ∗ in $x_1^{\ast }$ indicates that this integral was estimated from the data generated by the original model. Then, we pose the approximate model for $x_1$ as:

$${\dot{x}}_1=r_1{x}_1+a_{11}{x}_1^2+{\delta }_1{x}_1+{\delta }_2\mathcal{I}\left(t\right).$$

(50)

The posterior means in this case were ${\overline{\delta }}_1\approx -0.778$ and ${\overline{\delta }}_2\approx -0.449$. The different trajectories for this example are shown in Figure 4.

Example 6

(Approximate IS for $\beta =3$, $\alpha =1$). Here, we present the final example, with the same original model as Example 4 (approximate AS for $\beta =3$, $\alpha =1$), and the same approximate model as the previous example.

Here we find ${\overline{\delta }}_1\approx -0.452$ and ${\overline{\delta }}_2\approx -0.168$. The results are shown in Figure 5.

The approximate integral models also showed quite good agreement with the original model.

Note that each example here preserved the species correspondence: the species of interest was $x_1$, and that is precisely what was modeled. This would readily generalize to approximate models for $\alpha >1$. While the resulting differential equations now include either a derivative or an approximate integral, the numerical solution is rather easily obtained without introducing much numerical complexity. In [31], Morrison provided an expanded investigation into similar AS-type approximate models.

5. Discussion

In summary, we can reduce a GLV system exactly by introducing integral terms of the remaining species or by introducing their higher derivatives. In a sense, some variables may be eliminated in exchange for this extra information about the reduced set. A few interesting results emerged from this work: First, a reduction from $\beta$ to $\beta -1$ variables is always possible. Second, a reduction from $\beta$ to $\alpha$ is possible if the interaction terms in a specified set are zero. With $\lambda =\alpha /\beta$, the fraction of terms in this set is approximately ${(1-\lambda )}^2/2$. Third, we saw that the two different methods in the $\beta =2$ case yielded a single ODE for $x_1$, but either by maintaining the structure of the original ODE for $x_1$ or for $x_2$. Fourth, after reducing to $\beta -1$ variables, an entanglement of species prevents further (exact) reduction without breaking the original model in some other way.

In this paper, this break was achieved by setting the interaction terms to zero. Of course, in some systems, this may be a reasonable assumption. This work also serves as a motivation for or evidence of how to approximate the role of the eliminated species using only the reduced set. The preliminary results suggested reduced models that they are no longer exact, but still preserve the species correspondence and are computationally feasible. A true validation test of this method would determine the accuracy and precision of the predictions with respect to real ecological (agricultural, chemical, etc.) data. Specific implementations and their effectiveness are an open and rich area for future study.