A model of malaria population dynamics with migrants

1.   Introduction

Malaria disease has been a major cause of morbidity and mortality in many regions of the world, particularly in Africa where year after year we see millions of people suffering of malaria. The World Health Organization (WHO) sets targets for control of malaria and its elimination in a number of countries ^[1]. The relevant authorities need to plan the most efficient strategies for malaria control and elimination ^[2]. Understanding the drivers of malaria disease is crucial in this regard. Over the years mathematical modelling has been an important tool in quantifying the effects of the different drivers and informing the optimal combinations of interventions ^{[3,4,5,6,7,8,9,10]}. Human lifestyle and behaviour are crucial factors which, if changed suitably, may assist in curbing the spread of malaria. The movement of people between regions with very different levels of infection, could cause new infection cases that otherwise might not have been ^[11]. In this article we study the effect of infected migrants into a given population for the case of malaria, using a new compartmental model in ordinary differential equations (ODEs) designed for this purpose. An SEIR model was applied to a measles problem in a similar manner in ^[11]. Related studies can also be found in ^[12,13,14], but their models are structurally very different from that in the current paper.

In the literature up to now, a population with two distinguishable components have been approached by way of a 2-group model, see ^[12,15] for instance. Here we follow a much simpler single group model, which allows for migrants not only to enter into, but also to exit from a local population. To the best of the authors' knowledge, the only example of such a model is in the paper ^[11]. For vector-borne diseases such a model has not been published before. In conventional models with inflow of infected such as ^[13,16,17], the immigrants usually become part of the single population. The current model allows the recovered migrants to exit the local population. Our model of a local poulation with migrants is simpler than the two-group models such as those of ^[12,15] and others. This is possible since we do not study the entire group of migrants. Instead, we assume that we know the rate of influx of infected among the migrants, and then we are interested in quantifying its effect on disease transmission in the local population. Thus the model can be applied to test and quantify the effect of disease resulting from infected tourists, infected business travellers, traveling sports people, touring musicians or cultural groups, and similar traveling people from elsewhere, visiting a local population. This is important in particular when the disease in point is well under control in the local population.

We show how a steady flow of migrants into and out of a local population, keeps the number of infections in the population above a certain positive threshold level. The model also explicitly accommodates the indoor residual spraying (IRS). As an application we consider the malaria region of South Africa over the period 1971–2001. The region experienced a very significant upsurge of malaria case numbers towards the end of the said period ^[18]. Our model is utilized to determine the extent to which the upsurge is due to a combination of influx of malaria infected migrants into the region and the effect of IRS. We test different scenarios of infected immigrants, who leave the region subsequently to recovery, and we record the local infections caused by the infected immigrants. Even with limited data, only annual case numbers, we obtain meaningful results.

2.   The model

The basic form of the model of malaria disease dynamics that we introduce in this paper has been popularly used in the literature. In particular, models with inflow of infectives appear in ^[12,13] for instance. Our model has the distinct feature of allowing for also the outflow of migrants.

2.1.   The model description

The model comprises human compartments $S$, $E_h$, $I_h$ and $R$ to represent (resp.) susceptibles, latently infected, infectious and recovered. The vector population has the classes $V$, $E_v$ and $I_v$, being (resp.) the susceptibles, latently infected and infectious vectors. The flow diagram, Figure 1, depicts the movement of individuals into and out of the population, and also the transfer of individuals between the different compartments inside the population.

The parameters of the system are described in Table 1. Among these parameters we assume that:

$K_1, \ K_2,$ $K_3$, $\gamma_0$, $\gamma_1$ and $\gamma_2$ are non-negative,

all the other parameters are positive.

IRS is assumed to increase the mortality rates of the vector by increments $\gamma_i$ in the various compartments (cf. ^[8]). We assume that

This is based on the assumption that the bulk of vector bites that lead to transmission of malaria, will occur indoors, and indoors is where the vectors get exposed to the IRS.

We introduce the numbers $B_1$, $B_2$ and $K = K_0+K_1 +K_2- K_3$, with

$B_1 = K_1/(\rho_h+\mu) \ $ and $\ B_2 = (K_2+\rho_h B_1)/(\sigma+\mu +\delta)$.

We shall fix the value of $K_3$ to be

$K_3 = \sigma B_2$.

Motivation for this choice is given in Remark 5(c). The particular form and function of the numbers $B_1$, $B_2$ and $K_3$ constitute the main novelty of this model and pan out in Theorem 2 and in Remark 5.

We present the following system of ODEs.

