Limit Cycle Oscillations, response time and the time-dependent solution to the Lotka-Volterra Predator-Prey model

Limit-cycle oscillations occur in many areas of Biology and Physics [1, 2]. One popular model to describe these phenomena is the Lotka-Volterra model [1, 3] and its extensions. Lotka [4] and Volterra [5] derived it independently to describe the nonlinear interaction of predator and prey populations. Therefore, it is now widely known as the ‘Predator-Prey model’. In Physics, many nonlinear interactions can be described with this model. In Plasma Physics, on which we focus here, it describes the interaction between axisymmetric self-generated flows (zonal flows), which acts as the predator, and the microturbulence - the prey - that drives them [6]. In the plasma turbulence context, this model is a building-block for more extended models of the Low to High confinement (L-H) transition in fusion plasmas such as the Kim-Diamond model [7]. See e.g. Ref. [8] for a review. Refs. [9] and [10] analyzed such extended predator-prey models - which are not integrable - using bifurcation theory. The limit-cycles associated to the Lotka-Volterra model have one particular feature that is interesting: They exhibit multiple time-scales. Hence, such models may be useful to understand certain types of relaxation-oscillations, and intermittent transport in turbulent plasmas. The well-known Van-der-Pol oscillator [11] and the closely-related Rayleigh oscillator [12] are also examples of a system showing relaxation-oscillations. During the L-H transition, a dithering phase (also called Intermittent phase or I-phase) is often observed experimentally when the heating power is slowly ramped up [13, 14]. It can be understood as limit-cycle oscillations (LCO) between turbulence energy and zonal flow energy, as observed in gyrokinetic simulations [15, 16]. It is this feature that we focus on in this work. More precisely, we focus on the response-time between these two quantities. Refs. [17, 18, 19] showed - based on a different model - that this response-time is a key quantity to understand nonlinear interactions. In the simplified framework of Drift-wave turbulence, there are several candidates to explain the interaction between zonal density corrugations and the turbulence. One of the authors (M. Leconte) proposed a model based on the nonlinear modulation of the transport crossphase between density and potential perturbations [20, 21]. Another model, based on stochastic noise due to turbulence was proposed in Ref. [22]. The main results of this article are:
(i) a new analytical solution opens the possibility of directly fitting experimental data of LCO to extract its key predator-prey features, with only two fitting parameters, c.f. Eqs. (14,15,16,17).
(ii) the response time between turbulence energy and zonal energy increases as marginality is approached, in collision-less trapped electron mode turbulence simulations. A similar trend is observed in the Lotka-Volterra model.
The rest of this article is organized as follows: In Section 2, we describe the Lotka-Volterra Predator-Prey model, and we derive its time-dependent solution analytically. In Section 3, we apply the newly-found solutions to understand the predator-prey dynamics between zonal density perturbations and the turbulent particle flux observed in global gyrokinetic simulations of collision-less trapped-electron mode turbulence (CTEM) [23, 24, 25, 26]. Finally in Section 4, we discuss the results and give a conclusion.

We consider the following Lotka-Volterra Predator-Prey model:

Here, $x$ denotes the prey population, and $y$ denotes the predator population, and $\dot{u}=du/dt$, with $u=x,y$, denotes the time derivative. The parameter $\gamma$ denotes the prey growth-rate (birth-rate) in the absence of predator, and $\mu$ is the predator damping rate (death-rate) in the absence of prey. The coefficients $\alpha_{1},\alpha_{2}$ are positive constants. In applications to plasma turbulence, the predator is usually taken as zonal flow energy, and the prey as turbulence energy [6, 7, 8]. Here, we take the predator as zonal density corrugations driven by nonlinear modulation of the transport crossphase [20] and the prey as turbulence energy. In section 3, this model will be applied to gyrokinetic simulations of CTEM turbulence. There are several candidate mechanisms for the nonlinear generation of zonal density by the turbulence [20, 22] based on fluid models, although a specific application to CTEM has not been proposed yet. For this reason, we treat the Lotka-Volterra as a phenomenological model here. However, deriving a predator-prey like reduced-model directly from the bounce-averaged gyrokinetic equation would be an important task. One could use a similar method as in Ref. [19], where a derivation of the traffic-jam model from the nonlinear gyrokinetic equation is sketched.

