Model

The ecological model describes interaction dynamics between a predator and an age-structured prey composed of juvenile and adult developmental stages. The original model equations where the prey population is regulated through maturation and the predator feeds only on adult prey are as follows [2]: (1) where J, A, and P are state variables describing the size of juvenile, adult and prey populations respectively. Here bA is the linear fecundity function with b as reproduction rate of adults. The nonlinear function ϕJ/(1 + dJ²) describes maturation of juvenile population with maturation rate of ϕ/(1 + dJ²) that goes to zero as juvenile population approaches infinity. ϕ is the maximum maturation rate of juvenile at low density, and d is the strength of exploitation competition among juvenile population. Parameter n is the attack rate of predator on adult and parameter c is the predator conversion efficiency. Death rates are μ_J, μ_A, μ_P for juvenile, adult and predator populations respectively.

As suggested in [2] a scaled form of this model is used here to reduce the number of parameters: (2) where all rate parameters should now be interpreted as fractions of ϕ.

Steady states of the model

The steady states of the model are determined by setting all time evolutions to zero and simultaneously solving the resulting steady state equations. Deriving the nontrivial solutions of steady state equations is not always possible in higher dimensional dynamical systems. New methods have been introduced recently to approximate fixed points of higher order stochastic systems based on corresponding deterministic mean-field approximation [3]. Due to mathematical simplicity of P equation in our model, we are able to directly solve the corresponding steady state equations. From the third equation in Eq (2) we have P(Ac − μ_P) = 0 with P = 0 and A = μ_P/c both satisfying the equation. Here we ignore the P = 0 trivial solution since this corresponds to complete elimination of predator population. Instead we select the nontrivial solution A = μ_P/c which describes the equilibrium size of adult population. Using this equilibrium value in the second and third equations in Eq (2) we obtain: (3)

These equations can be rewritten as (4) (5) (6)

Eq (4) is a cubic polynomial in J, so its equilibrium values J° can be determined using the roots function in Matlab. Substitution in Eq (5) gives the corresponding value for P° and Eq (6) gives A°, the equilibrium values for predator and adult populations as a function of μ_P (death rate of predator population). Fig (1a)–(1c) shows the resulting steady state diagrams for each population. The juvenile and predator populations display multi-root regions while the adult density follows a linear trend over the full range of μ_P values. Two critical values of μ_P = B1, B2 on the border of multi-root region are indicated in Fig (1a) and (1b). Bifurcation points mark the locations where system dynamics undergoes a qualitative change. Eigenvalue analysis determines the stability properties and the types of bifurcation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Steady-state diagram and corresponding eigenvalues of the fisheries model as a function of predator death rate μ_P.

(a-c) steady states showing a S-bend shape for J and P variables with color-coded lower (blue), middle (green) and upper (red) branches. Population A shows a linear trend with μ_p. (d, e) Real and imaginary parts of linearized model. Saddle-node bifurcation points are marked B1 and B2 and indicated by vertical dashed grey lines. Model parameters are c = 1, b = 1, μ_J = 0.05, μ_A = 0.1. Eigenvalues are not relevant for since predator population goes negative.

https://doi.org/10.1371/journal.pone.0163003.g001

Linear stability analysis

We perform a full eigenvalue analysis of linearized system following standard methods described in [4]. We linearize the model by writing: (7) where is a small temporal perturbation, and Z ∈ {J, A, P}; Z° is the equilibrium point. It is informative to identify if , and perturbations grow or decay to locate unstable and stable points. Replacing model variables with their perturbed forms, and using Taylor series expansions, Eq (2) can be written (8) where the higher order terms are neglected since is small. Noting that f_1,2,3(J°, A°, P°) = 0, we obtain, (9) to describe the time evolution of the perturbations . Expressing in matrix form, (10) where (11) is the Jacobian matrix evaluated at the equilibrium point (J°, A°, P°) and f_i, i ∈ {1, 2, 3} are system functions defined in Eq 2. The exponential time-course for small perturbations away from steady state can be predicted from the eigenvalues of [5]. The equilibrium is stable when all eigenvalues have negative real parts, otherwise the equilibrium is unstable. Close to equilibrium the system dynamics is determined by the dominant eigenvalue, i.e., the eigenvalue whose real part is least negative. See Fig 1(d).

The lower branch of J steady-state diagram (top branch of P) is a stable focus-node since the system has one real negative eigenvalue and a pair of complex conjugate eigenvalues with negative real part. Note that the focus is dominant for μ_P ≤ B₃ but the node takes over for μ_P ≥ B₃. The top branch of J (bottom branch of P) forms stable focus-nodes for where the node aspect is dominant. In contrast, all points on the midbranch are unstable saddle-focus equilibria forming a separatrix between upper and lower branches.

We identify three critical μ_P values:

1. μ_B1 = 0.43525: a stable focus-node collides with an unstable saddle-focus. This is simply a saddle-node bifurcation at the left-hand turning point.
2. μ_B2 = 0.55281: saddle-node bifurcation at right-hand turning point.
3. μ_B3 = 0.53017: dominance is swapped between focus and node components. Stability is not modified, so a bifurcation does not happen here.

Having identified bifurcation points (B₁ and B₂) and special point (B₃), we run a series of stochastic numerical simulations, then validate against theoretical predictions.

Noise-induced fluctuations prior to bifurcation

We repeat the numerical experiments described in [1] but, rather than displaying coefficient of variation (standard deviation divided by mean), in Fig 2 we plot the unscaled fluctuation variance as a function of μ_p. The coefficient of variation is not ideal for tracking fluctuations since changes in the mean can confound changes in fluctuation amplitude; the latter is where the dynamic information resides. We first analyze the case when the noise is added to J population only as shown in Fig 2(a). From eigenvalue analysis we estimate the mortality rate at the right-hand saddle-node point as μ_P = μ_B2 ≈ 0.55280625013933332. We use geometrically spaced μ_P values to closely approach this catastrophic collapse point while noise is added to the death rate of J population in the same way as described by Boerlijst et al: (12) where a₁ is a coefficient to ensure that the noise amplitude is small, and ξ₁(t) represents a zero-mean, Gaussian-distributed white-noise process. We first calculate the steady state, and then use this as the initial value in an Euler-based numerical update of the differential equations using time steps of Δt = 0.01. The variance is extracted for each μ_P value and trends are plotted for 0.4 ≤ μ_p ≤ μ_B2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Numerically obtained fluctuation variance of model variables prior to saddle-node bifurcation.

(a–d) For each value of μ_P, starting with μ_P = 0.4, the model is simulated for 60,000 time units with resolution of 0.01 time units, from which the population variances are computed. μ_P is incremented geometrically towards the catastrophic collapse at μ_P = μ_B2. Death rates are perturbed every time unit using white noise with standard deviation of σ_noise = 2 × 10⁻⁶. (a) Noise added to the juvenile population (b) Noise added to the adult population (c). Independent noise added to all three populations. (d) Identical, fully correlated, noise added to all three populations. (e) An extra experiment where noise is added to P only population.

https://doi.org/10.1371/journal.pone.0163003.g002

Results are shown in Fig 2(a). Despite the general similarity with those reported in [1], there are subtle but important differences. We see that the fluctuation variance increases prior to catastrophe in all three populations, Fig 2(a). This increase is significant for juvenile and predator populations and very weak for adult population. We found that capturing the growth of fluctuation variance for adult and predator populations is numerically challenging, and requires:

1. Precise identification of control parameter value (μ_P) at phase transition.
2. Performing numerical experiments in close vicinity to catastrophe. We chose geometrically spaced μ_p values to closely approach the B₂ saddle-node, (13)
3. Application of considerably smaller amplitude white noise when performing the experiments close to phase transition. We used very small amplitude noise with standard deviation σ = 6.317 × 10⁻⁷ compared to σ = 0.005 used in [1], otherwise the noise-induced fluctuations are strong enough to cause a jump transition in the system, preventing us from dwelling close to bifurcation.

Fig 2 shows that fluctuation variance grows significantly prior to SN point for J and P population, and has a tiny or zero growth for A population regardless of the way noise is added to the system. The exception is when noise is added to the A population only (Fig 2(b)): for this case, none of the populations exhibit fluctuation growth. The results of an extra experiment where the noise is added to P only population is displayed in Fig 2(e).

To verify the accuracy of our simulations, we transform the stochastic model equations into Ornstein-Uhlenbeck (OU) form. Then we compare the theoretically-predicted OU statistics against the corresponding values extracted from numerical simulations.

