Recurrence quantification analysis of a three level trophic chain model.

Recurrence plots and recurrence quantification analysis based measures - recurrence rate, determinism, divergence, entropy, laminarity and trapping time, are used to detect transitions between periodic and chaotic states and also the laminar states of the complex three species food web model with improved growth rate function and predatory ability. The study of dynamical equation of such a complex biological system with transients in terms of these measures results in localization of bifurcation behavior of the system.

Introduction

In recent years, nonlinear system theory has been extensively used to understand the underlying complex processes involved in the behavior of biological, ecological, economic, physical, psychological, physiological dynamical systems [1], [2], [3], [4], [5], [6], [7], [8], [9]. Usually, the complexity of observed or simulated dynamical behavior is analyzed using a time series of a single dynamical variable of such a system. The parameters e.g., the spectrum of Lyapunov exponents, the correlation dimension, the entropy that characterizes the dynamics are computed on the assumption that the used steady state time series of the system is rather long, stationary and noise free. The use of recurrence plot (RP) and the recurrence quantification analysis (RQA) is free from such constraints and provide an alternative and meaningful investigation of various transitions of states of a dynamical system, in particular, the biological and ecological sciences, by incorporating the transients [10], [11], [12], [13], [14], [15], [16], [17], [18], [19].

In the present study, we apply the RP and RQA analysis [10], [18] to investigate the complexities of the dynamical system of three species food chain model [20]. This model incorporates an improved form of growth rate function based on resource utilization efficiency and also include the effect of reduced predatory ability caused by the presence of the top predator in the system. For specific values of the intrinsic growth rate factor (), considered as a control parameter of the system, the RPs based on the top predator time series, are shown to provide visually appealing structures that distinguish the regular and chaotic dynamics. In order to further unravel the complexity of the dynamical processes characterized by non-stationary drifts and state transitions [13], [14], [18], we simulate a dynamic time series for this model by varying the control parameter at each iteration. Recurrence quantification analysis acts as an empirical test to distinguish between random and relatively short time series non-stationary drifts and state transitions [33] with delays and embedding derived from quantification of recurrence plots [32]. The application of RQA to such a data results in measures e.g., recurrence rate (RR), determinism (DET), divergence (DIV), entropy (ENT), laminarity (LAM) and trapping time (TT), which are found to rightly ascertain the transition points. Unlike the other measures of the food chain, the measures LAM and TT refer to the vertical structures that appears in its RP and the observed distinct peak or maxima of these measures at specific values of the parameter identify the chaos to chaos transitions in the model.

The paper is organized under sections. In Section 2, we briefly outline the mathematical model and assumptions of considered food chain model along with some basic results on the coexistence equilibria and their stability. In Section 3, the bifurcation plot is obtained to identify different complex transition regions of the system on varying the intrinsic growth rate parameter . Section 3.2 describes the basic features of RPs and how its pattern is used to characterize the complexity of the food chain. In Section 3.3, the various RQA measures are introduced and the application of such measures to a dynamic time series of the food chain system reported in present work is shown to quantify not only periodic to chaotic/chaotic to periodic but also transitions type chaos to chaos.

Model

Mathematical model of three species food chain

A general class of three species model [21] consists of the predator population Z at the highest level preying on the population Y at the middle level that subsequently preys upon the species X at the bottom level. The model of such a chain could be represented as and the dynamics of such a system may be described by following set of three coupled differential equation, where represent a logistic growth function for the species X with and denoting the intrinsic growth rate and carrying capacity respectively. In recent years several forms of interaction between the species has been considered for the three level trophic chain model [22], [23]. In this work, the interaction between any two species is assumed to be of Holling II-type, such that the interaction terms becomes , , where and are constant parameters for the system. The factors and respectively represents the conversion rate of the prey to predator for the species Y and Z and the parameters and refers respectively to the dying rate for the species Y and Z.