Let us first consider the solution for the predator population $Y(t)$. Starting from the predator evolution Eq. (6). we write it in the form:

Here, the function is given by $g(X,Y)=-\frac{1}{(1-X)Y}$. Now, we use Eq. (8) to express the prey $X$ in terms of the predator $Y$ and total energy $E$:

where $W_{j}$ denotes the Lambert W function [29], and $j=0,-1$ is the associated branch.

Integrating both sides of Eq. (11) yields:

where we call $G_{LV}^{j}$ the ‘second Lotka-Volterra’ integral, defined as:

where $W_{j}=W_{0}$ or $W_{j}=W_{-1}$, depending on the branch of the Lambert $W$ function considered. The integrand of the Lotka-Volterra integral (14) for the branches j=0,-1 is shown v.s. $Y$ for different values of the parameter $\delta$, for an energy of $E=3-\ln 2$ [Fig. 3]. This value of energy corresponds to the initial conditions $(X_{0},Y_{0})=(2,1)$. Note the vertical asymptotes corresponding to the minima and maxima of the predator population $Y$. The integration constant $T_{Y_{\rm min}}$ is given by expression (B2) in Appendix.

Next, we invert Eq. (13) to obtain the time-dependent solution $Y(t,E)=\Lambda_{j}^{\rm preda}(t,E)$ for the predator population, where the function $\Lambda_{j}^{\rm preda}(t,E)$ is given by:

where the superscript ^-1 denotes the function inverse (not to be confused with the index $j=0,-1$ of the branch), and we use the shortcut notation $G_{LV}=G_{LV}^{j}$. Care must be taken when inverting Eq. (13), because the real-valued Lambert W function has two branches: $W_{0}(x)$ and $W_{-1}(x)$. As a second step, one applies a similar procedure for the prey population X, $dt=f(X,Y)dX$, with $f(X,Y)=\frac{1}{(\frac{1}{\delta}-Y)X}$, to obtain the time-dependent solution $X(t,E)=\Lambda_{j}^{\rm prey}(t,E)$ for the prey population:

where the superscript ^-1 denotes the function inverse, and we define the first Lotka-Volterra integral $F_{LV}$ as:

We call the analytical solutions $\Lambda_{j}^{\rm prey}(t,E)$, and $\Lambda_{j}^{\rm preda}(t,E)$ given by expressions (15) and (16) the ‘Lotka-Volterra functions’. The result is shown for a value $\delta=0.5$ of the parameter and an energy of $E\simeq 2$ [Fig. 4].