The initial conditions are:

The total $N(t)$ of the human population and the total $M(t)$ of the mosquito population satisfy the two identities, for all $t \geq 0$:

$S(t)+E_h(t)+I_h(t)+R(t) = N(t), \ \ V(t)+ E_v(t)+I_v(t) = M(t).$

In the case of influx only into the compartment of susceptibles, the model has a disease-free equilibrium,

We can calculate a threshold ${\cal R}$ for local stability of the disease-free equilibrium, following the next generation method ^[19]. This yields an invariant

This ${\cal R}$ has the property that if ${\cal R} < 1$, then the disease-free equilibrium is locally asymptotically stable, and if ${\cal R} > 1$, then the disease-free equilibrium is unstable. We immediately note that ${\cal R} < 1$ if and only if ${\cal R}^2 < 1$. We define the sets $G_0, G_1$ as below.

with

Theorem 1. (a) $B_1+B_2\leq (K_1+K_2-K_3)/\mu$.

Consider a fixed number $t_1 > 0$. Suppose that $X(t)$ is a local solution for which $X(t) \in G_0$ while $0 < t < t_1$.

(b) If $N(0) \le K/\mu$, then $N(t) \le K/\mu$ for all $0 < t\le t_1$.

(c) If $M(0) \le V_0$, then $M(t) \le V_0$ for all $0 < t\le t_1$.

Proof. (a) We prove that $K_1+K_2-K_3 -\mu(B_2+B_1)\geq 0.$

This proves (a).

Consider any local solution with $X(t)\in {\mathbb R}_+^7$ for all $0 < t\leq t_1$.

(b) We note that

Therefore $N(0) < K/\mu$ implies that $N(t) < K/\mu$ for all $0 < t\leq t_1$.

(c) Also,

By assumption, $\gamma_1-\gamma_0 > 0$ and $\gamma_2-\gamma_0 > 0$. Therefore,

Therefore, (c) follows.

The method of proof by contradiction for positivity of solutions that we follow in Theorem 2, is popularly used, for instance in ^[11,20,21].

Theorem 2 (Positive global solutions). The set $G_1$ is a positively invariant set.

Proof. Note that in view of Theorem 1, it suffices to prove that for paths with initial values in $G_1$, the path will stay inside $G_0$. Consider any point $y\in G_1$. Then there exists a solution $X(t)$ with initial value $X(0) = y$, and at least for some $t_0 > 0$ we shall have $X(t) \in G_0$ while $t\le t_0$. Suppose that, to the contrary of the statement of the theorem, there exists a time value $z$ with $X(z) \notin G_0$. Then the set

is non-empty and has (a finite) infimum, $t_1$. We shall show by contradiction that $t_1 = \infty$. We introduce the following two auxiliary variables:

$J_1(t) = E_h(t)- B_1$ and $J_2(t) = I_h(t)- B_2$.

Let us define the following function for $0\le t < t_1$.

Then each of the terms in the summation is a positive-valued function, and we know that for a constant $q > 0$,

In particular, if any of the variables $S, J_1, J_2, R, V, E_v$ or $I_v$ tends to 0 as $t\to t_1$, then $F_1 \to \infty$, i.e.,

We shall prove that the latter is impossible if $t_1$ is finite, and this will pose a contradiction.

We first note the following:

We calculate the derivative:

After cancellation of some terms which are obviously negative, we obtain an inequality:

$F_1'\leq abI_v + (\rho_h + \mu)+ (\delta +\sigma + \mu) + (\mu + \zeta) +acI_h+(\theta +\rho_v +\gamma_1)+(\theta +\gamma_2).$

Thus $F_1' \leq F_0$, with $F_0$ being the constant

Consequently for any $t\in [0, t_1)$, $0\leq F_1(t)\leq F_0 t_1$. This is the contradiction, completing the proof.

3.   Equilibria

For Theorem 3 below, the work on ^[5] is relevant, as there are no inflow of infections to consider. The proof that we present here is fairly direct and elementary.

Theorem 3. Consider the special case of our ODE model, with $K_1 = K_2 = 0$. If ${\cal R} < 1$, then the disease-free equilibrium is globally asymptotically stable.

Proof. In the proof we find it convenient to use the symbol

Let us suppose that ${\cal R} < 1$, which is equivalent to the inequality:

Then there exists a positive number $z < 1$, for which

We fix such a number $z$. We can also find (and take as fixed) positive numbers $g_3$ and $g_4$ which satisfy the inequality:

Now we proceed to introduce a set of numbers: $\{ g_0, g_1, g_2. ....g_6\}$ that will be used in order to define a Lyapunov function. The candidates $g_3$ and $g_4$ have already been chosen. Furthermore, we let

Now since $0 < z < 1$, we note that $g_1abK_0/\mu < g_6(\theta +\gamma_2)$. Thus we can find a positive number $g_0$ such that:

All of these numbers $g_i$ are positive. Now we define the following function

$F_2 = F_2(S, E_h, I_h, R, V, E_v, I_v)$.

$F_2 = g_0(K_0/\mu -S) +g_1E_h +g_2I_h +g_3R +g_4(L/(\theta + \gamma_0) -V)+g_5E_v+g_6I_v.$

We now show that $F_2$ is Lyapunov at the disease-free equilibrium point. Certainly $F_2$ is positive-definite. We can express $F_2'(t)$ as follows:

$\frac{dF_2}{dt} = Q_0(K_0/\mu-S) +Q_1E_h +Q_2I_h +Q_3R +Q_4(L/(\theta + \gamma_0)-V)+Q_5E_v+Q_6I_v,$

where

$Q_0 = -g_0\mu, \ \ Q_4 = -g_4(\theta+\gamma_0), \ \ Q_3 = -g_0\zeta-g_3(\mu+\zeta),$

$Q_1 = g_2\rho_h -g_1(\mu +\rho_h)$, $Q_5 = g_6\rho_v -g_5(\theta +\gamma_1 +\rho_v),$

$Q_2 = g_3\sigma+(g_4+g_5)acV -g_2(\mu +\delta+\sigma)$, $Q_6 = (g_0+g_1)abS - g_6(\theta +\gamma_2)$.

Clearly $Q_0 < 0$, $Q_3 < 0$ and $Q_4 < 0$. Considering the coefficient $Q_2$, we note that

and $g_2(\mu +\delta+\sigma) = z^3D$. Thus from inequality (1) it follows that $Q_2 < 0$. From the inequality (2) it follows that $Q_6 < 0$. From the choice of $g_6$ it follows that $Q_5 < 0$ and our choice of $g_2$ guarantees that $Q_1 < 0$. Therefore $F_2(t)$ is Lyapunov at the disease-free equilibrium point.

Thus, if ${\cal R} < 1$, then in order to secure elimination of malaria from the population, it is sufficient to stop infected immigrants from entering.

Theorem 4. A unique endemic equilibrium point exists when one of the following conditions (a) or (b) holds:

(a) $K_1+K_2 > 0$,

(b) $K_1+K_2 = 0$ and ${\cal R} > 1$.

Proof. In this proof we shall use the notation:

$\mu_1 = \mu +\delta +\sigma$, $\theta_1 = \theta +\rho_v +\gamma_1$.

The coordinates of the endemic equilibrium point, when it exists, will satisfy the following equations.

Let us write

Using the two different ways of expressing $V^*$, we derive the following identity.

Now we consider the two expressions for $S^*$, which we equate and we obtain an equation without $S^*$. From the latter equation we can eliminate all the variables other than $x$, ending up with a quadratic equation

the coefficients of which are as follows:

$Q_2 = a^2bc\rho_vL(\rho_0\mu_1-\zeta_0\sigma)+ac\rho_0 \mu_1\mu \theta_1(\theta+\gamma_2),$

$Q_1 = -a^2bc\rho_vL(K_0-\zeta_0K_3) - (\rho_0K_2+K_1)[\mu \theta_1(\theta+\gamma_2)ac+a^2bc\rho_vL]$

$+ \rho_0 \mu \mu_1\theta_1(\theta+\gamma_2)(\theta+\gamma_0),$

$Q_0 = - \mu \theta_1(\theta+\gamma_2)(\theta+\gamma_0) (\rho_0K_2+K_1)$.

Note that $\zeta_0 < 1 < \rho_0$ and $\sigma < \mu_1$, and thus, $\rho_0\mu_1- \zeta_0\sigma > 0$. Therefore $Q_2 > 0$.

(a) If $K_1+K_2 > 0$, then $Q_0 < 0$. Since $Q_2Q_0 < 0$, $P$ has a unique positive root which is $I^*_h$, and consequently there is a unique endemic equilibrium point.

(b) If $K_1+K_2 = 0$, then $Q_0 = 0$ and $Q_1 = \rho_0 \mu \mu_1\theta_1(\theta+\gamma_2)(\theta+\gamma_0) - a^2bc\rho_vLK_0$.

