Tomul LXVI, 2020, f. 2

Positively invariant semi-implicit discrete model for malaria propagation

Istv´an Farag´o · Rahele Mosleh

Abstract In this paper, we consider the malaria propagation process for infected pop- ulations for humans and mosquitoes with an extension of the classical Ross model. The numerical model is constructed by using the semi-implicitθ-method. We show that the numerical solution and the total populations of the extended Ross model are positively invariant in specific intervals for some step sizes. We demonstrate the validity of the the- oretical results of the semi-implicit method by numerical simulations for some examples.

The results are significant extension of the results of paper [9].

Keywords Ross model· Extended Ross model· Positively invariant · Positivity preservation·Semi-implicit method·Malaria propagation

Mathematics Subject Classification (2010) 65L05·34C60·92D30

1 Introduction

Malaria is an infectious and lethal host-vector disease that is transmitted by the bites of infectious female Anopheles mosquitoes as vectors to humans as hosts. It is an ancient disease and individuals have to deal with this disease since there is no effective vaccine. This phenomenon depends heavily on climatic factors, including temperature, altitude, precipitation, and raised humidity. This implies that the distribution and transmission of malaria is increased by high temperatures, rainy seasons and increased humidity, the reason it is spreaded out in tropical and subtropical regions such as Africa, Asia, Latin America and also some parts of Europe like Hungary, Austria, Italy [1], [13].

Some mathematical models such as Ross, Ross-Macdonald, delayed Ross- Macdonald, Anderson and May models and others have been built to gain a better insight into the spread of malaria and decrease its effects in the world. As an SIS-Type model, the Ross model is discussed for humans and mosquitoes in constant population sizes and the solutions are densities in

the interval [ 0,1] (c.f. [7]). In this analysis, we consider an expansion of the Ross model to take more effective malaria factors into account.

As the solutions of the extended Ross mode are the number of the indi- viduals, they should be positive. Analytically, it is shown that the extended Ross model is well-posed and the solutions are positively invariant in certain intervals (c.f. [6]). Moreover, there are two equilibrium points, disease-free and endemic, which under certain constraints, are globally stable (c.f. [2], [5], [8]). The numerical model obtained by using the explicit Euler method is investigated in [9] and the condition of the positivity preservation prop- erty is given. This approach allows us to use only relatively small step sizes.

Therefore in this paper we suggest to use implicite discretization method which enables us to select bigger step-size. The fully implicit method for the extended Ross model is rigid to prove the positivity preservation for the solutions, therefore a semi-implicit method is applied. We show that this kind of discretization allows us to select bigger step-size under which the postovoty is preserved. Moreover, in this paper we also investigate the upper bound preservation property of the solution.

As a brief outline, section 2 interprets the extended Ross model biologically and qualitatively. In section 3 we consider a semi-implicit method to prove the positively invariant solutions for humans and mosquitoes. Section 4 simulates the theoretical results of the semi-implicit method numerically with examples. Section 5 gives a summary of the results.

2 Mathematical models for malaria transmission

A diversity of mathematical models of malaria transmission have been de- veloped to explain the dynamics of malaria spread as epidemiological mod- els to some extent. Ronald Ross [12] introduced the first model of malaria transmission. In this model, the human and mosquito populations are split into two subclasses, infected and susceptible individuals, and the birth and mortality rates are the same.

2.1 Ross model

According to the Ross model, the dynamics of malaria transfer is described through a system of ordinary differential equations as below:

S˙_{h}(t) =r−αI_{m}(t)(1−I_{h}(t))−rS_{h}(t)
I˙h(t) =αIm(t)(1−Ih(t))−rIh(t)
S˙m(t) =µ−βI_{h}(t)(1−Im(t))−µSm(t)
I˙_{m}(t) =βI_{h}(t)(1−I_{m}(t))−µI_{m}(t)

(2.1)

with given nonnegative initial values satisfying the following conditions:

Sh(0) +Ih(0) = 1

Sm(0) +Im(0) = 1 (2.2)

Here S_{h}(t) denotes the density of host susceptible humans at time t, the
density of infected human at time t is defined by I_{h}(t), Sm(t) and Im(t),
respectively, signify the densities of susceptible and infected mosquitoes at
time t. Biologically, the other parametersα, β, r andµ are defined as the
proportion of bites that produce infection in humans, the proportion of bites
by which one susceptible mosquito becomes infectious, the rate of the death
for humans and the rate of the death for mosquitoes, respectively.

Since the solutions of the Ross models are densities, they should take values in the interval [0,1]. This property is proven in [7] analytically and numeri- cally. The Ross model (2.1) fails to introduce an adequate model for malaria transmission by the following reasons:

– There are no vital dynamics for the total populations. In other words, the birth and mortality rates are assumed to be equal, outcoming the total population sizes are constant for humans and mosquitoes.

– As malaria parasites have some latency period in humans and mosquitoes bodies, there is an intermediate state between the two susceptible and infected compartments known exposed. It means that the transition from the susceptible state to the infected state is not direct.

– Since the Ross model (2.1) is aSIS-type model, the recovery state,R(t), is not considered for humans.

2.2 Extended Ross model

In this section, we discuss an extended Ross model which is able to dispel the aforementioned deficits of the Ross model (2.1) to some extent (c.f. [10]). In this model we divide the humans population into four groups: susceptible humans, exposed humans, infectious humans and recovered humans.The population of mosquitoes is split into three subclasses, unlike the humans population. Due to the short life time, mosquitoes do not have enough time to meloriate and they die after infection. By some biological interpretations [10] the malaria propagation is sketched as below:

S˙_{h}(t) =Λ_{h}−^{bβ}_{1+}^{h}^{S}_{ν}^{h}_{h}^{(}^{t}_{I}^{)}_{m}^{I}^{m}_{(}_{t}^{(}_{)}^{t}^{)}−µ_{h}S_{h}(t) +ωR_{h}(t)
E˙_{h}(t) = ^{bβ}_{1+}^{h}^{S}_{ν}^{h}^{(}^{t}^{)}^{I}^{m}^{(}^{t}^{)}

hIm(t) −(α_{h}+µ_{h})E_{h}(t)
I˙_{h}(t) =α_{h}E_{h}(t)−(r+µ_{h}+δ_{h})I_{h}(t)
R˙_{h}(t) =rI_{h}(t)−(µ_{h}+ω)R_{h}(t)
S˙_{m}(t) =Λ_{m}−^{bβ}_{1+}^{m}^{S}_{ν}^{m}_{m}^{(}_{I}^{t}^{)}_{h}^{I}_{(}^{h}_{t}_{)}^{(}^{t}^{)}−µ_{m}S_{m}(t)
E˙m(t) = ^{bβ}_{1+}^{m}^{S}_{ν}^{m}^{(}^{t}^{)}^{I}^{h}^{(}^{t}^{)}

mIh(t) −(αm+µm)Em(t) I˙m(t) =αmEm(t)−(µm+δm)Im(t).

(2.3)

The following initial conditions are applied:

Sh(0) =S_{0}h, Ih(0) =I_{0}h, Eh(0) =E_{0}h, Rh(0) =R_{0}h,

Sm(0) =S0m, Em(0) =E0m, Im(0) =I0m. (2.4)

Here functionS_{h}(t) denotes the number of the host susceptible humans at
time t, the number of the host exposed humans to malaria infection at time
t is denoted byEh(t), the number of the infected and the recovered humans
at time t are characterized byIh(t) andRh(t) respectively.

In order to interpret a more accurate model, the model (2.3) takes into account population dynamics by considering various birth and death rates.

We assume that all the newly born children are healthy, so the recruitment
terms, (Λ_{h}, Λ_{m}), are the birth rates applied as inputs for the susceptible
class (S_{h}, Sm).

The biting rate of the mosquito is denoted by b and the parasites injected by
the mosquito into the human blood system with some probabilityβhspend
some a latent period in human body then the susceptible human moves
to the exposed class Eh(t). After the latency period the exposed human
is proceeded to the infected class with some α_{h} rate. The infected human
transfers to the recovery class with some rate r. The recovered stateR_{h}(t)
contains the recovered humans. The recovered human has some immunity
to the disease however, after a while the individual loses the immunity and
becomes susceptible again by ω pace. Each class of humans is decreased
by the natural death rate, µh, according to the model (2.3), and only the
infectious class is decreased by the death rate caused by disease,δh.
In a similar way, βm is the probability of the parasites transitioning into
the susceptible mosquitoS_{m}(t) through the bites of an infected human and
the susceptible mosquito moves to the exposed class E_{m}(t). After a given
time it becomes infectiousI_{m}. The mosquitoes population is decreased by
natural death rateµmor disease induced death rateδm.

The ratio 1

1 +ν_{h}I_{m}(t)denotes a saturating feature that prevents the force of
infection from infected mosquitoes to susceptible humans in whichν_{h}∈[0,1]

is the proportion of antibodies produced in humans bodies in response to the conflict of antigens produced by infected mosquitoes. By similar inter- pretation for mosquitoes,νm∈[0,1] is the rate of antibodies generated by infectious humans against the antigens contacted.

Since the malaria transmission occurs in an inharmonious population, the
epidemiological model (2.3) must partition the population into groups, in
which the members have similar characteristics such as mode of transmis-
sion, contact patterns, latent period, infectious period, genetic susceptibility
or resistance. Accordingly, the total humans population size at time t,V_{h}(t),
is defined as

Vh(t) =Sh(t) +Eh(t) +Ih(t) +Rh(t). (2.5)
The same definition holds for the total size of the mosquitoes population,
V_{m}(t), at time t

V_{m}(t) =S_{m}(t) +E_{m}(t) +I_{m}(t). (2.6)
Analytically, it is proven [6,10] that solutions of the extended Ross model
(2.3) are invariant in the intervals (0,^{Λ}_{µ}^{h}

h] and (0,^{Λ}_{µ}^{m}

m] at time t for humans and mosquitoes respectively. According to the aforementioned intervals, when

the rate of natural mortality for either humans or mosquitoes vanishes, the size of the corresponding total population is unbounded. Biologically, it means when there is no natural death, then the number of the total population blows up.

We note that 2.3-2.4 yields a Cauchy-problem for the system of nonlinear ordinary differential equations. Since the analytical solution cannot be de- fined, we have to use some numerical method, which results in the discrete model got the malaria propagation.

3 Semi-implicit method applied for the extended Ross model At the first step we apply the standard implicitθ- method for the extended Ross model to approximate the numerical solutions of the system (2.3) as below:

S_{h}^{i+1}−Sh^{i}

∆t = (1−θ)(Λ_{h}−^{bβ}_{1+}^{h}_{ν}^{S}_{h}^{h}^{i}_{I}^{I}_{m}^{i}^{m}^{i} −µ_{h}S_{h}^{i} +ωR^{i}_{h})
+θ(Λ_{h}−^{bβ}_{1+}^{h}^{S}_{ν}^{h}^{i+1}_{h}_{I}_{m}^{i+1}^{I}^{m}^{i+1} −µ_{h}S_{h}^{i}^{+1}+ωR^{i}_{h}^{+1})

E^{i+1}_{h} −Eh^{i}

∆t = (1−θ)(^{bβ}_{1+}^{h}_{ν}^{S}^{h}^{i}^{I}^{m}^{i}

hI_{m}^{i} −(α_{h}+µ_{h})E_{h}^{i})
+θ(^{bβ}^{h}^{S}^{h}^{i+1}^{I}^{i+1}^{m}

1+νhIm^{i+1} −(α_{h}+µ_{h})E_{h}^{i}^{+1})

I_{h}^{i+1}−Ih^{i}

∆t = (1−θ)(α_{h}E_{h}^{i} −(r+µ_{h}+δ_{h})I_{h}^{i})
+θ(α_{h}E_{h}^{i}^{+1}−(r+µ_{h}+δ_{h})I_{h}^{i}^{+1})

R^{i+1}_{h} −R^{i}_{h}

∆t = (1−θ)(rI_{h}^{i} −(µh+ω)R^{i}_{h})
+θ(rI_{h}^{i}^{+1}−(µh+ω)R^{i}_{h}^{+1})

S_{m}^{i+1}−S_{m}^{i}

∆t = (1−θ)(Λm−^{bβ}_{1+}^{m}_{ν}^{S}_{m}^{i}^{m}_{I}^{I}^{i}^{h}^{i}

h −µmS_{m}^{i} )
+θ(Λ_{m}−^{bβ}_{1+}^{m}^{S}_{ν}^{i+1}^{m}_{m}_{I}^{i+1}^{I}^{h}^{i+1}

h −µ_{m}S_{m}^{i}^{+1})

E^{i+1}_{m} −E_{m}^{i}

∆t = (1−θ)(^{bβ}_{1+}^{m}_{ν}^{S}^{m}^{i}^{I}^{h}^{i}

mI_{h}^{i} −(α_{m}+µ_{m})E_{m}^{i} )
+θ(^{bβ}^{m}^{S}^{m}^{i+1}^{I}^{h}^{i+1}

1+νmI^{i+1}_{h} −(αm+µm)E_{m}^{i}^{+1})

I_{m}^{i+1}−Im^{i}

∆t = (1−θ)(αmE_{m}^{i} −(µm+δm)I_{m}^{i} )
+θ(α_{m}E_{m}^{i}^{+1}−(µ_{m}+δ_{m})I_{m}^{i}^{+1}).

(3.1)

However, the system (3.1) is a system of nonlinear algebraic equations.

Therefore, there are some difficulties to prove the positivity preservation of the solutions by algebraic tools. Since the nonlinearity is a specific structure appearing in the force of infection terms only, we apply the following semi- implicit method a nonlocal discretization of the implicit method (3.1) :

S_{h}^{i+1}−Sh^{i}

∆t = (1−θ)(Λ_{h}−^{bβ}_{1+}^{h}^{S}_{ν}_{h}^{i}^{h}_{I}^{I}_{m}^{i}^{m}^{i} −µ_{h}S_{h}^{i} +ωR^{i}_{h})
+θ(Λh−^{bβ}_{1+}^{h}^{S}_{ν}^{h}^{i+1}_{h}_{I}_{m}^{i}^{I}^{m}^{i} −µhS_{h}^{i}^{+1}+ωR^{i}_{h})

E^{i+1}_{h} −E_{h}^{i}

∆t = (1−θ)(^{bβ}_{1+}^{h}^{S}_{ν}^{i}^{h}^{I}^{m}^{i}

hI_{m}^{i} −(αh+µh)E_{h}^{i})
+θ(^{γ}_{1+}^{h}^{S}^{i+1}^{h}_{ν} ^{I}^{m}^{i}

hI_{m}^{i} −(α_{h}+µ_{h})E_{h}^{i}^{+1})

I_{h}^{i+1}−I_{h}^{i}

∆t = (1−θ)(α_{h}E_{h}^{i} −(r+µ_{h}+δ_{h})I_{h}^{i})
+θ(α_{h}E_{h}^{i}^{+1}−(r+µ_{h}+δ_{h})I_{h}^{i}^{+1})

R^{i+1}_{h} −R^{i}_{h}

∆t = (1−θ)(rI_{h}^{i}−(µ_{h}+ω)R^{i}_{h})
+θ(rI_{h}^{i}^{+1}−µhR^{i}_{h}^{+1}−ωR_{h}^{i})

S_{m}^{i+1}−S_{m}^{i}

∆t = (1−θ)(Λ_{m}− ^{bβ}_{1+}^{m}_{ν}^{S}_{m}^{m}^{i}_{I}^{I}_{h}^{i}^{h}^{i} −µ_{m}S_{m}^{i} )
+θ(Λ_{m}− ^{bβ}_{1+}^{m}^{S}_{ν}^{m}^{i+1}_{m}_{I}^{i+1}^{I}^{h}^{i+1}

h −µ_{m}S^{i}_{m}^{+1})

E^{i+1}_{m} −Em^{i}

∆t = (1−θ)(^{bβ}_{1+}^{m}_{ν}^{S}^{m}^{i}^{I}^{h}^{i}

mI_{h}^{i} −(α_{m}+µ_{m})E_{m}^{i} )
+θ(^{bβ}^{m}^{S}^{m}^{i+1}^{I}^{h}^{i+1}

1+νmI_{h}^{i+1} −(αm+µm)E_{m}^{i}^{+1})

I_{m}^{i+1}−I_{m}^{i}

∆t = (1−θ)(αmE_{m}^{i} −(µm+δm)I_{m}^{i} )
+θ(α_{m}E^{i}_{m}^{+1}−(µ_{m}+δ_{m})I_{m}^{i}^{+1}).

(3.2)

The system (3.2) is a linear system for which we are able to define the solutions explicitly as below:

S_{h}^{i}^{+1}= S_{h}^{i}(1−∆t(1−θ)(_{1+}^{bβ}_{ν}^{h}^{I}^{m}^{i}

hI_{m}^{i} +µ_{h})) +∆t(Λ_{h}+ωR^{i}_{h})
1 +∆tθ(_{1+}^{bβ}_{ν}^{h}^{I}^{m}^{i}

hI_{m}^{i} +µ_{h}) , (3.3)
E_{h}^{i}^{+1}= E_{h}^{i}(1−∆t(1−θ)(α_{h}+µh)) +∆t(1−θ)^{bβ}_{1+}^{h}^{S}_{ν}^{h}^{i}^{I}^{m}^{i}

hI_{m}^{i} +∆tθ^{bβ}_{1+}^{h}^{S}_{ν}^{h}^{i+1}^{I}^{m}^{i}

hI_{m}^{i}

1 +∆tθ(α_{h}+µ_{h}) ,

(3.4)
I_{h}^{i}^{+1}=I_{h}^{i}(1−∆t(1−θ)(r+µ_{h}+δ_{h})) +∆t(1−θ)α_{h}E_{h}^{i} +∆tθα_{h}E_{h}^{i}^{+1}

1 +∆tθ(r+µ_{h}+δ_{h}) ,
(3.5)
R^{i}_{h}^{+1}= R^{i}_{h}(1−∆t((1−θ)µ_{h}+ω)) +∆t(1−θ)rI_{h}^{i} +∆tθrI_{h}^{i}^{+1}

1 +∆tθµh

, (3.6)
S_{m}^{i}^{+1}= S_{m}^{i} (1−∆t(1−θ)(_{1+}^{bβ}_{ν}^{m}^{I}^{h}^{i}

mI^{i}_{h}+µ_{m})) +∆tΛ_{m}
1 +∆tθ( ^{bβ}^{m}^{I}^{h}^{i+1}

1+νmI_{h}^{i+1} +µ_{m})

, (3.7)

E_{m}^{i}^{+1}=

E_{m}^{i} (1−∆t(1−θ)(αm+µm))+∆t(1−θ)^{bβ}_{1+}^{m}_{ν}^{S}^{m}^{i}^{I}^{i}^{h}

mI_{h}^{i} +∆tθ^{bβ}^{m}^{S}^{i+1}^{m} ^{I}^{h}^{i+1}

1+νmI_{h}^{i+1}

1+∆tθ(αm+µm) ,

(3.8)

I_{m}^{i}^{+1}= I_{m}^{i} (1−∆t(1−θ)(µ_{m}+δm)) +∆t(1−θ)α_{m}E_{m}^{i} +∆tθαmE_{m}^{i}^{+1}

1 +∆tθ(µ_{m}+δm) .

(3.9)
Since analytically we proved the solutions and total populations of the ex-
tended Ross model (2.3) are invariant on the intervals (0,^{Λ}_{µ}^{h}

h] and (0,^{Λ}_{µ}^{m}

m]
for humans and mosquitoes, respectively (c.f. [6], [8]), we shall consider this
property for the discrete model (3.2), too. For this aim, we need to guar-
antee that for any S_{h}^{i}, E_{h}^{i}, I_{h}^{i}, R^{i}_{h}, and V_{h}^{i} on the interval (0,^{Λ}_{µ}^{h}

h] and for
any S_{m}^{i} , E_{m}^{i} , I_{m}^{i} , andV_{m}^{i} on the interval (0,^{Λ}_{µ}^{m}

m] the discrete model (3.2)
results the solutionsS_{h}^{i}^{+1},E_{h}^{i}^{+1}, I_{h}^{i}^{+1}, R^{i}_{h}^{+1}, and V_{h}^{i}^{+1} on (0,^{Λ}_{µ}^{h}

h] and S_{m}^{i}^{+1},
E_{m}^{i}^{+1},I_{m}^{i}^{+1}, andV_{m}^{i}^{+1}on (0,^{Λ}_{µ}^{m}

m]. At the first step we consider the positivity preservation of the solutions.

Lemma 1.LetΛ_{h}, Λ_{m}>0 and the initial data of the system (3.2) are on
the intervals (0,^{Λ}_{µ}^{h}

h] and (0,^{Λ}_{µ}^{m}

m], respectively. If

∆t

(≤ _{(1}_{−θ}_{)(}_{bβ}_{h}^{µ}_{Λ}^{m}_{m}^{+}_{+}^{ν}_{µ}^{h}_{h}^{Λ}_{(}_{µ}^{m}_{m}_{+}_{ν}_{h}_{Λ}_{m}_{))} θ∈[0,1)

is any θ= 1, (3.10)

thenS_{h}^{i}^{+1} of the system (3.2) is positive.

Proof. From the first equation of the system (3.2) we have
S_{h}^{i}^{+1}(1+∆tθ( bβ_{h}I_{m}^{i}

1 +ν_{h}I_{m}^{i} +µ_{h})) =S_{h}^{i}+∆t(Λ_{h}−(1−θ)(bβhS_{h}^{i}I_{m}^{i}

1 +ν_{h}I_{m}^{i} +µ_{h}S_{h}^{i})+ωR^{i}_{h}).

(3.11)
At the first step, we assume thatθ does not equal one, i.e., 1−θ6= 0
If the right side of (3.11) is positive, then the required inequalityS_{h}^{i}^{+1}>0
holds. This means the condition

S_{h}^{i} > ∆t((1−θ)(bβhS_{h}^{i}I_{m}^{i}

1 +νhI_{m}^{i} +µ_{h}S_{h}^{i})−(Λ_{h}+ωR^{i}_{h})). (3.12)
When

(1−θ)(bβhS_{h}^{i}I_{m}^{i}

1 +ν_{h}I_{m}^{i} +µ_{h}S_{h}^{i})−(Λ_{h}+ωR^{i}_{h})≤0, (3.13)
then the inequality (3.12) is valid for any step size∆t. Otherwise, we have
the bound

∆t < S_{h}^{i}

(1−θ)(^{bβ}_{1+}^{h}_{ν}^{S}^{h}^{i}^{I}^{m}^{i}

hI_{m}^{i} +µhS_{h}^{i})−(Λh+ωR^{i}_{h}). (3.14)

Clearly,
S_{h}^{i}
(1−θ)(^{bβ}_{1+}^{h}^{S}_{ν}^{i}^{h}^{I}^{m}^{i}

hI_{m}^{i} +µhS_{h}^{i}) < S_{h}^{i}
(1−θ)(^{bβ}_{1+}^{h}^{S}_{ν}^{h}^{i}^{I}^{m}^{i}

hI_{m}^{i} +µhS_{h}^{i})−(Λh+ωR^{i}_{h}). (3.15)
Therefore, the condition

∆t≤ 1 +ν_{h}I_{m}^{i}

(1−θ)(bβhI_{m}^{i} +µh(1 +νhI_{m}^{i} )) (3.16)
is a sufficient condition of the positivity for S_{h}^{i}^{+1}. Since the right side of
(3.16) is a monotonically decreasing function ofI_{m}^{i} andI_{m}^{i} ∈(0,^{Λ}_{µ}^{m}

m], under the condition

∆t≤ µm+νhΛm

(1−θ)(bβhΛm+µh(µm+νhΛm)), (3.17)
the required positivity property for S_{h}^{i}^{+1} holds. For θ = 1, we have the
following equation

S_{h}^{i}^{+1}= S_{h}^{i} +∆t(Λ_{h}+ωR^{i}_{h})
1 +∆t(_{1+}^{bβ}_{ν}^{h}^{I}^{m}^{i}

hI_{m}^{i} +µ_{h}). (3.18)
Obviously, for this case,S_{h}^{i}^{+1} is positive unconditionally for any step size

∆t.

u t In the following, we analyze the positivity of the other classes for humans.

The second equation of the system (3.2) yields
E_{h}^{i}^{+1}(1+∆tθ(α_{h}+µ_{h})) =E_{h}^{i}+∆t((1−θ)(bβhS_{h}^{i}I_{m}^{i}

1+ν_{h}I_{m}^{i} −(α_{h}+µ_{h})E^{i}_{h})+θbβ_{h}S_{h}^{i}^{+1}I_{m}^{i}
1+ν_{h}I_{m}^{i} ).

(3.19)
The value of E_{h}^{i}^{+1} is positive provided that the right side of the equation
(3.19) is positive. Hence, when

(1−θ)(bβhS_{h}^{i}I_{m}^{i}

1 +νhI_{m}^{i} −(α_{h}+µ_{h})E_{h}^{i}) +θbβ_{h}S^{i}_{h}^{+1}I_{m}^{i}

1 +νhI_{m}^{i} ≥0, (3.20)
thenE_{h}^{i}^{+1}is positive for any step size∆t. Otherwise, for 1−θ6= 0 we have
the bound

∆t < E_{h}^{i}

(1−θ)((αh+µh)E_{h}^{i} −^{bβ}_{1+}^{h}^{S}_{ν}_{h}^{i}^{h}_{I}^{I}_{m}^{i}^{m}^{i} )−θ^{bβ}_{1+}^{h}^{S}_{ν}^{h}^{i+1}^{I}^{m}^{i}

hI_{m}^{i}

. (3.21)

Provided that (3.10) is satisfied, the following estimation is valid
E_{h}^{i}

(1−θ)(αh+µh)E_{h}^{i} ≤ E_{h}^{i}

(1−θ)((α_{h}+µ_{h})E_{h}^{i} −^{bβ}_{1+}^{h}_{ν}^{S}_{h}^{h}^{i}_{I}^{I}_{m}^{i}^{m}^{i} )−θ^{bβ}_{1+}^{h}^{S}_{ν}^{i+1}^{h} ^{I}^{m}^{i}

hI_{m}^{i}

. (3.22) Therefore,

∆t≤ 1

(1−θ)(α_{h}+µ_{h}) (3.23)

is a sufficient condition to attain the positiveE_{h}^{i}^{+1}>0. Forθ= 1 we get

E^{i}_{h}^{+1}= E_{h}^{i} +∆t^{bβ}_{1+}^{h}^{S}_{ν}^{i+1}^{h} ^{I}^{m}^{i}

hI_{m}^{i}

1 +∆t(α_{h}+µ_{h}) . (3.24)
Since in this caseS_{h}^{i}^{+1} (3.18) is positive unconditionally for any step size,
E_{h}^{i}^{+1}>0 for all step sizes. Similarly, we get I_{h}^{i}^{+1}>0 provided that (3.23)
is satisfied and

∆t≤ 1

(1−θ)(r+µ_{h}+δ_{h}), (3.25)
for 1−θ6= 0 and for all step sizes whenθ= 1. Moreover, we obtainR^{i}_{h}^{+1}>0
as long as (3.25) holds and

∆t≤ 1

(1−θ)µ_{h}+ω (3.26)

for 1−θ6= 0 and

∆t≤ 1

ω (3.27)

forθ= 1. Hence, we can summarize the results as follows.

Lemma 2.Assume thatΛ_{h}, Λm>0 and the initial data of the system (3.2)
are on the intervals (0,^{Λ}_{µ}^{h}

h] and (0,^{Λ}_{µ}^{m}

m], respectively. If

∆t ≤

min{_{(1}_{−}_{θ}_{)(}_{bβ}_{h}_{Λ}^{µ}^{m}_{m}^{+}_{+}^{ν}_{µ}^{h}_{h}^{Λ}_{(}_{µ}^{m}_{m}_{+}_{ν}_{h}_{Λ}_{m}_{))},_{(1}_{−}_{θ}_{)(}^{1}_{α}

h+µh),

1

(1−θ)(r+µh+δh),_{(1}_{−}_{θ}_{)}^{1}_{µ}

h+ω} θ∈[0,1)

1

ω θ= 1,

(3.28)

thenE_{h}^{i}^{+1}, I_{h}^{i}^{+1}, and R^{i}_{h}^{+1}of the system (3.2) are positive.

Similar results can be formulated for mosquitoes, as well.

Lemma 3.SupposeΛ_{h}, Λ_{m}>0 and the initial data of the system (3.2) are
on the intervals (0,^{Λ}_{µ}^{h}

h] and (0,^{Λ}_{µ}^{m}

m], respectively. If

∆t

(≤min{_{(1}_{−}_{θ}_{)(}_{r}_{+}^{1}_{µ}_{h}_{+}_{δ}_{h}_{)},_{(1}_{−}_{θ}_{)(}_{bβ} ^{µ}^{h}^{+}^{ν}^{m}^{Λ}^{h}

mΛh+µm(µh+νmΛh))} θ∈[0,1)

is any θ= 1, (3.29)

thenS_{m}^{i}^{+1} of the system (3.2) is positive.

Proof. From the fifth equation of the system (3.2) we have
S_{m}^{i}^{+1}(1+∆tθ( bβ_{m}I_{h}^{i}^{+1}

1 +νmI_{h}^{i}^{+1}+µ_{m})) =S^{i}_{m}+∆t(Λ_{m}−(1−θ)(bβmS_{m}^{i} I_{h}^{i}

1 +νmI_{h}^{i} +µ_{h}S_{m}^{i} )).

(3.30)
When 1−θ 6= 0 and (3.25) is satisfied, S_{m}^{i}^{+1} is positive provided that the
right side of the equation (3.30) is positive, i.e.,

S_{m}^{i} > ∆t((1−θ)(bβ_{m}S_{m}^{i} I_{h}^{i}

1 +νmI_{h}^{i} +µmS_{m}^{i} )−Λm). (3.31)
When

(1−θ)(bβmS_{m}^{i} I_{h}^{i}

1 +ν_{m}I_{h}^{i} +µ_{m}S_{m}^{i} )−Λ_{m}≤0, (3.32)
then the inequality (3.31) is valid for any step size∆t satisfied in the con-
dition (3.25). Otherwise, we have the bound

∆t≤ S_{m}^{i}

(1−θ)(^{bβ}_{1+}^{m}_{ν}^{S}^{m}^{i}^{I}^{h}^{i}

mI_{h}^{i} +µmS_{m}^{i} )−Λm

. (3.33)

Visibly,

S_{m}^{i}
(1−θ)(^{bβ}_{1+}^{m}_{ν}^{S}^{m}^{i}^{I}^{h}^{i}

mI^{i}_{h} +µmS_{m}^{i} ) < S_{m}^{i}
(1−θ)(^{bβ}_{1+}^{m}_{ν}^{S}^{i}^{m}^{I}^{h}^{i}

mI_{h}^{i} +µmS_{m}^{i} )−Λm

. (3.34) Accordingly, the condition

∆t≤ 1 +νmI_{h}^{i}

(1−θ)((1 +νmI_{h}^{i})µm+bβmI_{h}^{i}) (3.35)
is a sufficient condition of the positivity for S_{m}^{i}^{+1}. Since the right side of
(3.35) is a decreasing function ofI_{h}^{i} andI_{h}^{i} ∈(0,^{Λ}_{µ}^{h}

h], we attain

∆t≤ µh+νmΛh

(1−θ)(bβ_{m}Λh+µm(µ_{h}+νmΛh)). (3.36)
Supposeθ= 1, then we have the following expression

S_{m}^{i}^{+1}= S^{i}_{m}+∆tΛm

1 +∆t( ^{bβ}^{m}^{I}^{h}^{i+1}

1+νmI_{h}^{i+1} +µ_{m})

. (3.37)

As in this case I_{h}^{i}^{+1} is positive unconditionally for any step size, S_{m}^{i}^{+1} is
positive for any step size∆t.

u t The sixth equation of the system (3.2), yields

E_{m}^{i}^{+1}(1 +∆tθ(αm+µm)) =E_{m}^{i} +∆t((1−θ)(bβ_{m}S_{m}^{i} I_{h}^{i}

1 +νmI_{h}^{i} −(αm+µm)E_{m}^{i} )+

θbβmS_{m}^{i}^{+1}I_{h}^{i}^{+1}

1 +νmI_{h}^{i}^{+1} ). (3.38)
The value E_{m}^{i}^{+1} is positive as long as the right side of (3.38) is positive.

When

(1−θ)(bβ_{m}S_{m}^{i} I_{h}^{i}

1 +ν_{m}I_{h}^{i} −(α_{m}+µ_{m})E_{m}^{i} ) +θbβ_{m}S_{m}^{i}^{+1}I_{h}^{i}^{+1}

1 +νmI_{h}^{i}^{+1} ≥0, (3.39)
E_{m}^{i}^{+1} is positive for any step size∆t. Otherwise, for 1−θ6= 0 we have the
bound

∆t≤ E_{m}^{i}

(1−θ)((αm+µm)E_{m}^{i} −^{bβ}_{1+}^{m}_{ν}^{S}_{m}^{i}^{m}_{I}^{I}_{h}^{i}^{h}^{i})−θ^{bβ}^{m}^{S}^{m}^{i+1}^{I}^{i+1}^{h}

1+νmI_{h}^{i+1} )

. (3.40)

If (3.25) and (3.29) are satisfied, the following estimation is valid
E_{m}^{i}

(1−θ)(αm+µm)E_{m}^{i} ≤ E_{m}^{i}

(1−θ)((αm+µm)E^{i}_{m}−^{bβ}_{1+}^{m}_{ν}^{S}_{m}^{i}^{m}_{I}^{I}^{i}^{h}^{i}

h )−θ^{bβ}^{m}^{S}^{m}^{i+1}^{I}^{i+1}^{h}

1+νmI_{h}^{i+1} )
.
(3.41)
Hence,

∆t≤ 1

(1−θ)(αm+µm) (3.42)

is a sufficient condition to getE_{m}^{i}^{+1}>0. Ifθ= 1, we acquire

E_{m}^{i}^{+1}=

E^{i}_{m}+∆t^{bβ}^{m}^{S}^{m}^{i+1}^{I}^{h}^{i+1}

1+νmI_{h}^{i+1}

1 +∆t(αm+µm) . (3.43)

Since I_{h}^{i}^{+1} and S_{m}^{i}^{+1} are positive for all step sizes, E_{m}^{i}^{+1} is positive uncon-
ditionally for any step size. Similarly,I_{m}^{i}^{+1} is positive provided that (3.42)
holds and

∆t≤ 1

(1−θ)(δm+µm) (3.44)

for 1−θ6= 0 and without any restriction for any step size forθ= 1. There- fore, we can summarize the results as follows.

Lemma 4.SupposeΛ_{h}, Λ_{m}>0 and the initial data of the system (3.2) are
on the intervals (0,^{Λ}_{µ}^{h}

h] and (0,^{Λ}_{µ}^{h}

h], respectively. If

∆t

(≤min_{1}_{−θ}^{1} {_{(}_{bβ}_{m}_{Λ}_{h}^{µ}_{+}^{h}_{µ}^{+}_{m}^{ν}_{(}^{m}_{µ}^{Λ}_{h}_{+}^{h}_{ν}_{m}_{Λ}_{h}_{))},_{r}_{+}_{µ}^{1}

h+δh,_{α} ^{1}

m+µm,_{µ} ^{1}

m+δm} θ∈[0,1)

is any θ= 1,

(3.45)
then the solutionsE_{m}^{i}^{+1} andI_{m}^{i}^{+1}of the system (3.2) are positive.

Lemmas 1-4 give sufficient conditions for the positivity preservation prop- erty of the model (10). In the following we show a sharper property of the model.

Theorem 3.1 Assume that Λ_{h}, Λ_{m} >0 and the initial data of the system
(3.2) are on the intervals(0,^{Λ}_{µ}^{h}

h] and (0,^{Λ}_{µ}^{m}

m], respectively. If

h^{∗}=min{h_{1}, h_{2}}, (3.46)
where

h1 =min{ µ_{m}+ν_{h}Λ_{m}

(1−θ)(bβ_{h}Λm+µ_{h}(µm+ν_{h}Λm)), 1

(1−θ)(α_{h}+µ_{h}),
1

(1−θ)(r+µ_{h}+δ_{h}), 1

(1−θ)µ_{h}+ω} (3.47)
and

h_{2}=min 1

1−θ{ µ_{h}+νmΛ_{h}

(bβmΛ_{h}+µm(µ_{h}+νmΛ_{h})), 1

(αm+µm), 1 (µm+δm)},

(3.48)
for the semi-implicit model (3.2), then the solutions and total populations
are invariant on(0,^{Λ}_{µ}^{h}

h] and(0,^{Λ}_{µ}^{m}

m]for humans and mosquitoes respectively
for the discretization step size∆t ∈ (0, h^{∗}] for 1−θ 6= 0 and for step size

∆t∈(0,_{ω}^{1}]when θ= 1.

Proof. In this part, we shall prove invariance of the positivity preservation
of the total populations on the intervals (0,^{Λ}_{µ}^{h}

h] and (0,^{Λ}_{µ}^{m}

m] for humans and mosquitoes, respectively. For this aim, we should demonstrate if the number

of the total populations at time t_{i} are on the intervals V_{h}^{i} ∈ (0,^{Λ}_{µ}^{h}

h] and
V_{m}^{i} ∈ (0,^{Λ}_{µ}^{m}

m] for humans and mosquitoes respectively, the discrete model
(3.2) guarantees that V_{h}^{i}^{+1}∈ (0,^{Λ}_{µ}^{h}

h] andV_{m}^{i}^{+1}∈(0,^{Λ}_{µ}^{m}

m], too. The number of the total population in the discrete model (3.2) for humans at timeti+1

is defined as

V_{h}^{i}^{+1}=S_{h}^{i}^{+1}+E_{h}^{i}^{+1}+I_{h}^{i}^{+1}+R^{i}_{h}^{+1}. (3.49)
According to Lemma 1 and Lemma 2, if the initial data of the system (3.2)
are on the intervals (0,^{Λ}_{µ}^{h}

h] and (0,^{Λ}_{µ}^{m}

m], respectively and (3.28) satisfy, the
values ofS_{h}^{i}^{+1},E_{h}^{i}^{+1}, I_{h}^{i}^{+1}, andR^{i}_{h}^{+1} are positive, thereforeV_{h}^{i}^{+1}>0.

By adding the terms of the system (3.2) and the number of the total pop- ulation definition for humans (3.49) we attain

V_{h}^{i}^{+1}(1+∆tθµ_{h}) =V_{h}^{i}+∆t(Λ_{h}−(1−θ)µ_{h}V_{h}^{i}−(1−θ)δ_{h}I_{h}^{i}−θδ_{h}I_{h}^{i}^{+1})). (3.50)
If 1−θ6= 0 and (3.28) is satisfied, we have the bound