The collisionless trapped-electron mode (CTEM) is an instability due to the electron toroidal precession-drift resonance - a process similar to inverse Landau damping - in the low-collisionality regime [30]. Its sources of energy are the electron temperature gradient and the density gradient. CTEM instability is in the ion-scale range, with a typical poloidal wave-number $k_{\theta}\rho_{i}\sim 1$, where $\rho_{i}=\sqrt{m_{i}T_{i}}/eB$ is the ion-gyroradius, $T_{i}$ is the ion temperature and other notations are standard. From collisionless trapped-electron mode (CTEM) gyrokinetic simulations, zonal density perturbations $\delta n_{ez}$ and electron particle flux $\Gamma_{e}$ are obtained. Figure 5 shows the zonal density corrugations v.s. radius and time. One clearly observes the radially oscillating zonal pattern known as ‘staircase’ [17, 18, 19, 23, 24, 31, 32, 33, 34, 35, 36, 37, 38]. The time-trace of electron particle flux $\Gamma_{e}$ and zonal density energy $n_{ez}^{2}$ are shown [Fig 6a]. These quantities are spectral averages, e.g. $\Gamma_{e}(t)=\sum_{k_{r}\rho_{i}\in[0.4,1]}|\Gamma_{e}(t,k_{r})|$, around the radial wavenumber $k_{r}\rho_{i}=0.78$, which is a characteristic scale of the zonal staircase pattern. The simulations are performed with the gyrokinetic code gKPSP [39] which solves the nonlinear gyrokinetic equations for ions [40] and bounce-averaged kinetics for trapped electrons [41]. In the simulations, the equilibrium gradients are $R/L_{n}=2.2$, $R/L_{Ti}=2.2$, and $R/L_{Te}=4.0-12.0$ where $n,T_{e}$ and $T_{i}$ denote the electron density, electron temperature and ion temperature, and e.g. $L_{n}=-[\frac{1}{n}\frac{dn}{dr}]^{-1}$. At mid-radius, the inverse aspect ratio is $r/R=0.18$, safety factor $q=1.4$, magnetic shear $\hat{s}=q^{\prime}r/q=0.78$, $T_{i}=T_{e}$. Hydrogen is the main ion $m_{i}/m_{e}=1836$ and plasma elongation $\kappa=2$. The turbulent transport of this plasma is dominated by CTEM [24, 25]. A limit-cycle type of dynamics between $\Gamma_{e}$ and $n_{ez}^{2}$ is clearly observed in dynamical phase-space [Fig 6b], although its amplitude decreases with time, probably due to additional turbulent dissipation. This is probably the reason for the spiraling in Fig 6b.The Lotka-Volterra model does not take into account this additional dissipation, but the overall dynamics is similar to the model.

The response time is usually defined as the time-lag between maxima or minima of two signals. Here, for clarity, we define it as the time-lag between the first maximum of the two signals, as indicated by black arrows in Fig 6a . The response-time $\tau$ is shown v.s. the distance to threshold $R/L_{T}-R/L_{T}^{c}$ [26], related to the linear growth-rate $\gamma_{\rm TEM}$ via $\gamma_{\rm TEM}\sim\sqrt{\Big{|}\frac{R}{L_{T}}-\frac{R}{L_{T}^{c}}\Big{|}}$ [Fig. 7a]. Here, the critical gradient is $R/L_{T}^{c}=3$.

In the Lotka-Volterra model, the response time is given analytically by:

where $G_{LV}^{-1}(Y,E)$ is the second Lotka-Volterra integral (branch $j=-1$) given by expression (14).

For the Lotka-Volterra PP model, the response time is shown v.s. the parameter $1/\delta$ (i.e. $\gamma/\mu$), at fixed energy $E=3-\ln 2$ [Fig.7b]. The analytical result is shown (full symbols) and is compared to the numerical result (open symbols) obtained by solving numerically the Lotka-Volterra Eqs. (5,6). The agreement is reasonable.

One observes that both response times Fig. 7a and 7b increase with increasing drive. Conversely, both response time decrease with decreasing drive. This is consistent with a critical exponent behavior: the response time increases as marginal stability is approached, i.e. as $R/L_{T}\to R/L_{T}^{c}$ ($L_{T}=L_{Te}$) in the gyrokinetic simulation and $1/\delta\to 0$ in the Lotka-Volterra model. Writing $\tau=|R/L_{T}-R/L_{T}^{c}|^{-\alpha}$ and $\tau_{LV}=(\frac{1}{\delta})^{-\beta}$, the following scalings or obtained: $\alpha=0.67$ and $\beta=0.9$. This trend is also predicted in the traffic-jam model for avalanches of Ref. [19] for a 1D model of zonal ion temperature corrugations. Based on this model, the scaling of a characteristic response time was given as: $\frac{1}{2\tau}\sqrt{\frac{c_{0}^{2}\tau}{\chi_{2}}}\simeq\gamma_{\rm max}$, where $\gamma_{\rm max}\sim V_{E\times B}^{\prime}$, with $V_{E\times B}^{\prime}$ the zonal flow shearing rate, $c_{0}$ the initial avalanche speed and $\chi_{2}$ the heat diffusivity. Assuming that $\chi_{2}$ scales like $\chi_{2}\sim|R/L_{T}-R/L_{T}^{c}|^{\eta}$, as for stiff transport with $\eta>0$ a scaling exponent, this predicts that the response-time of Ref. [19] scales like $\tau\sim|R/L_{T}-R/L_{T}^{c}|^{-\eta}$.

