# Enhancement of Extreme Events through the Allee effect and its Mitigation
through Noise in a Three Species System

Deeptajyoti Sen  Indian Institute of Science Education and Research Mohali,
Knowledge City, Sector 81, Manauli PO 140306, India
Sudeshna Sinha  Indian Institute of Science Education and Research Mohali,
Knowledge City, Sector 81, Manauli PO 140306, India
<EMAIL_ADDRESS>
###### Abstract

We consider the dynamics of a three-species system incorporating the Allee
Effect, focussing on its influence on the emergence of extreme events in the
system. First we find that under Allee effect the regular periodic dynamics
changes to chaotic. Further, we find that the system exhibits unbounded growth
in the vegetation population after a critical value of the Allee parameter.
The most significant finding is the observation of a critical Allee parameter
beyond which the probability of obtaining extreme events becomes non-zero for
all three population densities. Though the emergence of extreme events in the
predator population is not affected much by the Allee effect, the prey
population shows a sharp increase in the probability of obtaining extreme
events after a threshold value of the Allee parameter, and the vegetation
population also yields extreme events for sufficiently strong Allee effect.
Lastly we consider the influence of additive noise on extreme events. First,
we find that noise tames the unbounded vegetation growth induced by Allee
effect. More interestingly, we demonstrate that stochasticity drastically
diminishes the probability of extreme events in all three populations. In fact
for sufficiently high noise, we do not observe any more extreme events in the
system. This suggests that noise can mitigate extreme events, and has
potentially important bearing on the observability of extreme events in
naturally occurring systems.

##  1 Introduction

The emergence of extreme events in the dynamical evolution of systems ranging
from weather [1] to power grids [2, 3] have catastrophic implications. So
understanding the underlying mechanisms that may trigger extreme events have
commanded considerable recent research interest [4]. An extreme event may be
defined as an event where one or more variables of a system, arising in nature
or in the laboratory, exhibits very large deviations from the mean value. So
the dynamics of the system is characterized by excursions to values that
differ significantly from the average. Further, though recurrent, these large
deviations are rare vis-a-vis the characteristic time scale of the system, and
their occurrences are aperiodic and uncorrelated in time. Without loss of
generality, an event is typically labelled ‘extreme’ if a state variable, in
the course of its temporal evolution, takes values that are several standard
deviations away from the average value, thereby signalling dynamical behaviour
beyond normal variability. Such extreme events have been observed in natural
systems such as rogue ocean waves [5], laboratory systems such as optical
systems [6], as well as financial phenomena like market crashes [7].
Interestingly, other definitions of extreme events, more appropriate to the
context, have also been employed. Notably, in the specific important problem
relating climate change to ecological dynamics, a synthetic definition of
extreme events involving both the driver and response system is proposed [8].
However in this work we will employ the most commonly used marker for extreme
events: namely, recurrent and aperiodic deviations larger than a prescribed
threshold from the mean value, are considered to be extreme events, with the
threshold typically taken to be 3-8 standard deviations from the mean.

A central direction in understanding extreme events is to find generic
mechanisms that can give rise such large fluctuations. Typical studies of
extreme events have involved stochastic models [9, 10], for instance random
walk models of transport on networks [11]. The emergence of extreme events in
deterministic dynamical systems, manifested as intermittent large amplitude
events, have also been investigated recently. Broadly speaking, the
statistical features of deterministic systems is an active research direction
that can lead to understanding extreme events that arise in the context of
deterministic dynamics [12, 13]. The search for dynamical systems that yield
extreme events, without the drive of external stochastic influences or
intrinsic random fluctuations, is a focus of much ongoing research effort from
the point of view of basic understanding of complex systems [14, 15, 16, 17].
Additionally, such probabilistic outcomes in dynamical systems are most
relevant in the applied context as well, such as in engineering sciences where
this direction of research leads to better assessment of risks [18].

In this broad direction, in this work we explore the dynamics of a vegetation-
prey-predator system, coupled through interactions of the Lotka-Volterra type.
Importantly, our model system also incorporates the biologically significant
Allee effect, which is one of the classic phenomena of population ecology. The
Allee effect reflects the beneficial effects on the growth of individuals
arising from conspecific interactions [19, 20, 21], and has bearing on the
long-term persistence of a population as a consequence of small size. Further
we investigate the effect of additive noise in the system, focusing on the
role of stochasticity in the emergence of extreme events. In a larger context,
this research direction also has bearing on the important general question of
the emergence of extremely large events in deterministic dynamical systems,
and the effect of noise on their sustained prevalence.

Our central findings in this work are as follows: the system under Allee
effect yields chaotic dynamics. Further, it leads to unbounded vegetation
growth for sufficiently strong Allee effect. Most significantly, Allee effect
aids the emergence of extreme events in this three-species chain.
Interestingly, we also observe that noise suppresses the unbounded blow-up of
vegetation induced by Allee effect. Lastly, sufficiently strong noise also
subdues the extreme events in vegetation, prey and predator populations, thus
suggesting a significant natural mechanism to mitigate extreme events in
population chains.

In Section 2 we present results arising in the model of three interacting
species, incorporating the Allee effect. In Section 3 we explore the effect of
additive noise in the system. We conclude with discussions of the scope of our
findings in Section 4.

##  2 Three-species food chain model incorporating the Allee Effect

Complex systems research in general, and theoretical ecology in particular,
has seen intense research activity in networks modelling interacting species,
often focusing on local and global stability properties [22, 23, 24, 25]. Here
we will focus on the emergence of extreme events and consider as our test-bed
the well-known model for the dynamics of the snowshoe hare and the Canadian
lynx populations, based on observed data [22]. Specifically, the system
incorporates a three species vertical food chain, consisting of vegetation
(denoted by $u$), prey (denoted by $v$) and predator (denoted by $w$).
Additionally, we will incorporate in this model a term to reflect the Allee
effect in the growth of the prey. The dynamics of this three species trophic
system is described by the following sets of coupled differential equations:

$\begin{split}\dot{u}&=f(u,v,w)\,=\,au\ -\ \alpha_{1}f_{1}(u,v),\\\
\dot{v}&=g(u,v,w)\,=\,\alpha_{1}f_{1}(u,v)\ A(v)\ -\ bv\ -\
\alpha_{2}f_{2}(v,w),\\\ \dot{w}&=h(u,v,w)\,=\,\alpha_{2}f_{2}(v,w)\ -\
c(w-w^{*}).\end{split}$ (1)

The interaction between the vegetation and prey populations is considered to
follow the type II functional response, described by the function
$f_{1}(u,v)$, where
$f_{1}(u,v)\ =\ \frac{uv}{1+ku}$. This well-known functional response is
characterized by a decelerating intake rate, stemming from the assumption that
the consumer is limited by its capacity to process food. The parameter $k$ is
the average time spent on processing a food item, which is termed the handling
time in literature.

The interaction of the predator population with the prey is considered to
follow the well-known Lotka-Volterra type interaction, described by the
function $f_{2}(v,w)$, where $f_{2}(v,w)\ =\ vw$. Here $\alpha_{1}$ denotes
the maximum growth rate of the prey, which is in general a product of
ingestion rate with a constant factor ($<1$), accounting the fact that not all
of the ingested resource (vegetation) converted into prey’s biomass. A similar
parameter for the predator is denoted by $\alpha_{2}$. The parameters $a$, $b$
and $c$ represent the intrinsic growth rates of the three species $u$, $v$ and
$w$ respectively. Further, the model allows the predator population to
maintain a equilibrium population $w^{*}$ when the prey concentration is very
low. In other words, predator can survive in the trophic system without
depending on prey.

Importantly, in contrast to work on similar systems [26], here we will also
explicitly consider the Allee effect. Specifically, the Allee effect is
considered in growth of the prey by introducing

$A(v)\,=\,\frac{v}{v+\theta},$ (2)

where $\theta$ is the Allee strength parameter representing the critical prey
density at which the probability of successful mating would be half. The Allee
effect here occurs due to difficulties in finding mates for sexual
reproduction and $A(v)$ describes the mating success at the low population
density [27, 28]. This kind of Allee effect in two and three dimensional
population models has drawn much research attention [27, 19, 20, 21, 29, 30].

In this work we consider the parameter values $a\,=\,1$, $b\,=\,1$,
$c\,=\,10$, $w^{*}\,=\,0.006$, $\alpha_{1}\,=\,0.5$, $\alpha_{2}\,=\,1$,
$k\,=\,0.05$ [22]. We explore the dynamics of the system under varying
$\theta$, through numerical simulations using the Runge-Kutta fourth order
algorithm. We have ascertained the stability and convergence of our results
with respect to decreasing step size.

Our main focus is to explore the dynamical consequences of the Allee effect in
the population of the prey, in this three-species model system. We will first
show that there is a sharp increase in the probability of obtaining unbounded
vegetation growth beyond a critical value of the Allee parameter. We will then
demonstrate that, interestingly, the system exhibits chaos as the Allee effect
becomes more significant. Then we will move on to the central focus of this
work, namely the influence of Allee effect on the emergence of extreme events.
We will explicitly demonstrate that there is a pronounced increase in the
propensity of extreme events under increasing Allee parameter $\theta$, and
also show that the extreme events occur with strikingly different
probabilities for the different species.

## Temporal evolution of the population densities

Our first observation is the emergence of explosive runaway growth in
vegetation when the Allee effect is too strong, i.e. when the Allee parameter
$\theta$ is sufficiently large, the vegetation grows in an unbounded manner
[31]. In order to quantify this blow-up, we estimate the probability of
unbounded vegetation growth from a large sample of random initial states,
followed over a long period of time. We also ascertain that the estimated
values are converged with respect to increasing sample size, and thus can be
considered to be robust numerically. The results thus obtained, for varying
Allee parameter $\theta$, are displayed in Fig. 1. It is clearly evident from
the figure that there exists a critical value of $\theta$, which we denote by
$\theta_{c}$, beyond which the vegetation has a probability of explosive
unbounded growth. So in the rest of this work, we will restrict our analysis
for the range $\theta\in[0,\theta_{c})$.

Figure 1: Probability of unbounded growth of the vegetation population ($u$)
with respect to $\theta$. Here an explosive blow-up is considered to have
occurred when the vegetation population exceeds a value of $10^{3}$. The
probability is estimated from a sample of $10^{3}$ initial states randomly
distributed in a hyper-cube ($u\in[0:4],v\in[0:2],w\in[0:5]$) in phase space.
Interestingly, the slight dip in the estimated probability of unbounded growth
in the narrow window around $\theta\sim 0.031$ arises from unbounded orbits
co-existing with a small set of initial states that evolve to bounded periodic
orbits.

