Polynomial Replicator Dynamics

• Polynomial Replicator Dynamics is a system defined on a simplex where the evolution of strategy frequencies follows a polynomial vector field.
• It extends classical replicator models by incorporating higher-order interactions, state-dependent payoffs, and tensor-valued contributions for richer dynamics.
• The framework supports integrable, conservative, and information-geometric formulations, bridging evolutionary game theory, ecology, and statistical inference.

Polynomial replicator dynamics are replicator systems on a simplex whose vector field is polynomial in the strategy frequencies. In the classical matrix case, the state $x \in \Delta^{n-1}$ evolves by

$\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$

so pairwise linear payoffs already produce a polynomial flow, while multi-player, tensor-valued, state-dependent, or spatially structured payoffs raise the polynomial degree without leaving the simplex. Recent work treats this polynomial structure not as a notational accident but as a unifying feature connecting stochastic imitation limits, higher-order ecological interactions, zero-sum recurrence, integrable competitive networks, information geometry, and canonical skew-symmetric representations of the same dynamics (Fontanari, 2024, Bratus et al., 7 Apr 2026, Yin et al., 29 Aug 2025).

1. Algebraic form and state-space structure

A general route to replicator equations starts from Kolmogorov-type equations for absolute abundances,

$\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$

and passes to relative frequencies

$$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$

When the interaction terms $g_i$ are homogeneous, the frequency dynamics are orbitally equivalent, after time reparametrization, to

$$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$

If $g_i(\mathbf u)$ is polynomial of degree $d$, then $\dot u_i$ is polynomial of degree $d+1$. The classical matrix replicator is the specialization $\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$0, giving a quadratic polynomial system on the simplex $\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$1 (Bratus et al., 7 Apr 2026).

In the two-strategy case, recent derivations often use the reduced scalar variable $\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$2, the frequency of strategy $\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$3, and obtain

$\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$4

Here $\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$5 and $\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$6 are expected payoffs and $\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$7 is only a time-rescaling constant. When payoffs come from a $\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$8 matrix

$\dot{x}_i = x_i\big((Ax)_i - x^\top A x\big),$9

the payoff difference $\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$0 is linear in $\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$1, while $\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$2 is quadratic, so the right-hand side is polynomial in $\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$3. In the multi-strategy formulation,

$\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$4

the simplex is invariant because $\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$5, boundary points are fixed whenever a coordinate vanishes, and interior equilibria satisfy payoff equalization conditions such as $\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$6 in the two-strategy reduction (Fontanari, 2024).

This algebraic viewpoint immediately separates polynomial replicator dynamics from arbitrary polynomial flows. The simplex constraint, the multiplicative factor $\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$7, and subtraction of mean fitness jointly enforce positivity, conservation of total mass, and the characteristic face structure of evolutionary state spaces (Fontanari, 2024, Bratus et al., 7 Apr 2026).

2. Microscopic origins in imitation dynamics

A central modern result is that standard replicator dynamics can be recovered as the deterministic limit of stochastic imitation based on instantaneous, rather than averaged, payoffs. In the well-mixed model analyzed in "Imitation dynamics and the replicator equation" (Fontanari, 2024), each time step chooses a focal individual and a model individual uniformly without replacement; if the model’s instantaneous payoff exceeds the focal’s, imitation occurs with probability

$\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$8

Under the assumptions of instantaneous payoffs, monotone imitation, linear dependence on payoff difference, well-mixed sampling, $\dot N_i(t) = N_i(t)\, g_i(\mathbf N),$9, and a large-population law-of-large-numbers limit, the macroscopic ODE is exactly the replicator equation after time rescaling. The same paper states that the derivation extends directly from pairwise games to $$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$0-player interactions, preserving the polynomial structure of the limiting dynamics (Fontanari, 2024).

This result corrects a common overstatement in the literature: the replicator equation does not require agents to compare long-run average payoffs or to imitate worse-performing strategies with positive probability. The 2024 proof shows that a simpler better-only rule based on one-shot realized payoffs already suffices. The same analysis also makes clear that if the imitation probability is replaced by a nonlinear function $$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$1, such as a Fermi rule, the limiting ODE generally ceases to be the standard replicator equation, even though it may remain polynomial for polynomial $$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$2 (Fontanari, 2024).

A related but distinct construction is the imitate-the-better-realization rule. Its mean dynamics depend only on ordinal comparisons of realized payoffs and yield quartic polynomial vector fields in the strategy frequencies. These ordinal imitative dynamics eliminate pure strategies iteratively strictly dominated by pure strategies and make strict equilibria locally stable, but they need not satisfy Nash stationarity or payoff monotonicity. In trivial two-payoff cases they reduce to the replicator equation up to a constant time change, and in symmetric Rock–Paper–Scissors they differ from replicator dynamics only by a positive polynomial speed factor; in other games they exhibit behaviors impossible under standard replicator dynamics (Loginov, 2019).

3. Higher-order interactions and degree escalation