Note that $Q_1 < 0$ if and only if ${\cal R} > 1$. Having ${\cal R} > 1$, we obtain $Q_2Q_1 < 0$. Consequently $P$ has a unique positive root which is $I^*_h$, and again there is a unique endemic equilibrium point.

Remark 5.. (a) The state $(E_{\rm mi}, I_{\rm mi}, R_{\rm mi})$ of the migrant sub-population is described by the following system of ODEs:

with initial conditions: $E_{\rm mi}(0) > B_1, \ I_{\rm mi}(0) > B_2, \ R_{\rm mi}(0) > 0.$

(b) Similarly as in Theorem 2 it can be proved that any solution

$(E_{\rm mi}(t), \ I_{\rm mi}(t), \ R_{\rm mi}(t))$

of the system in 5(a) above will be positive.

(c) If instead of taking $K_3 = \sigma B_2$ we venture to assign a higher value to $K_3$, then we find that there will not be an equilibrium point for the system, as it would render $R^*_{\rm mi} < 0$. This does not make sense biologically, and will certainly not be permissible.

If instead of taking $K_3 = \sigma B_2$ we choose a lower value to $K_3$, then $R^*_{\rm mi} > 0$. This means that some of the migrants will stay in the local population permanently. However, the aim of this model is to only consider short-term visitors rather than immigrants staying permanently in the local population. This is the motivation for our assumption: $K_3 = \sigma B_2$. Moreover, this results in a stable equilibrium point for the system of Remark 5(a), see Proposition 6.

The following proposition assures us that the outflow rate $K_3$ will keep the migrant population size stable.

Proposition 6. The equilibrium point $(B_1, B_2, 0)$ of the system in 5(a) is globally asymptotically stable.

Proof. Let us consider the function

$F_2(t) = \mu_1(E_{\rm mi}(t)-B_1)+(I_{\rm mi}(t)-B_2) +R_{\rm mi}(t).$

We can calculate its time derivative and simplify, eventually to obtain

$F_2'(t) = -\mu[E_{\rm mi}(t)-B_1]-(\mu_1+\delta)[I_{\rm mi}(t)-B_2] -\mu R_{\rm mi}(t)$.

Thus $F_2(t)$ is a Lyapunov function, and it follows that $(B_1, B_2, 0)$ is globally asymptotically stable.

4.   Applications

In this section we apply the theory to South Africa. Over the period 1971–1995, South Africa had an annual malaria prevalence which was at a level of well below 14000 cases p.a. on the average ^[18]. The annual prevalence shows a very clear spike in the period 1996–2001. Since malaria is closely correlated with certain climatic variables, such climatic anomalies could have been the reason for this upsurge of malaria, see ^[22] for instance. Another possible reason for the spike is that the IRS might have lost its effect due to mosquitos having become resistant to the insecticide ^[23]. Yet another possibility is that during the period of upsurge of malaria cases, there might have been a heavy influx of people who are symptomatically or asymptomatically infected with malaria. In this section of the paper we show how our model can assist in explaining the upsurge in terms of malaria imports and the waning of the effect of the insecticide. We utilise the model for estimating the average inflow over the period under the assumption that without the inflow, the relevant model parameters (all those other than the $K_i$-parameters) would have the same values as in the foregoing period. In the absence of data or sufficient information we tried a value for $L$ that we hoped would work well. We have worked on an average density of 15 vectors per human.

Note that the parameters $K_1$ and $K_2$ are not listed in Table 2, because computing values for them is part of the computational exercise that follows.

For the purposes of computation we use parameters in line with the paper ^[24]. This excludes of course the parameters pertaining to the population sizes. The essential computations cover two periods:

Period A: 01 January 1971 - 31 December 1995,

Period B: 01 January 1996 - 31 December 2001.

We assume that over the first period, Period A, the inflow of infected were insignificantly low whereas over the period B there was a significant inflow. During the 25 years of Period A, the annual number of cases were well below 14000. During Period B the number of malaria cases were well above 20000 p.a. and in one case even over 40000 p.a. Generally with compartmental models of ordinary differential equations it is assumed that the human population has a maximum size which is constant over time. For this reason we pick the end of 1995 as our reference date, and all the numbers referring to human population will be converted to the 1995 equivalents. We do these conversions before we calculate mean values to feed into the model. The South African population numbers over the entire period is obtained from Worldometer ^[29]. We assume that population growth over the malaria region was proportional to the national population growth of South Africa. In this manner, from ^[18] we calculate the mean number of malaria cases over Period A to be 5773 p.a. and over Period B it is 33451 p.a. The population at risk is calculated in a similar manner starting with its year-2000 value obtained from ^[25]. The 1995-converted value of the population at risk is 3.9942 million. The proportion of mortalities among the malaria cases is estimated from ^[23]. Using the data on Period A, through simulation we determine values of the parameters $a$ and $\delta$ and we calculate numerical values of an equilibrium point that corresponds to an annual number of cases and mortalities as per the data. This equilibrium point is:

$S^* = 3973600; \ E^*_h = 190; \ I^*_h = 2800; \ R^* = 10987;$

$V^* = 59120000; \ E^*_v = 80440; \ I^*_v = 201150,$

and it serves as the initial state for our computations over Period B.

Now we perform some computational analysis on Period B. We test the effect of varying the $\gamma_i$ values, which are the additional rate of deaths in the vector population, induced by the IRS. The graphs are shown in Figures 2 and 3. Note that these two graphs are only exploring the effect of IRS on the malaria case numbers. We do not consider the inflow of infected migrants yet. The baseline values of $\gamma_1$ and $\gamma_2$ were taken as 0.01 of the compartment per day, and resulted in 5772 cases per year. Reduction of the $\gamma_i$ values to 0.008, 0.005 or 0, respectively, yields 6203, 6940 or 8466 cases per year. The relevant $I_h$-trajectories appear in Figure 2. The $I_h$-trajectories over a longer period are shown in Figure 3, in order to give a view of how the curves will stabilize over time.

We now reach the highlight of this application. In what follows, we apply the model for estimating what was the average daily influx of infected migrants. Since the values of the $\gamma$-parameters are not known, we repeat this estimation process for different values of $\gamma_1$ and $\gamma_2$. In Figure 4 we show the results of combining a reduction in the $\gamma$ values with increasing influx of infected in each case to obtain an average annual number of cases equal to 33451 as was observed in the data. The graphs of the $I_h$-compartments are shown in Figure 4 and the combinations are as in Table 2. For the same parameters as for Figure 4, in Figure 5 the $I_h$-trajectories are shown over a period longer than Period B, in order to illustrate how the curves split and eventually stabilize over time. These graphs reveal the dramatic effect that inflow of infectives have on a region where malaria is reasonably well under control. Our attention is now focused on Figure 4 and comparison with the observed data.

Let us consider the first scenario presented in Table 2, having $\gamma_1 = \gamma_2 = 0.01$ per week. Then the model tells us that there must have been an average infux of less than 16 malaria infected migrants per day into the malaria region of South Africa. The precise number per year is $15.75\times 365 = 5749$. The average annual number of observed cases during Period A was 5773 and during Period B it was 33451. This means that the inflow of 5749 migrants resulted in an additional number of cases equal to $33451 - 5773 = 27 686$ p.a. on average.

Now we calculate $27686/5749 = 4.816$. This means that every infected immigrant caused 4 extra local infections (naively interpreted). This is valuable information revealed by the model.

For smaller values of $\gamma_1$ and $\gamma_2$, the corresponding influx of infected immigrants was also smaller. This is well expected, since the infected mosquitos now have a smaller mortality rate. Consequently the number of local cases per infected immigrant will now increase.

These graphs reveal the dramatic effect that inflow of infectives have on a region where malaria is reasonably well under control.

5.   Conclusions

The model that we present in this study is of a completely new class. The model accommodates both inflow of infected migrants and outflow of the surviving recovered migrants. Application of such a model is versatile and the application in this paper is relatively simple. We have illustrated with an example, the process of estimating the rate of influx that was responsible for an observed number of malaria cases in the South African malaria region during the late 1990's. It is likewise possible to compute the expected number of cases that will arise from a given influx of infected migrants. One can equally well make future projections of local cases if the rate of importation of cases is known. Already there is an application to measles ^[11] in this regard. Our application gives fairly clear results on the effects of infective immigrants on malaria epidemiology in South Africa, even though the input date was limited only to annual case numbers. Accuracy of the results can be improved if the initial state of the system is known. In this case we simply use equilibrium values of the foregoing period as our initial state. Further improvements in the accuracy can be achieved by using a climate-driven model such as ^[9,10]. However in the latter case the mathematics around inflow and outflow of migrants will be much more complex, and climatic data or predictions will be required.

Acknowledgements

We are grateful to the Reviewers for all the constructive comments, which resulted in a much improved paper.

Conflict of interest

The authors have no competing interests.