Ornstein-Uhlenbeck (OU) analysis

We include additive white noise in all three equations to transform model (2) to a stochastic form suitable for OU analysis, (14) where a_1,2,3 are scaling constants that ensure that the fluctuations are small, and ξ_1,2,3 are independent zero-mean, Gaussian-distributed white-noise sources [6], (15) where δ_mn is the dimensionless Kronecker delta, δ(⋅) is the dirac delta with a dimensionless total area under its curve equal to one, and the 〈⋅⋅⋅〉 represents the ensemble average over time. Noise samples (implemented by Matlab’s randn function) are scaled as (16) where describes a zero-mean, unit-variance Gaussian random number generator; the scaling by ensures that ξ(t) tends to an infinite-variance white noise in the limit Δt → 0 [6, 7, 8]. (Notice that the noise is not added to the mortality rates of populations as implemented by in Boerlijst et al [1] and shown in Eq (12)). We arrange the model equations in matrix form: (17) with as the Jacobian matrix defined in Eq (11). We recast into standard OU form, [9, 10]: (18) where is the drift matrix and D is a diagonal 3 × 3 diffusion matrix (19)

Based on well documented OU statistics [9, 10], we can immediately write down theoretical expressions for covariance matrix of the model. The stationary covariance matrix Σ of an OU process is defined by (20) where Σ is the 3 × 3 stationary covariance matrix (21) from which one can extract the theoretical variance predictions. Noting that Eq (20) is the continuous Lyapunov equation, one can derive Σ as: (22) where vec() represents the vectorised form of a matrix and ⊗ is the Kronecker product [11]. The vectorization operation is defined as the stacking of the columns of a matrix into a vector. This method is a generalization of the approach suggested by Gardiner which is limited to 2 × 2 matrices [10].

We use Eq (22) to compute the theoretical covariance matrix. Eq (21) shows that the diagonal elements of this matrix represent the individual fluctuation variances of each of the three populations. We perform similar experiments as shown in Fig 2(c) where independent noises are added to all three populations and plot the numerically obtained variances as solid lines in Fig 3(a)–3(c); theoretical predictions are superimposed as dashed lines. Good agreement between theory and experiment is evident in all three populations. Variance of J fluctuations increases towards the SN point as shown in Fig 3(a). Variances of A and P populations show a decreasing behaviour, but P variance recovers and increases significantly close to catastrophe as shown in (b), (c).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Fluctuation variance prior to saddle-node bifurcation with independent noises added to all three populations.

Experimental and theoretical variances are plotted for all three populations while approaching catastrophe. A fixed step Euler method with Δt = 0.01 is used for numerical simulations each for 6000 time units. Starting with μ_P = 0.4, it is incremented towards the saddle-node point at μ_B2. The distance from bifurcation point is geometrically reduced enabling more experiments close to bifurcation. Three independent white noise sources are added to each population as described in Eqs (14)–(16) with standard deviation of σ_noise = 2 × 10⁻⁶.

https://doi.org/10.1371/journal.pone.0163003.g003

One can perform similar analysis when the noise is added to only J or P population. The qualitative form of the results is not strongly sensitive to the number of independent noise sources, e.g., compare the three-noise case of Fig 3 with S1 and S2 Figs.

These results show that although a saddle-node bifurcation occurs in this three-dimensional dynamical system, not all system variables display early warning indicators. Considering the steady-state diagrams of Fig 1, one notices that J and P populations switch to a new state at bifurcation points while A variable does not. The linear steady state diagram for A population does not support bifurcation, so it is not expected to show critical fluctuations.

On the other hand one might expect to capture reliable early warning signs on both J and P fluctuations due to their similar multi-root steady-state diagrams, but the results do not support this expectation. While the J population demonstrates consistent growth towards the SN point, the P population does not. Instead, we only observe significant fluctuation growth very close to SN.

Boerjlist et al. [1] used the concept of dominant eigenvector to explain the strong growth in J fluctuations and apparent invisibility of P fluctuations. We tested this idea by extracting the dominant eigenvalue and corresponding eigenvector components for a broad range of μ_P values while approaching the B₂ saddle-node point. Fig 4 displays the dominant eigenvalue (previously displayed as part of Fig 1(d)) and the corresponding eigenvector decomposed into its J, P and A components. Close to state transition, the dominant eigenvector has a significant J component, a small P component, and an almost zero A component.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Tracking of dominant eigenvalue and its decomposed eigenvector towards saddle-node bifurcation point.