The polynomial degree of replicator dynamics is controlled by the order of interaction encoded in the payoff law. For two-strategy $$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$3-player games, expected payoffs become polynomials in the frequency $$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$4, and the reduced scalar replicator equation has degree up to $$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$5 (Fontanari, 2024). This already moves beyond the quadratic matrix case without leaving the replicator form.

A more explicit higher-order construction appears in "Higher Order Dynamics in the Replicator Equation Produce a Limit Cycle in Rock-Paper-Scissors" (Griffin et al., 2023). There the fitness of strategy $$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$6 is

$$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$7

where $$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$8 encodes pairwise interactions and the symmetric matrices $$u_i(t)=\frac{N_i(t)}{\sum_{k=1}^n N_k(t)}, \qquad \sum_i u_i=1.$$9 are slices of a rank-three tensor describing triadic effects. The resulting replicator system keeps the usual multiplicative structure

$g_i$0

but now $g_i$1 contains quadratic terms in $g_i$2, and the vector field contains cubic, and effectively higher-order, contributions on the simplex (Griffin et al., 2023).

The generalized Rock–Paper–Scissors example in that paper is important because it alters the catalogue of admissible phase portraits. Standard three-strategy matrix replicator dynamics cannot generate a nondegenerate Hopf bifurcation with a limit cycle. After adding triadic terms through a payoff tensor, the interior equilibrium undergoes a subcritical Hopf bifurcation, and an unstable limit cycle appears around the Nash point. The paper numerically characterizes the parameter regimes where this occurs, showing that higher-order polynomial terms are not merely quantitative perturbations of pairwise selection; they can create dynamical objects forbidden in the classical quadratic theory (Griffin et al., 2023).

The broader implication is structural. In matrix games, polynomiality comes from linear payoffs. In tensor or multi-player games, polynomiality comes from contracting the state against higher-rank tensors or against group-composition polynomials. The replicator form survives, but the degree, bifurcation structure, and admissible invariant sets change substantially (Griffin et al., 2023, Fontanari, 2024).

4. Qualitative dynamics: recurrence, periodicity, permanence, and boundary selection

In zero-sum evolutionary games, the polynomial structure interacts with antisymmetry to impose strong global constraints. "From Darwin to Poincaré and von Neumann" proves a complete dichotomy: for zero-sum evolutionary games, replicator dynamics exhibits Poincaré recurrence if and only if there is an interior Nash equilibrium; if no interior equilibrium exists, then every interior trajectory converges to the boundary, and every strategy not in the support of any equilibrium vanishes in the limit (Boone et al., 2019).

The same paper shows that two degrees of freedom already suffice for periodicity. In three-strategy zero-sum games with an interior equilibrium, every interior orbit is periodic, not merely recurrent. This recovers the classical closed-orbit geometry of Rock–Paper–Scissors and places it in a broader recurrence theorem for zero-sum polymatrix and evolutionary games (Boone et al., 2019).

Polynomial replicator dynamics also includes canonical non-zero-sum regimes with very different asymptotics. In the 2026 survey chapter, independent replication produces quadratic dynamics with monotone mean fitness and convergence to the type with maximal replication rate; autocatalytic replication produces cubic dynamics with an unstable interior equilibrium and asymptotically stable vertices; and first-order hypercyclic replication yields a quadratic cyclic system that is permanent and, for $g_i$3, carries a stable limit cycle (Bratus et al., 7 Apr 2026). The hypercycle therefore furnishes a polynomial replicator system in which coexistence and sustained oscillation arise from catalytic structure rather than from zero-sum antisymmetry.

These results separate three notions that are often conflated. Periodic orbits in three-strategy zero-sum matrix games are not isolated limit cycles but recurrent closed trajectories. Stable limit cycles can appear in higher-order polynomial replicator systems such as triadic Rock–Paper–Scissors. Permanence, by contrast, concerns uniform repulsion from the boundary and can hold in systems such as the hypercycle even when the phase portrait is not Hamiltonian or zero-sum (Boone et al., 2019, Griffin et al., 2023, Bratus et al., 7 Apr 2026).

5. Integrable and conservative formulations

Some polynomial replicator dynamics are not merely recurrent or permanent but completely integrable. "Completely Integrable Replicator Dynamics Associated to Competitive Networks" studies tournament replicators with skew-symmetric interaction matrices $g_i$4, for which the replicator equation reduces on the simplex to

$g_i$5

This is a quadratic polynomial system and simultaneously a Lotka–Volterra system with skew-symmetric interaction matrix. For circulant tournaments and for an infinite family obtained by embedding such tournaments into one another, the paper constructs explicit conserved quantities, a quadratic Poisson bracket, and enough integrals in involution to prove Liouville–Arnold integrability. As a consequence, these flows are quasi-periodic on invariant tori, and the five-species Allesina–Levine competitive network produces quasi-periodic dynamics rather than convergence to an equilibrium (Paik et al., 2022).