Next we investigate the temporal evolution of the population densities. Fig. 2
shows representative time series of vegetation, prey and predator and the
corresponding attractors in 3-dimensional phase space. To broadly illustrate
the influence of Allee effect, we present this for three values of $\theta$,
with increasing magnitude. When Allee effect is absent in the prey population
i.e. $\theta\,=\,0$, we observe that all populations fluctuate periodically
and are confined to a periodic orbit [see Fig. 2a and Fig. 2b]. For a larger
Allee parameter ($\theta\,=\,0.01776$), the populations of vegetation, prey
and predator all evolve in an aperiodic manner, as evident in Fig. 2c, with
the corresponding chaotic attractor shown in Fig. 2d. On further increasing
$\theta$ to $0.02475$, the size of the chaotic attractor increases, as evident
from Fig. 2e,f. Therefore increasing magnitude of the Allee effect parameter
drives the system into a chaotic state from the periodic state. To corroborate
this observation, we present the bifurcation diagram of the prey population,
with respect to the Allee parameter $\theta$, in Fig. 3. We observe the onset
of chaos through the usual period-doubling cascade, initiating at period-$4$
at $\theta=0$. Subsequently, we also observe narrow periodic windows, with
period-$3$ being the most prominent one. The most significant implication of
this bifurcation diagram is the emergence of chaotic dynamics in a very large
range of the Allee parameter $\theta$. That is, with no Allee effect or under
very weak Allee effect the system is periodic, while a strong Allee effect
typically induces chaos in this three-species system.

Figure 2:  Left panels display the time series for the vegetation ($u$), prey
($v$) and predator ($w$) populations in the system given by Eqns. 1, and the
right panels display the corresponding phase space attractor. The Allee effect
parameter is $\theta\,=\,0$ (a-b), $\theta\,=\,0.01776$ (c-d), and
$\theta\,=\,0.02475$ (e-f). The red dashed line shows the mean $\mu$, and the
black dashed line represents the threshold level of $7$ standard deviations
above the mean (i.e. $\mu+7\sigma$). Figure 3: Bifurcation diagram of prey
populations with respect to Allee parameter $\theta$. Here we display the
local maxima of the prey population. The parameter values in Eqn. 1 are as
mentioned in the text.

## Extreme events induced by Allee Effect

One of the most interesting observations from the time series presented in the
section above is the following: when the magnitude of the Allee parameter
$\theta$ is low, vegetation and prey densities are confined to low values.
However, the predator densities deviate very significantly away from their
mean. Now for very small $\theta$ the system is attracted to a periodic orbit,
and so the large deviations are completely correlated with time and occur
periodically. So they cannot be considered to be extreme events, as they are
neither aperiodic, nor rare. But for larger $\theta$, both predator and prey
densities can sometime shoot up over 7 standard deviations away from the mean
value. This is evident clearly in Fig. 2(c) & (e) where one can see that both
predator and prey populations exceed the $7\sigma$ threshold from time to
time. The instants at which prey and predator populations exceed the $7\sigma$
threshold are now completely uncorrelated with time. This is consistent with
the underlying chaotic dynamics that emerges under increasing Allee parameter
$\theta$.

Figure 4: Figure marking the time instances at which a population exceeds the
$7\sigma$ threshold, for different values of Allee parameter $\theta$, for the
case of (top to bottom) vegetation, prey and predator populations.

In order to illustrate this, we mark the time instances at which a population
exceeds the $7\sigma$ threshold, for different values of Allee parameter
$\theta$. Fig. 4 shows this for the vegetation, prey and predator populations.
The density of points signifying the occurrence of extreme events is clearly
the highest for the predator population. This indicates that the predator
population has the greatest propensity for large deviations. It is also clear
that vegetation has the least number of extreme events in the same time
window. The uncorrelated nature of the extreme events is also evident in the
scatter of these points, except in the small periodic windows that occur for
certain special ranges of $\theta$. The increasing density of these points
also illustrate the increasing probability of extreme events in the
populations with increasing Allee parameter $\theta$.

In order to understand the phenomena quantitatively, we first estimate the
maximum densities of vegetation, prey and predator populations (denoted by
$u_{max}$, $v_{max}$ and $w_{max}$ respectively) for varying the Allee
parameter $\theta$. To estimate this, we find the global maximum of the
populations sampled over a time interval $T\,=\,1000$, averaged over a large
set of random initial conditions.

Figure 5: Global maximum of vegetation $u_{max}$ (blue), prey $v_{max}$ (red)
and predator (black) populations, with respect to the Allee parameter
$\theta$, scaled by their values obtained for $\theta\,=\,0$. Clearly, when
Allee parameter $\theta$ is sufficiently large, the maximum prey and predator
populations are an order of magnitude larger than that obtained in systems
with no Allee effect.

Fig. 5 shows $u_{max}$, $v_{max}$ and $w_{max}$, for Allee parameter
$\theta\in[0,\theta_{c})$, scaled by their values at $\theta=0$. These scaled
maxima help us gauge the relative change in the maximum population densities
arising due to the Allee effect. It is evident from our simulation results
that the magnitude of the global maximum of vegetation does not change very
significantly for increasing Allee parameter $\theta$, with its magnitude
around $\theta_{c}$ being approximately 4 fold the value at $\theta\,=\,0$.
However, the magnitude of maximum prey and predator populations change very
significantly with respect to Allee parameter $\theta$ and exceeds over $10$
fold the value obtained for $\theta\,=\,0$.

We then go on to numerically calculate the probability density of the
vegetation, prey, and predator population densities, for increasing Allee
effect parameter $\theta$. The tail of this probability density function
reflects the influence of the Allee effect on the probability of obtaining
extreme events. To illustrate this, we show the probability density function
for the prey population in Fig. 6, for three different values of $\theta$.
Extreme events are confined to the tail of the distribution that lie beyond
the vertical red line, marking the $\mu+7\sigma$ value in the figure. So it is
clear from these probability distributions that the Allee effect in prey
population promotes the occurrence of extreme events as the tail of the
distribution is flatter and extends further with increasing Allee parameter
$\theta$.

Figure 6: Probability Density Function (PDF) of the prey population $v$, for
the system given by Eqs. (1), with increasing magnitude of $\theta$ with (a)
$\theta\,=\,0$, (b) $\theta\,=\,0.015$ and (c) $\theta\,=\,0.02$. The
threshold for extreme event $\mu+7\sigma$ is denoted by vertical red dashed
line.

Figure 7: (Left) Return Map of $\Delta t_{i+1}$ vs $\Delta t_{i}$, and (right)
Probability distribution of $\Delta t_{i}$ fitted with exponentially decaying
function, where $\Delta t_{i}$ is the $i^{th}$ interval between successive
extreme events, where an extreme event is defined at the instant when the prey
population crosses the $\mu+7\sigma$ line (cf. Fig 2). Here $\theta=0.024$.
Figure 8: Probability of obtaining extreme event in unit time ($P_{ext}$),
with respect to Allee parameter $\theta$, estimated by sampling a time series
of length $T=5000$, and averaging over $500$ random initial states. Here we
consider that an extreme event occurs when a population level crosses the
threshold $\mu+7\sigma$. $P_{ext}$ for vegetation, prey and predator are
displayed in blue, red and black colors respectively. Note that there exists a
narrow periodic window around $\theta\sim 0.02$ (cf. Fig. 9), and so the large
deviations in this window of Allee parameter are not associated with true
extreme events, as they occur periodically.

In order to ascertain that the extreme values are uncorrelated and aperiodic
we examine the time intervals between successive extreme events in the
population. Fig. 7 (left panel) shows representative results for the return
map of the intervals between extreme events in the prey population and it is
clearly shows no regularity. The probability distribution of the intervals is
also Poisson distributed and so the extreme population buildups are
uncorrelated aperiodic events, as clearly evident from the right panel of the
figure.

In order to further quantify how Allee effect influences extreme events, we
estimate the probability of obtaining large deviations, in a large sample of
initial states tracked over a long period of time. We denote this probability
by $P_{ext}$, and we calculate it by following a large set of random initial
conditions and recording the number of occurrences of the population crossing
the threshold value in a prescribed period of time, with this time window
being several orders of magnitude larger than the mean oscillation period.
This time-averaged and ensemble-averaged quantity yields a good estimate of
$P_{ext}$. With no loss of generality, we choose the threshold for determining
extreme events to be $\mu+7\sigma$, i.e. when the variable crosses the
$7\sigma$ level, it is labelled as extreme.

This probability, estimated for all three populations is shown in Fig. 8.
First, it is clear from Fig. 8, that the probability of the occurrence of
extreme events is the lowest for vegetation, and the highest for predator
populations, for any value of the Allee parameter $\theta\in[0,\theta_{c})$.
We also observe that, for values of the Allee parameter $\theta$ lower than a
critical value denoted by $\theta^{u}_{c}$ the probability of obtaining
extreme events in the vegetation population tends to zero. Beyond the critical
value $\theta^{u}_{c}$, the vegetation population starts to exhibit extreme
events. A similar trend emerges for the prey population. However, the critical
value of the Allee parameter $\theta$ necessary for the emergence of a finite
probability of extreme events, denoted by $\theta^{v}_{c}$, is much smaller
than $\theta^{u}_{c}$. So for the prey population, a weaker Allee effect can
induce extreme events.

Note that some mechanisms have been proposed for the generation of extreme
events in deterministic dynamical systems, which typically have been excitable
systems. These include interior crisis, Pomeau-Manneville intermittency, and
the breakdown of quasiperiodic motion. However the extreme events generated by
these mechanisms occur typically at very specific critical points in parameter
space, or narrow windows around it. The first important difference in our
system here is that the extreme events do not emerge only at some special
values alone. Rather, there is a broad range in Allee parameter space where
extreme events have a very significant presence. This makes our extreme event
phenomenon more robust, and thus increases its potential observability. This
also rules out the intermittency-induced mechanisms that have been proposed,
as is evident through the lack of sudden expansion in attractor size in our
bifurcation diagram (Fig. 3) in general.

However, interestingly, the system does have one parameter window where there
is attractor widening and this gives rise to a markedly enhanced extreme event
count. The peak observed in Fig. 8 can be directly correlated with a sudden
attractor widening leading to a marked increase of extreme event in a narrow
window of parameter space located near the crisis (see Fig. 9). Additionally,
for a narrow window around $\theta\sim 0.02$, the emergent dynamics is
periodic. So the large deviations are no longer uncorrelated, and so they are
not extreme events in the true sense.

Figure 9: Bifurcation diagram of prey populations with respect to Allee
parameter, in the range $\theta\in[0.0189:0.0191]$. Here we display the local
maxima of the prey population. The parameter values in Eqn. 1 are as mentioned
in the text.

Lastly we notice that the predator population shows extreme events for all
values of $\theta\in[0,\theta_{c})$. So the predator population is most prone
to experiencing unusually large deviations from the mean. We also observe that
the probability of occurrence of extreme events in the predator population is
not affected significantly by the Allee effect. This is in marked contrast to
the case of vegetation and prey, where the Allee effect crucially influences
the advent of extreme events. Also, for the predator population there is no
marked transition from zero to finite $P_{ext}$ under increasing Allee
parameter $\theta$, as evident for vegetation and prey populations.

##  3 Effect of Noise in the System: Mitigation of Blow-Ups and Extreme Events

Most realistic population models are not deterministic, as noise is ubiquitous
in nature. So stochasticity must be incorporated into the models, since there
are many external influences, such as migration, diversity and environmental
fluctuation, present in the real ecosystem. For instance, in important earlier
works the role of environmental fluctuation dependent fitness in population
dynamics, namely parametric noise, in extinction and persistence has been
studied [32].

In this work, we explore the interplay of stochasticity and extreme events, by
investigating the system under the influence of additive noise. It is of much
relevance to explore if noise has any significant effect on the dynamics, for
instance on the unbounded vegetation growth under increasing Allee effect and
on the emergent chaotic attractors. The other question of utmost interest is
the following: does noise mitigate or aid the emergence of extreme events.
This will be the focus of our investigation in this section.

Specifically we investigate the dynamics of the three species system (1) under
additive random noise $\xi(t)$, given by the following dynamical equations:

$\begin{split}\dot{u}&=f(u,v,w)+\xi_{1}(t),\\\
\dot{v}&=g(u,v,w)+\xi_{2}(t),\\\ \dot{w}&=h(u,v,w)+\xi_{3}(t),\end{split}$ (3)

where $f(u,v,w)$, $g(u,,w)$ and $h(u,v,w)$ have the functional forms given in
Eq. 1, and $\xi_{i}(t),\,i\,=\,1,2,3$ are Gaussian white noises with zero mean
and correlation function is given by
$<\xi_{i}(t),\xi_{j}(t^{{}^{\prime}})>\,=\,\sigma\delta(t-t^{{}^{\prime}})\delta_{ij}$
for $i,j\,=\,1,2,3$. Here $\sigma$ represents the strength of the noise. To
simulate the dynamics of this noise-driven system, we numerically solve the
stochastic differential system by the explicit Euler-Maruyama scheme.

First we investigate if the boundedness of the system (3) is affected by the
presence of additive noise. Recall that the system, without noise, blows up as
the magnitude of the Allee parameter $\theta$ increases beyond a threshold
(cf. Fig. 1). Therefore it is important to examine if noise suppresses or
enhances the probability of unbounded growth in the system under Allee effect.

Fig. 10 displays representative results for the probability of blow-ups in the
population of vegetation, estimated for varying noise strengths $\sigma$, for
the Allee parameter $\theta\,=\,0.1$. Note that the system without noise (i.e.
$\sigma=0$) had significant probability of unbounded vegetation growth for
this value of $\theta$ (see Fig. 1). It is clearly evident from the results in
Fig. 10 that the probability of blow-ups for vegetation rapidly decreases to
zero with increasing magnitude of noise strength. So the presence of noise
helps to keep the populations bounded, indicating the constructive role of
noise in the stability of this three species system.

Next we examine how noise influences the extreme events which were observed to
emerge in this system in the presence of Allee effect. We first examine the
temporal evolution of the population densities and their corresponding phase-
space attractors. In Fig. 11 we display illustrative results of the times
series and the phase-space attractors of the system governed by Eqn. (3), for
different noise strength $\sigma$. It is observed that when the strength is
very low ($\sigma\sim 10^{-4}$), all population densities fluctuate in an
aperiodical manner and settle down to a chaotic attractor, as shown in Fig.
11a. Also extreme events occur in all populations for very low noise
strengths, as is clear from the figure where the vegetation, prey and predator
populations can be seen to cross the $\mu+7\sigma$ threshold. However with
increasing strength of noise, these extreme events disappear from all
populations in the system. Additionally, the populations are seen to fluctuate
in a more regular almost-periodic manner (see Fig. 11b). On further increase
of the noise strength all populations settle down to a quasi-fixed state, as
evident from Fig. 11c. This suggests that noise transforms the chaotic
behaviour of the system to a noisy fixed point. Importantly, the very long
time intervals we sampled did not yield a single extreme event. That is, under
increased noise strengths there is no evidence of extreme events any more, in
either the vegetation, prey or predator populations. Thus we arrive at the
following important conclusion: Noise leads to quasi-fixed (non-zero)
populations and the suppression of extreme events in this three species
system.

Further, in order to quantify how noise influences the emergence of extreme
events, we again estimate the probability of obtaining extreme events,
$P_{ext}$, under varying noise strength $\sigma$. The results are exhibited in
Fig. 12. It is clear from the figure that the probability of obtaining extreme
events is the lowest for vegetation and the highest for the predator
population. This is consistent with the observations for the system without
noise (see Fig. 8). The new significant result here is that the probability of
obtaining extreme events decreases to zero for increasing the noise strength
$\sigma$. Therefore in presence of sufficiently strong additive noise, extreme
events are suppressed in vegetation, prey and predator populations. This
points to the novel finding that stochasticity can lead to the mitigation of
extreme events.

Figure 10: Probability of unbounded vegetation growth in presence of additive
noise, with respect to noise strength $\sigma$. As in Fig. 1, a blow up to be
considered to occur when the vegetation population exceeds $10^{3}$. Here the
Allee parameter $\theta\,=\,0.1$, and the other system parameters are the same
as in Fig. 1. The probability is estimated from a sample of $10^{3}$ initial
states randomly distributed in a hyper-cube ($u\in[0:4],v\in[0:2],w\in[0:5]$)
in phase space.

Figure 11: Time series and phase dynamics of the system (3) with different
$\sigma$. (a) $\sigma\,=\,10^{-4}$, (b) $\sigma\,=\,10^{-2}$ and (c)
$\sigma\,=\,10^{-1}$. We keep all other parameters values are same as before
except $\theta\,=\,0.024$. Figure 12: Probability of obtaining extreme events
in unit time ($P_{ext}$) of the system (3) with respect to the noise strength
$\sigma$. Here we also consider a extreme event occurs when population level
cross the $\mu+7\sigma$ threshold. $P_{ext}$ for vegetation, prey and predator
populations are shown by blue, red and black colour respectively.

##  4 Discussion

In summary, we explored the dynamics of a three-species trophic system
incorporating the Allee Effect in the prey population. Our focus is on the
emergence of extreme events in the system. In particular we address the
significant question of whether or not Allee effect suppresses or enhances
extreme events. Our key observations are as follows: First, under Allee effect
the regular periodic dynamics changes to chaotic, as evident from the
emergence of chaotic attractors for increasing Allee parameter $\theta$.
Further, we find that the system exhibits unbounded growth in the vegetation
population (a “blow-up”) after a critical value of the Allee parameter. The
most significant result is the observation of a critical Allee parameter
beyond which the probability of obtaining extreme events becomes non-zero for
all three population densities. Though the emergence of extreme events in the
predator population is not affected much by the Allee effect, the prey
population shows a sharp increase in the probability of obtaining extreme
events after a threshold value of the Allee parameter $\theta$, and the
vegetation population also yields extreme events for sufficiently strong Allee
effect. An interesting open problem in this context would be to check the
observation that the extreme events in the predator population are more
pronounced than in prey and vegetation across other models, in order to
establish the generality of this important trend in a larger class of models.

Lastly we consider the influence of additive noise on extreme events. First,
we find that noise tames the unbounded vegetation growth induced by Allee
effect. More interestingly, we demonstrate that stochasticity drastically
diminishes the probability of extreme events in all three populations. In fact
for sufficiently high noise, we do not observe any more extreme events in the
system. This indicates that noise can mitigate extreme events, and has
potentially important impact on the observability of extreme events in
naturally occurring systems.

Conflict of Interest

The authors declare that they have no conflict of interest.

Availability of data and material

Data will be made available on reasonable request.

Authors’ contributions

SS conceived the problem, DS did all the simulations, DS and SS analyzed the
results and wrote the manuscript together.

