Memory and mutualism in species sustainability: A time-fractional Lotka-Volterra model with harvesting.

We first present a predator-prey model for two species and then extend the model to three species where the two predator species engage in mutualistic predation. Constant effort harvesting and the impact of by-catch issue are also incorporated. Necessary sufficient conditions for the existence and stability of positive equilibrium points are examined. It is shown that harvesting is sustainable, and the memory concept of the fractional derivative damps out oscillations in the population numbers so that the system as a whole settles on an equilibrium quicker than it would with integer time derivatives. Finally, some possible physical explanations are given for the obtained results. It is shown that the stability requires the memory concept in the model.

Introduction

Ecological models are of great importance for environmental decision making because they provide stakeholders with a conceptual framework and a “laboratory” for studying the consequences of alternative policies and management scenarios [1]. One common method that scholars use to improve our understanding of environmental phenomena is the Lotka-Volterra (or predator-prey) model – an important and popular prototype model appearing in various fields of applied mathematics – due to its descriptive power, tractability and diverse applications [2], [3], [4], [5], [6], [7], [8], [9]. Therefore, many efforts have been made so far to propose more realistic models incorporating mutualism [10], [11], parasitism [12], [13], and the impact of harvesting [14], [15], [16], [17], [18], [19], [20].

Some scholars extend models formulated with fractional derivatives through applying fractional calculus (fc). In 2013, by subjecting the predators to harvesting, a modified fractional version of predator-prey model with a type II functional response was proposed [15]. Later, the fractionalised model was developed further by adding an economic interest equation to the model [21]. More recently, in addition to using the Caputo derivative to study the complex behaviour of the phenomena, scholars applied the newly formulated fractional version of the Adams-Bashforth method in their research [22], [23], [24], [25], [26], [27]. As another, for example, Kolade proposed a three-component time-fractional system representing the interaction among prey, intermediate-predators, and predators and examined the model behaviour under Caputo or the Atangana– Baleanu fractional derivatives [28].

In this study, we consider a new type of fractional model for three species (two predators and one prey) with type II mutualistic predation. An example of the type of interaction the model idealises is the mutualistic predation of spotted dolphins and yellowfin tuna upon schools of lanternfish [29]. Since the understanding of harvesting is an important issue for fishering, we assumed that all species have a market value, and also considered the impact of by-catch for dolphins when fishing for tuna. Making a comparison between the integer derivative model and fractional model under different scenarios, we examine the concept of “memory” on our model. We show that stability is more robust when the species exhibit “memory”. We examine the impact of harvesting on the system with and without the memory concept as well.

Preliminaries on fractional calculus and main models

Fractional calculus is a powerful tool which has been employed in different fields of science to model complex systems with non-local behaviour and long-term memory [30], [31], [32], [33], [34], [35]. This approach can lead to models capturing more of the phenomena under scrutiny while still keeping the model parameters to a minimum.

The adjective “fractional” in fc is a historical remnant and this calculus is a generalisation from integer order derivatives to arbitrary real-valued order, not merely rational order. Generally speaking, it can be formulated as follows. where and are the fractional derivative and fractional integral respectively [36]. Various definitions of fractional calculus have been proposed. All definitions coincide when the order is integer, however, this need not be the case for non-integer order. Therefore, different physical interpretations, known as the memory concept, are proposed in the fractional case [30].

We will use the Caputo definition (for other formula, see [30], [36], [37]). where , and is a weight function whose task is the storage of the system memory over time [30].

Larger values of α increase the weight on the integrand Df close to t, emphasizing the memory of nearby values of f. When alpha is close to zero, models formulated with the Caputo derivative (), will retain close to complete memory of the past starting at time a. It is expected that when the system maintains a near total memory of its past then the system resists changing over time. Based on this interpretation, therefore, we are expecting that oscillations in the population numbers of a species damp out in a system with the fractional time derivative of the order less than unity. We have examined this issue in the forthcoming sections and shown that this is the case with fractional prey-predator models.

Single predator model