Let us first discuss the analytical solution of the Lotka-Volterra model. To our knowledge, this is the first time that such a closed-form exact solution of the model was obtained. This analytical solution (14,15,16,17) may possibly be used to fit experimental data of limit cycle oscillations in fusion devices. It provides a simpler alternative compared to the method used in Ref. [16], which extracted directly predator-prey coefficients from gyrokinetic simulation data by solving the Lotka-Volterra system numerically. More precisely, our method extracts two fitting-parameters: the ratio of linear predator damping-rate to prey growth-rate $\delta=\mu/\gamma$, and the ‘energy’ $E$ associated to the limit-cycle. This can be used by the experimental community to know if a signal is ‘predator-prey’, and to extract its key parameters. We sketch briefly how one could proceed: From the experimental prey & predator signals $X(t)$ and $Y(t)$, one could minimize the squared-difference between $X(t)$ and $\Lambda_{j}^{\rm prey}(t)$, and between $Y(t)$ and $\Lambda_{j}^{\rm preda}(t)$ - i.e. a least-square fit - with respect to the fitting-parameters $\delta$ and $E$.
There are some limitations to our model. First, this model differs from the drift-wave zonal-flow model of Ref. [6] in that the self-damping term is neglected in the turbulence evolution equation. Although this term is important to obtain steady saturated states in the long time limit, and associated bifurcation between states, this parameter has no effect on the transient limit-cycle regime - apart from displacing its center - on which we focus here. Second, we acknowledge that we did not formally invert Eqs. (14) and (17) by solving e.g. Eq. (13) for $Y$. Instead, we inverted them graphically by plotting $X$ v.s. $F_{LV}$ and $Y$ v.s. $G_{LV}$. The inversion could be easily accomplished by writing the Lotka-Volterra integrals (14,17) as a Taylor series, and then inverting this series, using the Taylor series inverse-function formula. This is left for future work. Third, we verified our analytical solution only for one value of the energy $E$ and associated initial conditions $(X_{0},Y_{0})=(2,1)$. It would be useful to verify it for different initial conditions (different energies), but this is beyond the scope of this article. Let us now discuss the comparison of the response-time between the Lotka-Volterra model and gyrokinetic CTEM simulations. In Ref. [44], a stochastic version of the Kim-Diamond model was studied. The temporal cross-correlation between zonal flow energy and turbulence energy was computed and the associated response-time was obtained (Fig 1c in this reference). Ref. [19] predicted -in the framework of the ‘traffic-jam’ model- that the response-time scales like $\tau\sim|R/L_{T}-R/L_{T}^{c}|^{-\eta}$, i.e. the response time increases as marginal stability is approached. We showed that a similar trend, i.e. $\tau$ increasing when the prey drive $~{}1/\delta$ decreases is obtained when applying the Lotka-Volterra predator-prey model to the interaction between zonal density corrugations and turbulence. This is also consistent with Ref. [34], where the probability of finding large scale (and thus slowly evolving) zonal structures is maximal near-marginality (cf. Fig 4 in the latter Reference). More connection between the zonal staircase generation and predator-prey modeling would be interesting for future work.
In conclusion, we derived analytically the solutions to the well-known Lotka-Volterra Predator-Prey model. We applied the newly-found solutions to calculate analytically the Predator-Prey response time, and we compared its scaling with that of gyrokinetic simulations of CTEM turbulence. As the Lotka-Volterra model is one of the simplest model to describe self-organized systems, we believe that having an analytical insight into the dynamics of this model is the first step to allow a better understanding of more complicated models.

