Nonlocal Generalized Models of Predator-Prey Systems

The method of generalized modeling has been applied successfully in many different contexts, particularly in ecology and systems biology. It can be used to analyze the stability and bifurcations of steady-state solutions. Although many dynamical systems in mathematical biology exhibit steady-state behaviour one also wants to understand nonlocal dynamics beyond equilibrium points. In this paper we analyze predator-prey dynamical systems and extend the method of generalized models to periodic solutions. First, we adapt the equilibrium generalized modeling approach and compute the unique Floquet multiplier of the periodic solution which depends upon so-called generalized elasticity and scale functions. We prove that these functions also have to satisfy a flow on parameter (or moduli) space. Then we use Fourier analysis to provide computable conditions for stability and the moduli space flow. The final stability analysis reduces to two discrete convolutions which can be interpreted to understand when the predator-prey system is stable and what factors enhance or prohibit stable oscillatory behaviour. Finally, we provide a sampling algorithm for parameter space based on nonlinear optimization and the Fast Fourier Transform which enables us to gain a statistical understanding of the stability properties of periodic predator-prey dynamics.


Introduction
Predator-prey systems have been a cornerstone in mathematical biology for many decades [4]. Standard textbooks on dynamical systems, differential equations and ecology provide a plethora of models that aim at capturing the interaction between a predator population Y and a prey population X. Examples for modeling the situation by ordinary differential equations (ODEs) are [6,5] (LV) where p i , k i are parameters and (LV) are the Lotka-Volterra equations and (RM) is the Rosenzweig-MacArthur model [32]. These two models are the most common examples of a large class of different models of the form where α > 0 is a parameter describing biomass conversion efficiency and the functions S, G and M represent prey growth, predation, and predator mortality, respectively. Because the parameter α can always be removed by scaling the variable Y and re-labelling the functions we will always assume α = 1 from now on. Generalized models [20,33] directly work with the formulation (1) without specifying functional forms for S, G and M. Previous works on generalized models [14,15,12,52] focused on analyzing the dynamics close to stationary states. Beyond the structure of the equations (1) this analysis requires only the assumption that steady states exist in the class of models under consideration. The central idea of generalized modeling is to parametrize all possible Jacobians that can be encountered in steady states in the class of systems under consideration. Using a specific renormalization procedure, one can define parameters that are easily interpretable (and often also directly measurable) in the context of an application. Applications of generalized models to ecology can be found in [19,18,17,20,21,3,47,22,48,50]. Let us emphasize that we do not claim that stability results from generalized models have not been observed before in some specific models; in fact, the literature on stability of planarpredator prey systems is very large. For instance, questions of local and global stability have been investigated in various predator-prey systems [13,26,49,40,36].
In the present paper we go beyond the previous analysis and study nonstationary dynamics in the context of generalized modeling. We extend the theory of generalized models to arbitrary periodic solutions in the context of the predator-prey system (1). We show that this mathematical extension of generalized models yields several new phenomena in comparison to generalized models for steady states. For example, we can define time-periodic generalized parameters of the predator-prey model and we prove that these functions obey a system of ODEs (a flow on moduli space). Using Floquet theory [7] and Fourier analysis [28] we derive analytical conditions for the solvability of the moduli space flow and obtain an analytical stability formula. In this context, a main result is that the stability formula for periodic solutions only depends on two constants that can be calculated via a discrete convolution. Using this formula we can identify parameters and conditions that enhance the stability of predator-prey cycles. Furthermore, we develop an algorithmic approach to sample the function space of parameters by solving an auxiliary optimization problem, which will be instrumental for future applications to larger systems.
The paper is structured as follows: In Section 2, we recall the necessary tools from steadystate generalized models, Floquet theory, and Fourier analysis. In Section 3, we calculate the generalized vector field for non-equilibrium solutions. In Section 4, we derive the flow on moduli space. In Section 5, we compute the generalized scale and elasticity functions for several specific functional forms to gain a better understanding how generalized and specific models link up. In Section 6, we use tools from Fourier analysis to derive algebraic conditions from the moduli space flow. In Section 7, we provide an analytical stability analysis of periodic orbits in generalized predator-prey models. Using this result we can identify which situations increase or decrease stability and interpret the results in an ecological context. In Section 8, we develop a sampling technique for generalized scale and elasticity functions that is based on solving an auxiliary optimization problem and the Fast Fourier Transform. We also use this sampling-based approach to improve our understanding of stabilizing and destabilizing factors of the predator-prey system. In Section 9, we conclude with a brief summary and outline the large range of applications and theoretical challenges that can be found in non-equilibrium generalized models.

Background
In this section we introduce essential tools and techniques that will be used throughout this work. Further, we use the opportunity to fix the notation. Below, we denote a general ordinary differential equation ODE by and assume always that F is sufficiently smooth. In the following we are going to recall the necessary tools from steady state generalized models, Floquet theory and Fourier analysis.

