Higher-Order Lotka-Volterra Dynamics

• Higher-Order Lotka-Volterra dynamics are generalized ecological models where species growth depends on both pairwise and non-additive higher-order interactions.
• The models utilize polynomial formulations, tensor methods, and complementarity problems to analyze equilibria, stability, and dynamic regimes.
• Applications include predicting coexistence, oscillatory behaviors, and species-abundance distributions in complex ecosystems.

Higher-order Lotka-Volterra dynamics denotes a class of generalized Lotka-Volterra systems in which the per-capita growth of each species depends not only on additive pairwise interactions but also on non-additive higher-order interactions (HOIs), typically represented by triadic or higher multilinear terms. In this setting, the joint presence of multiple species can produce effects that are not expressible as a sum of pairwise contributions, including indirect interaction modification by third parties. Current work treats HO-LV models as polynomial dynamical systems, tensor equations, complementarity problems, random disordered ecosystems, and stochastic reaction-diffusion processes, and uses them to analyze coexistence, stability, multistability, oscillations, quasi-periodicity, chaos, extinction, and species-abundance structure (Cui et al., 2024, Sidhom et al., 2024, Kang et al., 29 Jul 2025).

1. Canonical formulations of HO-LV systems

A standard ecological HO-LV model assigns species densities $x_i(t)\ge 0$, intrinsic per-capita growth rates $r_i>0$, pairwise coefficients $a_{ij}\in\mathbb{R}$, and higher-order coefficients $b_{ijk}\in\mathbb{R}$. The quadratic-triadic form analyzed in recent work is

$$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$

with $a_{ii}\le 0$ and $b_{iii}\le 0$ allowed to model self-limitation. In tensor-vector form this becomes

$$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$