A broader conservative framework is given in "Conservative Replicator and Lotka-Volterra Equations in the context of Dirac\textbackslash big-isotropic Structures" (Alishah, 2018). Starting from a replicator equation with a formal equilibrium, the paper provides an algorithm for constructing constants of motion of logarithmic–polynomial type. After a nonlinear coordinate change and a time reparametrization, the transformed dynamics becomes Hamiltonian with respect to a Dirac or big-isotropic structure. Through the replicator–Lotka–Volterra equivalence, this enlarges the known class of conservative LV systems beyond the usual skew-symmetrizable cases and permits predator–prey models with interactions among different predators and among different preys (Alishah, 2018).

Taken together, these results show that polynomial replicator dynamics supports several distinct conservative mechanisms: Poisson-integrable skew-symmetric tournament systems, Dirac-geometric Hamiltonian formulations for broader matrix classes, and classical KL-type invariants in zero-sum games. Conservation is therefore not confined to antisymmetric matrix games, even though antisymmetry remains a privileged special case (Paik et al., 2022, Alishah, 2018).

6. Information geometry and inference-theoretic interpretations

A different line of work interprets replicator dynamics as a geometric inference process rather than as a purely game-theoretic selection law. "The Replicator Equation as an Inference Dynamic" establishes the formal analogy between the discrete replicator update

$g_i$6

and Bayes’ rule

$g_i$7

In this dictionary, the population state is a probability distribution over hypotheses, fitness plays the role of likelihood, and the normalization by mean fitness plays the role of marginal evidence (0911.1763).

In continuous time, the same paper emphasizes the equation

$g_i$8

as a natural-gradient flow on the simplex endowed with the Fisher information, or Shahshahani, metric. For an interior ESS $g_i$9, the Kullback–Leibler divergence

$$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$0

is a local Lyapunov function if and only if $$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$1 is evolutionarily stable. The argument identifies

$$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$2

so the ESS inequality is exactly the negativity condition for potential information to decrease along trajectories (0911.1763).

The same framework also expresses replicator solutions as exponential-family trajectories,

$$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$3

with $$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$4 and $$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$5. This interpretation is particularly natural for polynomial replicator systems because polynomial fitness landscapes on the simplex become information-geometric flows driven by polynomial likelihood-like terms. The information-theoretic viewpoint therefore does not replace the polynomial one; it recasts it in a differential-geometric language that explains why KL divergence, natural gradients, and exponential parameterizations repeatedly appear in stability and convergence analyses (0911.1763).

7. Representation problems and modern extensions

A major recent development concerns identifiability. "On Zero-sum Game Representation for Replicator Dynamics" shows that the payoff matrix cannot, in general, be inferred from observed strategy trajectories alone, because distinct payoff matrices can induce exactly the same replicator vector field. More strongly, for every polynomial replicator dynamics

$$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$6

with polynomial $$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$7 satisfying $$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$8 on the simplex, there exists a skew-symmetric polynomial matrix $$\dot u_i = u_i\Bigl(g_i(\mathbf u)-\sum_{k=1}^n u_k g_k(\mathbf u)\Bigr).$$9 such that

$g_i(\mathbf u)$0

Hence every polynomial replicator dynamics admits a zero-sum polynomial representative. The paper constructs this representation explicitly and shows that equivalence classes of payoff matrices are substantially larger in the state-dependent polynomial setting than in the constant-matrix case (Yin et al., 29 Aug 2025).

Structured population extensions preserve polynomiality but alter the meaning of the terms. In "Replicator dynamics with diffusion on multiplex networks", diffusion is linear in absolute agent counts but becomes nonlinear when expressed in fractions $g_i(\mathbf u)$1. The exact diffusion term contains quadratic products of layer and node frequencies, and the usual assumption of constant local population size induces a hidden selective pressure favoring faster-diffusing strategies. The resulting system is therefore a network-structured polynomial replicator dynamics rather than a simple addition of Laplacians to the classical equation (Requejo et al., 2016).

Continuous-trait and continuum-strategy models furnish another extension. "Mathematical description of continuous time and space replicator-mutator equations for quadratic fitness landscapes" studies a continuous-trait replicator–mutator PDE with quadratic fitness and Gaussian mutation. For Gaussian initial data, the population remains Gaussian, and the dynamics closes exactly on the mean and covariance, yielding Riccati-type moment equations and explicit analyses of the flying kite effect, survival of the flattest, and tracking of fixed, moving, oscillating, or randomly fluctuating optima (Pathiraja et al., 2024). A related continuum model with a nonsymmetric time-dependent payoff operator admits self-similar probability-density solutions that converge to a Dirac mass as $g_i(\mathbf u)$2, showing that polynomial-style replicator structure also survives in infinite-dimensional PDE settings (Papanicolaou et al., 2014).

These developments broaden the subject in two directions. First, polynomial replicator dynamics is now understood as an equivalence class of vector fields rather than as a unique payoff-matrix specification. Second, the same polynomial architecture persists under ordinal imitation, network diffusion, higher-order interaction tensors, and continuous-trait mutation–selection PDEs, making the topic a common language across evolutionary game theory, ecological dynamics, online learning, and statistical mechanics of adaptive populations (Yin et al., 29 Aug 2025, Requejo et al., 2016, Pathiraja et al., 2024).