Generalized Models
Let us start by reviewing generalized modeling [20] for ODEs with equilibrium points. A detailed mathematical approach to generalized models can be found in [33]. For the present discussion we restrict ourselves to review generalized models in the context of a planar predator-prey system [21]. Such systems describe the interaction of a population of prey X and a population of predators Y. The prey population grows at rate S(X), predation occurs at rate G(X, Y ) and natural mortality of the predator at rate M(Y ). Denoting the prey density as X and predator density as Y we capture the dynamics by where S, M ∈ C r (R + , R + ) and G ∈ C r (R + × R + , R + ) are sufficiently smooth functions. Generalized modeling assumes that (3) admits an equilibrium point ( We normalize the equilibrium defining new coordinates This transformation moves the equilibrium to (x, y) = (1, 1). The next step is to normalize the rate functions A direct substitution of (4)-(5) into (3) gives where the prefactors of the form S(X * )/X * , G(X * , Y * )/X * , etc. represent normalized fluxes in the steady state and are also called scale parameters Since (x, y) = (1, 1) is an equilibrium point we know that the following holds: Therefore (6) can be re-written as The Jacobian at the equilibrium (x, y) = (1, 1) is then given by =: where ∂ x , ∂ y denote partial derivatives and we refer to the constants as elasticities. The scale parameters and elasticities are also referred to as generalized parameters.
In the following we will use the insight that every power law function corresponds to an elasticity that is identical to the exponent of the power law. For example, if we assume that the mortality M(Y ) is a linear function M(Y ) = KY then we find Hence we can relate the growth properties of the unspecified functions forms to the elasticities. The stability of the equilibrium (x, y) = (1, 1) can be inferred from the eigenvalues of J(1, 1) and hence only depends on the generalized parameters. This admits a bifurcation analysis of all steady state models of the form (3) in generalized parameter space. Despite the large class of models that one treats simultaneously it is often easy to interpret scale parameters and elasticities in applications [20]. Thereby a generalized model enables us to draw conclusions about a whole class of differential equations, for further examples see [22,48,46,42,12].
We note that generalized modeling can also be applied to equilibria for delay equations [27,33], spatially homogeneous states for partial differential equations [3] and to stochastic differential equations [33].

Floquet Theory
For analyzing the stability of periodic solutions in GM we resort to the framework offered by Floquet Theory. Suppose (2) has a period orbit γ(t) = γ(t + T ) with minimal period T . Let Σ denote a suitable (N − 1)-dimensional transversal section to Γ and consider the associated Poincaré map P : Σ → Σ. This map has a fixed point X γ ⊂ Σ associated to the periodic orbit γ i.e. P (X γ ) = X γ . Recall [8,34] that the stability of γ is determined by the N − 1 eigenvalues (or characteristic/Floquet multipliers) λ 1 , . . . , λ N −1 of the matrix DP (X γ ). If |λ j | < 1 for all j ∈ {1, . . . , N − 1} then the periodic orbit is stable, if there exists λ j such that |λ j | > 1 then the orbit is unstable and eigenvalues with |λ j | = 1 signal bifurcations under parameter variation. We can study the stability of γ by considering the non-autonomous linear variational equation where A(t) is periodic. An N × N matrix M(t) that satisfies is called the fundamental matrix solution of (13). The constant matrix M(T ) is called the monodromy (or circuit) matrix. It has eigenvalues 1, λ 1 , λ 2 , . . . , λ N −1 where the trivial eigenvalue 1 is associated to the direction tangent to the periodic orbit that links the variational equation to the Poincaré map P . Furthermore, the Liouville formula holds. Floquet's theorem states that there exists a T -periodic coordinate change C(t) and a constant matrix R such that M(t) = C(t)e tR .
Since M(0) = Id it follows that C(0) = C(T ) = Id and we find that the monodromy matrix can be expressed as An elegant explicit formula for the Floquet multiplier from (15) is only available for N = 2.
In general the computation of Floquet multiplier thus requires numerical approaches, which typically start with computing the periodic solution with a suitable boundary value method such as collocation or finite differences [34,9]. The variational equation (14) is solved on suitable sub-intervals of the periodic orbit discretization as an initial value problem to obtain M(T ). The eigenvalues of M(T ) are then obtained yielding the Floquet multipliers. Although, in certain circumstances, such as large multipliers, the computation can be numerically problematic [10,37]. Let us point out that Floquet theory has not been widely applied in the context of ecology [30] although it is a standard tool in the mathematical theory of dynamical systems [8]. Klausmeier [30] suggests that "Floquet theory [is] a useful tool for studying the effects of temporal variability on ecological system". In the context of our approach, Floquet theory is not only a tool for a particular model but we will also show that it nicely extends to generalized models.

Fourier Series
Since we work with periodic solutions to ODEs and also other time-dependent periodic functions we briefly recall basic facts about Fourier series to fix normalization constants and notation. Assume that f : R → R is T -periodic so that we can identify the domain of f as the circle R/(T Z) ∼ = S 1 . We can formally write the complex Fourier series F [f ] of f as follows: where the Fourier coefficientsf (k) arê Observe thatf (k) =f (−k), where the overbar denotes complex conjugation. Further, f (0) = 1/T T 0 f (t)dt is the time average of the periodic function. The convergence question F [f ](t) → f (t) is extremely intricate depending on the properties of f [53,28]. In the following, all functions we are going to approximate by Fourier series will be in C r (S 1 , R) for some sufficiently large r or even r = ∞. In this case, uniform convergence is immediate. A very important practical result in this context is to control the Fourier coefficients.
Theorem 2.1 is a version of the Riemann-Lebesgue Lemma for smooth functions and can provide an extremely rapid decay of the Fourier coefficients. This justifies (for the smooth case!) dropping higher-order terms |k| > κ for some rather small suitable κ ∈ N. The remaining sum is expected to be a good approximation to the original periodic function f . We write We remark that it can be convenient to re-write the complex Fourier series (16) as a real Fourier series where the real Fourier coefficients relate to the complex ones bŷ for k ∈ N 0 . Another important tool in Fourier analysis we will need are convolutions. Recall that the discrete convolution of two periodic functions f and g is defined as Obviously the convolution operator ' * ' is associate, commutative and distributive.

Non-Equilibrium Planar Predator-Prey Systems
We return to the planar predator-prey system (3) from Section (2.1) given by Denote the vector field of (17) by F (X, Y ). The vector field is only considered on the first (positive) quadrant F : R + × R + → R 2 as predator-prey densities are assumed to be non-negative. We want to analyze the class of vector fields (17) under the assumption that it admits a non-equilibrium orbit that is bounded as |t| → ∞. From an ecological point of view the most interesting case are limit cycles, so-called predator-prey cycles. We assume that (17) has a periodic orbit γ(t) = (γ 1 (t), γ 2 (t)) with period T . The definition of the model implies that γ i > 0 for i ∈ {1, 2} and all t.
In the following, we are going to slightly extend the notation employed already in Section (2.1) by re-using names for variables and generalized parameters. As in the case of equilibria one can consider a normalizing coordinate change which maps the periodic orbit to the point (x, y) = (1, 1) =: 1. The ODEs (17) and the product rule imply Therefore the new equations can be written as In analogy to the equilibrium case we introduce normalized functions and define the scale parameters which are now time-dependent T -periodic scale functions. We will often suppress the timedependence in the notation and just write, for instance, β s instead of β s (t). Using (19)- (20) in (18) we find For applying Floquet theory we linearize (21) around the limit cycle which yields the matrix We can re-write A(1; t) in terms of the more familiar elasticities, leading to where the four time-dependent T -periodic elasticity functions are The periodicity and time-dependence becomes more apparent once we write out the detailed definitions, for example The previous calculations show that we can introduce replacements for the generalized parameters for equilibrium points in the context of periodic orbits. In particular, the scale parameters and elasticities become time-dependent and periodic. The term "generalized functions" is already used in a different context [51]. Therefore, we refer to elasticity functions and scale functions directly. To analyze the stability of the periodic solution we use Floquet theory (see Section 2.2). For planar systems the stability of the periodic orbit is determined by computing the only non-trivial Floquet multiplier λ. Liouville's formula implies that We can thus express the Floquet multiplier as a function depending on elasticity and scale functions. This is analogous to writing the eigenvalues of the Jacobian as functions of the generalized parameters in the equilibrium case.

The Moduli Space Flow
In analogy to the generalized exploration of local dynamics, the stability of the limit cycle can be studied by assuming plausible values for the generalized parameters (here, scale and elasticity functions). The value of generalized models lies in their ability to cover the whole range of possibilities that are plausible in the system. For an unbiased analysis it is essential that we consider only those values of parameters that are consistent with the set up of the system. For instance, in case of equilibrium generalized models we must demand that the parameter values which we assume do not preclude the existence of an equilibrium solution in the class of systems. Likewise, only those scale and elasticity functions should be considered which are mutually consistent and thus could arise in at least one example system in the class of models under consideration. To understand this problem we briefly go back to the equilibrium scenario (see Section 2.1). Suppose we just choose a set of generalized parameters where we assume that all parameters are positive. One natural question is if there exist specific functions S, G and M that lead to the generalized parameters (23).  (7) and (12) hold i.e. there exists a differential equation of the form (3) that has the given set of generalized parameters.
Using a slight modification of this approach we define G(X, Y ) = X g * x Y g * y and get g y = g * y as well as g y = g * y . We also must have β s = β 1 = β * 1 and β m = β 2 = β * 2 which translates into the conditions We can always choose p 1 and p 2 to satisfy (C2) and (C4). Then we can use X * and Y * to satisfy (C1) and (C3). The result follows.
Although there are certainly many other ways for constructing functions that are consistent with a given set of scale and elasticity parameters, already the existence of one such set of functions proves that the assumed parameter values could be encountered in the class of models under consideration. This observation is of central importance for sampling procedures, by which high dimensional generalized models are typically analyzed [46,22].
For non-equilibrium systems the situation is different since one has to ask whether a whole set of given functions can potentially arise from a system of the form (17).
Proof. We start by deriving the equation for β s . We know that β s = S(γ 1 )/γ 1 and direct differentiation with respect to time via the quotient and chain rules gives Noting that s x = γ 1 S ′ (γ 1 )/S(γ 1 ) and using the definition of β s the equation transforms to Since (γ 1 , γ 2 ) is a trajectory of (17) we must have γ ′ 1 = S(γ 1 ) − G(γ 1 , γ 2 ). This implies upon substitution into (25) that which is the first equation in (24). The calculation for β ′ m is similar. For β ′ 1 we find The calculation for β ′ 2 is similar to the one for β ′ 1 . The main conclusion is that the elasticities and scale functions which parametrize the ODE (21) satisfy an ODE themselves. Because one often uses the terms "parameters" and "moduli" interchangeably, Theorem 4.2 implies that the time-dependent parameters of generalized models generate a flow on moduli space. The following remark describes the relevance of this viewpoint in some other research areas.
Remark: The term "moduli space" is perhaps most commonly used in algebraic geometry which, broadly speaking, is the study of solutions of algebraic equations [24,25]. The solutions form algebraic varieties (e.g. curves). Often suitable parametrized families of algebraic varieties again have the structure of an algebraic variety, where the latter object is the moduli space of parametrized families. The study of the geometry of moduli spaces has also been transported into different branches of physics such as quantum field theory [1]. In dynamical systems theory, a classical moduli space argument is made in the renormalization analysis of parametrized families one-dimensional maps [23], where the renormalization transformation can be viewed as a map generating a dynamical system on moduli space. A very similar situation occurs for billiard dynamics where the so-called Teichmüller flow on the space of lattices appears [44,38].
We note that the positive quadrant is an invariant set for (24) which means that this property lifts from the predator-prey family of vector fields to the moduli space. From Theorem 4.2 we can immediately infer a condition for the existence of a generalized model with given elasticities. Note that the existence of periodic solutions in Corollary 4.3 is only a necessary condition for the existence of a generalized model. We also observe that for the equilibrium case the conditions β 1 = β s and β 2 = β m give β ′ s = β ′ m = β ′ 1 = β ′ 2 = 0 consistent with steady state generalized modeling. It is also interesting to ask what happens if we do not specify the elasticities.
Taking the idea of deriving a differential equation one step further we consider Applying similar substitutions as in the proof of Theorem 4.2 we obtain Similar calculations can be carried out for g ′ x , g ′ y and m ′ y . These suggest that specifying a suitable scaled version of second partial derivatives of S, G and M will provide a system of eight ODEs. This procedure could be continued iteratively. It is interesting to note that closing a system of ODEs at a given order is a problem that also occurs in the context of moment closure for networks [29,16] and for moment equations of stochastic differential equations [45,11]. To get a better understanding of the stability of non-equilibrium generalized models and the flow on moduli space we proceed to consider a few typical specific functions S, G and M that appear in predator-prey models.

Specific Functions
In this section we calculate the generalized elasticity and scale functions for several wellknown predator-prey models. All the model parameters k l (for l ∈ N) we are going to use below are positive due to modeling considerations. We start with the growth of the prey S(X). Typical choices are (growth with strong Allee effect), 0 < k 2 < k 3 .
We start by looking at linear growth. We find For estimating the impact of linear growth on stability we consider the formula (22) and view it as a product of exponentials. The term involving β s and s x is Therefore, a linear prey growth does not contribute to the non-trivial Floquet multiplier because s x = 1 and exp( T 0 0dt) = 1. Different types of polynomial growth with a single term can be treated analogously since for S(X) = k 1 X p we find where the elasticity function coincides with the result for equilibrium generalized models. This allows us to write . Considering (26) we find that increasing p increases the Floquet multiplier and therefore has always a destabilizing effect, whereas decreasing p has a stabilizing effect. For logistic growth we obtain This implies Considering (26) again we find where the integral is positive because k 2 > 0. This means that increasing k 2 or increasing T 0 γ 1 dt will promote stability as the Floquet multiplier will move closer to 0. For logistic growth increasing k 2 corresponds to decreasing the carrying capacity k 1 /k 2 of the population. This can be interpreted as a manifestation of the paradox of enrichment [43,17] which captures the observation that increasing the carrying capacity generally has a destabilizing effect on attractors observed in ecological models. Furthermore, the expression obtained for logistic grows permits discussion of the contribution of the shape of the limit cycle to stability. For t ∈ [δ 1 , T − δ 2 ], where δ 1,2 > 0 are small, we find that which implies that the integral T 0 γ 1 (t)dt is small as well. Therefore limit cycles where the number of prey is extremely small for long times are not expected to be a basis for a stable ecosystem. For the Allee effect we find Considering the contribution of this term to the Floquet multiplier yields Increasing k 2 and/or k 3 will decrease stability. This is natural as these parameters represent the threshold to growth and the carrying capacity, providing another example for the paradox of enrichment. Note that the shape of the limit cycle can influence stability. In particular, the same conclusion to assumption (28) holds. In the case of the Allee effect the de-stabilization effect for long periods of low prey density even enters quadratically in the term T 0 k 1 2γ 2 1 dt. This confirms the intuitive conclusion that imposing a threshold to growth is a de-stabilizing factor for non-equilibrium systems when the prey density is small.
We proceed to consider the mortality of the predator. A very common functional form used in a large number of models is so-called density independent (linear) mortality Terms occurring in the Floquet multiplier . . . (Holling type I) Therefore, linear predator mortality has no effect on the stability of the periodic solution.
The interaction term between prey and predator is usually the most complicated and debated choice for the model. Some common choices are considered in Table 1. The observation that β 2 [g y − 1] vanishes for all functions considered in Table 1, can be directly linked to the ecological assumption that predators hunt independently of each other. The functions that are therefore used in practice are generally linear in the density of predators and the impact of predator dependence on stability vanishes. The same assumption cannot generally be made for prey dependence of predation, leading to more complex expressions for the impact on stability.
Therefore, we are going to make the assumptions g y = 1 and m y = 1 (30) from now. Regarding the Floquet multiplier formula (22) the assumptions (30) simplify the situation to investigating The influence of g x and β 1 on stability is not obvious since there is a non-trivial interaction with the shape of the limit cycle. The flow on moduli space given by (24) simplifies to where we can view β m as a parameter and simply drop the last equation.  Figure 1: Dynamics in a specific example. (a) Stable periodic orbit γ(t) of (33) (solid black) and two other trajectories (dashed magenta) with initial conditions marked by stars; the parameters are given in (34). Five points (black dots) are shown on the limit cycle for orientation purposes which are equally space over one period. (b) Scale functions in moduli space (black) for γ solving (32); a trajectory (solid magenta) with slightly perturbed initial conditions is also shown where the same elasticities as for the periodic orbit were used for numerical integration. (c) Time series of β s (t) for part (b).(d) Time series γ 1 (solid black) and γ 2 (dashed black). (e) Scale functions associated to γ: β s (t) (red), β m (t) (green), β 1 (t) (blue) and β 2 (t) (cyan). (f) Elasticity function associated to γ: s x (t) (red), g x (t) (blue) and g y (t) (green); note that m y = 1 = g y .
Example 5.1. For gaining an intuitive understanding one can consider the flow on the moduli space in a specific example. The combination of logistic prey growth, Holling-type-II interaction and linear predator mortality gives us the Rosenzweig-MacArthur predator-prey model that can produce periodic solutions where we use the parameters Figure 1 shows that integrating a slightly perturbed initial condition trajectory does seem to diverge from the exact periodic solution in moduli space. Furthermore, even for a classical planar predator-prey system, the scale and elasticity functions are quite complicated for non-equilibrium solutions. In fact, prescribing the elasticities is much more difficult than just picking a set of fixed parameters for equilibrium generalized models.
To verify the necessary condition from Corollary 4.3 for periodic solutions we must ask for solvability of the boundary value problem (BVP) where β(t) := (β s (t), β 1 (t), β 2 (t)). It is well-known that BVPs can have one, many or no solutions [2]. Furthermore determining solvability conditions is usually not easy and even using numerical methods may be dangerous; for example, if a numerical algorithm fails to provide a solution to (35) this may just be due to the numerical problems that can arise when solving BVPs [2].

Fourier Decomposition
The previous discussion of specific functions and the Rosenzweig-MacArthur model motivate the need for a more concrete version of the moduli space conditions (35) and of the Floquet multiplier (31). The natural step is to use a decomposition of the periodic functions into Fourier series; see Section 2.3. Using discrete convolution we can easily re-write the problem (35) on moduli space.
for all k ∈ Z where we have also used the notation1(0) = 1 and1(k) = 0 for k = 0 and employed the obvious definition for addition of infinite sequences.
Proof. To complete the proof we only have to recall another basic fact from Fourier analysis. For two T -periodic functions we f , g we have This formula for Fourier coefficients of products of functions yields the right-hand side of equation (36) as a direct consequence of Theorem 4.2. The left-hand side of equation (36) follows from direct differentiation which is allowed since all our periodic functions are assumed to be sufficiently smooth; see Section 2.3.  Figure 2: Absolute value of the first nine Fourier coefficients (|k| ≤ 4) associated to the stable periodic orbit γ(t) of (33); the parameters are given in (34). The coefficients of the phase space coordinates as well as the generalized elasticity and scale functions are shown.
Since (36) is an infinite set of algebraic equations it may look like we have not considerably simplified the problem of finding scale functions that are consistent with prescribed elasticities. However, the rapid decay of Fourier coefficients provided by Theorem 2.1 allows us to approximate the solution of (36) by focusing on the first few harmonics with |k| ≤ κ ≪ ∞. Example 6.2 (Example 5.1 continued). Just for illustration purposes we look a the Fourier coefficients of generalized scale and elasticity functions in an example. Figure 2 shows the results for the Rosenzweig-MacArthur model from Section 5 with κ = 4. We can clearly see that the Fourier coefficients decay very rapidly; it is also interesting to observe thatŝ x (k) is bimodal for logistic growth whereas the other coefficients show a uni-modal distribution for the first few harmonics. For the Rosenzweig-MacArthur model the algebraic relations (36) on the Fourier coefficients become  (37); parameter values used are given in (34). The black coefficients (lines shifted slightly left) are the coefficients of the derivatives β ′ s , β ′ 1 and β ′ 2 and the green coefficients (lines shifted slightly right) are associated to the periodic functions on the right-hand side of (37). The agreement of the two sets of coefficients is clearly visible. Figure 3 shows the values of the Fourier coefficients for the Rosenzweig-MacArthur example where we see that the algebraic conditions (37) are satisfied as proven in Proposition 6.1. Furthermore, it is evident that due to the convolution a wider support κ M is necessary i.e. the algebraic equations (36) must be satisfied for |k| ≤ κ M where κ M > κ and κ is our truncation for the Fourier coefficients of the phase space periodic orbit.

Stability Analysis
To get a better understanding of stability we can also use the Fourier series approach to re-express the Floquet multiplier (31) given by The next results shows how the different Fourier coefficients enter in formula (38).
Theorem 7.1. For the non-equilibrium generalized predator-prey model with g y = 1 = m y the single Floquet multiplier of a T -periodic orbit is given by i.e. whether |λ| > 1 or |λ| < 1 depends only on the difference of two zeroth-order Fourier coefficients C 1 and C 2 that arise from two discrete convolutions.
Proof. We start by looking at the first summand in integral in (38) which gives where the last step follows from the fact that T 0 e 2πikt/T dt = 0 for k = 0. From this calculation we find C 1 and in a similar way also C 2 . Using the two factors C 1,2 and m y = 1 = g x in (38) the result (39) follows. Since T > 0 the modulus |λ| only depends on the difference of C 1 and C 2 ; if C 1 − C 2 > 0 then |λ| > 1 and if C 1 − C 2 < 0 we obtain |λ| < 1.
Before we consider in more detail the dependency of stability on C 1 and C 2 we briefly investigate the influence of the period T . Although T does not effect the stability of a periodic orbit directly it does have an interesting biological interpretation. If T ≫ 1 then the period amplifies stability and instability. For example, when C 1 − C 2 > 0 then a long period moves the multiplier even further away from |λ| = 1 and trajectories near the unstable periodic orbit will escape vary quickly. On the other hand, if C 1 − C 2 < 0 and |λ| < 1 then a very large period T moves the multiplier even closer to the super-attracting case 0 ≤ |λ| ≪ 1. A very short period 0 < T ≪ 1 has the effect of moving the multiplier very close to |λ| ≈ 1. This means that when the periodic orbit is unstable, it will take a very long time to escape from it. The last effect can be interpreted as inducing meta-stability i.e. when the period of the predator-prey cycle is short then the predator-prey system stays near a metastable state for a long time although it is eventually unstable. This could lead to the conjecture that fast oscillations could be beneficial to survival for predator-prey populations during periods when external parameters entering C 1 and C 2 drive the system, potentially only temporarily, to a state when |λ| > 1.
As a next step, we want to understand better how the Fourier coefficients of β s , β 1 , s x and g x influence C 1 and C 2 .
In practice, we never use the infinite sum formulas from Proposition 7.2 but truncate them at a finite order. Using the explicit formulas for C 1 and C 2 we can directly draw several conclusions regarding periodic solutions depending on generalized scale and elasticity functions (recall: we still use g y = 1 = m y ). If all Fourier coefficients of higher-order k ≥ 1 are small, then stability of periodic solutions is dominated by the terms and Since the scale functions are always positive the time averages of the elasticity functionŝ s x (0) andĝ x (0) determine the signs of C 1 and C 2 . Therefore, average sub-linear elasticitŷ s x (0) < 1 and average super-linear conversionĝ x (0) > 1 enhance stability. In ecological terms 0 ≤ŝ x (0) < 1 means that, on average, the prey growth should be limited by external factors andŝ x (0) > 1 means that, on average, the predation rate should be sensitive to prey abundance; see also [21] for an interpretation of the generalized parameters for the equilibrium case. Both conditions make intuitive sense: if the prey grows without external limitation then solutions may be expected to diverge from a periodic solution and become unbounded while insensitivity of predation to prey growth could potentially drive a system to extinction. Of course, also the inverse relationships hold so thatŝ x (0) > 1 andĝ x (0) < 1 act towards de-stabilization. In this context, the scale functions act as amplifiers. For example, if C 1 < 0 and C 2 > 0 then a large average growth rateβ s (0) ≫ 1 and a large conversion ratê β 1 (0) ≫ 1 will enhance stability even more since the Floquet multiplier moves closer to the super-attracting regime λ ≈ 0. In this case, initial conditions will be attracted much quicker to a stable limit cycle. If the leading-order terms between growth and predation balance and the stability properties are dominated by higher-order harmonics. The leading-order terms also become irrelevant for elasticity functions which average close to onê In this scenario we have to focus on the relations between the higher-order Fourier coefficients of β s and s x as well as β 1 and g x . Let us assume for simplicity thatŝ x (0) = 1 =ĝ x (0) so that we can focus on the higher harmonics. Then stability enhancing conditions are Figure 4 depicts several different situations in the complex plane for the first two higher-order harmonics ofβ s (k) andŝ x (k) (k = 1, 2). In Figure 4(a) the first two higher-harmonics are in "anti-phase" so that the angles between the coefficients are separated by π. This means that In such a situation, we expect that C 1 < 0 by disregarding higher orders so that stability is enhanced. Figure 4(b) shows the situation where there is only a small phase difference between the coefficients ("in-phase") which gives There is also a possible situation where a competition between the different order harmonics arises as illustrated in Figure 4(c). We can now also give an ecological interpretation of these conditions. Stability C 1 < 0 is enhanced if s x (t) and β s (t) oscillate with a phase separation near π which means that a period of high sensitivity of prey abundance should coincide with a period of low prey growth and vice versa. Note that these conditions also make sense intuitively and suggest that prey growth is most efficient if there is a small number of prey and there are no limiting factors from the environment. Similar considerations also apply to the stability enhancing condition C 2 > 0. A small phase separation between β 1 (t) and g x (t) increases stability of the predator-prey limit cycle. Observe that g x (t) can be interpreted as the dependence of predation on prey abundance and β 1 (t) as a predation rate (normalized by the total number of prey) [21]; the stability conditions mean that a high predation rate should coincide with a high dependence of predation on prey abundance. In other words, if the dominating factor to prey abundance is predation then it is good for the predator to hunt a lot to increase stability of the limit cycle.
Note that although the conclusions stated above seem to be "obvious" in an ecological context, it is by no means clear how to prove them. That they can be obtained by an analysis of nonlocal generalized models underlines the applicability of the approach.

Sampling
Recall that due to Proposition 4.1 it was straightforward for equilibrium generalized models to choose a set of generalized parameters, just random sampling produces a set of parameters that is consistent with at least one specific model. Random sampling of generalized parameters has been exploited to correlate different aspects of the dynamical system to stability [42,46]. For non-equilibrium systems we must certainly check the necessary condition from Corollary 4.3. One possibility is the following algorithm which allows sampling of elasticity and scale functions: (A1) Prescribe a set of T -periodic elasticity functions by their Fourier coefficients. For simplicity we will always choose T = 1 and assume g y = 1 = m y .
(A2) Choose a truncation order κ M for the algebraic system (37) so that the necessary condition reads for |k| ≤ κ M . where we view c s , c 1 and c 2 as vectors of dimension 2κ M + 1 and real and imaginary parts are applied component-wise.
(A5) Observe that F (X 0 ) = 0 if and only if the Fourier coefficients encoded in X 0 satisfy the algebraic equations in (A2). Therefore we can attempt to solve the following optimization problem with a random initial condition, say x = x l .
Solving the optimization problem for different random initial conditions is expected to yield different values for X m that solve the algebraic constraint in (A2). This means that we get a set of Fourier coefficients {X m (l)} L l=1 where L denotes the sample size and the index l ∈ N indicates the dependence on the initial condition.
The main technical difficulty of the algorithm (A1)-(A5) is that it involves the solution of the optimization problem (42). This is computationally much more expensive than the direct random sampling for equilibrium generalized models. It is known [41] that the main computational cost in optimization is often given by the difficulty of the function evaluations of F (x). For our case, this seems to be the case since we have to compute several discrete convolutions to evaluate F (x). However, the convolution computation is inexpensive due to the Fast Fourier Transform [31].
Now we want to demonstrate that the algorithm can be used for a sampling analysis of stability similar to the one used in [22]. Let us point out that we do not attempt a full detailed statistical analysis here but that we only aim at a proof-of-principle. We solved (42) for 110000 initial conditions for κ M = 2 using a standard algorithm for nonlinear optimization [39,35]. Each sequence of Fourier coefficients in the initial condition consists of five real numbers e.g. β s (0), Re(β s (1)), Im(β s (1)), Re(β s (2)), Im(β s (2)), (43) which were sampled uniformly and independently from the interval [0.5, 1.5]. We discarded all solutions of the optimization algorithm that did not satisfy the positivity condition which is required by the definition of the scale functions and the invariance of the positive quadrant for the moduli space flow. The 63587 remaining solutions x m (l) satisfied the the optimization problem (and therefore the moduli space flow) at least up to a tolerance of 10 −4 i.e. |x m (l)| < 10 −4 for all l; the average value was E[x m (l)] ≈ 1.73 · 10 −6 . We have also calculated the single Floquet multiplier λ l associated to each solution using Proposition 7.2. Re(β s (2)) Im(β s (1)) Im(β s (2)) Figure 5: Histogram of the 5·63587 Fourier coefficients for β s obtained from the optimization of (42) with uniformly sampled initial conditions (43). The columns show the five different real numbers with their observed number on the vertical axes. The first row shows coefficients associated to a stable Floquet multiplier and the second row those with an unstable Floquet multiplier. Observe that the number of stable coefficients is substantially larger than the number of unstable ones. Figure 5 shows some of the output of the computation. We plot the Fourier coefficients associated to the scale function β s . The top row in Figure 5 corresponds to coefficients with stable periodic orbit (|λ| < 1) and the bottom row to coefficients with an unstable periodic orbit (|λ| > 1). We see that, despite the initial uniform sampling, the results for each coefficient of β s closely resemble normal distributions. The same observation also applies for the other scale and elasticity functions. In total we find that 37873 solutions associated to a stable multiplier and 25714 unstable ones. From this discrepancy one may either conjecture that the moduli space flow constraint could bias ecosystem towards stability or that our choice of initial uniform random sampling over a particular region in parameter space causes the bias towards stability.
In Table 2 we list mean and variance for each coefficient. Several observations can be made based on Table 2. The scale functions β 2 and β m have a much bigger variance than β s and β 1 . This could indicate that the prey growth rate and the prey-per-capita predation rate have to obey much smaller ranges in ecosystems compared to the predator-per-capita rates describing consumption and mortality. It is also interesting that the mortality rate β m allows for much larger amplitude higher-order harmonics whereas e.g. |β s (2)| is always comparatively small. The elasticities show no consistent variance decay towards higherorder harmonics although the coefficients themselves seem to decay. From the ecological perspective this suggest that predator-prey systems may exhibit a wide diversity in terms of sensitivities s x and g x . β s(k)βs(0) Re(βs(1)) Im(βs (1)) Re(βs (2)) Im(βs (2) (1)) Im(βm (1)) Re(βm (2)) Im(βm (2) (1)) Im(ŝx (1)) Re(ŝx (2)) Im(ŝx (2) (1)) Im(ĝx (1)) Re(ĝx (2)) Im(ĝx (2)  To understand how the different coefficients relate to stability we calculate the Pearson correlation coefficient. For two vectors of observations {a l } and {b l } it is defined as . Figure 6 shows r(a, λ l ) where a is a sequence of real or imaginary parts of the Fourier coefficients e.g. {a l } = {Re(β s,l (k))}. One important conclusion to draw from the correlation coefficients is that although a Fourier coefficient does not appear in the stability formula for the Floquet multiplier it may still correlate positively or negatively with stability. For example,β 2 (1) andβ 2 (2) show a negative correlation with Floquet multiplier. This effect can be caused by the fact that the scale and elasticity functions are not independent i.e. they are related via the moduli space flow.
It is very important to observe that we can recover conclusions, which we found already analytically in Section 7, from the statistical analysis. For example, the coefficientĝ x (0) correlates negatively with stability which means that decreasing it increases the Floquet Reβs (1) Imβs (1) Reβs (2) Imβs(2) Figure 6: Pearson correlation coefficient r = r(a, b) between stability and the different harmonics; positive correlation is indicated in blue and negative correlation in black. Each panel represents correlation for the Fourier coefficients from left to right. For example, the top left panel shows the values r(λ l ,β s,l (0)), r(λ l , Re(β s,l (1))), r(λ l , Im(β s,l (1))), r(λ l , Re(β s,l (2))), r(λ l , Im(β s,l (2))) from left to right where λ l is the Floquet multiplier with index l. multiplier and acts towards destabilization. This is precisely the result we have already obtained analytically in Section 7. Let us point out again that the basic statistical analysis we have provided is incomplete but that it definitely does show that the proposed sampling techniques based on the FFT, optimization and correlations can help to understand stability of periodic solutions.

Outlook
In this paper we have extended the method of generalized modeling from equilibrium to non-equilibrium systems. This extension has been achieved in the context of a classical predator-prey system with periodic solutions. The main re-scalings and definitions from the equilibrium case can be carried over to periodic orbits. However, the resulting generalized ODEs differ from the steady state case in several respects. The algebraic form is different due to the time dependent re-scalings and also the generalized parameters become timedependent elasticity and scale functions. The Jacobian A(t) of the system has to be analyzed using Floquet theory that describes the stability of periodic orbits. For planar vector fields we have been able to use Liouville's formula λ = exp T 0 T r(A(t))dt which facilitated several analytical calculations. We have discovered that the generalized elasticity and scale functions have to satisfy a flow moduli space. Then we used Fourier analysis to find computable conditions from the moduli flow. Discrete convolutions turned out to be the key to stability analysis providing explicit interpretable stability results. In the last part of the paper, we suggested a sampling algorithm that uses optimization methods to find elasticity and scale functions that satisfy the (algebraic) moduli space flow. During our analysis we have also obtained several ecological conclusions about arbitrary predator-prey models that can be written in the generalized form (3).
In principle, we can extend the theory described here without any technical problems to limit cycles in N-dimensional systems for N > 2; see also [33] regarding generalizations to R N in the equilibrium context. The main difference in R N will be that we have to compute several the Floquet multipliers numerically since Liouville's formula only provides the product of the eigenvalues.
One can also consider a generalization to non-equilibrium system beyond periodic orbits. For example, the generalized model (21), as well as Theorem 4.2 on the moduli space flow, carry over directly to other situations such as homoclinic trajectories [34] or chaotic dynamics [23]. For instance, instead of posing periodic boundary conditions of the form i.e. that we have asymptotic limits of the generalized elasticity and scale functions to their value at a saddle-point equilibrium. For chaotic dynamics one must search for aperiodic bounded trajectories in moduli space. Note that this raises interesting mathematical as well as application questions. For example, the moduli space flow may provide new insights when a dynamical system may be chaotic. Finally, also our sampling analysis can obviously be extended. Beyond a more detailed statistical validation, we could consider higher-dimensional food webs [22] which leads to a problem in R N .