This model can be extended to include a new growth rate term of the form instead of a logistic growth term for the species X [20]. It is to be noted that the new form of growth factor incorporate the resource utilization efficiency which is high in case the factor and low in case . The latter situation restricts the population growth quickly [24]. Further a function of the form may be introduced to account for reduction in the predatory ability of the species Y in the presence of the predator Z. As a result the nonlinear dynamical equations for the three species food chain could be written as, Using the following transformations: eq. (2) can be reduced to the where we have defined For co-existence of the three species, the equilibrium points, of the system depend on the values of various parameters e.g., , , , and and could be determined as, where , . For positive equilibrium points i.e., , and , we have chosen biologically relevant parameters values as in [20] i.e., , , , , , , and .

The stability of equilibrium point may be studied using the Jacobian matrix where , and .

The characteristic equation obtained for the Jacobian suggests asymptotic stability of the equilibrium point S as discussed earlier in [20]. The stability analysis of other trivial equilibrium points and of the system eq. (3) similarly suggest these points to be saddle points. We may further note that the absence of the top predator () results in a subsystem comprising intermediate predator y which may survive on its prey x. In such a situation the positive equilibrium point could be written as , , provided and . The conditions for asymptotic stability of such equilibrium points are described in [20]. The complexity of a nonlinear dynamical system gets revealed through its bifurcation diagram and the following section describe it when the growth rate parameter is varied.

Analysis

Bifurcation diagram

The temporal behavior of the population x, y and z is obtained, as a result of numerical integration of the system eq. (3) with initial condition as in [20], using Runge-Kutta method with step size of 0.01 for different values of the intrinsic growth parameter . Values of other biological parameter, considered above, are retained. The time series of the top predator z, for and 6.0, are shown in Fig. 1(a, b).

Figure 1

Time series of the top predator z in three species food chain model with k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1, b5 = 0.3 and (a) r0 = 2.3, (b) r0 = 6.0.

Power spectrum analysis may be used as a preliminary tool to characterize a time series to be regular or chaotic/noisy. For the intrinsic growth rate parameter , the time series (Fig. 1a) of the predator z reveals periodic dynamics as characterized by the occurrence of sharp lines in its Fourier spectrum (Fig. 2a). However, the temporal behavior of top predator population time series (Fig. 1b), for , is observed to be very complex. Its power spectrum shows the presence of multi-frequency splitting of lines (Fig. 2b), a characteristic of chaotic/noisy dynamics.

Figure 2

Power spectra of the top predator z time series in three species food chain model with k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1, b5 = 0.3 and (a) r0 = 2.3, (b) r0 = 6.0.

For different values of the control parameter of the three level trophic chain, the simulated time series of the various species in the steady state may be searched for their respective successive maxima to construct the bifurcation diagram for each species. As in [20], the details of the bifurcation diagram (Fig. 3) for the species z are observed to be similar to one for the species x and y. We therefore consider subsequently the bifurcation diagram (Fig. 3) of the system for the top predator z only with control parameter in the range . Fig. 4(a-c) show the details of the Fig. 3 in the range , and respectively. For instance in Fig. 4(a), we observe the dynamics of three level trophic chain to remain doubly periodic until around when the system makes a transition to period-4 state. With increase in , we observe the system to exhibit period-8 behavior (), and chaos in the region (Fig. 4b). Interestingly, the region in this bifurcation diagram shows complex dynamics interspersed with fine regular windows. In Fig. 4(c), we observe the chaotic dynamics of the system changing to period-3 behavior at and subsequently a wide chaotic band showing the coexistence of expanding and collapsing dynamics appears. In addition a few windows where transition from chaos → regular and regular → chaos also appears, say around and beyond .

Figure 3

Bifurcation diagram of the three species food chain model [eq. (3)] with k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Figure 4

Details of bifurcation diagram of Fig. 3 for (a) 2 ≤ r0 ≤ 4, (b) 4 ≤ r0 ≤ 6, (c) 6 ≤ r0 ≤ 8, with k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Recurrence plot analysis (RPA)