V_{h}^{i}^{+1}(1 +∆tθµ_{h})≤V_{h}^{i}(1−∆t(1−θ)µ_{h}) +∆tΛ_{h}. (3.51)
LetV_{h}^{i}≤ ^{Λ}µh^{h}, from the right side of the inequality (3.51) we obtain

V_{h}^{i}(1−∆t(1−θ)µh) +∆tΛh< Λ_{h}

µ_{h}(1−∆t(1−θ)µh+∆tµh). (3.52)
According to (3.51)-(3.52) we get

V_{h}^{i}^{+1}(1 +∆tθµ_{h})≤ Λ_{h}

µ_{h}(1 +∆tθµ_{h}). (3.53)
Consequently,

V_{h}^{i}^{+1}< Λh

µh

. (3.54)

Whenθ= 1, the value ofI_{h}^{i}^{+1}is positive for any step size. Therefore from
(3.50) we have the following bound

V_{h}^{i}^{+1}(1 +∆tµ_{h})≤V_{h}^{i}+∆tΛ_{h} (3.55)
IfV_{h}^{i}≤ ^{Λ}_{µ}_{h}^{h}, we attain

V_{h}^{i}^{+1}< Λh

µh

. (3.56)

The relationV_{h}^{i}≤ ^{Λ}_{µ}_{h}^{h} is valid for any i. Hence, for the initial total popula-
tions for humans we haveV_{h}^{0}≤ ^{Λ}_{µ}^{h}_{h} implying the invariancy of the number of
the total populations for humans in the interval (0,^{Λ}_{µ}^{h}

h]. As a result, the so-
lutions of the system (3.2) for humans are bounded and positively invariant
in (0,^{Λ}_{µ}^{h}

h]. Similarly, the solutions for mosquitoes are bounded and positively
invariant in (0,^{Λ}_{µ}^{m}

m].

u t