In this section, we considered a fractional model of predator and prey. Using constant harvest quota and , we assume that either both species have market value or that one species is caught as by-catch. By doing so, we tried to answer such questions as how does the harvesting of the species affect the natural equilibrium of the ecology? How heavily can a species be harvested and still be sustainable [20], [38]? where X and Y are the population densities of prey and predator respectively, T is time, r is the prey growth rate, K is the environmental carrying capacity for the prey, a is the feeding rate of predators, σ is the predator growth saturation factor and k is the predator death rate. All parameters are positive reals.

After substituting the rescalings , and in to eq. (3), we arrived at the following dimensionless system where , , , and .

Two predators model

In this section, the system (3) is extended to model interactions amongst three species: one prey and two predators. The predators are not treated as isolated hunters. Rather, we consider the predator species to be cooperative. The model of the three species with type II mutualism [10], [39] functional response for the predators is where X is the prey population density (lanternfish), Y and Z are population densities of distinct predators (tuna fish and dolphin respectively), T is time, r is the prey growth rate, K is the environmental carrying capacity for the prey, a, b, and ξ are the feeding rate of predators, and are the predator growth saturation factor and and are the predator death rate. As before, all parameters are positive reals and the terms for are the harvesting functions.

The number of parameters in system (5) can be reduced by considering the following transformations thus we arrive the following dimensionless system: where , , , , , , , and .

Model analysis

Existence and uniqueness of the non-negative solution

For and , the Mittag-Leffler function is defined by the following series: where Γ is the gamma function.

Assume that and where . Then for any initial conditions , all the solutions of model (4) are non-negative and unique for all .

(Proof by contradiction): let , then . Also assume If , then after taking Laplacian transform (Table C1& C2 of [30]) from the first equation of system (4), we have This leads to a contradiction. In the similar way, when , we can obtain the contradiction. Therefore, and as a result, will be positive for all .

We now show that the system (4) satisfies the locally Lipschitz condition [40] needed to establish the existence and uniqueness of solutions to system (4).

Consider a mapping with norm such that We show that where . where and . That is, with initial condition , an unique solution exists for system (4). □

Assume that and where . Then for any initial conditions , all the solutions of model (6) are non-negative and unique for all .

(Proof by contradiction): let , then . Also assume If , then from the first equation of system (6), we have This leads to a contradiction. In the similar way, for and , we can obtain the contradiction. Therefore, and as a result, will be positive for all .

As for the existence and uniqueness, after following the same procedure as Theorem 3.1 and considering a mapping function (Eq. (9)) as , we have found that under the following , equation (10) holds for all and . Hence, with initial condition , an unique solution exists for system (6). □

Stability analysis

In this section, we have examined the stability of equilibrium solutions of the proposed models. All the calculated stability conditions are summarised into Table 1 and Table 2.

Table 1

Asymptotic stability analysis of system (4). The notation ✓ indicates that the relevant condition is satisfied.

System(4)
Ei: Equilibrium point	Stability conditions	Fig. 2
E_1: Extinction point	ρ<ε_1α≤1	…
E_2: Predator-free point	ψ<(1+ε_2)[ρ(ρ−ε_1)+ϕ]α≤1	…
E_3: Coexistence points	ω<1−ε_1ρb^2−4c<0α<2π|arg⁡(−b2±i4c−b^22)|	0.13<1−0.4→✓(0.53)^2−3.7<0→✓α∈(0,1.178)→✓

Table 2

Asymptotic stability analysis of system (6). The notation ✓ indicates that the relevant condition is satisfied.

System(6)
Ei: Equilibrium point	Stability conditions	Fig. 3
E_1: Extinction point(0,0,0)	ρ<ε_1α≤1	…
E_2: Predator-free point(x^⋆,0,0)	ε_1<ρε_2>ψ(ρ−ε_1)/[ρ+ϕ(ρ−ε_1)]ε_3>βη(ρ−ε_1)/[ρ+ϕ_1(ρ−ε_1)]α≤1	…
E_3: Partial coexistence point(x^⋆,y^⋆,0)	ε_1<ρ(1−ω)ε_2<ψ/ϕε_3>βω(η+γ)/[1+ϕ_1ω]ρ<2(1+ϕω)γωψα∈(0,α_1)	…
E_4: Partial coexistence point(x^⋆,0,z^⋆)	ε_1<ρ(1−ω)ε_2>ψω(γ+1)/[η(1+ϕω)]ε_3<ηβ/ϕ_1ρ<2(1+ϕ_1ω)γωηβα∈(0,α_1)	…
E_5: Coexistence point(x^⋆,y^⋆,z^⋆)	ε_1<ρ+η−γγ_1ε_2>ψω_1/[1+ϕω_1]ε_3>ηβω_1/[1+ϕ_1ω_1]ψε32+ηβε22>ψε32γ_1+βε22γ−ψωγγ_1βε_1ε_2α∈(0,1)	0.12<0.62→✓0.43>0.3812→✓0.06>0.0563→✓0.016>−879.9→✓0.96∈(0,1)→✓

To obtain the stability of dynamic models with integer order, one very common way is to take advantage of the Routh-Hurwitz (rh) conditions. where J and are the Jacobian matrix and equilibrium point of the dynamic system respectively.

If , then will be a saddle point (regardless of the sign of ),

If , then will be asymptotically stable when ,

If , then is unstable when .

In addition to rh conditions, in fractional calculus, one needs to examine another condition (Lemma 3.1) [41].

[37, page 158]

Consider function with . An equilibrium point () is asymptotically stable if all eigenvalues (, ) of the associated Jacobian matrix in satisfy the following condition (Fig. 1)

Figure 1

Stability region of the fractional-order system.

What stands out from the Fig. 1 is that the stability domain of the fractional order system is larger than the corresponding domain for integer systems when . Therefore it is expected that a non-integer derivative will increase the stability of the system.

For , implies that all the roots are real and implies that there is only one real root and one complex and its complex conjugate. Also

Consider the following polynomial the discriminant of the polynomial is defined by where is the derivative of f and where and is an determinant.

[41]

Consider characteristic equation (12), then the conditions that satisfies (11) are as follows.

i) For , the condition is .

ii) For , the conditions are either rh conditions or iii) For , if , , then (11) is satisfied for all .

iv) If , , then (11) is satisfied for all .

v) For general n, , is a necessary condition for (11).

Single predator model

The system (4) has three stationary points: extinction of both species; predator only extinction; and predator/prey co-existence, where

According to the stationary points, non-zero harvesting clearly shifts the prey-only equilibrium, , to lower densities of prey. For the co-existence point, , the prey population density is unaffected by harvesting of the prey itself but is shifted to higher densities as the rate of predator harvesting is increased. The harvesting of prey is effectively competition for the predators reducing their equilibrium population in concert with any direct harvesting. Also where is the Jacobian matrix of the system (4) in the stationary points .

From (13), it can easily be deduced that in extinction point (). According to Lemma 3.1 therefore, is asymptotically stable for if and only if . That is, if the rate of harvesting of the prey outstrips its growth rate, the system is driven to extinction.

Also, it can be concluded that the linearised system has two eigenvalues of at , and thus according to Lemma 3.1, is asymptotically stable for all only when which corresponds to the coexistence point being biologically irrelevant. For fixed point , we have Thus, the trace and determinant of the Jacobian evaluated at Eq. (14) are as follows: According to rh, the sufficient condition for asymptotically stability of the system (4) is when .

To satisfy the Lemma 3.1 as the necessary condition, we must have (Fig. 2)

Figure 2

The damping of population oscillations produced by system (4) with initial densities (0.2,0.25) and ψ = 19, ρ = 1, ψ = 15, ϕ = 2, ε1 = 0.4 and ε2 = 1. (a and b) Numerical values of the lanternfish and tuna versus time respectively. (c): Combined predator-prey population densities with respect to time and (d): Phase portrait of the system.

Two-predator model

There are five stationary points for model (6): where is total population extinction, is prey only, and are equilibria of partial co-existence (the prey with one of the predators), and is the co-existence of all three species with .

, ;

s.t ;

s.t ; and

s.t .

Linearising the system (6) about the stationary points , we can determine each point's linear stability by considering the eigenvalues of the resulting Jacobian matrix. From (16) and Lemma 3.1, it can easily be deduced that is asymptotically stable for if only if . At also, we have three following eigenvalues thus according to Lemma 3.1, is asymptotically stable for all only when

As for , we have the following characteristic equation. Therefore, according to the Lemma 3.1 and rh, the necessary and sufficient conditions for asymptotically stability of the system (6) at for all is where , and .

Regarding , according to (16), we will also have the following characteristic equation. Using the rh condition (ii) and the Lemma 3.1, E4 is asymptotically stable for all only if where , , and .

Also, the following characteristic equation is given for . thus according to Lemma 3.2, under the following conditions are asymptotically stable for all . where , , , and (Fig. 3d).

Figure 3

The damping of population oscillations produced by system (6) with integer (a) and fractional order (b). (c and d): Phase portrait of the integer and fractional system respectively. The initial densities (0.2,0.3,0.3), (ρ,ψ,β)=(0.61,1,7), (ϕ,ϕ1)=(1.4,0.02), (ε1,ε2,ε3)=(0.12,0.43,0.06) and η = 0.01.

Physical interpretation of the models

In this section, we have numerically [42] simulated both systems (4) and (6) under the stability conditions found in the previous section. Also, we presented some possible physical explanations for the obtained results.

From the simulations presented in Figure 2, Figure 3, the effect of reducing the order of the time derivative can be seen. As the fractional order α is decreased, the system (with Caputo derivative) stabilises faster. That is the higher “memory” the system has of past states, the greater the damping of oscillations in system dynamics. The simulations demonstrate that, even with quite moderate reductions in α, the amplitude of population density oscillations is strongly retarded.

From the numerical simulations we see that the harvesting appears to enhance the stability of the system. When no harvesting is conducted, oscillations in the population numbers are immense (Fig. 4a). Once at least one of the harvesting parameters is non-zero these oscillations are damped out (Fig. 4b), aligning with the earlier results presented in the literature [15], [20].

Figure 4

The impact of harvesting upon the system (4) with fractional order (0 < α ≤ 1), ρ = 1,ψ = 15,ϕ = 2 and (x(0),y(0))=(0.2,0.1).

We emphasize the novel relationship between memory concept and harvesting. The numerical integration of system (4) (blue lines) indicates that when interaction among species is immense, such that the population numbers fluctuate wildly, harvesting could be a solution to stabilise the system. For some parameter values, harvesting works less well as a stabilising mechanism and could destabilise the coexistence fixed point of the system (Fig. 4c). In other words, the modelling assumption of the “largeness” of the population, allowing the use of differential equations, may very well breakdown. In these situations, harvesting reduces the range of values through which the population numbers fluctuate (thus reducing the time for the system to settle at equilibrium) but also shifts the range to lower values. The stability of the fixed point as per the stability condition is likely incorrect because the population returns to healthy levels after having dropped to number better interpreted as having gone extinct (Fig. 4c). To smooth out this inconsistency, our findings suggest that another parameter, defined as species' memory, should be incorporated in the model (Fig. 4d). A possible explanation for this could that when species are heavily harvested, they rely much more on their memory to prevent the species extinction. For example, if they are heavily harvested during the spring season, they produce more offspring intrinsically during that period to make up their lack of population. In other words, they might have a variable growth rate to compensate for the different physical circumstances from one breeding season to another; they might find a safer place to reproduce so that they reduce the death-rate of infants, or even migrate to another place where they can decrease mortality rates of offspring through less exposure to predators. Therefore, we believe that such scenarios could be explained by the memory concept.

Discussion

The mutualistic interactions amongst species is of great importance in the field of conservation ecology, so gaining an understanding of such interactions can make a noteworthy contribution to species maintenance. With this in mind, we introduced a modified Lotka-Volterra model to study interaction amongst three species with mutualistic predation. The motivation is based on the observed feeding behaviour of spotted dolphins and yellowfin tuna upon schools of lanternfish. Along with obtaining the stability conditions for the model, we also investigated the impact of “memory effects” on the species interaction via fractional calculus. Our analysis reveals that the fractional system dampens out induced oscillating inherent in predator-prey models and reaches the local stable point sooner than the integer model does. In other words, stability is more robust when the species exhibit “memory”.

Moreover, assuming that either both species have market value or that one species is caught as by-catch, we have investigated the effect of constant rate harvesting within the proposed models. We have also discussed the local stability behaviour of all the equilibrium states of the system. The output of all models shows that stability and extinction of the ecosystem are affected by economic interest/harvesting. In fact, applying constant harvest at a rate below the threshold of the stability condition causes the system to stabilise faster. In other words, exploitation of a species can be regarded as beneficial to the ecosystem as a whole as the system will reach steady-state sooner. Further, the greater the collaboration among species (large mutualism coefficients), the more harvesting that can occur before putting the populations at risk of extinction (over-fished).

An assumption of a differential equation based population models is that the quantities under consideration are sufficiently large such that the addition or loss of an individual to that population may be considered an infinitesimal change. That is, the population may take on a continuum of values rather than being strictly discrete. Further, it is well known that many such models can exhibit, for certain parameter values and initial conditions, large fluctuations — particularly those with Hopf bifurcations — whereby the assumption of the “largeness” of the population may very well breakdown Fig. 5. In practical terms the population has either already gone extinct, even though the model predicts a return to healthy levels, or its reduced to a level vulnerable to external perturbations such that it future viability should be discussed in probabilistic terms. These issues imply the model should incorporate stochasticity at low population levels. Alternatively, assuming a memory like behaviour can be attributed to the population, a non-integer order time derivative arguably extends the validity of the model by preventing the wild swings in population numbers that represent ‘essentially extinct’ to ‘thriving’ (Fig. 5a).

Figure 5

Phase plane of system (4) with integer order derivative (5a) and fractional order derivative (5b). All initial densities are (x,y)=(0.2,0.15) with parameter values (ρ,ψ,ϕ,ε1,ε2)=(1,50,0.02,0,0). For these parameter values criterion (15) holds.

It has been argued that ecological systems, while being stable for moderate numbers of species interactions or moderate strength of the inter-species connections, will become unstable to small perturbations once some threshold value for interaction number or interaction strength is breached [43], [44]. In 2012, similar findings were presented, claiming that any feasible co-existence of the system species will be unstable when the pair-wise competitive interactions are sufficiently strong [45]. Further, Goh presented findings that the continuum of globally stable Lotka-Volterra models of mutualism among three or more species is smaller than the continuum of globally stable Lotka-Volterra models of competition among the same number of species [46], [46]. Thus, if diversity has an adverse effect on stability in competitive system it has even more so in mutualistic ones.

More recently, however, Butler and O'Dwyer have to some extent overturned this understanding through proposing a consumer-producer-resource model. They demonstrated, for a model of N bacteria species consuming N abiotic resources that the stability is guaranteed for all feasible equilibria. For an extension of the model where the bacteria also produce some or all of the resources to mutual benefit of all consumers, stability of all feasible solutions can be guaranteed provided that mutualistic interactions are symmetric [47], [47].

Our system, on the other hand, is one with moderate to strong interactions between mutualistic predators and a prey (biotic resource) species. Our formulation allows weak to strong memory effects via the fractional derivative order. The three-species model suggests the stability at the co-existence is not guaranteed by symmetric mutualism. Rather, the model needs to have a memory to satisfy the stability condition. As future work we will investigate how stability of feasible solutions in the model changes as the number of predators and prey increase.

Declarations

Author contribution statement

M.M. Amirian: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

I.N. Towers, Z. Jovanoski, A.J. Irwin: Analyzed and interpreted the data; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.