Complex systems such as the three species food web system [eq. (3)] with given values of associated biological parameters could be categorized by using the nonlinear dynamical tools e.g., Lyapunov Characteristic exponent (LCE), correlation dimension, Kolmogorov-Sinai (KS) entropy, etc. However, such analysis assumes the systems temporal behavior to be stationary and noise free. The recurrent behavior of natural systems involves periodicities and also irregular cyclicities. It is in this regard the idea of recurrence plot (RP) was introduced [11] to visualize the behavior of a trajectory in phase space of a complex nonlinear dynamical system. The method was subsequently extended to provide various measures to quantify the RP's [11], [12], [13], [14], [15] to unravel the complexities of the observed dynamical system.

Consider a discrete time series P with N data points such that Here the values to corresponds to the values of the time series P sampled at times . In case the true dimension m of the dynamical system is known, a sequence of following vectors may be constructed from this time series as, where τ refers to the time lag or delay parameter and m the embedding dimension of the phase space.

Recurrence plots (RP) could be formed by considering the distances between the m dimensional reconstructed points. In fact RP is a symmetrical array, wherein a ‘dot’ is marked at a point , if a point is close to another point . Mathematically, RP is expressed as a matrix [11], [12], [13], [14], [15], Here is a Heavyside function, ϵ is a pre-defined threshold distance between neighboring points and is an Euclidean norm [15], [25], [26].

RPs display characteristic patterns which may be used to characterize typical dynamical behavior. The RP method can be used to characterize both stationary and non-stationary data of any dynamical system. It is to be noted that for deterministic time series, RP exhibit both horizontal and vertical lines since a dynamical system has the property of being recurrent. For the time series of a stochastic dynamical system, such lines in RP are of very small size and appear by chance. The distribution of recurrent points in RP, in such a case, is observed to be almost homogeneous. However, in case of periodic system, the RP comprises parallel longer diagonal lines.

Here, in order to construct the RPs of the three species food chain model [eq. (3)], we consider the population time series of the species z corresponding to a given value of intrinsic growth parameter . The delay parameter τ is found by observing the first minima of the mutual information curve [27]. Similarly the embedding dimension, m, is determined using the method of false nearest neighbor [28], [35]. Knowing the delay parameter τ and the embedding dimension m, eq. (7) allows one to reconstruct the state space to estimate characteristic properties of the dynamical system [18]. The reconstructed state space further facilitate the formation of RP using eq. (8). The phase space attractor along with RP constructed using the z-time series with intrinsic growth rate and the threshold parameter [29], [30] are shown in Fig. 5(a, b). The RP for clearly shows parallel diagonal lines characteristics of a periodic signal comprising two almost irrationally related frequencies [31]. Similarly the 3-D attractor and its RP in Fig. 6(a, b) demonstrate the transition of the dynamical behavior of the system to period-4 stable limit cycle when the intrinsic growth parameter is increased to 3.9. The phase space chaotic attractor and the RP, shown in Fig. 7(a, b) display the complex chaotic dynamics of the system for . It is also observed that the regular dynamics set in the system as is further increased to, say where the period-3 stable limit cycle develops (cf. Fig. 8(a, b)). However, we find the dynamics once again becomes chaotic at and the RP clearly illustrate the unstable periodic orbits of the chaotic attractor (Fig. 9(a, b)) [15]. It is interesting to observe that complexity of the dynamics of the three species food chain as revealed by RPs for a specific value of confirms the results obtained in section 3 using the bifurcation diagram. The other measures characterizing the dynamical complexities using the RPs are discussed in section 3.3.

Figure 5

(a) Phase space period-2 attractor in three species food chain model, (b) Recurrence plot of period-2 behavior of the z-time series with r0 = 2.3, k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Figure 6

(a) Phase space period-4 attractor in three species food chain model, (b) Recurrence plot of period-4 behavior of the z-time series with r0 = 3.9, k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Figure 7

(a) Phase space chaotic attractor in three species food chain model, (b) Recurrence plot of chaotic behavior of the z-time series with r0 = 6.2, k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Figure 8

(a) Phase space period-3 attractor in three species food chain model, (b) Recurrence plot period-3 behavior of the z-time series with r0 = 6.45, k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Figure 9