The authors would like to thank P.H. Diamond, Y. Kosuga, Raghvendra Singh, X. Garbet, T. Kobayashi, Min-Jun Choi, Jae-Min Kwon and Sumin Yi for helpful discussions. We also thank the anonymous Referees for valuable comments. This work was supported by R & D Program through Korean Institute of Fusion Energy (KFE) funded by the Ministry of Science and ICT of the Republic of Korea (KFE-EN2141-7). The data that support the findings of this study are available from the corresponding author upon reasonable request.

First note that harmonic functions $\sin(\theta)$ and $\cos(\theta)$ can be defined as the inverse of integrals, via:

respectively. Jacobi elliptic functions arise from a generalization of these formulas, when the trigonometric circle is generalized to an ellipse (although historically, Jacobi discovered these functions using a different approach).

One can make an analogy between the solutions $\Lambda_{j}^{\rm preda}(t,E)$, $\Lambda_{j}^{\rm prey}(t,E)$, with $j=0,-1$ and the Jacobi elliptic functions sn(t,k) and cn(t,k), respectively. Note that Jacobi elliptic functions appear naturally in describing the non-circular bounce/transit periodic motion of guiding centers in tokamak geometry (e.g. the ‘banana’ orbits, etc …) and its applications to residual zonal flow problem [43]. One way to define the Jacobi elliptic functions is to view them as extended trigonometric functions [42]. Consider an ellipse for which:

Defining normalized variables $X=x/b$ and $Y=y/b$, one obtains, after some algebra:

Here $\kappa=a/b=E/b^{2}$ is the elongation, where $E=ab$ is the energy, proportional to the surface area of the ellipse, i.e. the action invariant, $I=\oint pdq=\pi ab$, where $q=X$ and p=Y are the usual position and momentum conjugate variables. Hence elongation can also be viewed as ‘normalized energy’. One clearly sees the analogy between Eq. (A3) - which represents limit cycles of normalized energy $E/b^{2}=\kappa$ - and Eq. (8) for the Lotka-Volterra system. Both represent energy conservation. The only difference is the topology of the limit-cycle, and hence the form of the energy integral. Note also that in both cases, the shape of the limit-cycle depends on the value of the energy.

Dividing by $\kappa$ and defining the elliptic modulus $k=\sqrt{1-\frac{1}{\kappa^{2}}}$, this can be written:

It is well-known that Eq. (A4) describes phase-space contours associated to the solutions $X/\kappa={\rm cn}(t,k)$ and $Y={\rm sn}(t,k)$, where ${\rm cn}(t,k)$ and ${\rm sn}(t,k)$ are Jacobi elliptic functions, and $\kappa=1/\sqrt{1-k^{2}}$.

The quantity $T_{Y_{\rm min}}$ takes the form [27]:

where $T_{Y_{\rm lower}}=\int_{Y_{\rm min}}^{Y_{\rm max}}\frac{dY}{Y\Big{[}1+W_{-1}(-Y^{-1/\delta})e^{Y-E}\Big{]}}$ is the lower semi-period for the predator population. Similarly, the quantity $T_{X_{\rm min}}$ reads:

where $T_{X_{\rm upper}}=-\int_{X_{\rm min}}^{X_{\rm max}}\frac{dX}{X/\delta\Big{[}1+W_{-1}(-\delta X^{-\delta})e^{\delta(X-E)}\Big{]}}$ is the upper semi-period for the prey population. Note that similar formulas are given in Ref. [27] for the quantities $T_{Y_{\rm upper}}$ and $T_{X_{\rm lower}}$. The exact period $T_{\rm LCO}$ of the orbit (or limit cycle) is then given by the sum of the two semi-periods: