Aphids, Ants and Ladybirds: a mathematical model predicting their population dynamics

TL;DR

This paper introduces a new mathematical model using differential equations and finite difference schemes to predict the population dynamics of aphids, ants, and ladybirds, highlighting the impact of mutualism and competition.

Contribution

It develops a novel mathematical model combining existing ecological equations with a discretization scheme to analyze population stability and interactions among aphids, ants, and ladybirds.

Findings

Mutualism increases population peaks of aphids and ants.

Aphids and ladybirds populations collapse asymptotically.

Ant populations stabilize at a constant value.

Abstract

The interaction between aphids, ants and ladybirds has been investigated from an ecological point of view since many decades, while there are no attempts to describe it from a mathematical point of view. This paper introduces a new mathematical model to describe the within-season population dynamics in an ecological patch of a system composed by aphids, ants and ladybirds, through a set of four differential equations. The proposed model is based on the Kindlmann and Dixon set of differential equations, focused on the prediction of the aphids-ladybirds population densities, that share a prey-predator relationship. The population of ants, in mutualistic relationship with aphids and in interspecific competition with ladybirds, is described according to the Holland and De Angelis mathematical model, in which the authors faced the problem of mutualistic interactions in general terms. The set of differential equations proposed here is discretized by means the Nonstandard Finite Difference scheme, successfully applied by Gabbriellini to the mutualistic model. The constructed finite-difference scheme is positivity-preserving and characterized by four nonhyperbolic steady-states, as highlighted by the phase-space and time-series analyses. Particular attention is dedicated to the steady-state most interesting from an ecological point of view, whose asymptotic stability is demonstrated via the Centre Manifold Theory. The model allows to numerically confirm that mutualistic relationship effectively influences the population dynamic, by increasing the peaks of the aphids and ants population densities. Nonetheless, it is showed that the asymptotical populations of aphids and ladybirds collapse for any initial condition, unlike that of ants that, after the peak, settle on a constant asymptotic value.

Figures

Tables

Table 1. Table 1: Steady-states and eigenvalues of the ( 20 ).

Steady-state	(h,x,y,z)	Eigenvalues
${\Gamma }^0$	(0,0,0,0)	(1, 0.999996, 0.999934, 1.000030)
${\Gamma }^1$	(0.599525,0,0,0)	(1, 0.999935 $\pm$ 0.000023 $i$, 1.000030)
${\Gamma }^2$	(1.796894, 0, 30, 0)	(1, 0.999897 $\pm$ 0.000022 $i$, 0.999970)
${\Gamma }^3$	(0, 0, 30, 0)	(1, 1.005961, 0.991550, 0.997031)

Equations

$$\frac{dh(t)}{dt}=ax(t)$$

$$\frac{dx(t)}{dt}=(r_1-h(t))x(t)-\frac{\nu px(t)y(t)}{b+px(t)+y(t)}+{\alpha }_{12}\frac{x(t)y(t)}{{ϵ}_1+y(t)}$$

$$\frac{dy(t)}{dt}=r_{2}y(t)+{\alpha }_{21}\frac{x(t)y(t)}{{ϵ}_2+x(t)}-d_2{y}^2(t)$$

$$\frac{dz(t)}{dt}=-\frac{\nu z^2(t)}{b+px(t)+z(t)}-ky(t)z(t),$$

$$\frac{d\xi (t)}{dt}=f(\xi (t),K),$$

$$\frac{{\xi }_{n+1}-{\xi }_n}{\Delta t}=F({\xi }_n,K),$$

$$\Delta t⟶\varphi (\Delta t)+O(\Delta t^2),$$

$$\varphi (\Delta t)=\frac{1-e^{-q\Delta t}}{q}.$$

$$j_{ij}=\frac{\partial {\Phi }_i}{\partial {\xi }_j},$$

$$q\geq\max_{\Omega}\Big{\{}\frac{\lambda^{2}}{2|\Re(\lambda)|}\Big{\}}\,\,\,\rm if\,\,\,\Re(\lambda)\neq 0\,\,for\,\,\,\lambda\in\Omega.$$

$$f(\xi ,K)=(-1{)}^p{\xi }^{m}g(\xi ,K),$$

$$\displaystyle F(\xi_{n},K)=\frac{(-1)^{p}}{2}\bigg{[}\big{(}1-(-1)^{p+1}\big{)}\xi_{n}^{m}+\big{(}1-(-1)^{p}\big{)}\xi_{n}^{m-1}\xi_{n+1}\bigg{]}g_{n}(\xi_{n},K).$$

$$\displaystyle\frac{\xi_{n+1}-\xi_{n}}{\phi(\Delta t)}=\frac{(-1)^{p}}{2}\bigg{[}\big{(}1-(-1)^{p+1}\big{)}\xi_{n}^{m}+\big{(}1-(-1)^{p}\big{)}\xi_{n}^{m-1}\xi_{n+1}\bigg{]}g_{n}(\xi_{n},K),$$

$$\displaystyle\xi_{n+1}=\frac{\frac{(-1)^{p}}{2}g_{n}(\xi_{n},K)\phi(\Delta t)\big{(}1-(-1)^{p+1}\big{)}\xi_{n}^{m}+\xi_{n}}{1-\frac{(-1)^{p}}{2}\xi_{n}^{m-1}\big{(}1-(-1)^{p}\big{)}g_{n}(\xi_{n},K)\phi(\Delta t)}.$$

$${\xi }_{n+1}(p=0)=g_n({\xi }_n,K)\varphi (\Delta t){\xi }_n^m+{\xi }_n,$$

$${\xi }_{n+1}(p=1)=\frac{{\xi }_n}{1+g_n({\xi }_n,K)\varphi (\Delta t){\xi }_n^{m-1}}.$$

$$0\le {\xi }_{n+1}(p=1)\le {\xi }_{n+1}(p=0).$$

$$\frac{dx(t)}{dt}=-\frac{x^3}{x^2+\theta }.$$

$$\displaystyle F(x_{n},\theta)=-\frac{1}{2}\bigg{[}0+2x_{n}^{2}x_{n+1}\bigg{]}\frac{1}{x_{n}^{2}+\theta}=-\frac{x_{n}^{2}x_{n+1}}{x_{n}^{2}+\theta}.$$

$$\frac{x_{n+1}-x_n}{\varphi (\Delta t)}=-\frac{x_n^2{x}_{n+1}}{x_n^2+\theta },$$

$$\displaystyle x_{n+1}=\frac{x_{n}(x_{n}^{2}+\theta)}{x_{n}^{2}\big{(}1+\phi(\Delta t)\big{)}+\theta},$$

$$\left\{\begin{array}[]{l}\frac{h_{n+1}-h_{n}}{\phi(\Delta t)}=ax_{n}\\ \frac{x_{n+1}-x_{n}}{\phi(\Delta t)}=r_{1}x_{n}-h_{n}x_{n+1}-\frac{\nu px_{n+1}z_{n}}{b+px_{n}+z_{n}}+\alpha_{12}\frac{x_{n}y_{n}}{\epsilon_{1}+y_{n}}\\ \frac{y_{n+1}-y_{n}}{\phi(\Delta t)}=r_{2}y_{n}+\alpha_{21}\frac{x_{n}y_{n}}{\epsilon_{2}+x_{n}}-d_{2}y_{n}y_{n+1}\\ \frac{z_{n+1}-z_{n}}{\phi(\Delta t)}=-\frac{\nu z_{n}z_{n+1}}{b+px_{n}+z_{n}}-ky_{n}z_{n+1},\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}h_{n+1}=\Phi_{1}(h_{n},x_{n},y_{n},z_{n})\\ x_{n+1}=\Phi_{2}(h_{n},x_{n},y_{n},z_{n})\\ y_{n+1}=\Phi_{3}(h_{n},x_{n},y_{n},z_{n})\\ z_{n+1}=\Phi_{4}(h_{n},x_{n},y_{n},z_{n})\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\Phi_{1}(h_{n},x_{n},y_{n},z_{n})=ax_{n}\phi(\Delta t)+h_{n}\\ \Phi_{2}(h_{n},x_{n},y_{n},z_{n})=\frac{x_{n}(b+px_{n}+z_{n})\big{(}\left(\epsilon_{1}+y_{n}\right)\left(r_{1}\phi(\Delta t)+1\right)+\alpha_{12}y_{n}\phi(\Delta t)\big{)}}{\left(\epsilon_{1}+y_{n}\right)\big{(}\phi(\Delta t)(h_{n}(b+px_{n}+z_{n})+\nu pz_{n})+b+px_{n}+z_{n}\big{)}}\\ \Phi_{3}(h_{n},x_{n},y_{n},z_{n})=\frac{y_{n}\left(\epsilon_{2}+x_{n}\right)\left(r_{2}\phi+1\right)+\alpha_{21}x_{n}y_{n}\phi(\Delta t)}{\left(\epsilon_{2}+x_{n}\right)\left(d_{2}y_{n}\phi(\Delta t)+1\right)}\\ \Phi_{4}(h_{n},x_{n},y_{n},z_{n})=\frac{z_{n}(b+px_{n}+z_{n})}{\phi(\Delta t)\big{(}ky_{n}(b+px_{n}+z_{n})+\nu z_{n}\big{)}+b+px_{n}+z_{n}}.\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\Phi_{1}(h_{n},x_{n},y_{n},z_{n})=h_{n}\\ \Phi_{2}(h_{n},x_{n},y_{n},z_{n})=x_{n}\\ \Phi_{3}(h_{n},x_{n},y_{n},z_{n})=y_{n}\\ \Phi_{4}(h_{n},x_{n},y_{n},z_{n})=z_{n}\\ \end{array}\right.$$

$$\displaystyle\Gamma^{1}=(\tilde{h}_{1},0,0,0)\,\,\,and\,\,\,\Gamma^{2}=\bigg{(}\tilde{h}_{2},0,\frac{r_{2}}{d_{2}},0\bigg{)}.$$

$$X_{n+1}=A{X}_n+F(X_n),$$

$$Aq=\lambda q,pA=\lambda p,pq=1.$$

$$X=qu+v.$$

$$v=G(u)=C_2{u}^2+C_3{u}^3+\dots ,$$

$$F(qu+G(u))=F_2{u}^2+F_3{u}^3+\dots ,$$

$$\left\{\begin{array}[]{l}C_{2}=(\lambda^{2}I-A)^{-1}(I-qp)F_{2},\\ C_{3}=(\lambda^{3}I-A)^{-1}\left((I-qp)F_{3}-2\lambda pF_{2}C_{2}\right),\\ \vdots\end{array}\right.$$

Full text

1 Introduction

It is well known that aphids and ants share a symbiotic relationship: ants collect honeydew produced by aphids using it as food and the ants reciprocate by protecting aphids against predators [18, 16]. In such cases the symbiotic relationship is defined mutualism [1]. Holland and De Angelis in 2010 built a mathematical model that links the consumer functional responses of a mutualistic species with resources supplied by another; through the phase-plane analysis they shown that their set of differential equations correctly predicts the enhanced population growth rates of both species sharing the mutualistic relationship. Although in this paper the relationship between aphids and ants was hypothesized as mutualistic, an alternative and interesting vision is given in [20].

One of the most aggressive predators threatening the aphids colonies is the ladybird [18], in both larval and adult stages. The voracity of the ladybirds has left researchers thinking to use them as a biological control agent against the proliferation of the aphids, since they constitute a problem for crops [23]. Although there is a large scientific community which considers the contribution of the ladybirds to be effective in eradication of aphids, several authors demonstrated through mathematical models supported by experimental tests that they are ineffective [5, 7]. For a complete discussion and the full bibliography available on this topic refer to [7]. The mechanism regulating the population dynamics of ladybirds and aphids has been described by Kindlmann and Dixon in 2003 [5]. The authors proposed a population dynamical model that incorporates an optimization of egg distribution, offering an explanation as to why ladybirds have little effect on the aphids population dynamics.

The relationship between aphids, ants and ladybirds has been known and studied from a an ecological [19] and chemicals point of view [21], while a mathematical model that describes this dynamical system it seems not to have been developed yet.

This paper introduces a new set of differential equations describing the within-season population dynamics of aphids, ladybirds and ants. It includes the model proposed by Kindlmann in 2003 to reproduce the within-season aphids-ladybirds interaction and extends it by adding the mutualistic aphids-ants relationship, following the model proposed by Holland and De Angelis in [8].

Starting from the continuous-time model, a discrete-time version was proposed by applying the Non Standard Finite Difference (NSFD) scheme [10]. This mathematical framework allows to numerically integrate nonlinear differential equations of interest for a large variety of scientific fields, such as equations modeling the stellar structure, the dynamics of HIV transmission, heat transport equations and many other topics [13]. The high reliability of the NSFD respect to standard finite difference schemes proved in several works [10, 11, 13, 14]. Following the NSFD rules, the continuous-time model proposed by Holland and De Angelis in 2010 has been converted into a set of difference equations, proving that the NSFD performs better respect to a standard finite difference approach during the transient dynamics, especially with high time steps [4]; since the mutualistic interaction has a fundamental role in the model proposed in this paper, the discretization measures that have been successfully applied in [4] have been adopted also in this study.

2 Aphids, ants and ladybirds: the continuous-time model

In this section a set of differential equations describing the within-season population dynamics of a system composed by aphids, ants and ladybirds living in an ecological patch is proposed. The system, that for conciseness reasons can be cited later as AAL, is the following:

[TABLE]

with initial conditions $h(0)=0,\,x(0)=x_{0}\geq 0,\,y(0)=y_{0}\geq 0\,z(0)=z_{0}\geq 0$. The (1d) is a system of four non-linear differential equations regulating $h(t),x(t),y(t),z(t)$, respectively the cumulative density of the aphids, densities of aphids, ants and ladybirds at time $t$. Meaningful values of $h(t),x(t),y(t),z(t)$ are non-negative.

The parameter $a$ is the scaling constant relating aphids cumulative density to its own dynamics; $r_{1}$ and $r_{2}$ are respectively the maximum potential growth rate of the aphids and ants; $\nu$ is the ladybirds voracity; $p$ is the ladybirds preference for aphids; $b$ is a parameter of the functional response of the aphids; $\alpha_{12}$ is a positive term quantifying the advantages that the presence of ants induces on the growth of aphids; $\alpha_{21}$ is a positive term quantifying the advantages that the presence of aphids induces on the growth of ants; $\epsilon_{1}$ and $\epsilon_{2}$ are the half-saturation terms for aphids and ants respectively; $k$ is the coefficient of interspecific competition between ants and ladybirds; $d_{2}$ is the self-limiting term of ants population. The initial condition $z_{0}$ is defined by the number of eggs laid there by adults [5]. In order to clarify how the (1d) works, each equation is described below:

• •

The formulation given in (1a) expresses the cumulative density $h(t)$ as a regulatory term for aphids, instead of the instantaneous density [5].

• •

The (1b) describes changes in aphids density through the sum of three contributes: the first is an auto-regulatory term allowing aphids population density to decline even in the absence of natural enemies [5]; the second expresses the decrease of the aphids population, related to predator voracity $\nu$ and predator’s preference for prey $p$ [5]; the last term, modulated by the coefficient $\alpha_{12}$, quantifies the advantages that the presence of ants induces on the growth of the aphids, according to the research of Holland and De Angelis about the mutualistic relationships [8].

• •

The (1c) reproduces the ants population density, in which the first and last terms represent respectively the linear effect due to the intrinsic population growth rate and a quadratic one to modify the growth rate with density dependent self–limitation; the second term, regulated by the coefficient $\alpha_{21}$, quantifies the advantages that the presence of aphids induces on the growth of the ants.

• •

The (1d) describes the decreases of the ladybirds population caused by two terms: a first reproducing the cannibalistic effect [5] and a second one expressing the competition between ladybirds and ants [9].

The analysis of the continuous-time model was not reported in the paper due to the achievement of partial results. The complexity of the proposed set of differential equations makes the effort to carry out the stability analysis of the steady-states remarkable, and suitable to be dealt in a separate work.

3 Discrete-time model

To carry out this step, the continuous variable $t\in[0,\infty)$ must be replaced by the discrete variable $n\in\,\mathbb{N}$ and the variables $h(t),x(t),y(t)$ and $z(t)$ must take discrete values $h_{n},x_{n},y_{n}$ and $z_{n}$. The result is a set of difference equations.

The numerical integration of ordinary differential equations using traditional methods could produce different solutions from those of the original ODE [2, 14]. In particular, using a discretization step-size larger than some relevant time scale, is possible to obtain solutions that may not reflect the dynamics of the original system. To overcome this problem, Ronald Mickens, in 1989, suggested what is known as the Nonstandard Finite Difference (NSFD) method [10].

As introduced, in the AAL model the mutualistic interaction has an important role. Since in the paper [4] the NSFD scheme best has performed respect to a standard finite difference scheme, in the present research the author has followed the same way, extending the approach also for the other terms of the AAL dynamical system.

3.1 Non-Standard Finite Difference schemes (NSFD)

Let $f:\mathbb{R}^{n}\longrightarrow\mathbb{R}^{n}$ sufficiently smoothed, and $\xi(t):[0,+\infty)\longrightarrow\mathbb{R}^{n}$ the coordinates. Given the differential equation

[TABLE]

with initial condition $\xi_{0}=\xi(t=t_{0})\in\mathbb{R}^{n}$ and the system parameters identified by $K=(K_{1},K_{2},\ldots)$. We also suppose that

[TABLE]

is the difference equation corresponding to (2). In order to construct a NSFD scheme, the following rules have to be respected [14]:

1. I.

The order of the discrete derivative should be equal to the order of the corresponding derivatives of the differential equation. 2. II.

Denominator functions for the discrete representation must be nontrivial. The following replacement is then required:

[TABLE]

where $\phi(\Delta t)$ is such that $0<\phi(\Delta t)<1$, $\forall\Delta t>0$. The explicit form of $\phi(\Delta t)$ is given by the following expression [12]:

[TABLE]

Let the Jacobian $J_{\Phi}(\xi)=(j_{ij})_{4\times 4}$, with $\xi=(h_{n},x_{n},y_{n},z_{n})$ and

[TABLE]

let $\Omega$ the set of the eigenvalues of the Jacobian, the optimal value of $q$ must respect the condition

[TABLE] 3. III.

Nonlinear terms could be replaced by nonlocal discrete representations. 4. IV.

Special conditions that hold for the solutions of the differential equations should also hold for the solutions of the finite difference scheme. 5. V.

The scheme should not introduces spurious solutions.

Definition 3.1** (PESN condition).**

A NSFD scheme preserving the solution positivity is called Positive and Elementary Stable Nonstandard (PESN) method.

The positivity condition is particularly advantageous to discretize dynamical systems describing a population evolution, being a population density never negative.

3.1.1 Nonlocal representation and positivity

In order to construct a positivity-preserving scheme, it is required that the nonlocal representation of the nonlinear terms (condition III) is chosen carefully. To the author’s knowledge, there are no attempts to formalize the way in which the nonlocal approximation are constructed, leaving this fundamental step to the experience. In this paper, a simple rule was proposed.

Theorem 3.1**.**

Let the differential equation (2), $g(\xi,K)$ a function depending in nonpolynomial way on $\xi$, always positive for all $\xi\in\mathbb{R}^{n}$, and $K_{1},K_{2},\ldots\in\mathbb{R}^{+}$. If

[TABLE]

with $m\in\mathbb{N}^{+}$ and $p=0,1$. Then, a sufficient condition to have a NSFD respecting the positivity condition is:

[TABLE]

Proof.

By substituting the (9) in (3),

[TABLE]

[TABLE]

It is straightforward to show that

[TABLE]

and

[TABLE]

Being $\xi_{n+1}(p=0)\geq\xi_{n}$ and $\xi_{n+1}(p=1)\leq\xi_{n}$, it follows that

[TABLE]

The (14) shows the positivity of the solution. ∎

Remark*.*

The rule (9) allows to build a minimalist positivity-preserving finite difference scheme. More advanced nonlocal representations may be required [14].

Example 3.1**.**

Let the differential equation

[TABLE]

with $\theta\in\mathbb{R}^{+}$. By comparing (15) with (8), it follows that $g(x,\theta)=\frac{1}{x^{2}+\theta}$, $p=1$ and $m=3$. Following the (9), it result

[TABLE]

The NSFD scheme is:

[TABLE]

Then, with a simple algebraic manipulation

[TABLE]

that is a positivity-preserving scheme.

3.2 AAL: the discrete-time model

The finite difference scheme of the AAL dynamical system was built by extending the work carried out in [4], in which a stable NSFD scheme for a system of two populations in mutualistic relationship has been proposed. In that work, the proposed NSFD scheme respected the rules listed in the Subsection 3.1, including the nonlocal approximation to ensure the positivity of the solutions (although not yet formalized). For the AAL, the same approach was followed, by extending the NSFD rules also to the non-mutualistic terms. The set of difference equations proposed is the following:

[TABLE]

with the conditions $h_{0}=0,\,x_{0}\geq 0,\,y_{0}\geq 0\,z_{0}\geq 0$. Each right term of the equations 19 are obtained by taking advantage from (9). The (19) were explicated respect to $h_{n+1},x_{n+1},y_{n+1},z_{n+1}$ as follows:

[TABLE]

in which

[TABLE]

All the difference equations of (20) are positive, being positive all the parameters and variables.

3.3 Steady-state analysis

To find the steady-states of (21) the following condition has to be satisfied:

[TABLE]

By analytically solving the (22), the following steady-states $\Gamma=(\tilde{h},\tilde{x},\tilde{y},\tilde{z})$ were found:

[TABLE]

The first equation in (22) implies that only $\tilde{x}_{n}=0$ is admissible for $\phi(\Delta t)>0$, while it gives no information about the asymptotic values $\tilde{h}_{1}$ and $\tilde{h}_{2}$. In order to estimate them, and complete the knowledge about $\Gamma^{1}$ and $\Gamma^{2}$, a numerical simulation was carried out by assigning the following values to the parameters: $r_{1}=r_{2}=0.3$, $a=0.000005$, $k=0.005$, $\nu=1$, $b=0$ and $p=1$ according to [22]; $\alpha_{12}=\alpha_{21}=0.6$, $d_{2}=0.01$ and $\epsilon_{1}=\epsilon_{2}=0.3$ according to [8, 4]; the discretization parameters are $dt=10^{-5}$ and $q=1.5$. After $10^{7}$ iterations, by assuming the initial values $(h_{0},x_{0},y_{0},z_{0})=(0,5,0,1)$, it results $\tilde{h}_{1}\simeq 0.599525$; let $(h_{0},x_{0},y_{0},z_{0})=(0,100,40,40)$, it results $\tilde{h}_{2}\simeq 1.796894$. Also, by varying the initial conditions to $(h_{0},x_{0},y_{0},z_{0})=(0,0,0,1)$ a new fixed point $\Gamma^{0}=(0,0,0,0)$ arises. Finally, with $(h_{0},x_{0},y_{0},z_{0})=(0,0,40,40)$, a fourth steady-state $\Gamma^{3}=(0,0,30,0)$ was obtained.

In order to know the nature of the four steady-states $\Gamma^{0}-\Gamma^{3}$, the eigenvalues of the Jacobian matrix need to be calculated.

Definition 3.2**.**

Given $\sigma(J_{\Phi})$ the set of the Jacobian eigenvalues, a steady-state $\tilde{\xi}\in\mathbb{R}^{n}$ is defined:

• •

hyperbolic if $\nexists\lambda\in\sigma(J_{\Phi})$ so that $|\lambda|=1$,

• •

nonhyperbolic if $\exists_{1}\lambda\in\sigma(J_{\Phi})$ so that $|\lambda|=1$.

For each steady-state the eigenvalues are summarized in Table 1 and represented in the complex plane of Figure 1: the eigenvalues characterized by the same symbol are relative to the same steady-state. Since all the four steady-states have a unitary eigenvalue, following the Definition 3.2 they are nonhyperbolic.

4 Stability analysis

For nonhyperbolic fixed points the Hartman–Grobman is not applicable and the effort to infer about their stability is significantly greater. The centre manifold theory helps when some of the eigenvalues of the Jacobian have unitary absolute value and the others have absolute value less than one. It allows reducing the dimensionality of a multi-dimensional dynamical system around a nonhyperbolic fixed point and determine its stability properties by studying the centre manifold [24]. As visible in Figure 1 and Table 1, the steady-states $\Gamma^{0},\Gamma^{1}$ and $\Gamma^{3}$ have one eigenvalue with absolute value falling outside the unitary circle, then the centre manifold theory is not applicable. In the present research the stability of these fixed points is not analytically shown but inferred by means the phase-portrait analysis (see the Subsection 5.1). $\Gamma^{2}$ has one eigenvalue with unitary module and the others less than one, then the centre manifold theory is applicable. Nonetheless, from a strictly ecological point of view, the initial conditions leading to $\Gamma^{2}$ are the most realistic, then the effort to determinate the stability of $\Gamma^{2}$ is required.

4.1 Centre Manifold Theory

Following the approach of [17], consider the $m$-dimensional system

[TABLE]

where $X_{n}\in\mathbb{R}^{m}$ and $F(X_{n})=O(||X||^{2})$. Let $\lambda$ an eigenvalue of $A$, $p\in\mathbb{R}^{1\times m}$ and $q\in\mathbb{R}^{m}$ the eigenvectors which satisfy

[TABLE]

Let $u=pX\in\mathbb{R}$ and $v=X-qu\in\mathbb{R}^{m}$, then is simple to show that $X$ can be decomposed as

[TABLE]

If $A$ has only one eigenvalue such that $|\lambda|=1$ and the other less than one, there exists a function $v=G(u)$ such that $G(0)=0$ and $G^{\prime}(0)=0$. Let the center manifold

[TABLE]

and

[TABLE]

where $C_{2}$ and $C_{3}$ are given by:

[TABLE]

Substituting (27) in (26):

[TABLE]

then,

[TABLE]

in which $F_{2}$ and $F_{3}$ can be evaluated by comparison of (24) and (28). Given the centre manifold described by (31), in order to evaluate the stability of the nonhyperbolic fixed point, it is possible to make use of the following theorem:

Theorem 4.1** (See [3]).**

Let $\tilde{\xi}$ be a fixed point of a map $f$ such that $f^{\prime}(\tilde{\xi})=1$. If $f^{\prime}(\tilde{\xi})$, $f^{\prime\prime}(\tilde{\xi})$, and $f^{\prime\prime\prime}(\tilde{\xi})$ are continuous at $\tilde{\xi}$, then the following statements hold:

• A.

if $f^{\prime\prime}(\tilde{\xi})\neq 0$, then $\tilde{\xi}$ is unstable;

• B.

if $f^{\prime\prime}(\tilde{\xi})=0$ and $f^{\prime\prime\prime}(\tilde{\xi})>0$ then $\tilde{\xi}$ is unstable;

• C.

if $f^{\prime\prime}(\tilde{\xi})=0$ and $f^{\prime\prime\prime}(\tilde{\xi})<0$ then $\tilde{\xi}$ is asymptotically stable.

4.2 Stability of the steady-state $\Gamma^{2}$

In this part of the study only the most important steps were reported. The eigenvectors left $p$ and right $q$ associated to the unitary eigenvalue of $\Gamma^{2}$ are:

[TABLE]

The (20) was first rewritten by evaluating the Taylor series around $\Gamma^{2}$, in order to simplify the equations. Making the change of variables

[TABLE]

the fixed point was shifted to the origin. The set of difference equations can be written

[TABLE]

in which $\hat{\Phi}_{i}(\hat{h}_{n},\hat{x}_{n},\hat{y}_{n},\hat{z}_{n})$, $i=1,2,3,4$, expresses the function $\Phi(h_{n},x_{n},y_{n},z_{n})$ after both Taylor series expansion and the change of variable. The system can be rewritten

[TABLE]

in which $A$ is a $4\times 4$ matrix containing the coefficients of the linear terms and $F_{i}(\hat{h}_{n},\hat{x}_{n},\hat{y}_{n},\hat{z}_{n})$ represent the nonlinear ones. By using the same parameters values adopted in Subsection 3.3, the matrix A is:

[TABLE]

With regard to $F(\hat{h}_{n},\hat{x}_{n},\hat{y}_{n},\hat{z}_{n})$, by neglecting the terms with coefficients less than $10^{-7}$, it is:

[TABLE]

Taking into account the first equation of (34),

[TABLE]

it is possible to show that the centre manifold equation is:

[TABLE]

that respects the conditions $f^{\prime\prime}(\tilde{u})=0,\,\,\,f^{\prime\prime\prime}(\tilde{u})<0$, implying that $\Gamma^{2}$ is asymptotically stable.

A necessary, but not sufficient, condition for bifurcation of a fixed point is that the fixed point is nonhyperbolic [24]; in the present research the bifurcation problem was not addressed due to the complexity of the equations.

5 Numerical simulations

5.1 Phase-space analysis

The phase-space analysis of the system (20) was carried out and, to facilitate the visualization, three phase-portraits, containing the $x$-$y$, $x$-$z$ and $y$-$z$ populations, are represented respectively in figures 2b,c and d. Each phase-portrait was calculated by considering the initial conditions displayed in Figure 2a: $h_{0}=0$, $0\leq x,y,z\leq 60$. The initial conditions were reported also in the three phase-portraits, however, in plots b and c, the $x$-axis is in $10^{4}$ units due to the important range of variation of the aphids population, then the range of variation of $x$ initial conditions is not recognizable. Each initial condition point of Figure 2a is plotted with a symbol indicating the steady-state it leads to (i.e. the points characterized by $x=0,y=0,0\leq z\leq 60$ are marked by an empty circle and lead to the steady-state $\Gamma^{0}$, according to the legend).

Note that where the steady-states are superimposed each other it means that the variables considered in that phase-portrait do not make the difference, e.g., in $x$-$z$ plane represented in Figure 2b, $\Gamma^{2}$ and $\Gamma^{3}$ are coincident since is the variable $h$ (not shown) that makes the difference.

About the phase-space analysis, it is possible to deduce that:

• •

the steady-state $\Gamma^{0}$ is reached if initial populations are characterized by $x_{0}=y_{0}=0$; $\Gamma^{1}$ is reached if $y_{0}=0$, regardless of $x_{0}$ and $z_{0}$;

• •

$\Gamma^{3}$ is reached if $x_{0}=0$, regardless of $y_{0}$ and $z_{0}$;

• •

the system evolves toward $\Gamma^{2}$ for any other initial condition. In this case, the phase-portraits show that the trajectories experience a peak of the $x$ and $y$ populations and a convergence towards $\Gamma^{2}$.

Numerical values of the maximum populations will be given in the next Section.

5.2 Time series analysis

The trajectories produced by the system (20) were depicted in Figure 3 by using the initial conditions $h_{0}=0,x_{0}=y_{0}=60,z_{0}=20$. The trajectories for the couples $h,x$ and $y,z$ were represented respectively in Figure 3a and b. As we can see, the system evolves towards the steady-state $\Gamma^{2}$. The aphid density grows up to a maximum of $8058$ individuals at time $t=9.96$ and drops to 0 individuals for $t\gtrsim 20$. As a consequence of the aphids-ants mutualistic relationship, the ants increase their density population up to 90 individuals in the period of greatest abundance of aphids, although they must face a sharp decline and a stabilization on $30$ individuals after the density population of aphids collapses. The ladybird population density results monotonic decreasing, becoming zero together with the aphids.

Nonetheless, the shape of their population curve is sensitively influenced by the $k$ coefficient of interspecific competition between ants and ladybirds. In Figure 4b the $z$ population was calculated by assuming $k$ in range $10^{-5}$–$10^{-2}$. It is possible to see that lower $k$ values make the ladybird population supported by the aphid colony size. In particular, with $k=10^{-5}$–$10^{-6}$ the ladybird colony is almost constant during the period of abundance of the aphids, while greater values imply an increasingly rapid decline of the ladybirds population size. The maximum of the aphids population density is also affected by the magnitude of $k$, although there is a variation of only 100 individuals about between $k=10^{-5}$ and $10^{-2}$, see Figure 4a.

The influence of the parameter $\alpha_{12}$, quantifying the advantages that the presence of ants induces on the growth of aphids, was also studied. From the curves represented in Figure 4c it is possible to conclude that the maximum of the aphids abundance sensitively increases with $\alpha_{12}$, while this parameter scarcely influences the ants density. Conversely, the parameter $\alpha_{21}$, quantifying the advantages that the presence of aphids induces on the growth of ants, sensitively influences the growth of the ant colony, determining only a small percentage growth of the aphid colony, as visible in Figure 4d.

6 Conclusions

By taking advantage of the mathematical model proposed by Kindlmann & Dixon in 2003 [5] to describe the ladybirds-aphids (prey-predator) population density and the model built by Holland & De Angelis in 2010 [8] describing the population dynamics of two mutualistic species, in this research a new set of differential equations was proposed to describe the within-season population dynamics of an ecological patch hosting a community of ladybirds, aphids and ants. This framework describes the population dynamics in presence of both prey-predator and mutualistic relationships, based on the following facts:

• •

between aphids and ants exists a well documented mutualistic relationship;

• •

ladybirds behaves as predator of aphids;

• •

ladybirds and ants may experience an interspecific competition.

The proposed model, first formulated in continuous-time, then in discrete-time domain by taking advantage of the NSFD scheme, is characterized by four nonhyperbolic steady-states $\Gamma^{0}$-$\Gamma^{3}$. The phase-space analysis highlighted that the initial conditions most realistic from an ecological point of view lead to the steady-state $\Gamma^{2}$, characterized by a collapsed colony of aphids and ladybirds and a population of ants equal to $r_{2}/d_{2}$ ($r_{2}$ is their maximum potential growth rate and $d_{2}$ their self-limiting term). The asymptotic stability of $\Gamma^{2}$ was also shown via the centre manifold theory. It is an feature that the mutualistic relationship does not influence any of the asymptotic population densities.

The time-series analysis shown that, with sufficiently large initial populations, the populations of aphids and ladybirds reach a maximum, then definitely collapse according to the $\Gamma^{2}$ population densities. Then, the presence of aphids and ladybirds constitutes a transitory phenomenon. Nonetheless, the presence of the mutualistic relationship determine a sensitive grows of the aphids and ants population densities peaks.

Since the four steady-states are nonhyperbolic, future research activities should be carried out to investigate the presence of bifurcations. Furthermore, a validation with experimental data, actually not publicly available to the author’s knowledge, should be carried out in order to validate the reliability of the AAL model in predicting the populations evolution.

Acknowledgments

I would like to thank Professors W. Nazeer, M. Imran and to an anonymous reviewer for giving me precious suggestions addressed to improve the quality of the paper.