(a) Phase space chaotic attractor in three species food chain model, (b) Recurrence plot of chaotic behavior of the z-time series with r0 = 7.0, k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Recurrence quantification analysis (RQA)

The analysis of observed horizontal and vertical lines in a RP of several known complex dynamical systems by [15], [18], [34], [36] led to development of several complexity measures characterizing the complex behavior of nonlinear dynamical systems. In the following, we briefly outline some of these complexity measures which further characterize the state transition observed in the dynamic time series of the top predator z (Fig. 10).

Figure 10

Dynamic time series of the three-species food web model [eq. (3)], with k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Recurrence rate percentage (RR)

RR measure defined as, provides the density of recurrence points in a RP. The measure, in fact provides the probability of finding similar states of a dynamical system and therefore could be used as an indicator to observe changes in system dynamics.

Determinism percentage (DET)

DET measure may be defined as, where corresponds to number of diagonal lines of length l.

DET is the ratio of recurrence points forming diagonal structures of length to all recurrent points in a RP. It is a measure of predictability or determinism and therefore for a deterministic dynamical system DET is high (RP consists of longer diagonal structures and a very few minute diagonals or none) while for chaotic or stochastic dynamical system, it is very low (RP shows predominance of short diagonal structures and hardly any longer diagonal lines).

Maximal length () and divergence (DIV)

and DIV measures may be defined as, where equals the total number of diagonal lines.

The measure DIV is related to largest positive Lyapunov exponent [13], [14], [16], [17], [26], [34]. In RP of periodic system, the diagonal lines are longer whereas they are very short in RP of chaotic dynamical system.

As a result, for regular dynamics DIV values tend to be smaller while for the system exhibiting irregular/chaotic behavior, DIV values are higher.

Shanon entropy (ENT)

ENT defines the Shanon entropy as, where is the probability of finding a diagonal line of length l. It corresponds to the complexity of deterministic structures in RP. For a system with higher complexity, ENT value is high.

Laminarity (LAM)

LAM is defined as: where corresponds to the frequency distribution of vertical lines of length ν in RP [15], [36].

The definition of LAM, above is analogous to the determinism, and corresponds to ratio of recurrent vertical structures to the entire set of recurrence points.

Trapping time (TT)

TT measure is defined as: TT is termed as trapping time and refers to the information in RP about the extent and length of the vertical structures. It provides the estimates of the average time the system spends in a given state of the dynamical system [36].

Now, we attempt to characterize the dynamical transitions that takes place as result of varying the intrinsic growth rate parameter of the system in terms of the measures described above. We note here that the bifurcation characteristics of the system (Figs. 3 and 4) are investigated presuming that the system is in steady state regime after the removal of transients. Complexity of the biological system such as, eq. (3), is often characterized by non-stationary drifts and changing states. RQA enables us to investigate the bifurcation of the model by retaining the transients while dynamically varying the control parameter .

The RQA analysis is carried on the dynamic time series of the species z obtained as a result of numerical integration of the system [eq. (3)] using fourth order Runge-Kutta method with time step of 0.20 [13], [14], [16]. At every integration step the intrinsic growth parameter get incremented by 0.000175 and the range covered is while holding other parameters as constant. Fig. 10 illustrates the resulting time series of length for the species z and incorporates transient effects at each steps. Using the delay parameter , embedding dimension , the RQA analysis is repeatedly applied to each window of data size 1000 with 20% overlap.

Fig. 11 shows the variation of the RQA measures RR, DET and DIV with increase in intrinsic growth parameter . The RR exhibit clearly the deterministic structure for . For the system to make a transition from periodicity to chaos, we observe that the RR, DET decreases while DIV measure increases. In the present model, unlike the measure RR, we note a drop in the value of the measure DET at which is indicative of system dynamics changing from periodic-2 to period-4 state. Again a significant decrease in the measure RR in the region is suggestive of the transition from periodic to chaotic state (cf. Fig. 4b). This is further supported by the peak observed in the plot of the measure DIV at (Fig. 11). It is to be noted that in [20], the largest Lyapunov exponent () has been computed for different values of for the present three species food chain model and it is shown that is less than zero in the region . Therefore regularity of the dynamics in this region is clearly identified in the plot of the DIV measure. Further, the observed sharp peaks in DIV plot at various values of localizes the more chaotic state of the system.