## References

  * [1] Lubchenco, J. & Karl, T. R.  Extreme weather events.  _Phys. Today_ 65, 31 (2012).
  * [2] Kinney, R., Crucitti, P., Albert, R. & Latora, V.  Modeling cascading failures in the north american power grid.  _The European Physical Journal B-Condensed Matter and Complex Systems_ 46, 101–107 (2005).
  * [3] Strogatz, S.  How the blackout came to life.  _The New York Times_ (2003).
  * [4] Albeverio, S., Jentsch, V. & Kantz, H.  _Extreme events in nature and society_ (Springer Science & Business Media, 2006).
  * [5] Dysthe, K., Krogstad, H. E. & Müller, P.  Oceanic rogue waves.  _Annu. Rev. Fluid Mech._ 40, 287–310 (2008).
  * [6] Solli, D. R., Ropers, C., Koonath, P. & Jalali, B.  Optical rogue waves.  _Nature_ 450, 1054–1057 (2007).
  * [7] Lillo, F. & Mantegna, R. N.  Power-law relaxation in a complex system: Omori law after a financial market crash.  _Physical Review E_ 68, 016119 (2003).
  * [8] Smith, M.  An ecological perspective on extreme climatic events: a synthetic definition and framework to guide future research.  _Journal of Ecology_ 99, 656–663 (2011).
  * [9] Majumdar, S. N. & Ziff, R. M.  Universal record statistics of random walks and lévy flights.  _Physical review letters_ 101, 050601 (2008).
  * [10] Schehr, G. & Majumdar, S. N.  Universal order statistics of random walks.  _Physical review letters_ 108, 040601 (2012).
  * [11] Kishore, V., Santhanam, M. & Amritkar, R.  Extreme events on complex networks.  _Physical review letters_ 106, 188701 (2011).
  * [12] Balakrishnan, V., Nicolis, C. & Nicolis, G.  Extreme value distributions in chaotic dynamics.  _Journal of Statistical Physics_ 80, 307–336 (1995).
  * [13] Nicolis, C., Balakrishnan, V. & Nicolis, G.  Extreme events in deterministic dynamical systems.  _Physical Review Letters_ 97, 210602 (2006).
  * [14] Ansmann, G., Karnatak, R., Lehnertz, K. & Feudel, U.  Extreme events in excitable systems and mechanisms of their generation.  _Physical Review E_ 88, 052911 (2013).
  * [15] Karnatak, R., Ansmann, G., Feudel, U. & Lehnertz, K.  Route to extreme events in excitable systems.  _Physical Review E_ 90, 022917 (2014).
  * [16] Kingston, S. L., Thamilmaran, K., Pal, P., Feudel, U. & Dana, S. K.  Extreme events in the forced liénard system.  _Physical Review E_ 96, 052204 (2017).
  * [17] Moitra, P. & Sinha, S.  Emergence of extreme events in networks of parametrically coupled chaotic populations.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 29, 023131 (2019).
  * [18] Premraj, D. _et al._ Dragon-king extreme events as precursors for catastrophic transition.  _Europhysics Letters_ 134, 34006 (2021).
  * [19] Courchamp, F., Berec, L. & Gascoigne, J.  _Allee effects in ecology and conservation_ (Oxford University Press, 2008).
  * [20] Sen, D., Ghorai, S., Sharma, S. & Banerjee, M.  Allee effect in prey’s growth reduces the dynamical complexity in prey-predator model with generalist predator.  _Applied Mathematical Modelling_ 91, 768–790 (2021).
  * [21] Dennis, B.  Allee effects: population growth, critical density, and the chance of extinction.  _Natural Resource Modeling_ 3, 481–538 (1989).
  * [22] Blasius, B., Huppert, A. & Stone, L.  Complex dynamics and phase synchronization in spatially extended ecological systems.  _Nature_ 399, 354–359 (1999).
  * [23] Sinha, S. & Sinha, S.  Evidence of universality for the may-wigner stability theorem for random networks with local dynamics.  _Physical Review E (Rapid Communication)_ 71, 020902 (2005).
  * [24] Sinha, S. & Sinha, S.  Robust emergent activity in dynamical networks.  _Physical Review E_ 74, 066117 (2006).
  * [25] Choudhary, A. & Sinha, S.  Balance of interactions determines optimal survival in multi-species communities.  _PLoS ONE_ (2015).
  * [26] Chaurasia, S. S., Verma, U. K. & Sinha, S.  Advent of extreme events in predator populations.  _Scientific Reports_ 10, 1–10 (2020).
  * [27] Boukal, D. S. & Berec, L.  Single-species models of the allee effect: extinction boundaries, sex ratios and mate encounters.  _Journal of Theoretical Biology_ 218, 375–394 (2002).
  * [28] Rowe, S., Hutchings, J. A., Bekkevold, D. & Rakitin, A.  Depensation, probability of fertilization, and the mating system of atlantic cod (gadus morhua l.).  _ICES Journal of Marine Science_ 61, 1144–1150 (2004).
  * [29] McCarthy, M.  The allee effect, finding mates and theoretical models.  _Ecological Modelling_ 103, 99–102 (1997).
  * [30] Scheuring, I.  Allee effect increases the dynamical stability of populations.  _Journal of Theoretical Biology_ 199, 407–414 (1999).
  * [31] Choudhary, A., Kohar, V. & Sinha, S.  Taming explosive growth through dynamic random links.  _Scientific Reports (Nature)_ 4308 (2014).
  * [32] Roth, G. & J. Schreiber, J.  Pushed beyond the brink: Allee effects, environmental stochasticity, and extinction.  _Journal of Biological Dynamics_ 8, 187–205 (2014).


∎

11institutetext: Eric Okyere 22institutetext: Department of Basic Sciences,
University of Health and Allied Sciences, PMB 31, Ho, Volta Region, Ghana
Tel.: +233-508389028
22email<EMAIL_ADDRESS>33institutetext: F. T. Oduro, S. K. Amponsah, I.
K. Dontwi, N. K. Frempong44institutetext: Department of Mathematics, Kwame
Nkrumah University of Science and Technology, Kumasi, Ghana

# Fractional Order Malaria Model With Temporary Immunity

E. Okyere  F. T. Oduro  S. K. Amponsah  I. K. Dontwi  N. K. Frempong

(March, 2016)

###### Abstract

In this paper, we propose and analyze a new fractional order malaria model
with temporary immunity. We extend the existing malaria model to include
fractional derivatives. Basic reproduction number is computed using the next
generation matrix method. Model equilibria are determined and their local
asymptotic stability analysis are considered. An efficient Adams-type
predictor-corrector method is used to solve the initial value problem.
Finally, phase plane portraits are generated to describe the dynamics between
susceptible fractions of human hosts and infected fractions of mosquito
vectors.

###### Keywords:

Fractional calculus Predictor-corrector method Malaria model

##  1 Introduction

Malaria is one of the deadly vector-borne diseases affecting the developing
part of the world and mostly in Africa. The reports by the World Health
Organization (see, e.g, WHO, 2005, 2012, 2014) indicates that, malaria is a
serious threat to human life and remains a dangerous infectious disease. Its
economic burden on the affected countries are extremely huge and needs serious
public health attention. Since its discovery, many researchers have studied
the disease dynamics mathematically to capture and understand the complex
interaction between infected mosquitoes and the susceptible human population.

Mathematical modeling of infectious diseases has played a key role in
understanding the dynamical processes of disease transmission and control
strategies. The first malaria epidemiological model was constructed and
formulated by Ross (1911). He used deterministic compartmental models to
describe malaria infection dynamics between human and mosquito populations.
Ngwa and Shu (2000) proposed a malaria model and considered global stability
analysis for the disease free equilibrium point. Mandal et al. (2011)
considered an extensive survey and review on modeling of malaria infection.
Backward bifurcation analysis is very important in mathematical modeling of
infectious diseases and has been investigated by many authors in malaria
infection models (see, e.g, Keegan and Dushoff, 2013, Ghosh et al., 2013,
Chiyaka et al., 2008, Chitnis, 2005, Chitnis et al., 2006, Ngonghala et al.,
2012). Águas et al. (2012) developed an age structure model to study the
effects of relapse in malaria parasites. A periodic model based on the Ross-
MacDonald model was considered by Gao et al. (2014). They studied the effects
of spatial heterogeneity and temporal on disease dynamics. Their proposed
model incorporated periodic behavior in mosquito ecology and seasonal human
movement to capture variations of malaria spread among different regions.

Tumwiine et al. (2010) proposed and analyzed a malaria infection model with
infective immigrants. Tumwiine et al. (2007b) formulated an epidemic model
with standard incidence for the dynamics of malaria infection. They applied
additive compound matrices approach to show global stability analysis of the
endemic equilibrium. Using the model in Tumwiine et al. (2007b), the same
authors (Tumwiine et al., 2007a) studied the oscillatory behavior of malaria
disease in a population with temporary immunity. Their numerical results shows
that the endemic equilibrium converges to a steady state. Okosun and Makinde
(2011) developed a deterministic compartmental model to investigate the impact
of drug resistance in malaria infection and the same authors (Okosun and
Makinde, 2014) formulated an optimal control model for malaria and cholera co-
infection.

Abdullahi et al. (2013) have assessed the effectiveness of drugs in a malaria
transmission model. Barley et al. (2012), Mukandavire et al. (2009) and Abu-
Raddad et al. (2006) formulated a model to study infectious disease co-
infection of HIV and malaria. Okosun and Makinde (2013) studied the spread of
malaria infection with an optimal control model. Pongtavornpinyo et al.
(2008), they studied the spread of anti-malaria resistance. Lawi et al. (2011)
studied co-infection model for malaria and meningitis among children. They
observed that, when the threshold parameter $R_{0}<1$, the disease-free
equilibrium might exhibit some instabilities in the global sense. They further
applied centre manifold theorem to investigate asymptotic stability of the
endemic equilibrium. Antia et al. (2008) studied the dynamics of acute malaria
transmission using an age-structured models. A mathematical model for malaria
infection which incorporates weather has been formulated and studied by Hoshen
and Morse (2004). Koella and Anita (2003) developed an epidemic model to
explore the incidence of anti-malarial resistance.

Mathematical modeling of infectious diseases dynamics and other important
areas of studies such as economics, finance and engineering has extensively
been explored using the theory and applications of the classical differential
equations. But in recent times, the theory and applications of fractional
calculus has become extremely useful and important in modeling of biological
processes and other areas of studies due to the memory property of fractional
derivatives. Many authors have contributed significantly in compartmental
modeling of infectious disease dynamics using fractional differential
equations. Sardar et al. (2015) have developed a mosquito-transmitted disease
model with fractional differential equations. Pinto and Machado (2013)
described fractional order models for malaria infection. They modified the
integer order compartmental model formulated by Chiyaka et al. (2008) to
numerically study control strategies in malaria dynamics. Fractional order
models for HIV and TB co-infection have been analysed by Pinto and Carvalho
(2014). Fractional order models for HIV infection have been considered and
analysed by several authors (see, e.g, Ye and Ding, 2009, Gokdogan et al.,
2011, Liu et al., 2014, Arafa et al., 2014, Ding and Ye, 2009, Ding et al.,
2012, Liu and Lu, 2014). Liu and Chen (2015) constructed fractional order
competition model for love triangle. Ahmad and El-Khazali (2007) proposed a
mathematical model for love using fractional order systems. El-Saka (2015)
investigated backward bifurcation in a fractional order vaccination model.

Javidi and Ahmad (2014) developed fractional order model for cholera
infection. Arafa et al. (2012) applied homotopy analysis method and fourth
order Runge-Kutta method to generate analytical and numerical results for
fractional order childhood disease models. An epidemic model with fractional
order for influenza A have been formulated and analysed (see, e.g, El-Shahed
and Alsaedi, 2011, González-Parra et al., 2014). Rihan et al. (2014) studied
deterministic fractional order SIRC model with Salmonella bacteria infection.
Ahmed and El-Saka (2010) presented and discussed fractional order systems for
Hepatitis C virus. The deadly ebola disease which killed a lot of people in
some part of Africa has been modelled by a system of fractional order
derivatives (Area et al., 2015). An endemic model with constant population has
been analytically and numerically studied by Okyere et al. (2016) using Caputo
fractional derivatives.

In this work, we formulate a fractional order model for malaria infection with
temporary immunity. We consider the malaria disease model studied and analysed
by Tumwiine et al. (2007a, b). The authors developed their model using
deterministic integer order differential equations. Our new model will be
constructed using fractional order derivatives. We are motivated by this
method of mathematical formulation in epidemiological modeling due to the
effective nature of fractional derivatives.

The rest of the article is as follows. In section 2, we formulate the
fractional order malaria model. Section 3, deals with non-negativity of model
solutions. In section 4, we compute model equilibria and the basic
reproduction number and then investigate local asymptotic stability analysis.
We consider numerical simulations in section 5. In section 6, we discuss
simulation results and finally conclude the paper in section 7.

##  2 Model Derivation

In this section, we extend and modify the integer order malaria model by
Tumwiine et al. (2007a, b) to become fractional order malaria model. The
theory and applications of fractional calculus has become extremely useful and
important in modeling of biological processes and other areas of studies due
to the memory property of fractional derivatives. Many authors have
contributed significantly in compartmental modeling of infectious disease
dynamics using fractional differential equations. There are several
definitions of fractional derivatives and integrals in fractional calculus
(see, e.g, Petras, 2011, Podlubny, 1999), but in this work, our model
construction will be based on Caputo derivative. Our motivation for this type
of fractional derivative is based on the fact that, it has an advantage on
initial value problems.

###### Definition 1

(Podlubny, 1999, Petras, 2011) Fractional integral of order $\alpha$ is
defined as

$I^{\alpha}g(t)=\frac{1}{\Gamma(\alpha)}\int_{0}^{t}\\!\frac{g(x)}{(t-x)^{1-\alpha}}\,\mathrm{d}x$

for $0<\alpha<1,\ t>0.$

###### Definition 2

(Podlubny, 1999, Petras, 2011) Caputo fractional derivative is defined as

$D^{\alpha}g(t)=\frac{1}{\Gamma(k-\alpha)}\int_{0}^{t}\\!\frac{g^{k}(x)}{(t-x)^{\alpha+1-k}}\,\mathrm{d}x.$

for $k-1<\alpha<k.$

In this article, we consider the scaled malaria model model presented by
Tumwiine et al. (2007a, b). The authors, transformed their compartmental
malaria model ($S_{H}I_{H}R_{H}-S_{V}I_{V}$) into systems of proportions
($s_{h}i_{h}r_{h}-s_{v}i_{v}$) for the human and mosquito populations describe
by systems of ordinary differential equations.

$\begin{array}[]{lllll}\displaystyle\frac{ds_{h}}{dt}&=&\displaystyle{\lambda_{h}(1-s_{h})-abms_{h}i_{v}+\nu
i_{h}+\gamma r_{h}+\delta s_{h}i_{h}},\\\\[15.0pt]
\displaystyle\frac{di_{h}}{dt}&=&\displaystyle{abms_{h}i_{v}-(\nu+r+\lambda_{h}+\delta)i_{h}+\delta
i^{2}_{h}},\\\\[15.0pt]
\displaystyle\frac{d{r_{h}}}{dt}&=&\displaystyle{ri_{h}-(\gamma+\lambda_{h})r_{h}+\delta
i_{h}r_{h}},\\\\[15.0pt]
\displaystyle\frac{ds_{h}}{dt}&=&\displaystyle{\lambda_{v}(1-s_{v})-aci_{h}s_{v}},\\\\[15.0pt]
\displaystyle\frac{d{i_{v}}}{dt}&=&\displaystyle{acs_{v}i_{h}-\lambda_{v}i_{v}},\end{array}$
(1)

with $s_{h}(t)+i_{h}(t)+r_{h}(t)=1$ and $s_{v}(t)+i_{v}(t)=1$,

where

•

$s_{h}$: proportion of susceptible human hosts at time $t$

•

$i_{h}$: proportion of infected human hosts at time $t$

•

$r_{h}$: proportion of immune human hosts at time $t$

•

$s_{v}$: proportion susceptible mosquito vectors at time $t$

•

$i_{v}$: proportion of infected mosquito vectors at time $t$

•

$a$: average daily biting rate on man by a single mosquito (infection rate)

•

$b$: proportion of bites on man that produce an infection

•

$\nu$ recovery rate of human hosts from the disease

•

$c$: probability that a mosquito becomes infectious

•

$\gamma$: per capita rate of loss of immunity in human hosts

•

$r$: rate at which human hosts acquire immunity

•

$\delta$: per capita death rate of infected human hosts due to the disease

•

$m=\frac{N_{V}}{N_{H}}$: number of female mosquitoes per human host

•

$N_{H}$: total human population

•

$N_{V}$: total mosquito population

•

$\mu_{h}$: per capita natural death rate of humans

•

$\mu_{v}$: per capita natural death rate of mosquitoes

•

$\lambda_{h}$: per capita natural birth rate of humans

•

$\lambda_{v}$: per capita natural birth rate of mosquitoes

The new model we consider in this paper is the system of non-linear fractional
differential equations in the sense of Caputo fractional derivatives.

$\begin{array}[]{lllll}\displaystyle
D^{\alpha}s_{h}&=&\displaystyle{\lambda^{\alpha}_{h}(1-s_{h})-a^{\alpha}bms_{h}i_{v}+\nu^{\alpha}i_{h}+\gamma^{\alpha}r_{h}+\delta^{\alpha}s_{h}i_{h}},\\\\[15.0pt]
\displaystyle
D^{\alpha}i_{h}&=&\displaystyle{a^{\alpha}bms_{h}i_{v}-(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})i_{h}+\delta^{\alpha}i^{2}_{h}},\\\\[15.0pt]
\displaystyle
D^{\alpha}r_{h}&=&\displaystyle{r^{\alpha}i_{h}-(\gamma^{\alpha}+\lambda^{\alpha}_{h})r_{h}+\delta^{\alpha}i_{h}r_{h}},\\\\[15.0pt]
\displaystyle
D^{\alpha}s_{v}&=&\displaystyle{\lambda^{\alpha}_{v}(1-s_{v})-a^{\alpha}ci_{h}s_{v}},\\\\[15.0pt]
\displaystyle
D^{\alpha}i_{v}&=&\displaystyle{a^{\alpha}cs_{v}i_{h}-\lambda^{\alpha}_{v}i_{v}},\end{array}$
(2)

with $s_{h}(t)+i_{h}(t)+r_{h}(t)=1$ and $s_{v}(t)+i_{v}(t)=1$.

It is important to remark that when the fractional order $\alpha\rightarrow
1$, the fractional order malaria model (2) with the model restrictions
$s_{h}(t)+i_{h}(t)+r_{h}(t)=1$ and $s_{v}(t)+i_{v}(t)=1$, becomes the integer
order malaria model (1) studied by Tumwiine et al. (2007a, b).

##  3 Non-negative Solutions

Let $\mathbf{R^{5}_{+}}=\\{X\in\mathbf{R^{5}}\mathrel{\mathop{\mathchar
58\relax}}X\geq 0\\}$, where $X=(s_{h},i_{h},r_{h},s_{v},i_{v})^{T}$. We apply
the following Lemma in (Odibat and Shawagfeh, 2007) to show the theorem about
the non-negative solutions of the model.

###### Lemma 1

(Odibat and Shawagfeh, 2007) (Generalized Mean Value Theorem).
Suppose that $w(x)\in C[a,b]$ and $D^{\alpha}w(x)\in C(a,b]$ for $0<\alpha\leq
1$, then we have

$w(x)=w(a)+\frac{1}{\Gamma(\alpha)}D^{\alpha}w(\xi)(x-a)^{\alpha}$

with $a\leq\xi\leq x,\ \forall x\in(a,b].$

###### Remark 1

Assume that $w(x)\in C[a,b]$ and $D^{\alpha}w(x)\in C(a,b]$, for $0<\alpha\leq
1$. It follows from Lemma 1 that if $D^{\alpha}w(x)\geq 0,\forall x\in(a,b)$
then $w(x)$ is non-decreasing $\forall x\in[a,b]$ and if $D^{\alpha}w(x)\leq
0,\forall x\in(a,b)$, then $w(x)$ is non-increasing $\forall x\in[a,b].$

###### Theorem 3.1

The malaria disease model (2) has a unique solution and it remains in
$\mathbf{R^{5}_{+}}$.

###### Proof

By applying Theorem $3.1$ and Remark $3.2$ in (Lin, 2007), the existence and
uniqueness of the solution of equation (2) in $(0,\infty)$ follows. The domain
$\mathbf{R^{5}_{+}}$ for the model problem is positively invariant, since

$\begin{array}[]{lllll}\displaystyle
D^{\alpha}s_{h}|_{s_{h}=0}&=&\displaystyle{\lambda^{\alpha}_{h}+\nu^{\alpha}i_{h}+\gamma^{\alpha}r_{h}\geq
0},\\\\[15.0pt] \displaystyle
D^{\alpha}i_{h}|_{i_{h}=0}&=&\displaystyle{a^{\alpha}bms_{h}i_{v}\geq
0},\\\\[15.0pt] \displaystyle
D^{\alpha}r_{h}|_{r_{h}=0}&=&\displaystyle{r^{\alpha}i_{h}\geq 0},\\\\[15.0pt]
\displaystyle
D^{\alpha}s_{v}|_{s_{v}=0}&=&\displaystyle{\lambda^{\alpha}_{v}\geq
0},\\\\[15.0pt] \displaystyle
D^{\alpha}i_{v}|_{i_{v}=0}&=&\displaystyle{a^{\alpha}cs_{v}i_{h}\geq
0},\end{array}$ (3)

on each hyperplane bounding the non-negative orthant, the vector filed points
into $\mathbf{R^{5}_{+}}$.

##  4 Analysis of Model Equilibria

In order to determine the model equilibria, we proceed this way.

Let

$\begin{cases}D^{\alpha}s_{h}=0&\\\ D^{\alpha}i_{h}=0&\\\
D^{\alpha}r_{h}=0&\\\ D^{\alpha}s_{v}=0&\\\ D^{\alpha}i_{v}=0\end{cases}$ (4)

By solving equation (4) and ignoring the steps involved, the disease-free
equilibrium is given by $F_{o}=(1,0,0,1,0)$ and the endemic equilibrium as
$F_{*}=(s^{*}_{h},i^{*}_{h},r^{*}_{h},s^{*}_{v},i^{*}_{v})$, where
$s^{*}_{h},r^{*}_{h},s^{*}_{v}$ and $i^{*}_{v}$ are expressed in terms
$i^{*}_{h}$ as follows:

$\begin{array}[]{lllll}\displaystyle
s^{*}_{h}&=&\displaystyle{\dfrac{(\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h})[(\lambda^{\alpha}_{h}+\gamma^{\alpha}-\delta^{\alpha}i^{*}_{h})(\lambda^{\alpha}_{h}+\nu
i^{*}_{h})+\gamma^{\alpha}r^{\alpha}i^{*}_{h}]}{(\lambda^{\alpha}_{h}+\gamma^{\alpha}-\delta^{\alpha}i^{*}_{h})[(\lambda^{\alpha}_{h}-\delta^{\alpha}i^{*}_{h})(\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h})+a^{2\alpha}bmci^{*}_{h}]}},\\\\[15.0pt]
\displaystyle
r^{*}_{h}&=&\displaystyle{\dfrac{r^{\alpha}i^{*}_{h}}{(\lambda^{\alpha}_{h}+\gamma^{\alpha}-\delta^{\alpha}i^{*}_{h})}},\\\\[15.0pt]
\displaystyle
s^{*}_{v}&=&\displaystyle{\dfrac{\lambda^{\alpha}_{v}}{(\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h})}},\\\\[15.0pt]
\displaystyle
i^{*}_{v}&=&\displaystyle{\dfrac{a^{\alpha}ci^{*}_{h}}{(\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h})}},\end{array}$
(5)

Next, we compute the basic reproduction number $R_{o}$ using the next
generation method on the fractional order malaria model (2). By using the
method and the notations by Driessche and Watmough (2002), the special
matrices $F$( for new infections terms) and $V$ (for the remaining transitions
terms) related to the model problem (2) are respectfully given by

$F=\begin{bmatrix}0&&&a^{\alpha}bm\\\ \\\ a^{\alpha}c&&&0\par\end{bmatrix}$
(6)

and

$V=\begin{bmatrix}\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha}&&&0\\\
\\\ 0&&&\lambda^{\alpha}_{v}\par\end{bmatrix}$ (7)

The basic reproduction number is the spectral radius of the generation matrix
$F{V^{-1}}$, where

$F{V^{-1}}=\begin{bmatrix}0&&&\frac{a^{\alpha}bm}{\lambda^{\alpha}_{v}}\\\ \\\
\dfrac{a^{\alpha}c}{\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha}}&&&0\par\end{bmatrix}$
(8)

Therefore we have

$R_{o}=\sqrt{\dfrac{a^{2\alpha}bmc}{{\lambda^{\alpha}_{v}}(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})}}$

###  4.1 Local Asymptotic Stability Analysis of Disease-Free Equilibrium

###### Theorem 4.1

The equilibrium $F_{o}$ of the malaria disease model (2) is locally
asymptotically stable if $R_{o}<1$ and unstable if $R_{o}>1.$

###### Proof

The Jacobian matrix $J(F_{o})$ for the malaria disease model (2) evaluated at
$F_{o}$ is given by

$J(F_{o})=\begin{bmatrix}-\lambda^{\alpha}_{h}&\nu^{\alpha}+\delta^{\alpha}&\gamma^{\alpha}&0&-a^{\alpha}bm\\\
0&-(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})&0&0&a^{\alpha}bm\\\
0&r^{\alpha}&-(\lambda^{\alpha}_{h}+\gamma^{\alpha})&0&0\\\
0&-a^{\alpha}c&0&-\lambda{\alpha}_{v}&0\\\
0&a^{\alpha}c&0&0&-\lambda^{\alpha}_{v}\\\ \end{bmatrix}$ (9)

The disease-free equilibrium point $F_{o}$ is locally asymptotically if all
the eigenvalues $\lambda_{i},i=1,2,3,4,5$ of $J(F_{o})$ satisfy the following
condition [Ahmed et al. (2007), Matignon (1996)]: $\mathinner{\\!\left\lvert
arg(\lambda_{i})\right\rvert}>\frac{\alpha\pi}{2}$.

From the square matrix $J(F_{o})$, it is obvious that, three of the
eigenvalues are given by
$-\lambda^{\alpha}_{h},-(\lambda^{\alpha}_{h}+\gamma^{\alpha})$ and
$-\lambda^{\alpha}_{v}$. It can be seen that all the three eigenvalues are
negative and hence satisfy the stability condition: $\mathinner{\\!\left\lvert
arg(\lambda_{i})\right\rvert}>\frac{\alpha\pi}{2}$.

The rest of the eigenvalues can determined from the sub-matrix given as:

$B=\begin{bmatrix}-(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})&&&a^{\alpha}bm\\\
\\\ a^{\alpha}c&&&-\lambda^{\alpha}_{v}\par\end{bmatrix}$ (10)

The trace of matrix $B$ is given by

$tr(B)=-(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha}+\lambda^{\alpha}_{v})<0$

The determinant of $B$ yields :

$\displaystyle detB$
$\displaystyle=\lambda^{\alpha}_{v}(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})-a^{2\alpha}bmc$
$\displaystyle=\lambda^{\alpha}_{v}(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})\left[1-\dfrac{a^{2\alpha}bmc}{\lambda^{\alpha}_{v}(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})}\
\right]$
$\displaystyle=\lambda^{\alpha}_{v}(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})\left[1-R^{2}_{o}\right]$

Therefore $detB>0$ if $R_{0}<1$.

Since $detB>0$ if $R_{0}<1$ and $tr(B)<0$, it follows that the two remaining
eigenvalues will have negative real parts and we can therefore conclude that,
the disease equilibrium is locally asymptotically stable, otherwise unstable.

###  4.2 Local Asymptotic Stability Analysis of Endemic Equilibrium

In this subsection, we discuss the local stability of the endemic equilibrium
point. Our analysis is based on the construction by Ahmed et al. (2006, 2007).
From the model constraints $s_{v}+i_{v}=1$ and $s_{h}+i_{h}+r_{h}=1$, we can
reduce the dimension of the malaria model in order to simplify and reduce the
complexity of the analysis.

As in Tumwiine et al. (2007a, b), let $s_{v}=1-i_{v}$ and
$r_{h}=1-s_{h}-i_{h}$ then the 5-dimensional fractional order model (2)
reduces to a 3-dimensional problem with Jacobian matrix evaluated the at
endemic equilibrium $F_{*}$ as:

$J(F_{*})=\begin{bmatrix}-(\lambda^{\alpha}_{h}+\gamma^{\alpha}+a^{\alpha}bmi^{*}_{v}-\delta^{\alpha}i^{*}_{h})&\nu^{\alpha}+\delta^{\alpha}s^{*}_{h}-\gamma^{\alpha}&-a^{\alpha}bms^{*}_{h}\\\
\\\
a^{\alpha}bmi^{*}_{v}&-(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha}-2\delta^{\alpha}i^{*}_{h})&a^{\alpha}bms^{*}_{h}\\\
\\\
0&a^{\alpha}c(1-i^{*}_{v})&-(\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h})\par\end{bmatrix}$

The characteristics equation of the Jacobian matrix, $J(F_{*})$ is

$\lambda^{3}+b_{1}\lambda^{2}+b_{2}\lambda+b3=0$ (11)

where

$\begin{split}b_{1}=&\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h}+\lambda^{\alpha}_{h}+\gamma^{\alpha}_{v}+a^{\alpha}bmi^{*}_{v}-\delta^{\alpha}i^{*}_{h}+\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha}-2\delta^{\alpha}i^{*}_{h}\\\
\\\
b_{2}=&(\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h})(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha}-2\delta^{\alpha}i^{*}_{h})-a^{\alpha}bmc(1-i^{*}_{v})s^{*}_{h}+(\gamma^{\alpha}-\nu^{\alpha}-\delta^{\alpha}s^{*}_{h})a^{\alpha}bmi^{*}_{v}\\\
&+(\lambda^{\alpha}_{h}+\gamma^{\alpha}+a^{\alpha}bmi^{*}_{v}-\delta^{\alpha}i^{*}_{h})(\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h}+\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha}-2\delta^{\alpha}i^{*}_{h})\\\
\\\
b_{3}=&(\lambda^{\alpha}_{h}+\gamma^{\alpha}+a^{\alpha}bmi^{*}_{v}-\delta^{\alpha}i^{*}_{h})[(\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h})(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha}-2\delta^{\alpha}i^{*}_{h})-a^{2\alpha}bmc(1-i^{*}_{v})s^{*}_{h}]\\\
&+a^{3\alpha}b^{2}mc(1-i^{*}_{v})s^{*}_{h}+a^{\alpha}bmi^{*}_{v}(\lambda^{\alpha}_{v}+a^{\alpha}ci^{*}_{h})(\gamma^{\alpha}-\nu^{\alpha}-\delta^{\alpha}s^{*}_{h})\end{split}$

Let $D(g)$ represent the discriminant of the polynomial function

$g(x)=x^{3}+b_{1}x^{2}+b_{2}x+b3$

where

$D(g)=-\left|\begin{array}[]{cccccccc}1&&b_{1}&&b_{2}&b_{3}&&0\\\
0&&1&&b_{1}&b_{2}&&b_{3}\\\ 3&&2b_{1}&&b_{2}&0&&0\\\ 0&&3&&2b_{1}&b_{2}&&0\\\
0&&0&&3&2b_{1}&&b_{2}\end{array}\right|=18b_{1}b_{2}b_{3}+(b_{1}b_{2})^{2}-4b_{3}b^{3}_{1}-4b^{3}_{2}-27b^{2}_{3}$

By following the results in (Ahmed et al., 2006, 2007) and similar
applications of their method in (El-Shahed and Alsaedi, 2011, Ozalp and
Demirci, 2011, Al-Sulami et al., 2014), we obtain the following proposition.

###### Proposition 1

One assume that $F_{*}$ exist in $\mathbf{R^{3}_{+}}$.

(i) If the discriminant of $g(x)$, $D(g)$ is positive and Routh-Hurwitz
conditions are satisfied , i.e. $D(g)>0,\ b_{1}>0,\ b_{3}>0$, then the endemic
equilibrium $F_{*}$ locally asymptotically stable.

(ii) If $D(g)<0,\ b_{1}\geq 0,\ b_{2}\geq 0,\ b_{3}>0,\ \alpha<2/3$, then the
endemic equilibrium $F_{*}$ locally asymptotically stable.

(iii)If $D(g)<0,\ b_{1}>0,\ b_{2}>0,\ b_{1}b_{2}=b_{3},\ \alpha\in(0,1),$ then
the endemic equilibrium $F_{*}$ locally asymptotically stable.

(iv) If $D(g)<0,\ b_{1}<0,\ b_{2}<0,\ \alpha>2/3$, then the endemic
equilibrium $F_{*}$ is unstable.

##  5 Numerical Simulations

Analytical solutions of nonlinear coupled fractional order differential
equations are as difficult as the classical integer order system of nonlinear
equations. Therefore, there is need to consider numerical discretizations and
approximations for such systems or models. Over the past years, several
numerical schemes have been developed to solve fractional order models. One of
the powerful approximation methods that has been considered by many authors
(see, e.g, El-Sayed et al., 2007, Ozalp and Demirci, 2011, Al-Sulami et al.,
2014, Cho et al., 2015, Okyere et al., 2016) is the efficient Adams-type
predictor-corrector method developed by Diethelm et al. (2002, 2004). Li and
Tao (2009) have considered a rigorous error analysis of this numerical scheme.

For our numerical approximations, we re-define the variables in the model
problem (2) as:

Let $s_{h}=X,\ i_{h}=Y,\ r_{h}=Z,\ s_{v}=S,\ i_{v}=I.$

Therefore following the approach in (Diethelm et al., 2002, 2004, Cho et al.,
2015, Al-Sulami et al., 2014) and ignoring the details in the derivation of
the scheme, we obtain the discretized system corresponding to the fractional
order malaria model (2)as follows:

$\begin{split}X_{k+1}=&X_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}a_{{}_{j,k+1}}\biggl{[}\lambda^{\alpha}_{h}(1-X_{j})-a^{\alpha}bmX_{j}I_{j}+\nu^{\alpha}Y_{j}+\gamma^{\alpha}Z_{j}+\delta^{\alpha}X_{j}Y_{j}\biggr{]}\\\
+&\frac{1}{\Gamma(\alpha)}\biggl{[}\lambda^{\alpha}_{h}(1-X^{p}_{k+1})-a^{\alpha}bmX^{p}_{k+1}I^{p}_{k+1}+\nu^{\alpha}Y^{p}_{k+1}+\gamma^{\alpha}Z^{p}_{k+1}+\delta^{\alpha}X^{p}_{k+1}Y^{p}_{k+1}\biggr{]}\\\
\\\
Y_{k+1}=&Y_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}a_{{}_{j,k+1}}\biggl{[}a^{\alpha}bmX_{j}I_{j}-(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})Y_{j}-\delta^{\alpha}Y^{2}_{j}\biggr{]}\\\
+&\frac{1}{\Gamma(\alpha)}\biggl{[}a^{\alpha}bmX^{p}_{k+1}I^{p}_{k+1}-(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})Y^{p}_{k+1}+\delta^{\alpha}(Y^{P}_{K+1})^{2}\biggr{]}\\\
\\\
Z_{k+1}=&Z_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}a_{{}_{j,k+1}}\biggl{[}rY_{j}-(\gamma^{\alpha}+\lambda^{\alpha}_{h})Z_{j}+\delta
Y_{j}Z_{j}\biggr{]}+\frac{1}{\Gamma(\alpha)}\biggl{[}rY^{p}_{k+1}-(\gamma^{\alpha}+\lambda^{\alpha}_{h})Z^{p}_{k+1}+\delta
Y^{p}_{k+1}Z^{p}_{k+1}\biggr{]}\\\ \\\
S_{k+1}=&S_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}a_{{}_{j,k+1}}\biggl{[}\lambda^{\alpha}_{v}(1-S_{J})-a^{\alpha}cY_{j}S_{j}\biggr{]}+\frac{1}{\Gamma(\alpha)}\biggl{[}\lambda^{\alpha}_{v}(1-S^{p}_{k+1})-a^{\alpha}cY^{p}_{k+1}S^{p}_{k+1}\biggl{]}\\\
\\\
I_{k+1}=&I_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}a_{{}_{j,k+1}}\biggl{[}a^{\alpha}cS_{j}Y_{j}-\lambda^{\alpha}_{v}I_{j}\biggr{]}+\frac{1}{\Gamma(\alpha)}\biggl{[}a^{\alpha}cS^{p}_{k+1}Y^{p}_{k+1}-\lambda^{\alpha}_{v}I^{p}_{k+1}\biggl{]}\end{split}$

where

$\begin{split}X^{p}_{k+1}=&X_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}b_{{}_{j,k+1}}\biggl{[}\lambda^{\alpha}_{h}(1-X_{j})-a^{\alpha}bmX_{j}I_{j}+\nu^{\alpha}Y_{j}+\gamma^{\alpha}Z_{j}+\delta^{\alpha}X_{j}Y_{j}\biggr{]}\\\
\\\
Y^{p}_{k+1}=&Y_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}b_{{}_{j,k+1}}\biggl{[}a^{\alpha}bmX_{j}I_{j}-(\nu^{\alpha}+r^{\alpha}+\lambda^{\alpha}_{h}+\delta^{\alpha})Y_{j}-\delta^{\alpha}Y^{2}_{j}\biggr{]}\\\
\\\
Z^{p}_{k+1}=&Z_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}b_{{}_{j,k+1}}\biggl{[}rY_{j}-(\gamma^{\alpha}+\lambda^{\alpha}_{h})Z_{j}+\delta
Y_{j}Z_{j}\biggr{]}\\\ \\\
S^{p}_{k+1}=&S_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}b_{{}_{j,k+1}}\biggl{[}\lambda^{\alpha}_{v}(1-S_{J})-a^{\alpha}cY_{j}S_{j}\biggr{]}\\\
\\\
I^{p}_{k+1}=&I_{0}+\frac{1}{\Gamma(\alpha)}\sum_{j=0}^{k}b_{{}_{j,k+1}}\biggl{[}a^{\alpha}cS_{j}Y_{j}-\lambda^{\alpha}_{v}I_{j}\biggr{]}\end{split}$

$a_{{}_{j,k+1}}=\dfrac{h^{\alpha}}{\alpha(\alpha+1)}\begin{cases}k^{\alpha+1}-(k-\alpha)(k+1)^{\alpha}&\quad\text{if
}j=0\\\ (k-j+2)^{\alpha+1}+(k-j)^{\alpha+1}-2(k-j+1)^{\alpha+1}&\quad\text{if
}1\leq j\leq k\\\ 1&\quad\text{if }j=k+1\end{cases}$ and

$b_{{}_{j,k+1}}=\frac{h^{\alpha}}{\alpha}\left[(k-j+1)^{\alpha}-(k-j)^{\alpha}\right],\
\ \ 0\leq j\leq k$

We then implement our numerical scheme in Matlab and the simulation results
are discussed in the next section.

##  6 Results and Discussion

In this section, we discuss the endemic trajectories of our numerical
solutions with model parameter values adapted from the work by Tumwiine et al.
(2007a) and basic reproduction number $R_{o}=1.5$. The simulated solutions for
the scaled human population $s_{h}(t),\ i_{h}(t)$ and $r_{h}(t)$ are shown in
Figure 1 and that of the scaled mosquito populations $s_{v}(t),\ i_{v}(t)$ are
shown in Figure 2. In Figures 3 and 4, we considered four different values of
the fractional $\alpha$, $\alpha=1,\ 0.99,\ 0.95\ 0.90$. The oscillatory
behavior of the fractional order malaria model (2) are interesting and can be
compared to the research work be Tumwiine et al. (2007a). It is important to
observe that when $\alpha=1$, the model problem (2) is equivalent to the
classical initial value problem (2) and that as the fractional order $\alpha$
increases the behaviour of the fractional order model solutions approaches
that of the integer order model. Phase plane portrait is displayed in Figure 5
to describe the dynamics between susceptible fractions of human hosts and
infected fractions of mosquito vectors.

Figure 1: Solutions of fractional order model with $\alpha=1$.

Figure 2: Solutions of fractional order model with $\alpha=1$.

Figure 3: Solutions of fractional order model with $\alpha=1,0.99,0.95,0.90$.

Figure 4: Solutions of fractional order model with $\alpha=1$. Figure 5: Phase
portrait for fractional order model with $\alpha=1,0.99,0.95,0.90$.

##  7 Conclusion

In this research, we have proposed and studied fractional order model for the
integer order model developed and studied by Tumwiine et al. (2007a, b). We
determined the basic reproduction number $R_{o}$ using the method and the
notations by Driessche and Watmough (2002). We have investigated local
asymptotic stability analysis of the model equilibria. Furthermore, we have
numerically describe the trajectories of the model using Adams-type predictor-
corrector method. We have extended the existing malaria model propose by
Tumwiine et al. (2007a, b) to include fractional derivatives which is very
interesting and important in the field of mathematical biology.

## References

  * Abdullahi et al. (2013) M. B. Abdullahi, Y. A. Hasan, and F. A. Abdullah.  A mathematical model of malaria and the effectiveness of drugs.  _Applied Mathematical sciences_ , 7:3079–3095, 2013.
  * Abu-Raddad et al. (2006) L. J. Abu-Raddad, P. Patnaik, and j. G. Kublin.  Dual infection with hiv and malaria fuels the spread of both diseases in sub-saharan africa.  _Science_ , 314:1603–1606, 2006.
  * Águas et al. (2012) R. Águas, M. U. Ferreira, and M. G. M. Gomes.  Modeling the effects of relapse in the transmission dynamics of malaria parasites.  _Journal of Parasitology Research_ , 2012, Article ID 921715:8, 2012.
  * Ahmad and El-Khazali (2007) W. M. Ahmad and R El-Khazali.  Fractional-order dynamical models of love.  _Chaos, Solitons and Fractals_ , 33:1376–1375, 2007.
  * Ahmed and El-Saka (2010) E. Ahmed and H. A. El-Saka.  On fractional order models for hepatitis c.  _Nonlinear Biomedical Physics._ , 4:1–3, 2010.
  * Ahmed et al. (2006) E. Ahmed, A. M. A. El-Sayed, and H. A. A. El-Saka.  On some routh–hurwitz conditions for fractional order differential equations and their applications in lorenz, rössler, chua and chen systems.  _Physics Letters A_ , 358:1–4, 2006.
  * Ahmed et al. (2007) E. Ahmed, A. M. A. El-Sayed, and H. A. A. El-Saka.  Equilibrium points, stability and numerical solutions of fractional-order predator–prey and rabies models.  _Journal of Mathematical Analysis and Applications_ , 325:542– 553, 2007.
  * Al-Sulami et al. (2014) H. Al-Sulami, M. El-Shahed, J. J. Nieto, and W. Shammakh.  On fractional order dengue epidemic model.  _Mathematical Problems in Engineering_ , 2014, Article ID 456537:6, 2014.
  * Antia et al. (2008) R. Antia, A. Yates, and J. C. de Roode.  The dynamics of acute malaria infections. i. effect of the parasite’s red blood cell preference.  _Proceeding of the Royal Society. B_ , 275:1449–1458, 2008\.
  * Arafa et al. (2012) A. A. M. Arafa, S. Z. Rida, and M. Khalil.  Solution of fractional order model of childhood diseases with constant vaccination strategy.  _Mathematical Sciences Letters_ , 1:17–23, 2012.
  * Arafa et al. (2014) A. A. M. Arafa, S. Z. Rida, and M. Khalil.  A fractional-order model of hiv infection with drug therapy effect.  _Journal of the Egyptian Mathematical Society._ , 22:538–543, 2014.
  * Area et al. (2015) I. Area, H. Batarfi, J. Losada, J. J. Nieto, W. Shammakh, and Á. Torres.  On a fractional order ebola epidemic model.  _Advances in Difference Equations 2015:278_ , 2015.
  * Barley et al. (2012) K. Barley, D. Murillo, S. Roudenko, A. M. Tameru, and S. Tatum.  A mathematical model of hiv and malaria co-infection in sub-saharan africa.  _Journal of AIDS Clinical Research_ , 3:173. doi: 10.4172/2155-6113.1000173, 2012.
  * Chitnis (2005) N. Chitnis.  _Using mathematical models in controlling the spread of malaria_.  PhD thesis, University of Arizona, Tucson, Arizona, USA, 2005.
  * Chitnis et al. (2006) N. Chitnis, J. M. Cushing, and J. M. Hyman.  Bifurcation analysis of a mathematical model for malaria transmission.  _SIAM J. Appl. Math_ , 67:24–45, 2006.
  * Chiyaka et al. (2008) C. Chiyaka, J. M. Tchuenche, W. Garira, and S. Dube.  A mathematical analysis of the effects of control strategies on the transmission dynamics of malaria.  _Applied Mathematics and Computation_ , 195:1641–662, 2008\.
  * Cho et al. (2015) Y. Cho, I. Kim, and D. Sheen.  Afractional-order model for minmod millennium.  _Mathematical Biosciences_ , 262:36–45, 2015.
  * Diethelm et al. (2002) K. Diethelm, N. J. Ford, and A. D. Freed.  A predictor-corrector approach for the numerical solution of fractional differential equations.  _Nonlinear Dyamics_ , 29:3–22, 2002.
  * Diethelm et al. (2004) K. Diethelm, N. J. Ford, and A. D. Freed.  Detailed error analysis for a fractional adams method.  _Nonlinear Dyamics_ , 36:31–52, 2004.
  * Ding and Ye (2009) Y. Ding and H. Ye.  A fractional-order differential equation model of hiv infection of cd4+ t-cells.  _Mathematical and Computer Modelling_ , 50:386–392, 2009\.
  * Ding et al. (2012) Y. Ding, Z. Wang, and H. Ye.  Optimal control of a fractional-order hiv-immune system with memory.  _IEEE Trans. Contr. Sys. Techn._ , 20:763–769, 2012.
  * Driessche and Watmough (2002) P. V. D. Driessche and J. Watmough.  Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission.  _Mathematical Biosciences_ , 180:29–48, 2002.
  * El-Saka (2015) H. A. A. El-Saka.  Backward bifurcations in fractional-order vaccination models.  _Journal of the Egyptian Mathematical Society._ , 23:49–55, 2015.
  * El-Sayed et al. (2007) A. M. A. El-Sayed, A. E. M. El-Mesiry, and H. A. A. El-Saka.  On the fractional-order logistic equation.  _Applied Mathematics Letters 20_ , 20:817–823, 2007.
  * El-Shahed and Alsaedi (2011) M. El-Shahed and A. Alsaedi.  The fractional sirc model and influenza a.  _Mathematical Problems in Engineering_ , 2011, Article ID 480378:9, 2011.
  * Gao et al. (2014) D. Gao, Y. Lou, and S. Ruan.  A periodic Ross-Macdonald model in a patchy environment.  _Discrete Continuous Dyn Syst Ser B_ , 19:3133–3145, 2014\.
  * Ghosh et al. (2013) M. Ghosh, A. A. Lashari, and Xue-Zhi Li.  Biological control of malaria: A mathematical model.  _Applied Mathematics and Computation_ , 219:923–7939, 2013\.
  * Gokdogan et al. (2011) A. Gokdogan, A. Yildirim, and M. Merdan.  Solving a fractional order model of hiv infection of cd4+ t cells.  _Mathematical Problems in Engineering_ , 54:2132–2138, 2011\.
  * González-Parra et al. (2014) G. González-Parra, A. J. Arenas, and B. Chen-Charpentier.  A fractional order epidemic model for the simulation of outbreaks of influenza A(H1N1).  _Mathematical Methods in the Applied Sciences_ , 37:2218–2226, 2014.
  * Hoshen and Morse (2004) M. B. Hoshen and A. P. Morse.  Spread of anti-malarial drug resistance: Mathematical model with implications for act drug policies.  _Malaria Journal_ , 3:32, 2004.
  * Javidi and Ahmad (2014) M. Javidi and B. Ahmad.  A study of fractional-order cholera model.  _Appl.Math.Inf.Sci_ , 8:2195–2206, 2014.
  * Keegan and Dushoff (2013) L. T. Keegan and J. Dushoff.  Population-level effects of clinical immunity to malaria.  _BMC Infectious Diseases 2013_ , 13:428, 2013.
  * Koella and Anita (2003) J. C. Koella and R. Anita.  Epidemiological models for the spread of anti-malarial resistance.  _Malaria Journal_ , 2:3, 2003.
  * Lawi et al. (2011) O. Lawi, J. Y. T. Mugisha, and N. Omolo-Ongati.  Mathematical model for malaria and meningitisco-infection among children.  _Applied Mathematical Sciences_ , 5:2337–2359, 2011.
  * Li and Tao (2009) C. Li and C. Tao.  On the fractional adams method.  _Computers and Mathematics with Applications_ , 58:1573–1588, 2009.
  * Lin (2007) W. Lin.  Global existence theory and chaos control of fractional differential equations.  _Journal of Mathematical Analysis and Applications_ , 332:709–726, 2007.
  * Liu and Chen (2015) W. Liu and W. Chen.  Chaotic behavior in a new fractional-order love triangle system with competition.  _Journal of Applied Analysis and Computation_ , 5:103–113, 2015.
  * Liu et al. (2014) Y. Liu, P. Lu, and I. Szanton.  Numerical analysis for a fractional differential time-delay model of hiv infection of cd4+ t-cell proliferation under antiretroviral therapy.  _Abstract and Applied Analysis_ , 2014, Article ID 291614:13, 2014.
  * Liu and Lu (2014) Z. Liu and P. Lu.  Stability analysis for hiv infection of cd4+ t-cells by a fractional differential time-delay model with cure rate.  _Advances in Difference Equations 2014_ , 298, 2014.
  * Mandal et al. (2011) S. Mandal, R. R. Sarkar, and S. Sinha.  Mathematical models of malaria - a review.  _Malaria Journal_ , 10:202, 2011.
  * Matignon (1996) D. Matignon.  Stability results for fractional differential equations with applications to control processing.  _Computational Engineering in Systems Applications_ , 2:963–968, 1996.
  * Mukandavire et al. (2009) Z. Mukandavire, A. B. Gumel, W. Garira, and J. M. Tchuenche.  Mathematical analysis of a model for hiv-malaria co-infection.  _Math Biosci Eng_ , 6:333–362, 2009.
  * Ngonghala et al. (2012) C. N. Ngonghala, G. A. Ngwa, and M. I. Teboh-Ewungkem.  Periodic oscillations and backward bifurcation in a model for the dynamics of malaria transmission.  _Mathematical Biosciences_ , 240:45–62, 2012.
  * Ngwa and Shu (2000) G. A. Ngwa and W. S. Shu.  A mathematical model for endemic malaria with variable human and mosquito populations.  _Mathematical and Computer Modelling_ , 32:747–763, 2000\.
  * Odibat and Shawagfeh (2007) Z. M. Odibat and N. T. Shawagfeh.  Generalized taylor’s formula.  _Applied Mathematics and Computation_ , 186:286–293, 2007\.
  * Okosun and Makinde (2011) K. O. Okosun and O. D. Makinde.  Modelling the impact of drug resistance in malaria transmission and its optimal control analysis.  _International Journal of the Physical Sciences_ , 6:6479–6487, 2011.
  * Okosun and Makinde (2013) K. O. Okosun and O. D. Makinde.  Optimal control analysis of malaria in the presence of non-linear incidence rate.  _Appl. Comput. Math._ , 12:20–32, 2013.
  * Okosun and Makinde (2014) K. O. Okosun and O. D. Makinde.  A co-infection model of malaria and cholera diseases with optimal control.  _Mathematical Biosciences_ , 258:19–32, 2014.
  * Okyere et al. (2016) E. Okyere, F. T. Oduro, S. K. Amponsah, I. K. Dontwi, and N. K. Frempong.  Fractional order sir model with constant population.  _British Journal of Mathematics and Computer Science_ , 14:1–12, 2016.
  * Ozalp and Demirci (2011) N. Ozalp and E. Demirci.  A fractional order seir model with vertical transmission.  _Mathematical and Computer Modelling_ , 54:1–6, 2011.
  * Petras (2011) I. Petras.  _Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation_.  Springer, 2011.
  * Pinto and Carvalho (2014) C. M. A. Pinto and A. R. M. Carvalho.  New findings on the dynamics of hiv and tb coinfection modelsl.  _Applied mathematics and computation_ , 242:36–46, 2014\.
  * Pinto and Machado (2013) C. M. A. Pinto and J. A. T. Machado.  Fractional model for malaria transmission under control strategies.  _Computers and Mathematics with Applications._ , 66:908–916, 2013.
  * Podlubny (1999) I. Podlubny.  _Fractional Differential Equations_.  Academic Press, New York, 1999.
  * Pongtavornpinyo et al. (2008) W. Pongtavornpinyo, S. Yeung, I. M. Hastings, A. M. Dondorp, N. P. J. Day, and N. J. White.  Spread of anti-malarial drug resistance: Mathematical model with implications for act drug policies.  _Malaria Journal_ , 7 :229, 2008.
  * Rihan et al. (2014) F. A. Rihan, D. Baleanu, S. Lakshmanan, and R. Rakkiyappan.  On fractional sirc model with salmonella bacterial infection.  _Abstract and Applied Analysis_ , 2014, Article ID 136263:9, 2014.
  * Ross (1911) R. Ross.  The prevention of malaria.  _John Murray, London_ , 1911.
  * Sardar et al. (2015) T. Sardar, S. Rana, S. Bhattacharya, K. Al-Khaled, and J. Chattopadhyay.  A generic model for a single strain mosquito-transmitted disease with memory on the host and the vector.  _Mathematical biosciences_ , 263:18–38, 2015.
  * Tumwiine et al. (2007a) J. Tumwiine, J. Y. T. Mugisha, and L. S. Luboobi.  On oscillatory pattern of malaria dynamics in a population with temporary immunity.  _Computational and Mathematical Methods in Medicine_ , 8:191–203, 2007a.
  * Tumwiine et al. (2007b) J. Tumwiine, J. Y. T. Mugisha, and L. S. Luboobi.  A mathematical model for the dynamics of malaria in a human host and mosquito vector with temporary immunity.  _Applied Mathematics and Computation_ , 189:1953–1965, 2007b.
  * Tumwiine et al. (2010) J. Tumwiine, J. Y. T. Mugisha, and L. S. Luboobi.  A host-vector model for malaria with infective immigrants.  _J. Math. Anal. Appl_ , 361:139–149, 2010.
  * WHO (2005) WHO.  World malaria report 2005.  Technical Report, World Health Organization/UNICEF, Geneva, Switzerland., 2005.
  * WHO (2012) WHO.  World malaria report 2012.  Technical Report, World Health Organization/UNICEF, Geneva, Switzerland., 2012.
  * WHO (2014) WHO.  World malaria report 2014.  Technical Report, World Health Organization/UNICEF, Geneva, Switzerland., 2014.
  * Ye and Ding (2009) H. Ye and Y. Ding.  Nonlinear dynamics and chaos in a fractional-order hiv model.  _Mathematical Problems in Engineering_ , 2009, Article ID 378614:9, 2009.