(a) The real part of dominant eigenvalue and (b) the corresponding decomposed eigenvector as a function of predator mortality rate μ_P.

https://doi.org/10.1371/journal.pone.0163003.g004

Although behaviour of the eigenvalue components matches the variance trends—and Boerlijst et al [1] reported similar findings—one may still ask the important question: Why do the P fluctuations not increase all the way towards SN point? Since SN bifurcation can occur in one-dimensional systems [12], perhaps it is the case that J is the only variable in control of SN bifurcation in this model. We will demonstrate that observability analysis can provide insights for this argument.

Observability analysis

Inferring whole system dynamics based on measuring a single time series of the system is a challenging problem. Observability analysis reveals which variable (or combination of variables) best represents internal dynamics. See S1 Appendix for a general description of observability analysis. Here we compute the observability coefficients for the system when monitored through J, P, and A populations to quantify the system observability from individual variables. The observability coefficient is derived from the observability matrix constructed from either Lie derivatives or a coordinate transformation that maps from the system to a nominated single variable [13]. Both procedures are mathematically the same, resulting in the same observability matrix. Here we construct the observability matrices from Lie derivatives. The same observability matrices can be obtained using coordinate transformation as shown in S2 Appendix. The Lie derivatives of the system are with f_i, i ∈ {1, 2, 3} previously defined in Eq 2. The next step is to construct matrix by taking derivatives of and with respect to J, A and P variables

We are interested in observability analysis via each individual system variable, so define measurement vectors, corresponding to J, A and P variables respectively. Then the state-dependent observability matrix of each variable is (23)

These matrices can be constructed for every point at space space, but we choose to compute them only at stable steady states in order to confine our study to ecologically significant stable states. A system is defined to be observable via variable x if the observability matrix is full rank, otherwise not observable [14]. Equivalently if is nonsingular (does not have a zero eigenvalue), then the system is observable, otherwise not. In practice, an observable system may gradually lose observability due to a varying system parameter. Aguirre et al. [13] developed a method to quantify degree of observability using observability coefficient. Following their method to compute the observability coefficient through individual system variables one should first extract the state dependent observability measure as (24) where is the maximum eigenvalue of matrix estimated at point (J(t), A(t), P(t)) (likewise for λ_min). According to this definition 0 ≤ δ_x ≤ 1 and the lower bound is reached when the system is unobservable for variable x. The next step is to compute the observability coefficient (which is a constant) as the time average of δ_x(J(t), A(t), P(t)) along the trajectory (J(t), A(t), P(t)), t ∈ {t_initial, t_final}. Confining our analysis to stable steady states, we extract locally without a need for time averaging, giving a local measure of stability: (25) where observability matrices are computed at steady states. We use this coefficient to determine the observability of the fisheries model while gradually approaching the B₂ saddle-node point by increasing the mortality rate of the prey population.

We extract the observability coefficients using observability matrices obtained by Lie derivatives. The results are shown in Fig 5. and cross at μ_P ≈ 0.475 showing that for lower values of the predator mortality rate it is best to measure P whereas for higher rates of μ_P it is best to measure J if one is to monitor the internal dynamics. Focusing at the close vicinity of SN point, while the system is significantly observable from J population (), the observability coefficients of P and A populations are and . These values indicate that the system is mainly observable from J prior to SN point. The coefficient trends over the μ_P range show a decrease in system observability via and but an increase for . Notice that recovers slightly just prior to SN point as shown in the inset. These trends are in agreement with variance of fluctuations as displayed in Figs (2a), (2c) and (2e) and 3. The small observability coefficients of P and A variables explain the challenging task of capturing signs of slowing down on these variables in numerical simulations. This implies that signs of slowing down would be very hard to observe on predator or adult populations in field or even in controlled laboratory setups.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Observability coefficients of fisheries model.

Observability coefficient of system variables , and as a function of μ_P calculated using the definition in Eq (24). Lie derivatives are used to construct the observability matrix from which the observability coefficients are calculated. See text for details. A geometrically spaced vector of predator mortality rate is used to cover the range 0.4 ≤ μ_P ≤ μ_B2.

https://doi.org/10.1371/journal.pone.0163003.g005