Figure 11

RQA statistics RR, DET and DIV computed for the dynamic time series of the top predator, z with k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Fig. 12 illustrates the variation of other RQA statistics viz. ENT, LAM and TT, with the parameter . As in the case of the measures RR, DET and DIV the computation of the measure ENT is based on the distribution of the diagonal lines of the RPs. Various low values observed in the entropy of the diagonal line distribution corresponds to lower dynamical complexity in the model. For instance, it is readily observed that the sharp dip in the ENT measure at is identified with transition from chaotic to regular period-3 behavior of the system. Various windows exhibiting higher values of the measure ENT clearly demonstrates the complex chaotic dynamics that may prevail in the present system at different values of the parameter . In addition to the bifurcation diagram, the observed variation in ENT measure further identifies the domain of variation in wherein the system behavior is more complex than the other. It is to be noted that measures such as LAM and TT, however, are based on the vertical structures in RP and may be used to identify chaos → chaos transitions (cf. [15]). The occurrence of vertical/horizontal line structure in RP corresponds to a time length during which the state of the system changes slowly or remain almost unchanged. As a result, the state becomes laminar state or intermittent as the system is trapped for a finite time. The sharp decrease in LAM and TT at and 5.25 exhibit the periodic-chaotic and chaotic-periodic transitions. However, the transition between various states of the system in the vicinity of the intrinsic growth parameter in the bifurcation diagram (Fig. 3) appears to indicate chaos to chaos transitions. Consequently vertical/horizontal lines tends to appear frequently at such chaos to chaos transitions than in other chaotic regions [cf. [15], [36]]. Therefore the observed maxima of LAM and TT in the region may point to the onset of laminar behavior of the system. The appearance of the vertical structures in RP of the model at (cf. Fig. 9b) further confirms the onset of laminar behavior as revealed in the plot of LAM and TT.

Figure 12

RQA statistics ENT, LAM and TT computed for the dynamic time series of the top predator, z with k0 = 3, k1 = 4.5, b1 = 4.0, b2 = 1, b3 = 2, b4 = 1 and b5 = 0.3.

Conclusion

The RP and RQA framework is used to investigate the complex dynamics of the food chain, incorporating the effects of resource utilization and the reduced predatory ability in the presence of top predator z, on varying the intrinsic growth rate parameter . The distinctive patterns of RP, constructed for the food chain at different values of the control parameter , are shown to characterize the periodic, quasi-periodic and chaotic regimes of the model. RPs also enable localization of the bifurcation sequence of the model. Various values of at which the RP of the system exhibit period-doubling transitions, are found to be in agreement with the values obtained from its bifurcation diagram. The RQA statistics - RR, DET, DIV, ENT, LAM and TT are further used to characterize the complexity of the model by incorporating the dynamical transient and the non-stationary behavior in its time evolution. In this work, RQA is therefore applied to the time series of the top predator z of the model generated by incrementing the value of at each iterations. The dynamical transition from periodic to chaotic state is identified with decrease in value of RR and DET while the value of DIV increases. Various peaks in DIV plot are found to represent the chaotic state of the system. The ENT measure shows several windows of higher complexity in the system and peaks/dips in the plot precisely identifies the dynamical regimes of the system like the bifurcation diagram. The LAM and TT measures, unlike the foregoing measures, are based on the vertical structures of the RP. These measures not only identifies the periodic to chaotic or chaotic to periodic transition but also the chaos to chaos transition in the food chain model, say at .

Declarations

Author contribution statement

Rashmi Bhardwaj: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data.

Saureesh Das: Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Funding statement

Authors are thankful to GGSIP University for providing financial support.

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.