where $(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$. The same framework extends naturally to arbitrary polynomial order $m$ by replacing $r_i>0$0 with higher multilinear tensor contractions (Cui et al., 2024).

A related generalized Lotka-Volterra formulation adds a small species-specific dispersal term: $r_i>0$1 Here $r_i>0$2 is the pairwise interaction matrix, $r_i>0$3 is the third-order interaction tensor, and $r_i>0$4 is a small dispersal or immigration rate that prevents extinction in the deterministic limit and facilitates interior coexistence equilibria. In this formulation, $r_i>0$5 and $r_i>0$6 encode competition-like effects, $r_i>0$7 and $r_i>0$8 encode cooperation or facilitation, and mixed signs capture realistic heterogeneity. No explicit saturating Holling forms are used; the functional response is polynomial through the HOI terms (Kang et al., 29 Jul 2025).

The random-community literature generalizes this further to arbitrary interaction order $r_i>0$9. One representative form is

$a_{ij}\in\mathbb{R}$0

with quenched Gaussian higher-order couplings specified by order-dependent mean $a_{ij}\in\mathbb{R}$1, variance $a_{ij}\in\mathbb{R}$2, and symmetry parameter $a_{ij}\in\mathbb{R}$3. This setting isolates the effect of interaction order on stability and coexistence in the thermodynamic limit (Sidhom et al., 2024).

Cyclic three-agent models provide a more specialized HO-LV class. In a well-mixed rock-paper-scissors-type system with three-agent reactions, the mean-field equations are cubic in the densities because three-body interactions enter via mass action; monomials such as $a_{ij}\in\mathbb{R}$4, $a_{ij}\in\mathbb{R}$5, and $a_{ij}\in\mathbb{R}$6 replace the quadratic terms of pairwise cyclic LV models. This makes higher-order interaction order explicit at the level of the reaction scheme, not merely as a tensor generalization of pairwise GLV (Palombi et al., 2019).

2. Equilibria, tensor methods, and complementarity structure

Because the HO-LV systems above are positive, every equilibrium satisfies either $a_{ij}\in\mathbb{R}$7 or vanishing of the per-capita growth factor. For a strictly positive equilibrium of the quadratic-triadic model,

$a_{ij}\in\mathbb{R}$8

This is a non-homogeneous polynomial system, and one of the principal analytical developments is to study it through tensor theory and complementarity methods rather than through direct phase-plane arguments (Cui et al., 2024).

The central tensor notion is the S-tensor, defined by the existence of some $a_{ij}\in\mathbb{R}$9 such that $b_{ijk}\in\mathbb{R}$0 componentwise. Every $b_{ijk}\in\mathbb{R}$1-tensor is an S-tensor. For the non-homogeneous tensor equation

$b_{ijk}\in\mathbb{R}$2

Theorem 3.1 states that if $b_{ijk}\in\mathbb{R}$3 are S-tensors that share a common positive vector $b_{ijk}\in\mathbb{R}$4, then the system has a unique positive solution. The proof proceeds by rewriting each tensor into positive and nonnegative parts, constructing an increasing fixed-point map $b_{ijk}\in\mathbb{R}$5, building an order interval $b_{ijk}\in\mathbb{R}$6, applying an ordered Banach fixed-point lemma, and then establishing uniqueness by a scaling-based contradiction. This result is specifically presented as an existence and uniqueness theorem for the non-homogeneous polynomial systems arising from HO-LV equilibria (Cui et al., 2024).

The same equilibrium problem admits a polynomial complementarity reformulation. With

$b_{ijk}\in\mathbb{R}$7

positive equilibria solve $b_{ijk}\in\mathbb{R}$8, while general equilibria correspond to the polynomial complementarity problem

$b_{ijk}\in\mathbb{R}$9

The HO-LV interpretation is two-sided: any positive solution of the non-homogeneous polynomial system yields a PCP solution, and PCP solutions encode both interior and boundary equilibria via complementarity conditions. For generalized row strictly diagonally dominant tensors and matrices, explicit lower and upper bounds on $$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$0 are derived for the quadratic PCP, providing a priori localization of equilibrium magnitudes (Cui et al., 2024).

These developments matter because the equilibrium problem is intrinsically nonlinear and non-homogeneous. The tensor and PCP formulations convert coexistence analysis into algebraic conditions on S-tensors, $$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$1-tensors, diagonal dominance, and solvability properties such as compactness or the GUS-property, rather than requiring direct construction of Lyapunov functions for each coefficient configuration. This suggests a systematic algebraic route to existence and uniqueness in models where higher-order coefficients proliferate.

3. Monotonicity, Jacobians, and extensions of classical LV structure

For the quadratic HO-LV system, the Jacobian can be written explicitly. In particular,

$$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$2

and

$$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$3

where $$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$4 is the bracketed per-capita factor in the dynamics. This formula underlies the monotone-systems analysis of cooperative and competitive HO-LV models (Cui et al., 2024).

In the fully cooperative model, with nonnegative inter-specific coefficients and nonpositive intra-specific coefficients, the Jacobian is an irreducible Metzler matrix, implying irreducible monotonicity in $$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$5. The origin is always unstable because the Jacobian at $$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$6 is $$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$7 with positive eigenvalues. Boundary equilibria are unstable, and local stability of an interior equilibrium follows from Hurwitzness of the Jacobian; a sufficient criterion is the existence of $$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$8 with $$\dot{x}_i(t)=r_i\,x_i(t)\left(1+\sum_{j=1}^n a_{ij}\,x_j(t)+\sum_{j,k=1}^n b_{ijk}\,x_j(t)\,x_k(t)\right),\quad i=1,\ldots,n,$$9 for an irreducible Metzler $a_{ii}\le 0$0. Under diagonal dominance-like bounds on the Jacobian, the unique positive equilibrium is globally asymptotically stable. When higher-order terms vanish, the model recovers the classical cooperative LV results associated with Goh and Smith (Cui et al., 2024).

The same paper extends analogous classical results to two-faction competition and pure competition. In the two-faction case, within-faction interactions are cooperative, cross-faction interactions are competitive, and the Jacobian can be permuted by $a_{ii}\le 0$1 into an irreducible Metzler matrix. The resulting taxonomy includes unstable origins, locally stable one-faction-wins-all equilibria under strong loser-faction competitive impacts, locally stable mixed boundary equilibria when winners’ subsystems admit stable interior equilibria, unstable strict-subset-winner states within a single faction, and locally stable all-species coexistence when competitive and HOI terms except self-competition are sufficiently small. The paper identifies bistability as a biological interpretation of alternative community assemblies (Cui et al., 2024).

In pure competition,

$a_{ii}\le 0$2

every positive orbit is upper bounded, the interior equilibrium exists and is unique under S-tensor assumptions on $a_{ii}\le 0$3 and $a_{ii}\le 0$4, and stronger $a_{ii}\le 0$5 conditions imply global asymptotic stability. The paper also gives higher-order extensions of classical winner-take-all and winner-share-all results, including multistability of boundary equilibria under strong competitive effects on losers and internal stability of winners’ subsystems. A common misconception is that higher-order terms necessarily destroy the classical Lotka-Volterra analytical structure; the results here are presented instead as natural extensions of many analogous classical theorems (Cui et al., 2024).

4. Random HO-LV communities and dynamic mean-field theory

For quenched random higher-order interactions, generating-functional analysis yields a dynamic mean-field description in the large-$a_{ii}\le 0$6 limit. The effective single-species process has colored noise, a retarded memory kernel, and self-consistent order parameters $a_{ii}\le 0$7, $a_{ii}\le 0$8, and $a_{ii}\le 0$9. Under a fixed-point ansatz, the stationary abundance satisfies

$b_{iii}\le 0$0

where $b_{iii}\le 0$1, $b_{iii}\le 0$2, and $b_{iii}\le 0$3 are summed effective parameters. The abundance distribution is therefore a clipped Gaussian: a Gaussian clipped at $b_{iii}\le 0$4, with the negative mass accumulating into a delta peak at $b_{iii}\le 0$5 for extinct species (Sidhom et al., 2024).

This DMFT yields closed self-consistent expressions for the mean abundance $b_{iii}\le 0$6, second moment $b_{iii}\le 0$7, integrated response $b_{iii}\le 0$8, and survivor fraction $b_{iii}\le 0$9. Linear instability of the fixed point occurs at

$$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$0

which reduces to $$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$1 for a single interaction order $$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$2. The analysis shows that increasing heterogeneity $$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$3 generally reduces stability and diversity, while more antisymmetric interactions stabilize the fixed point. It also identifies a new mathematical condition for a sudden onset of diverging abundances: for $$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$4, the self-consistent $$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$5 can attain a finite maximum $$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$6, beyond which no fixed-point solution exists even though $$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$7 remains finite at the boundary (Sidhom et al., 2024).

The higher-order random model behaves similarly to the pairwise model in some respects, but several differences are specific to $$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$8. More competition, meaning $$\dot{x} = R\,\mathrm{diag}(x)\big(\mathbf{1}+Ax+Bx^2\big),$$9, can increase linear stability and the diversity in the community, an effect not seen in the pairwise model. Fully antisymmetric pairwise interactions do not reach linear instability at finite $(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$0, whereas higher-order fully antisymmetric interactions can. The phase diagram is correspondingly more complex than in pairwise GLV, with stable, unstable bounded, and divergent regimes separated by instability and divergence boundaries (Sidhom et al., 2024).

A complementary random-ensemble approach with immigration reaches related conclusions from direct numerical integration rather than DMFT. In that setting, HOIs consistently shrink the collapse region across $(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$1 random ensembles, and larger pairwise variability $(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$2 lowers the mean HOI strength $(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$3 needed for coexistence. The effective community matrix

$(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$4

is state-dependent, so negative $(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$5 contracts spectra and can counteract destabilizing pairwise heterogeneity. This suggests that HO-LV stability is controlled not only by disorder statistics but also by the abundance-dependent reweighting of those statistics at coexistence (Kang et al., 29 Jul 2025).

5. Oscillations, quasi-periodicity, chaos, extinction, and spatial organization

In HOI-GLV systems with immigration, the coexistence equilibrium satisfies

$(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$6

and the Jacobian simplifies at equilibrium to

$(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$7

The diagonal term $(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$8 is always stabilizing, while the HOI term can move eigenvalues leftward or generate oscillatory regimes when complex conjugate eigenpairs cross the imaginary axis. Simulations exhibit stable equilibria, limit cycles, quasi-periodic torus trajectories, and intermittent chaos. Representative cases include a stable equilibrium for $(Bx^2)_i=\sum_{j,k=1}^n b_{ijk}x_jx_k$9, $m$0, $m$1, $m$2, $m$3; a limit cycle when the HOI mean is increased to $m$4; quasi-periodicity in HOI-only systems with $m$5, $m$6, $m$7, and $m$8 means near $m$9; and chaos at $r_i>0$00 or $r_i>0$01 in $r_i>0$02 communities, confirmed by period-$r_i>0$03 orbits, Poincaré maps, and strong sensitivity to initial conditions (Kang et al., 29 Jul 2025).

The same model reproduces empirical species-abundance distributions across birds, bats, bees, insects, plankton, and bacteria. Rank-abundance curves are S-shaped, Shannon entropies closely match empirical values, and K–S tests fail to reject similarity with $r_i>0$04 across the reported datasets, whereas pairwise-only GLV gives K–S $r_i>0$05-values below $r_i>0$06. The simulation protocol uses MATLAB R2018b ODE45 over $r_i>0$07 time units, initial conditions $r_i>0$08, an extinction threshold $r_i>0$09, and $r_i>0$10 random draws $r_i>0$11 per parameter pixel for coexistence fractions. In this setting, HOIs rescue coexistence, suppress collapse, and broaden the range of dynamical regimes available to the deterministic GLV framework (Kang et al., 29 Jul 2025).

Three-agent cyclic HO-LV models expose the same issues in a more microscopic setting. In the homogeneous symmetric case, the interior coexistence point is neutrally stable when equilibrating and polarizing three-agent reaction rates coincide, and the linearization yields a degenerate Hopf bifurcation. Finite-size effects then turn neutrally stable orbits into erratically diverging trajectories. A Kramers-Moyal expansion gives a Fokker-Planck approximation with diffusion constant $r_i>0$12, smaller than in the original cyclic Lotka-Volterra model by a factor of $r_i>0$13, and the Laplace transform of the extinction probability takes the form

$r_i>0$14

Extinction probabilities derived in this approximation agree well with Gillespie simulations, and typical extinction times scale with $r_i>0$15 (Palombi et al., 2019).

Spatial locality qualitatively changes the phase portrait. On a two-dimensional lattice with strictly local three-agent sampling, all bifurcations disappear and the reactive fixed point is stable for every choice of reaction rates. In metapopulation models with patch capacity $r_i>0$16 and pair-exchange mobility, rotating spiral waves reappear, the wavelength and front speed follow a complex Ginzburg-Landau reduction near onset, and three-agent chase reactions induce nonlinear diffusion terms such as

$r_i>0$17

The reported propagation instabilities, far-field breakup, and rectification patterns show that higher-order interactions affect not only reaction kinetics but also effective transport (Palombi et al., 2019).

6. Scope, terminological variants, and current limitations

The term higher-order is used in several related but nonidentical senses. In the explicit HOI literature, it refers to non-additive multispecies interaction terms such as $r_i>0$18 or more general tensors $r_i>0$19. In the nonlinear-feedback literature, higher-order effects arise because the per-capita growth depends on a saturating function of the interaction field,

$r_i>0$20

so a Taylor expansion of the Hill response generates quadratic and higher polynomials in abundances. For finite saturation level $r_i>0$21, unbounded growth is not possible, the stable fixed-point phase extends to much larger $r_i>0$22, decreasing $r_i>0$23 has a stabilizing effect similar to decreasing the symmetry parameter $r_i>0$24, and co-operation is no longer detrimental to stability (Sidhom et al., 2019).

A distinct literature uses higher-order to refer to non-integer derivative order rather than higher interaction order. In the two-dimensional fractional-order Lotka-Volterra system

$r_i>0$25

the coexistence equilibrium $r_i>0$26 is locally asymptotically stable iff $r_i>0$27 when $r_i>0$28, but the system can never have a stable focus, and for $r_i>0$29 both $r_i>0$30 and $r_i>0$31 are unstable for all $r_i>0$32. Here the higher-order feature is long-memory via Caputo derivatives, not explicit triadic interaction coefficients (Merrikh-Bayat, 2013).

Other usages are broader still. In the cosmological LV mapping for interacting dark energy, dark matter, matter, and radiation, the “higher-order” aspect refers primarily to higher dimensionality, with $r_i>0$33 and $r_i>0$34 bilinear LV systems and optional self-quadratic terms induced by a quadratic equation of state. The paper states explicitly that no triadic interaction terms of the form $r_i>0$35 are introduced. This makes the terminology nonuniform across fields, even when the Lotka-Volterra label is shared (Aydiner, 2016).

Across the ecological HOI literature, several limitations recur. Deterministic tensor and monotone-systems results focus on ODEs without environmental stochasticity or spatial structure; random-interaction models often assume dense random couplings with connectance approximately $r_i>0$36; extending monotone-system methods to more than two factions is nontrivial because permutation to Metzler form fails for at least three factions; and HOI coefficients are numerous, with $r_i>0$37 parameters for triads, which raises identifiability challenges from limited data. The proposed remedies in the literature are regularization, sparsity assumptions, experimental designs targeting triadic interactions, and equilibrium bounds from complementarity theory to guide calibration (Cui et al., 2024, Kang et al., 29 Jul 2025).

These distinctions clarify a frequent misconception. Not every nonlinear or higher-dimensional LV model is a higher-order interaction model in the strict HOI sense, but the broader family of higher-order LV dynamics now includes explicit HOI tensors, nonlinear response functions that induce effective higher-order terms, stochastic and random-coupling generalizations, and fractional-order memory extensions. What unifies them is that the classical pairwise additive LV closure is no longer sufficient to describe the interaction structure or the resulting dynamical regimes.