An Epidemiological Model of Malaria Accounting for Asymptomatic Carriers

Asymptomatic individuals in the context of malarial disease are subjects who carry a parasite load, but do not show clinical symptoms. A correct understanding of the influence of asymptomatic individuals on transmission dynamics will provide a comprehensive description of the complex interplay between the definitive host (female Anopheles mosquito), intermediate host (human), and agent (Plasmodium parasite). The goal of this article is to conduct a rigorous mathematical analysis of a new compartmentalized malaria model accounting for asymptomatic human hosts for the purpose of calculating the basic reproductive number (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_0$$\end{document}R0) and determining the bifurcations that might occur at the onset of disease-free equilibrium. A point of departure of this model from others appearing in the literature is that the asymptomatic compartment is decomposed into two mutually disjoint sub-compartments by making use of the naturally acquired immunity of the population under consideration. After deriving the model, a qualitative analysis is carried out to classify the stability of the equilibria of the system. Our results show that the dynamical system is locally asymptotically stable provided that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_0<1$$\end{document}R0<1. However, this stability is not global, owning to the occurrence of a sub-critical bifurcation in which additional non-trivial sub-threshold equilibrium solutions appear in response to a specified parameter being perturbed. To ensure that the model does not undergo a backward bifurcation, we demand an auxiliary parameter denoted \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varLambda <1$$\end{document}Λ<1 in addition to the threshold constraint \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_0<1$$\end{document}R0<1. The authors hope that this qualitative analysis will fill in the gaps of what is currently known about asymptomatic malaria and aid in designing strategies that assist the further development of malaria control and eradication efforts. Electronic supplementary material The online version of this article (10.1007/s11538-020-00717-y) contains supplementary material, which is available to authorized users.


Introduction
Malaria is one of the most lethal and complex parasitic diseases in the world World Health Organization 2019). Throughout human history, malaria has burdened most regions of our planet and has had a profound impact on human history and evolution; for example, it has been credited for contributing to the decline of the Roman Empire (Sallares 2002, p. 14). In 2018, the World Health Organization (WHO) reported the occurrence of approximately 228 million new cases of malaria (range 206-258 million), which resulted in an estimated 405 thousand disease-induced deaths (World Health Organization 2019). It was estimated that 67% of these fatalities were experienced by children under the age of five.
The life cycle of the Plasmodium parasite can be broken down into two separate subcycles: the asexual cycle, occurring in humans (intermediate host) and the sexual cycle in mosquitoes (definitive host), in which maturity is reached. The sexual cycle begins when a susceptible mosquito feeds on the blood of an infectious human, ingesting sexual forms of the Plasmodium parasite previously developed in the human body, known as gametocytes. While in the midgut lumen of the mosquito, these gametocytes fuse to form diploid zygotes, which grow into elongated ookinetes. The motile ookinetes burrow into the outer membrane of the mosquito midgut and form ellipsoid shaped oocysts. Eventually, the oocysts rupture releasing thousands of haploid forms called sporozoites (Rosenberg and Rungsiwongse 1991). These sporozoites accumulate in the salivary glands of the mosquito, causing it to become infectious.
The asexual cycle begins when an infectious mosquito bites the host and injects saliva with anticoagulant agents that keep the wound open, thus allowing a blood meal and simultaneously injecting sporozoites into the skin (Cowman et al. 2012). The sporozoites travel through the blood vascular system to the liver where they invade the cells of the liver, known as hepatocytes. Inside the human, a Plasmodium infection goes through two cycles: a initial liver (hepatic, or exo-erythrocytic) stage lasting a few days, followed by a blood (erythrocytic) stage that lasts until the host clears naturally the infections, receives treatment, or dies.
The hepatic stage begins in the hepatocytes, where a proportion of the sporozoites undergo a process called pre-erythrocytic or hepatic schizogony, in which they multiply asexually to produce thousands of haploid daughter cells, known as merozoites. During this process, schizonts are formed, causing the hepatocytes to rupture. This allows the merozoites to enter the bloodstream.
The erythrocytic stage begins when free-floating merozoites invade erythrocytes in a matter of minutes. Inside the erythrocyte, parasites enter the ring stage in which some mutate into an enlarged ring-shaped form called trophozoites that mature into schizonts, causing the cell to burst and releasing more merozoites into the bloodstream. In the case of P. falciparum, this process of invasion and rupture of RBCs occurs synchronously every 48 h (Mideo et al. 2013). At this point in the sub-cycle, a small portion of the merozoites, for reasons incompletely understood, develop into gametocytes, or sexual forms; however, no sexual reproduction occurs inside the host (John et al. 2006;Yan et al. 2015). The infected human host is now ready to infect new susceptible mosquitoes, thus completing the Plasmodium life cycle.
The blood stage parasites are responsible for most clinical symptoms associated with the disease Yazdani et al. 2006). The periodic rupturing of the RBCs results in the release of various debri and waste products, which in turn activate the immune system and cause symptoms such as chills, fatigue, pain, and fever. The average duration for the infection of an RBC is dependent on the Plasmodium species. P. falciparum has the interesting pathological effect of sequestration, which occurs when infected RBCs containing mature forms of the parasite, i.e., trophozoites and schizonts, adhere to the walls of small diameter blood vessels, e.g., the endothelium of capillaries and venules (David et al. 1983). As a result of sequestration, the microcirculation is reduced and in some cases inflammatory processes take place. One of the common complications of this sequestration is cerebral malaria, which might cause patients to sustain brain injury, resulting in long-term neuro-cognitive impairment (MacPherson et al. 1985).
A human host is called asymptomatic when it is a carrier for the Plasmodium parasite, but displays no clinical symptoms. Asymptomatic carriers contribute to gametocyte circulation by providing a hidden reservoir for the parasite to take refuge. As pointed out by Laishram et al. (2012b), asymptomatic infections often go undetected, resulting in a major source of gametocytes for local mosquito vectors. Accordingly, asymptomatic carriers contribute to the persistence of malaria transmission within their localized populations (Bousema et al. 2004). Frequent exposure to the Plasmodium parasites leads to naturally acquired immunity to the symptoms of the disease, but not necessarily to the parasite, and as a result, it creates asymptomatic carriers in a given population (Staalsoe and Hviid 1998).
Asymptomatic malaria infections have been reported in various high and intermediate transmission areas such as Kenya and Nigeria (Bousema et al. 2004;Eke et al. 2006). Recently, asymptomatic infections have been reported in relatively low endemic areas such as Colombia, Ecuador, and the Amazonian region of Brazil (Cucunubá et al. 2008;Coura et al. 2006;Sáenz et al. 2017;Lopez-Perez et al. 2015). There is much evidence that asymptomatic malaria infections play a fundamental role in malaria transmission (Lindblade et al. 2013). Disease transmission dynamics may be affected by the amount of asymptomatic carriers in a given population over a specified time interval. Indeed, a positive correlation between high transmission and high asymptomatic prevalence has been reported in Nigeria, Senegal, Gabon, and the Amazonian regions of Brazil (Alves et al. 2002;de Andrade et al. 1995;Dal-Bianco et al. 2007;Eke et al. 2006).
In accordance with Bruce-Chwatt et al. (1980), we define malaria immunity as the state of resistance to the infection brought about by all processes which are involved in destroying the Plasmodia or limiting their multiplication. Natural innate immunity is an intrinsic property of the host. This type of immunity is characterized by an immediate inhibitory response to the introduction of the parasite which is independent of any previous infection. Their are two types of acquired immunity, namely active acquired immunity and passive acquired immunity. Active acquired immunity is defined as an enhancement of the hosts defense mechanism due to previous contact with the pathogen. Passive acquired immunity is characterized by either the mother to child transfer of protective antibodies in the pre-or post-natal developmental periods or by the injection of such antibodies. In this work, we are focused on active acquired immunity, more specifically a type of immunity that is acquired through means of exposure.
Humans experience various kinds of active acquired immunity which provide different kinds of protection. Adopting the definitions by Doolan et al. (2009), we define protection to be objective evidence of a lower risk of clinical disease, indicated by the absence of fever, that is, the oral temperature does not exceed the threshold 37 • C (Clark and Kruse 1990). Anti-disease immunity is conferred protection against clinical disease, which affects the overall risk and extent of morbidity associated with a given parasite density. Anti-parasite immunity is conferred protection against parasitemia, which affects the parasite density. Premunition provides protection against new infections by maintaining a generally asymptomatic parasitemia (Koch 1900;Sergent and Parrot 1935). In this article, we make use of a kind of premunition called naturally acquired immunity (NAI). As reported by Doolan et al. (2009), in holoendemic regions across sub-Saharan Africa, most people are continuously infected by P. falciparum, while the majority of infected adults rarely experience observable disease. This valid protection against infection is NAI corresponding to P. falciparum.
Mathematical models of malaria transmission have been studied by various authors; for a survey, we refer the reader to Mandal et al. (2011). Model formulations where partially immune humans are allowed to be infective have been studied by Ngwa and Shu (2000) and Roop-O et al. (2015). Both of these models accounted for partially immune humans in the recovered compartment. Asymptomatic malaria has also been previously modeled and studied. Of recent, an asymptomatic malaria model was introduced by Filipe et al. (2007). This model depends on a state-invariant control parameter φ, which stands for the proportion of human infections that develop disease. After letting 1/h denote the mean latent period in humans, the progression rates from the exposed to the symptomatic and asymptomatic classes were defined to be the products hφ and h(1 − φ), respectively. Additionally, asymptomatic humans were included in the recovered compartment.
In this article, we depart from the previous models of asymptomatic malaria by creating explicitly different compartments for symptomatic (Y ) and asymptomatic ( A) subjects, which in addition to susceptibles (S), exposed (E), and recovered (R) yields the acronym SEYAR. Another important point of departure with respect to previous models of asymptomatic malaria is that in the SEYAR model (4), the progression rates corresponding to the symptomatic and asymptomatic human classes are nonlinear functions of the time-dependent exposed proportion. Therefore, the SEYAR model does not fall into a sub-class of such models currently appearing in the literature. Moreover, we do not include the asymptomatic humans in the recovered compartment. This allows an effective isolation of the effect that asymptomatic carriers have on the disease transmission dynamics. Unlike other models appearing in the literature, the recovered human compartment studied in this article does not have an associated transmission probability, since it does not contribute to mosquito infection.
The derivation of the SEYAR model (4) hinges on a specific decomposition of the infected human compartment into two mutually disjoint sub-compartments accounting for asymptomatic and symptomatic carriers. This decomposition is accomplished by making use of a nonlinear exposure-dependent NAI function, which is the solution to the initial-value problem (1) derived in Sect. 2. Although there is much in the literature, the proper inclusion of asymptomatic carriers into the epidemiological modeling of malaria warrants a formal mathematical understanding.
This manuscript provides a new malarial model accounting for asymptomatic human hosts in terms of the NAI of the population under consideration and is organized as follows: Sect. 2 presents the model formulation, Sect. 3.1 covers the issue of well-posedness of the initial-value problem and provides an analysis of the total population dynamics, Sect. 3.2 contains a rigorous study of the local asymptotic stability of disease-free equilibrium (DFE) for the model with a mathematical and epidemiological interpretation of the reproductive threshold, Sect. 4 introduces some modifications of the SEYAR model along with their corresponding reproductive thresholds and addresses the impact of the asymptomatic class on the reproductive threshold of the original model, Sect. 5 is focused on nonlinear stability analysis and provides a classification parameter in which its size determines the type of bifurcation undergone by the dynamical system, Sect. 6 incorporates control measures into the dynamical system, Sect. 7.1 is focused on a sensitivity analysis of the reproductive threshold arising from the model, Sect. 7.2 consists of numerical results corresponding to the following three high transmission sites: Kaduna in Nigeria, Namawala in Tanzania, and Butelgut in Papua New Guinea, Sect. 8 consists of a summary of the results contained in Sects. 2-7.2 and a discussion regarding future direction and extensions, "Appendix A" contains formal proofs of the lemmas and theorems contained in Sects. 2-6, "Summary of Stability Theorems" of appendix is a summary of the main stability theorems used in the investigation of the local asymptotic stability of the equilibrium solutions studied in Sects. 3.2 and 5, and "Parameter Values" of appendix contains tables of numerical rates corresponding to the high transmission sites studied in Sect. 7.2.

Methods: Model Formulation
The formulation of the SEYAR model for the spread of malaria in the human and mosquito populations begins with dividing the total host-vector population into two compartments, denoted by N H (t) and N M (t), which stand for the total population sizes of the humans and mosquitoes, respectively, at a given time t. From this point on, whenever implied by the context of the discussion, the time t dependency is suppressed. Assuming a homogeneously mixed host population, we further decompose the compartments into the following five epidemiological classes: susceptible human S, exposed human E, symptomatic human Y , asymptomatic human A, and recovered human R, so that N H = S + E + Y + A + R. For simplicity of exposition, the state variable is identified with its corresponding class. For example, when we are considering a human from class A, it is understood that A is a function and not a class, in general.
A point of departure from the usual SEIR models, as studied by d' Onofrio (2002), Li and Muldowney (1995), Li et al. (1999), Roop-O et al. (2015) and Smith et al. (2001), resides in the mutually disjoint partitioning of the infected compartment into two sub-compartments, labeled asymptomatic A and symptomatic Y , i.e., For the mosquito population, we have the following three classes: susceptible mosquito M S , exposed mosquito M E , and infected mosquito M I . Accordingly, the total mosquito population is given by N M = M S + M E + M I .
As mentioned by Filipe et al. (2007), it is known that the infection rates between human and mosquito populations depend on numerous factors including the human biting rate of the mosquito σ (which is the number of bites per mosquito), transmission probabilities (to be later defined), and the number of infectious and susceptible of each species involved. Furthermore, we assume that the average number of mosquito bites suffered by humans depends on the total sizes of their respective populations in the community. As a result, the number of bites per human is σ N M N H . Therefore, the force of infection from mosquitoes to humans λ SE is defined to be the product of the number of bites per human, the transmission probability β M from a mosquito in the class M I to a human in class S, and the probability that a mosquito is infectious M I N M , i.e., Even though asymptomatic carriers might not get clinically ill, they still could harbor low levels of gametocytes in their bloodstreams and could to pass the infection onto mosquitoes (Vinetz and Gilman 2002). When a mosquito from class M S bites a human from class Y , the force of infection ω Y is defined as the product of the number of bites per mosquito σ , the transmission probability β Y from a human in Y to a mosquito in M S , and the probability that a human is in the symptomatic class Y N H . When a mosquito from the class M S bites a human from class A, the corresponding force of infection ω A is the product of the number of bites per mosquito σ , the transmission probability β A from a human in A to a mosquito in M S , and the probability that a human is in the asymptomatic class A N H . As pointed out by Laishram et al. (2012a), the parasites carried by asymptomatic hosts can be more infectious than those of symptomatic hosts. One could assume that a typical asymptomatic carrier has a higher NAI level than a symptomatic, so that β A ≤ β Y . Accordingly, the force of infection from humans to mosquitoes V SE is defined to be the sum of the forces of infection corresponding to the humans in classes Y and A, i.e., Let ν E I = τ , where τ is the reciprocal of the mean duration of the definitive host latent period. At a given time t ∈ R + , an individual's experience of malaria is dependent upon the degree of naturally acquired immunity that he or she has gained. Effective antiparasitic immunity is achieved only after many frequent infections (Carter and Mendis 2002;James et al. 1920;Macdonald et al. 1957). This important epidemiological observation, in combination with the discussion regarding naturally acquired immunity (NAI) in Sect. 1, implies that the rate of progression λ E A from the exposed class E to asymptomatic class A depends on the proportion of human individuals receiving sufficient protection from the average NAI accumulated in the population with respect to natural exposure.
Let u(t) denote the proportion of the human population fully protected by NAI. Since naturally acquired immunity to the Plasmodium parasite is acquired and accumulates over time in response to frequent exposure, the rate that this proportion of protected individuals changes depends on the rate that the human population is being exposed, up to a threshold value. To uncover this exposure dependency, firstly let the lower and upper protected proportion thresholds be given by u low and u high , respectively. It should be noted that 0 ≤ u low < u high < 1.
Let ε := E N H , upon assuming that the initial NAI protected proportion is given by the lower threshold u(0) := u low , these epidemiological principles lead to the following initial-value problem (IVP) being posed (1) Notice that the time derivative on both sides of the differential equation (1) is an unusual formulation. It is explained by the fact that the rate of the human population fully protected by NAIu is dependent upon the rate that the population is exposed to the pathogenε. For example, ifε increases or decreases, then it stands to epidemiological reason thatu should increase or decrease.
By making use of the integrating factor L = e t 0ε (s)ds := e ε−ε 0 , it follows thaṫ Upon rearranging terms and invoking a slight abuse of notation, to emphasize the exposure ε dependency of u, the solution is represented by the following equation: In the above, the symbol ε 0 := ε(0) stands for the initial exposed proportion of the human host population. The function N H is positive for all t ∈ R + ; thus, it directly follows that ε ∈ C 1 b (R + ). In general, care should be taken to ensure that the progression rate is mathematically well defined and epidemiologically sensible, i.e., a singularity should not arise and it should be non-negative. Provided that x is such that ε ∈ C 1 (R + ), then clearly the progression rate will not experience a singularity.
Besides the trivial singularity issue covered above, care must be taken to ensure the non-negativity of the nonlinear progression rate u. As a result, the solution to the IVP (1) implicitly imposes an additional constraint upon the initial data of the SEYAR model (4). The constraint is listed below in Lemma 1, which is proven in "Appendix." Theorem 1 (Initial data constraint for the SEYAR model) Let ϑ be defined as follows: To ensure the non-negativity of u, the initial data are required to satisfy the inequality Let the eight-dimensional vector of functions x = (S, E, Y , A, R, M S , M E , M I ) T be such that ε ∈ C 1 (R + ) and ϑ be defined as in Theorem 1. The nonlinear progression rate u is well defined in a mathematical and epidemiological sense.
Define λ E A = γ u(ε) where γ is the reciprocal of the mean duration of the human latent period. It is a direct consequence that λ EY = γ (1−u(ε)), so that λ E A +λ EY = γ . Since the naturally acquired immune proportion will grow in response to exposure, the rate of progression from the exposed class E to asymptomatic class A should increase, warranting the choice of λ E A . Furthermore, as the exposure rate increases, λ E A will eventually be maximized. This is consistent with the observation that the average amount of asymptomatic human hosts in a population should increase after frequent exposure over a sufficient time period. When the exposed proportion is equal to zero over a prescribed time interval, it follows that u = e ε 0 (u low − u high ) + u high over the interval. This quantity is a sum consisting of the upper threshold and a negative scaled difference of the lower and upper thresholds. If this infimum is achieved, then the progression rate λ E A will be minimal. This is due to the fact that if there is little exposure, then there is little NAI developed in the population, so that the rate of progression from E to Y will be maximal. Moreover, the progression rate from E to Y should decrease as the exposure rate increases.
If the population under consideration is free of malaria, then no immunity should be present and it is possible that u high = 0. In this scenario, we set u low = 0. As a result, u = 0 for all time and it follows that λ E A = 0. Thus, the progression rate from the exposed human class to the asymptomatic class is nullified. However, as the purpose of this article is to model asymptomatic malaria, it is assumed that malaria is present in the population.
These assumptions give rise to the following SEYAR model IVP (4), depicted in Fig. 1, describing the dynamics of malaria disease transmission in the human and mosquito populations: This figure is a schematic diagram of a malaria model including an asymptomatic compartment. The solid lines represent progression from one compartment to the next, while the dotted stand for the human-mosquito interaction. Humans enter the susceptible compartment either through birth or migration and then progress through each additional compartment subject to the rates described above (Color figure online) where u and ϑ are defined by Eqs. (2) and (3), respectively. For convenience, the model parameters are summarized in Table 1. All of the parameters are strictly positive, with the exception for the disease-induced death rate δ, which is allowed to be non-negative. The naturally acquired immune proportion will increase or decrease, depending on the rate that the human population is being exposed. Hypothetically, as the rate of exposure increases or decreases, this proportion should grow or shrink up to a threshold value. This epidemiological behavior is quantified by the solution to the IVP (1). Assume u low = 0.5, as above, and that it is possible for at most ninety percent of the population to acquire sufficient protection through means of natural exposure, i.e., u high = 0.9. Then, the initial exposed proportion ε 0 can be assumed to take any value in the interval In other words, one can assume at most 81% of the human population to be initially exposed. On the other hand, if one modifies the above assumptions so that the initial naturally acquired immune proportion is u low = 0.1, then ε 0 ∈ [0, ln 9 8 ] ≈ [0, 0.12], so that at most 12% of the human population can be assumed to be initially exposed.

Well-Posedness and Feasible Region
Although assuming that Φ ∈ C 1 provides sufficient regularity to ensure that system (5) is well posed, this work is primarily concerned with the stability of the system near equilibrium points. This requires additional regularity assumptions on the vector field Φ in order to invoke a variation of the center manifold theorem, proven by Castillo-Chavez and Song (2004). Consequently, from now on, it is necessary to assume that Φ ∈ C 2 ⊂ C 1 , that is, it is at least twice continuously differentiable. Moreover, to be reasonable in an epidemiological sense, the functions under consideration should posses a bounded first derivative, that is, they should be members of the class C 1 b (R + ). From this point on, the model is studied in the more regular (smaller) function space . In light of the mathematical and epidemiological well-posedness of the IVP (1) and structure of the underlying vector field, the additional C 2 regularity will be inherited by the solution x.
To this end, let and we rewrite (4) in the following compact form: Theorem 2 (Existence theory of the SEYAR Model) There exists a sufficiently regular unique solution x to the SEYAR model IVP (5) that can be continued to a maximal time interval. Additionally, x depends continuously on the initial data x 0 and model parameters involved.
For the proof of the above theorem, the reader is referred to "Proof of Theorem 2 (See Page 11)" of appendix. The dynamics of the total population are given by the following decoupled system: Due to the non-homogeneous term δY , the asymptotic behavior of the human population is more delicate. For the human population, we have the following theorem which provides a tighter lower bound on the attracting region for the model. In the absence of infection, the long-time behavior of N H is trivial, that is, the total human population converges to the equilibrium population density, as in the mosquito population. In the case of an infectious disease with disease-induced death rate as in Eq. (6), one would expect the population to converge to a smaller quantity, as there is a disease-induced death rate δ adding to the natural death rate ξ H . The long-time population size will be smaller, as it has to account for the additional disease-induced deaths suffered by symptomatic individuals. Listed below is the theorem concerning the feasible region of model (6).
Theorem 3 (Feasible region of the SEYAR model) Let (N H , N M ) be the solution of system (6) emanating from Theorem 2, with corresponding initial data (N H (0), N M (0)) ∈ R 2 + . Define the following compact sub-space Then, Γ is a forward invariant attractor for system (6).
For the proof of the above theorem, the reader is referred to "Proof of Theorem 3 (See Page 12)" of appendix. Mathematically speaking, Theorem 3 reduces the complexity of the analysis involved regarding the long-term dynamics of the system by allowing the replacement of a potentially unbounded (epidemiologically unreasonable) space with the smaller (epidemiologically reasonable) compact sub-space Γ . If we let L t = e ξ H t denote an exponential multiplier, then the solution for the total dynamics of the human population is given by The instantaneous rate of occurrence of death, i.e., the force of mortality ξ H , is assumed to be constant. As a result, the probabilities of living and dying up to t days are given by e −ξ H t and 1 − e −ξ H t , respectively. Assuming that N H (0) = 0, the above solution says that the total population N H (t) at time t is given by the weighted product of the carrying capacity Ω H ξ H and the distribution of humans that are left after those that have died due to natural causes 1 − e −ξ H t , minus a weighted average of humans that have died due to symptomatic infections. The later quantity is captured by the non-homogeneous forcing term − δL −t * Y (t), given by a convolution with the inverse multiplier. Since convolution is a smoothing operation, this emphasizes the fact that we are subtracting a "smoothing average" over past time of the humans which have died from symptomatic infections. A straightforward calculation yields the following differential inequality: From an epidemiological view point, the above inequalities imply that if the total population N H breaches its carrying capacity, then the weighted average of fatal symptomatic infections must increase to stabilize the population back to a healthy level.
In an epidemiological setting, one can define the term δ Y ∞ to be the potential maximum disease-related death of the human population. If Ω H > δ Y ∞ , then it directly follows that α = Ω H −δ Y ∞ ξ H . In this case, the above theorem provides a tighter lower bound for the feasible region corresponding to the dynamical system (6). In practice, the disease-induced death rate for humans δ is sufficiently small (and in some cases negligible), so that this inequality is satisfied.
In the previous variants of malaria models appearing in the literature, e.g., SIR, SEIR, SEIRS, etc., the quantity 0 is listed as the lower bound; however, this is unreasonable since the populations under consideration usually do not go extinct, unless Ω H = δ Y ∞ . If the maximum impact the disease is capable of having on the population is less than the recruitment rate, then their will always be accumulation over long time. For other members of the SIR model class, the lower bound would be the same except with Y replaced by I . In the case of an infectious disease, there will be additional disease-induced deaths, so that the total human population will not converge to the equilibrium population density.

Reproductive Threshold and Disease-Free Equilibrium
This section is focused on deriving a threshold value that characterizes the local asymptotic stability of the underlying dynamical system (4). This value, called the reproductive threshold, provides a way to estimate the reduction in transmission intensity required to eliminate malaria through vector-based control (Smith et al. 2007). The basic reproductive number R 0 corresponding to a given model is a threshold value which represents the average amount of new infections produced by a typical infectious individual in a completely susceptible population, at a disease-free equilibrium. This quantity is equal to the reproductive threshold for a class of simplified population models.
Disease-free equilibrium (DFE) points are solutions of a dynamical system corresponding to the case where no disease is present in the population. Define the diseased classes to be E, Y , A, M E , and M I . Notice that R is not considered to be a diseased class, as the asymptomatic class A has been effectively removed, cf. Roop-O et al. (2015). As a result, individuals in the R compartment are considered to be temporarily immune, but not infectious. After determining the DFE of system (4), this threshold value is used to address its local asymptotic stability. Upon equating the right-hand side of (4) to zero and solving, we arrive at the following unique DFE: Lemma 1 (Local asymptotic stability of the DFE for the SEYAR model) Define the following quantity: where U low := e ε 0 (u low − u high ) + u high . Then, the DFE x d f e for the SEYAR model (4) is locally asymptotically stable provided that R 0 < 1 and unstable if R 0 > 1 in Γ .
For the proof of Lemma 1, the reader is referred to "Proof of Corollary 1 (See Page 18)" of appendix. A verification of the reproductive threshold R 0 (7) is provided in the electronic supplementary material. Lemma 1 is proven by utilizing the next-generation method developed by Van den Driessche and Watmough (2002). The threshold value (7) has major epidemiological implications on the underlying dynamical system (4). To gain a deeper insight into the qualitative information encoded in this important quantity, we decompose it in the form of an epidemiological meaningful product in order to analyze each factors contribution: Due to the above lemma and in an epidemiological setting, it is desirable to have the reproductive threshold below unity. Listed below is the size contribution and biological description for each of the factors. The first factor is σ , which stands for the human biting rate. This factor is usually much less than unity due to the fact that female anopheline mosquitoes generally transmit fewer than 100 sporozoites per bite (Ponnudurai et al. 1991). As malaria is a mosquito borne disease, the agent Plasmodium will spread at a much slower rate, provided less vectors are introducing it into human hosts. Owning to the monotonicity of the square root function, it is sufficient focus on the size of each r i .
i The term r 1 = Ω M Ω H is the ratio of the mosquito and human recruitment rates. The population density of anopheles mosquitoes is a function of proximity to breeding sites. A study in the urban settings of Dakhar in 1998 (Trape et al. 1992) reports an average mosquito population density of 10.4 mosquito/room in general, and specifically 1.6 Anopheles/room during the rainy season; these counts correspond to mosquitoes collected after indoor pyrethrum spray and do not account for outdoor mosquito activity. These numbers seem to indicate that the population of Anopheles mosquitoes is larger than the population of humans in that urban setting. Likewise, a study in a rural area of Sudan in a 4 km 2 area inhabited by approximately 800 people reported a mosquito population in the range 135,000-330,000 mosquitoes (Costantini et al. 1996). Combined, these studies suggest that in an endemic site, the quotient r 1 is greater than unity. As the human and mosquito recruitment rates rank high on the sensitivity hierarchy of many epidemic models appearing in the literature, it is no surprise that this term is problematic with respect to the overall size of the threshold. ii The term r 2 = τ τ +ξ M is the reciprocal of the mosquito latent period τ divided by itself plus the mosquito mortality rate ξ M . This quantity is bounded above by one and is monotonically decreasing with respect to ξ M . It follows that the larger the mosquito death rate is, the smaller r 2 will be. iii The first fully human-dependent term r 3 = γ γ +ξ H is comprised of the reciprocal of the human latent γ period divided by itself plus the human mortality rate ξ H . As in the case of r 2 , this quantity is always less than one and monotonically decreases with respect to the human mortality rate. This is consistent with the fact that the fewer hosts there are for the parasite to invade, the less infections will arise. However, since increasing ξ H is not practical, this terms offers no control over R 0 . iv The term r 4 = ξ H ξ M is the ratio of the human and mosquito death rates. This particular quantity is always less than one as the mosquito death rate is much higher than the human death rate. v The term r 5 = β M ξ M is the ratio of the mosquito-to-human transmission probability β M and the mortality rate of the mosquito population ξ M . This quantity will be less than one provided β M < ξ M . In the case of a population with a relatively high vector transmission probability, the vector death rate must be large enough to make r 5 less than one. This implies a restriction on the size of β M . It will be shown in Sect. 5 that if β M breaches a certain threshold, then sub-threshold endemic equilibria can emerge. vi The second fully human-dependent term is given by the following equation: The quantity r 6 is a difference of ratios consisting of asymptomatic and symptomatic vital dynamics, along with transmission, recovery, and disease-induced death rates, weighted with a distribution consisting of the lower and upper NAIrate threshold of the human population scaled by an exponentiation of the initial exposure rate. For simplicity of exposition, we assume that the initial exposed proportion is zero, i.e., ε 0 = 0, so that U low = u low . If the initial exposed proportion is not equal to zero, then initially, U low < u low and the reproductive threshold will be larger under parameter configurations to be specified shortly. Under such a configuration, upon initial exposure, there will be more symptomatic individuals, but as exposure increases, less humans will die as naturally acquired immunity will begin to develop in the overall population. However, the following discussion is unaffected by this minor detail. As a result, the human-dependent factor is taken to be: Let the low and high thresholds be such that u low > 0 and u high < 1, respectively, and define the following compact sub-interval T u := [u low , u high ] ⊂ (0, 1) consisting of various NAI protected proportions corresponding to a given population. When subjected to a certain parameter restriction, the sizes of the quantities u low = u(0) and r 6 are inversely related, i.e., the larger u low is, the smaller r 6 will be, resulting in a relatively smaller R 0 . Hence, an additional way to control the size of R 0 arises, provided the parameters are such that this inequality restriction, to be mentioned below, holds. However, as we will see, if the parameters are such that this inequality reversed, then they are directly related.
λ Y R +ξ H +δ and define T ⊂ R + ∪ {+∞} to be an ordered subset of the nonnegative extended real numbers. In addition, let u solve Eq. (2) and denote {u(t) ∈ T u |t ∈ T } to be the set of NAI protected proportions experienced by a population over a prescribed, possibly infinite, time interval indexed with t and consider the following function: Subject to a given parameter configuration, the type of monotonicity obeyed by R 0 (u(t)) is dependent on the sign of the nonzero combination of parameters C 1 − C 2 . These observations motivate the following definition which classifies the configuration space of the model based on how the asymptotic dynamics of the corresponding population responds with respect to the NAI accumulated in response to the rate of exposure.
Upon the trivial case that C 1 − C 2 = 0, the system is said to have a null-configuration. Additionally, we refer to a system possessing such a configuration as either A-, Y -, or null-configured. A given human population is called A-,Y -, or null-dominant, provided that its corresponding configuration is.
Although r 6 is always positive, the combination of parameters C 1 −C 2 need not be. The positiveness of r 6 is a result of how the transmission probabilities are weighted by the NAI protected proportion. If the system is null-configured, then R 0 (u(t)) := C 0 √ C 2 and it is locally asymptotically stable, provided that C 0 < 1 and β Y < λ Y R + ξ H + δ, i.e., a high vector death rate, and the symptomatics are recovering or dying out at a faster rate than they are transmitting. Consider the case of an A-dominant configuration, i.e., As mentioned in Sect. 2, in the formulation of the SEYAR model (4), we do not assume an ordering on the asymptomatic and symptomatic transmission probabilities β A and β Y , respectively. However, provided that asymptomatic carriers transmit at a lower rate than that of symptomatic, the left-hand side of the above inequality is less than unity. Moreover, if asymptomatic individuals recover faster than symptomatic, provided the disease-induced death rate δ and recovery rates are such that λ AR > λ Y R + δ, then the right-hand side of the above inequality is greater than unity and inequality (8) is satisfied. In an epidemiological setting, such a configuration corresponds to holoendemic regions across sub-Saharan Africa. In such regions, the majority of people are continuously infected by P. falciparum, but only a small proportion display clinical symptoms (Águas et al. 2008). The high level of naturally acquired immunity present in the population allows them to live their daily lives feeling healthy despite a relatively high blood-parasite density (Doolan et al. 2009).
Analytically speaking, in the case of an A-dominant configuration, provided that the mosquito mortality rate ξ M can be made large enough so that the term C 0 compromised of fractional multipliers is sufficiently small, then R 0 (u t ) achieves its maximum R 0 (u low ) at the low NAI threshold u low and decreases to its infimum as the high NAI threshold u high is approached. Moreover, since the ordered set T is a subset of a separable metric space, we can extract an ordered countable subset and form the partition {u low = u(t 1 ) < · · · < u(t n ) < u(t n+1 ) < · · · } of the compact subinterval T u . In this case, the corresponding values of R 0 (u(t)) obey the following descending order {R 0 (u low ) > · · · > R 0 (u(t n )) > R 0 (u(t n+1 )) > · · · }. This is consistent with the fact that as a given population acquires natural immunity through exposure, the disease will start to spread at a slower rate. Conversely, if the system has a Y -dominant configuration, i.e., C 1 − C 2 > 0, then the monotonicity is reversed.
In mathematical terminology, provided σ < 1, one can always find a large enough ξ M to make C 0 < 1. In this case, the factor that will cause R 0 to breach unity is r 6 . From an epidemiological standpoint, regardless of the size of vector transmission probability β M or human biting rate σ , if enough mosquitoes are dying to significantly slow the disease transmission dynamics, then C 0 < 1 and, as a result, the size of the reproductive threshold R 0 will be determined by r 6 which depends on the human immune systems response to the Plasmodium parasite. This attests to the fact that in such vector transmitted diseases, the reproductive threshold will be lower in response to vector elimination up to a point and the factor allowing the disease to persist under This figure is a schematic diagram of an asymptomatic malaria model without a recovered compartment. The solid lines represent progression from one compartment to the next, while the dotted stand for the human-mosquito interaction. Humans enter the susceptible compartment either through birth or migration and then progress through each additional compartment subject to the rates described above (Color figure online) such low vector activity lies in the intricate relationship between the parasite and host. This delicate relationship is captured by the term r 6 .

Variations of the SEYAR Model and the Impact of the Asymptomatic Class
This section is primarily concerned with introducing some modifications of the SEYAR model and their corresponding reproductive thresholds.

Elimination of the R Compartment in the SEYAR Model
As mentioned earlier, premunition refers to instances when infected hosts carry low parasitemias that confer immunity to the symptoms of the disease. When these individuals are infectious, they are able to transmit the parasite to the mosquitoes; thus, they belong to the A compartment. In addition to showing no symptoms, some individuals exhibiting premunition are unable to transmit the parasite; these hosts are considered part of the R compartment. We consider next the scenario when sterilizing immunity via premunition is negligible, thus eliminating the need for an R compartment.
Corollary 1 (Local asymptotic stability of the DFE for the SEYAR model without an R compartment) The reproductive threshold of the dynamical system corresponding to Fig. 2 is structurally the same as that of dynamical system corresponding to Fig. 1.
The only difference resides in the relabeling of the terms This figure is a schematic diagram of the SEY AR malaria model including relapse rates. The solid lines represent progression from one compartment to the next, while the dotted stand for the humanmosquito interaction. Humans enter the susceptible compartment through recovery, birth, or migration and then progress through each additional compartment subject to the rates described above. The new solid lines have been highlighted to emphasize the model modification (Color figure online) A verification of Corollary 1 is provided in the electronic supplementary material. For the proof, the reader is referred to "Proof of Corollary 1 (See Page 18)" of appendix.

Including Relapse Rates in the SEY AR Model
During the course of P. vivax and P. ovale human infections, a number of sporozoites stay in the hepatocytes in a dormant stage for a variable amount of time ranging from weeks to years. This stage of infection is known as the hypnozoite stage and is not observed in P. falciparum, P. malariae, or P. knowlesi. It is unclear what causes a hypnozoite to become active, but when it does, it causes a relapse of the disease. In order to properly account for such biological behavior, relapse rates are included into the S EYAR model. The asymptomatic and symptomatic human relapse rates are denoted by λ R A and λ RY , respectively. This slight adjustment yields the P. vivax and P. ovale version of the SEYAR model shown in Fig. 3.
The model parameters and initial data restrictions are the same as in the original model. Moreover, the results of the mathematical analysis concerning the issues of well-posedness and nonlinear stability are similar to those of the original model. For example, the results regarding the LAS of the DFE for the above model variation are exactly the same as in the original model.
Corollary 2 (Local asymptotic stability of the DFE for the SEYAR model including relapse rates) The dynamical system corresponding to Fig. 3 has the same reproductive threshold as the original dynamical system corresponding to Fig. 1. A verification of Corollary 2 is provided in the electronic supplementary material. For the proof, the reader is referred to "Proof of Corollary 2 (See Page 20)" of appendix.
Therefore, the modified system has the same reproductive threshold as the original. This is due to the fact that the sub-matrices of the Jacobian F and V , resulting from the block matrix partitioning technique utilized in the next-generation method, are unaffected by this minor modification. The results regarding the LAS of the DFE corresponding to the new system are inherited from the original.
The above modification attests to the delicate relationship between the compartmentalized modeling of infectious diseases and the next-generation approach. As we have seen in the above example, multiple compartmentalized models result in the same reproductive threshold.

The Impact of the Asymptomatic Class on the Reproductive Threshold
The goal of this section is to investigate how the reproductive threshold arising from the SEYAR model behaves in the case where asymptomatic carriers are not transmitting the disease. The natural parameter space of the SEYAR model corresponding to the reproductive threshold is The set Θ consists of all possible positive, epidemiologically reasonable, parameter values for the SEYAR model in which the reproductive threshold depends upon. Although the disease-induced death rate δ is allowed to be non-negative, the analysis presented in this section is unaffected by this minor detail. Under this formalization, neglecting the disease transmission and recovery rates of asymptomatic human hosts on the reproductive threshold formally corresponds to restricting the model to the following A-nullified parameter configuration spaceΘ, defined as follows: Consider a typical element Θ 0 ∈ Θ listed below: In a similar fashion, the dual elementΘ 0 ∈Θ corresponding to the symptomatic class Y is given as follows: The dual elementΘ 0 is comprised of the same fixed parameter configuration as described by Θ 0 , but with the asymptomatic progression rates specified above set equal to zero. To emphasize the asymptotic dynamic influence of the asymptomatic class A on the reproductive threshold R 0 (7) arising from the SEYAR model (4), subject to a given fixed parameter configuration, the following notation is employed , the size relationship between the two quantities R A and R Y is captured below in the following theorem.
Theorem 4 (Impact of the asymptomatic class on the reproductive threshold) Let R 0 be the threshold quantity (7) arising from Lemma 1 and consider the following fixed parameter configuration vectors Θ 0 ∈ Θ andΘ 0 ∈Θ, corresponding to the SEYAR model (4), defined as above. Denote This theorem means that R 0 increases when there are asymptomatic individuals capable of transmitting the disease. The proof of the above theorem is located in "Proof of Theorem 4 (See Page 21)" of appendix. The above inequality is strict. Therefore, neglecting to account for asymptomatic carriers results in an underestimation of the reproductive threshold. To demonstrate the theoretical estimate provided in Theorem 4, the numerical values of R A and R Y , along with the entomological inoculation rate (E I R) and parameter configuration space classifications, introduced via definition (1) in Sect. 3.2 are presented in Table 3 of Sect. 7.2. These numerical values correspond to the following three high transmission sites: Kaduna in Nigeria, Namawala in Tanzania, and Butelgut in Papua New Guinea. The parameter values associated with each site are listed in "Parameter Values" of appendix.
Previously, it was shown that R 0 can be written in the following form: where each term √ r i is defined as in Sect. 3.2. Care needs to be taken when arbitrarily substituting zero parameter values into the reproductive threshold R 0 (9). Clearly, if the human biting rate σ = 0 or mosquito-to-human transmission probability β M = 0, then it follows that R 0 = 0. These quantities effectively nullifying the reproductive threshold stand to epidemiological reason and correspond to the following scenarios, respectively: (1) no mosquitoes are biting humans, and (2) mosquitoes are not transmitting the disease. Furthermore, if both the asymptomatic β A and symptomatic β Y transmission probabilities are equal to zero, then the reproductive threshold will be identically zero, as infected humans are not transmitting the disease to susceptible mosquitoes.
However, it is implied that while one (or possibly all) of such parameter values may be equal to zero, the other parameter values in which the threshold depends on will be such that it is well defined. This problem is primarily due to the inclusion of vital dynamics for the human and mosquito populations into the dynamical system (4). For example, if one lets the human recruitment rate Ω H = 0 or the mosquito mortality rate ξ M = 0, then a singularity occurs and R 0 ceases to be well defined. Although the mosquito mortality rate being equal to zero is unreasonable in an epidemiological sense, as all mosquitoes experience death, it is informative to study how this parameters size effects R 0 , as in general, it is desirable to increase such a parameter when introducing control measures into the system. This attests to the subtle fact that after including vital dynamics into a given compartmentalized infectious disease model, one cannot expect to obtain qualitative information about the reproductive threshold in the absence of vital dynamics by simply setting the corresponding terms equal to zero. To properly study the asymptotic behavior of such models without vital dynamics included, one would have to go back to the original model derivation and not include them from the beginning, then proceed to calculate the threshold arising from the modified system.
i Consider the following fixed parameter configuration vector Θ 1 ∈ R 13 >0 defined as follows: Then, by Eq. (9) and the fact that the square root function is uniformly continuous on [0, +∞) (so that the limit can be taken inside), it follows that where the superscripts appearing above each appropriate ith term denote the fact that these quantities are fixed and thus invariant under the limit operation.
There exists a natural number ξ M n 0 ∈ N such that for all ξ M n ≥ ξ M n 0 , it follows that R 0 ∈ [0, 1). This qualitative observation can be interpreted as follows: Provided a scenario such that all of the associated model parameters are fixed epidemiologically reasonable quantities, if control measures sufficient to increase the vector death rate are introduced into the model, then the corresponding reproductive threshold arising from the model will be lowered and eventually fall below unity. This is consistent with the discussion in Sect. 3.2.
ii In a similar fashion, consider the following fixed vectors (Θ 2 , Θ 3 ) ∈ R 13 >0 × R 13 >0 , defined as follows: Then, it follows that The above analysis demonstrates how R 0 behaves with respect to the human and mosquito recruitment rates. As we will see in Sect. 7.1, if a given population, e.g., the three sites from which the parameter values are taken for this work, has a relatively small human recruitment rate and relatively large mosquito recruitment rate, then the factor r 1 will dramatically contribute to the size of R 0 . For example, in the case of the Kaduna site r 1 ≈ 14,762,941.18. The fractional multipliers r i for i = 2, 3, 4 will reduce the resulting size of this quantity, as they are all strictly less than unity. Further reduction in the size of the remaining quantity in the decomposition of R 0 depends on the terms σ , r 5 , and r 6 . These remaining factors will lower the resulting threshold provided that the mosquitoes are biting a relatively small amount of humans per unit time and transmitting less than they are dying off. Additionally, the asymptomatic and symptomatic human hosts need to be transmitting at a sufficiently low rate. For this reason, we must introduce control measures which both reduce the various disease transmission probabilities involved and increase the vector death rate.
Setting the parameters τ and γ equal to zero obviously results in R 0 = 0. However, these scenarios are not considered as the mosquito and human latent periods τ and γ are considered to be intrinsic properties of the vector and host, respectively. Moreover, these specific terms only appear in the factors √ r i for i = 2, 3 which are strictly bounded above by unity. Additionally, we do not consider the bizarre cases that the human mortality rate ξ H = 0 or mosquito recruitment rate Ω M = 0. Although both cases result in effectively nullifying the threshold quantity R 0 , humans are not immortal and the mosquito recruitment rate is usually relatively large. Letting this particular parameter value be equal to zero would imply the epidemiologically unreasonable scenario that there are no mosquitoes present in the region being considered. Neither do we consider the qualitative behavior of R 0 if ξ H is sufficiently large, as introducing a human transmission-blocking control measure which also increases the human death rate of a given population is not an ethical control method. It is important to note that emphasis is made on the terms which are related to utilized control measures, i.e., the measures which have an effect on the mosquito death rate and the various transmission probabilities involved. Control measures will be formally introduced and covered in Sect. 6.

Endemic Equilibria and Bifurcation Analysis
An endemic equilibrium occurs when disease persists in the population. For this reason, endemic equilibrium (EE) points are equilibrium points where some of the state variables corresponding to the infected classes are positive. Most epidemic models exhibit a dichotomy in terms of bifurcations that occur at the threshold R 0 = 1, namely super-critical (forward) and sub-critical (backward). These have drastically different epidemiological implications. A forward bifurcation happens when R 0 crosses unity from below and, as a result, a small positive asymptotically stable super-threshold equilibrium appears and the disease-free equilibrium losses its stability. Backward bifurcation happens when R 0 < 1 and a small positive unstable sub-threshold equilibrium appears, while the disease-free equilibrium and a larger positive equilibrium are locally asymptotically stable.
From an epidemiological viewpoint, a forward bifurcation is more desirable as it results in the reproductive number being below unity to be sufficient to ensure that an epidemic does not occur. In the presence of a backward bifurcation, the reproductive number being below unity is no longer sufficient, as sub-threshold endemic equilibria can arise in response to perturbations of specific parameters.
Due to the presence of the term e −ε = ∞ n=0 (−1) n n! E N H n , one cannot obtain a closed-form expression for the endemic equilibrium solutions of system (4). However, we turn to a variant of the center manifold theorem, introduced by Castillo- Chavez and Song (2004), to show the existence of non-trivial equilibrium solutions of the SEYAR model (4) near the DFE. This section is focused on the nonlinear stability analysis corresponding to the SEYAR model. More precisely, the following theorem concerning its bifurcation behavior is proven.
Theorem 5 (Bifurcation analysis for the SEYAR model) Let R 0 = 1 and the positive quantities η 1 and η 2 be defined as follows: where the terms labeled K i , Q i , and Z i , for 1 = 1, . . . , 7, are defined in "Appendix".
If the parameter Λ is defined as then the SEYAR model (4) exhibits a sub-critical bifurcation, provided that Λ > 1, and super-critical bifurcation, provided Λ < 1.
The formal proof of the above theorem can be found in "Proof of Theorem 5 (See Page 25)" of appendix. Moreover, a verification of the entries of the Jacobian and Hessian evaluated at the DFE is provided in the electronic supplementary material. As previously mentioned, Theorem 5 is proven by making use of an application of the center manifold theorem (Castillo-Chavez and Song 2004), adapted to the case of nonlinear dynamical systems. As in the case of most malaria models appearing throughout scientific literature, the type of bifurcation experienced by the system is completely determined by the sign of the a-term (18), listed in the appendix as Theorem 7. In the case of the SEYAR model (4) a ∝ (η 2 − η 1 ), resulting in a size constraint on Λ. In an epidemiological setting, it is desirable for the bifurcation, if it exists, to be super-critical, i.e., forward. Theorem 5 tells us that if one wants to avoid the case of a sub-critical bifurcation from occurring, we must demand the quantity Λ < 1 in addition to R 0 < 1.

Incorporating Control Measures into the SEYAR Model
The entomological inoculation rate (EIR) is a meaningful epidemiologic predictor that serves as a good measure of malaria intensity in a given region (Killeen et al. 2000). In 2007, Smith et al. estimated the reproductive number for 121 African populations. These estimates can be found in Figure 2 (Smith et al. 2007, p. 0534), where two numerical plots are displayed in which the reproductive number estimates are compared with the entomological inoculation rate of the populations under consideration. One plot corresponds to heterogeneous biting and transmission-blocking immunity taken into account in the parameter estimates and the other without. In both cases, the quantities were shown to be directly proportional, that is, regions with a relatively large (small) EIR also have a relatively large (small) reproductive number. In other words, regions with relatively large EIR values of each region also have relatively large R 0 values. In areas with large R 0 , it is unlikely that one single control measure will be sufficient to stop the disease expansion (Smith et al. 2007).
In practice, vaccine-conveyed immunity is not one hundred percent effective. This fact is accounted for by the vaccine efficacy V f , which denotes the percentage of protection each vaccinated individual has. If V p denotes the proportion of the population that has been vaccinated, i.e., the vaccine coverage, the product V f V p stands for the fraction of the population under consideration that is protected, so that the remaining proportion (1 − V f V p ) is not directly protected, with respect to vaccine-conveyed immunity. As a result, vaccination controls are incorporated into the model by defining the weightv := (1 − V f V p ). Therefore, the control-modified progression rates are given by the following equations: Insecticide-treated nets (ITNs) are the most prominent malaria preventive measure for large-scale deployment in highly endemic areas such as sub-Saharan Africa (Lengeler 2004). ITNs are nets coated with synthetic pyrethroid insecticides. Many studies have shown them to both kill and repel mosquitoes. In a recent study, a regression analysis of the protective efficacy on the transmission intensity, as measured by the EIR, was performed at the following four different endemic regions of Africa: Burkina Faso, The Gambia, Ghana, and Kenya. It was noted that the protective efficacy was lower in areas with a higher EIR, which was consistent with the original hypothesis that relative impact is lower in areas with higher entomological inoculation rates (Lengeler 2004). Moreover, in the case of homogeneous biting, 99.95% ITN coverage was predicted to be necessary (Smith et al. 2007).
Regarding the ITN coverage, let the symbol ρ f denote the protective efficacy, i.e., the percentage reduction in malaria episodes due to bed net usage. Upon letting ρ p be the proportion of ITN usage, i.e., the percentage decrease in transmission due to the employment of ITNs, then the reduction in mosquito-to-human transmission is captured by the multiplier (1 − ρ f ρ p ). Additionally, let ξ ITN denote the maximum ITN-induced death rate for the mosquito population. Following (Agusto et al. 2013), it is assumed that ITN usage reduces the effective human-to-mosquito effective contact probabilities β A and β Y and increases the mosquito mortality rate ξ M . Thus, the effects of ITN usage on the disease transmission dynamics of the SEYAR model (4) are accounted for by the following modification: Therefore, the resulting control-modified variant of the SEYAR model (4) expanded in the original model parameters is Assuming that the initial exposed proportion is zero, i.e., ε 0 = 0, the corresponding vaccination control-modified reproductive threshold Rv 0 is given by the following formula:

Sensitivity Analysis
Uncertainty is usually present in data collection and presumed parameter values. In this section, a sensitivity analysis is applied to classify the parameters which have the highest impact on the reproductive threshold R 0 . This provides a way to determine which parameter values should be targeted by intervention strategies. A parameter with a relatively large sensitivity index should be estimated with precision, while a parameter with a relatively small sensitivity index does not require as much effort. Let Θ be defined as in Sect. 4.3, then it is of trivial consequence that R 0 ∈ C 1 (Θ). Due to this fact and that we have an explicit expression for R 0 of the SEYAR model (4), we arrive at the following definition.
Definition 2 (Sensitivity index of the reproductive threshold (Chitnis et al. 2008)) Consider the reproductive threshold R 0 ∈ C 1 (Θ) given by Eq. (7) in Sect. 3.2, and let {e i : 1 ≤ i ≤ 14} be the canonical basis in R 14 . Forρ ∈ Θ define ρ i := e i ,ρ , where ·, · denotes the inner product in 14-dimensional Euclidean space, then the normalized forward sensitivity index of R 0 with respect to the parameter ρ i is defined by the following differential equation: The normalized forward sensitivity index provides a way to quantify the relative change in the given expression when the parameter changes. The sensitivity index is well defined, provided that R 0 is at least in C 1 with respect to each parameter ρ i .
The analytic formulas for the sensitivity indices are complex and do not offer much qualitative insight; as a result, we evaluate the indices at the parameter values corresponding to each location. Table 2 contains the sensitivity indices of R 0 for the SEYAR model (4) evaluated at the parameter values given in "Parameter Values" of appendix. The parameters are ordered from the most sensitive to the least.
A verification of the numerical entries in Table 2 is provided in the electronic supplementary material. The sensitivity analysis conducted on the reproductive threshold of the model above shows that the most sensitive parameters are the mosquito mortality rate ξ M and the human biting rate σ . Conversely, the least sensitive is the diseaseinduced death rate for humans δ. The sensitivity indices listed in the above tables can be viewed as growth measurements of the reproductive threshold with respect to the parameter under consideration. Furthermore, the result of the sensitivity analysis shows that mosquito parameters are most important, which is consistent with the classic model of Ross-Macdonald; however, the magnitude of the sensitivity indices varies due to the reconfiguration of the compartments in our model, and the addition of new parameters.
Without loss of generality, attention is turned to the Kaduna location. Concerning Kaduna, an increase in ξ M by 10% will result in a decrease in R 0 by 12.75%. Similarly, an increase in σ by 10% will cause a 10% increase in R 0 . Also, it is worth noting the asymptomatic recovery and effective contact probabilities λ AR and β A are relatively less sensitive than the symptomatic human-related terms λ Y R and β Y . Furthermore, an increase in Ω M by 10% results in an increase of R 0 by approximately 5%, and The hierarchy of sensitivity that the SEYAR model parameters obey is common among many compartmentalized homogeneous population malaria models appearing in the literature, e.g., Chitnis et al. (2008) and Roop-O et al. (2015). The two most sensitive parameters correspond to the vector population. These parameters have the property that one is directly proportional to the reproductive number, while the other is inversely proportional.
Increasing the mosquito death rate will also reduce the human biting rate, as the average mosquito life span is shortened. This is beneficial from a practical standpoint, as the parameter which is desirable to increase has the additional effect of decreasing the parameter that is desirable to decrease. Theoretically, this aids in the designing of programs for disease control, as it isolates the parameters that should be targeted for reduction by intervention strategies. Insecticide-treated bed nets and indoor residual spraying are among the most common methods used for such purposes. In Sect. 6, these control measures are incorporated into the model.

Numerical Results
Displayed below is a table containing numerical values for the following: (i) the parameter configuration space classification C 1 − C 2 , introduced via Definition 1 in Sect. 3.2, (ii) the EIR corresponding to each location, (iii) the reproductive threshold accounting for asymptomatic carriers R 0 , given by Eq. (7), (iv) the reproductive threshold neglecting asymptomatic carriers R Y , as discussed in Sect. 4.3.
These numerical values correspond to the following three high transmission sites: Kaduna in Nigeria, Namawala in Tanzania, and Butelgut in Papua New Guinea. The parameter values associated with each site are listed in "Parameter Values" of appendix. A verification of the numerical threshold quantities presented in Table 3 is provided in the electronic supplementary material. One should observe that the sizes of the reproductive thresholds for the three sites under consideration are consistent with the sizes of the corresponding entomological inoculation rate values. Regions with a relatively large EIR value also have relatively large R 0 values. As mentioned in Sect. 6, in areas with large reproductive thresholds, it is unlikely that one single control measure will be sufficient to stop the disease expansion.
Additionally, the threshold quantities R 0 = R A and R Y obey the theoretical estimate provided in Theorem 4. Therefore, neglecting to account for asymptomatic carriers results in an underestimation of the reproductive threshold corresponding to each location. The theoretical estimate provided in Theorem 4 holds for all possible positive epidemiologically reasonable quantities. Furthermore, the numerical entries displayed in Table 2 of Sect. 7.1 and those listed in Table 3 are point-wise evaluations.
There have been reports which indicate that asymptomatic carriers in a given population may transmit the disease at a higher rate than the symptomatic (Vallejo et al.  Service (1965) c This value was taken from Smith et al. (1993) d This value was calculated by Killeen et al. (2000), using data obtained from the following sources (Burkot et al. 1988;Graves et al. 1990) e These quantities were calculated from Eq. (7) assuming that the initial exposed proportion is zero, i.e., ε 0 = 0 f These quantities were calculated by setting the asymptomatic progression rates equal to zero 2016). In this case, the sensitivity hierarchy will possess a different ordering and the size differences in the threshold quantities displayed above will be larger.

Conclusions and Discussions
This study clearly shows that the existence of asymptomatic individuals results in a strict underestimation of R 0 and provides the means to quantify this influence. It also provides the means to study NAI as the factor that drives asymptomaticity. As mentioned by Doolan et al. (2009), the exploration of NAI is key to the rational development and deployment of vaccines and other malaria control methods corresponding to any given population at risk. Therefore, it is a necessary foundation upon to build strategies of eradication by any means. The SEYAR model (4) accounts for the impact that the exposure-dependent naturally acquired immune proportion has on asymptomatic carriers and malaria disease transmission dynamics. Through making use of the IVP (1), the infected compartment I is effectively decomposed into two mutually disjoint sub-compartments accounting for symptomatic and asymptomatic individuals. This results in a model which does not fall into a sub-class of the type studied by Filipe et al. (2007).
Current asymptomatic models appearing in the literature are formed by inserting a state-invariant constant control parameter or a sum of transcendental expressions involving various state-invariant immunity acquisition rates. The SEYAR model is derived by a separation through means of the NAI proportion of a population which depends on exposure through Eq. (2). After deriving the model and addressing the issues of well-posedness and stability analysis, a nonlinear stability analysis is performed in which the bifurcation behavior of the model is characterized. A sensitivity analysis is carried out, and generalized control measures are introduced in the model. Moreover, numerical values of various quantities discussed throughout this work are provided for the following three high transmission sites: Kaduna in Nigeria, Namawala in Tanzania, and Butelgut in Papua New Guinea. A brief summary of highlights drawn from the conclusions of this work is presented in the form of a list below: 1. In Sect. 3.1, it was shown that the SEYAR model (4) is mathematically and epidemiologically well posed, provided the initial data satisfied suitable regularity assumptions. Additionally, in Theorem 3, a mathematically precise and epidemiologically reasonable lower bound for the feasible region of the model was provided. 2. In Lemma 1 of Sect. 3.2, the SEYAR model (4) was shown to satisfy the threshold condition. Moreover, R 0 was decomposed into a product to properly analyze the size contribution of each factor involved. Provided that the mosquito mortality rate ξ M can be made large enough so that the term C 0 compromised of fractional multipliers is sufficiently small, i.e., C 0 < 1, then the size of R 0 is completely determined by the human-dependent factor r 6 , defined in Sect. 3.2, which consists of a weighted difference of vital dynamics with the NAI proportion, recovery, and death rates. Motivated by the monotonic behavior of R 0 , a formal characterization of the parameter configuration space was introduced via Definition 1. 3. In Sect. 4.1, a modification was made to the SEYAR model which accounts for the case when sterilizing immunity via premunition is negligible, thus eliminating the need for an R compartment. Moreover, in Sect. 4.2, relapse rates are introduced into the model and the corresponding reproductive threshold is calculated. In Sect. 4.3, an estimate is provided which characterizes the impact that the asymptomatic class has on the reproductive threshold. More precisely, it was shown that neglecting asymptomatic carriers results in an underestimation of the threshold. 4. In Sect. 5 use is made of Theorem 7 to show the existence of non-trivial subthreshold equilibrium solutions near the DFE. More precisely, it was shown that the bifurcation experienced by the SEYAR model (4) is forward or backward depending on the size of an auxiliary threshold parameter Λ defined by (10). As it is desirable for no endemic equilibrium states to arise while R 0 < 1, we impose the additional requirement that Λ < 1. 5. In Sect. 6, control measures are incorporated into the SEYAR model (4). 6. In Sect. 7.1, a sensitivity analysis is conducted on the reproductive number of the model for parameter configurations corresponding to high transmission settings. Additionally, a table was provided which contains the sensitivity indices of R 0 with respect to each parameter. An ordering of the parameters from the most sensitive to least revealed that the most sensitive parameters were the mosquito mortality rate ξ M and human biting rate σ . The least sensitive was the diseaseinduced human death rate δ. Additionally, the asymptomatic recovery rate λ AR and the transmission probability β A were shown to be relatively less sensitive than the symptomatic human-related terms λ Y R and β Y . 7. In Sect. 7.2, numerical results corresponding to the three high transmission sites mentioned above are provided. The numerical values of the threshold quantities were shown to be consistent with the theory presented in Sect. 4.3. As the sites of interest are high transmission areas with relatively high entomological inoculation rates, the reproductive thresholds were shown to be comparably large.
As directions of future research, it will be interesting to apply the method of decomposing the infected compartment by means of the related rate IVP (1) in Sect. 2 to other infectious diseases where asymptomatic individuals play a fundamental role in the disease dynamics. Additionally, it will be informative to consider extensions of the SEYAR model formed by incorporating the other kinds of immunity mentioned in the introduction. Furthermore, it will be beneficial to introduce additional control measures into the models such as aerial fogging and a time-dependent treatment rate of symptomatic carriers. An important direction for future exploration is the study of P. vivax and P. falciparum co-infection.
In conclusion, the SEYAR model (4) has provided us with a precise mathematical understanding of the relationship between the exposure-dependent nature of NAI and asymptomatic malaria disease transmission dynamics.

A Appendix
This section of the appendix contains formal proofs of the lemmas and theorems presented in Sect. 3. Listed below is an overview of the notation and mathematical framework utilized in the proofs. For an in-depth look into the function classes utilized in this work, and many other closely related topics in the field of harmonic analysis, the reader is referred to Fabec and Ólafsson (2014) and Grafakos (2008).

A.1 Notation and Mathematical Framework
Firstly, let R + := {x ∈ R : x ≥ 0} be the space of nonnegative real numbers and dt denote the Lebesgue measure for any complex-valued measurable function ϕ. The L 1 (R + ) norm of ϕ is defined as Extensive use is made of the space of essentially bounded functions L ∞ (R + ) characterized by the following norm: Analogous to the case of the essential supremum, the essential infimum is defined as ess inf Some of the proofs require at least a bounded first derivative; however, a general ϕ ∈ C 1 (R + ) is not bounded. Due to this fact, we need to restrict our analysis to a smaller space possessing a higher degree of regularity. The most natural space to turn to is the sub-space denoted C 1 b (R + ) ⊂ C 1 (R + ). This is the Banach space of bounded continuous functions whose first derivative is also bounded, endowed with the norm The convolution of two functions is defined as ϕ 1 * ϕ 2 (t) = t 0 ϕ 1 (τ )ϕ 2 (t − τ )dτ . Whenever ϕ 1 * ϕ 2 (t) ∈ L 1 (R + ), the integral is finite and, as a result, well defined.

A.2.1 Proof of Theorem 1 (See Page 8)
Proof The continuous composition of functions u(t) = (u • ε)(t) solves the IVP (1) for all ε ∈ C 2 (R + ) in the solution space. By the non-negativity of ε, if the initial data ε 0 are such that e ε 0 (u low − u high ) + u high ≥ 0, then the solution u ∈ C 2 (R + ) is non-negative for all t ∈ R + . The above inequality is satisfied, provided ε 0 ≤ ϑ for ϑ defined by Eq. (3) in Sect. 2.

A.2.2 Proof of Theorem 2 (See Page 11)
Proof The functions under consideration represent bounded, continuous, smoothly varying, nonnegative biological quantities. Consequently, it is reasonable to assume that x ∈ C 1 b (R + ). This reasoning combined with the fact that the vector field Φ ∈ C 1 is locally Lipschitz implies the existence of a unique solution x. Moreover, the solution x depends continuously on the initial data and model parameters and can be continued to a maximal time interval (Smale et al. 2003).

A.2.3 Proof of Theorem 3 (See Page 12)
Proof The case for the mosquito population is trivial, viz. For the human population, let L t = e t 0 ξ H ds = e ξ H t , so thaṫ Therefore, the explicit solution is given by After passing to the limit, it follows that In the absence of symptomatic infection, i.e., Y = 0, it follows that L −t * Y (t) = 0, so that As expected, N H asymptotes to the equilibrium population density of the human population. If symptomatic infection occurs in the population, by Hölder's inequality it follows that where the right-hand side of the above inequality is finite due to the assumption that . By letting t → +∞, we arrive at the following estimate for the forcing term The above estimate in combination with (12) yields the following lower bound Due to the restriction N H ∈ R + for all t ∈ R + , it follows that where α is defined in the statement of the theorem. Furthermore, for system (6) we have thatṄ Combining the above observation with a similar argument for N M yields Inequalities (13) imply that if the solution leaves the region Γ , then its derivative will instantaneously become negative, forcing it back to Γ . Moreover, if x i (0) = 0 for any 1 ≤ i ≤ 8 in system (5), then it directly follows thatẋ i ≥ 0. Therefore, all trajectories tend to Γ and are forward invariant. Due to this fact, it is sufficient to study the dynamics of the system on the smaller compact sub-space Γ .

A.2.4 Proof of Lemma 1 (See Page 14)
Proof Firstly, we order the compartments so that the first five correspond to infected individuals and denote w = (E, Y , A, M E , M I , R, S, M S ) T . The corresponding DFE is Following the next-generation method, system (4) is rewritten in the following form: where F := (F 1 , . . . , F 8 ) T and V := (V 1 , . . . , V 8 ) T , or more explicitly In addition, the matrix V admits the decomposition V = V − − V + , where the component-wise definition is inherited. Biologically speaking: F i is the rate of appearance of new infections in compartment i, V + i stands for the rate of transfer of individuals into compartment i by any other means and V − i is the rate of transfer of individuals out of compartment i. It is easy to see that F, V − , V + satisfy assumptions (i)-(v) in Theorem 6. As mentioned in the beginning of Sect. 3.1, to study the stability of the equilibrium points, we assume that each of the above vector fields is at least twice continuously differentiable. Define U low := e ε 0 (u low − u high ) + u high , and let F and V be the following sub-matrices of the Jacobian of the above system, evaluated at the solution w d f e and F V −1 is given by the following matrix: Let I denote the 5 × 5 identity matrix, so that the characteristic polynomial P(λ) of the matrix F V −1 is given by The solution set {λ i } 1≤i≤5 is given by Therefore, the reproductive threshold for the SEYAR model (4) is given by The proof of the lemma regarding the local asymptotic stability (LAS) of the DFE w d f e corresponding to the SEYAR model (4) is now complete after invoking Theorem (6) in "Summary of Stability Theorems" of appendix.

A.2.5 Proof of Corollary 1 (See Page 18)
Proof Let the DFE v d f e be defined as v d f e = 0, 0, 0, 0, 0, A straightforward calculation shows that the sub-matrices F and V of the Jacobian evaluated at the DFE v d f e corresponding to Fig. 2 are given by Therefore, the proof directly follows as in Corollary 2.

A.2.6 Proof of Corollary 2 (See Page 20)
Proof A straightforward calculation shows that the sub-matrices F and V of the Jacobian evaluated at the DFE w d f e corresponding to Fig. 3 are given by The reproductive threshold R 0 := ρ F V −1 for any given compartmentalized infectious disease model is completely determined by the matrices F and V . Therefore, it is of trivial consequence that the models corresponding to Figs. 1 and 3 possess identical reproductive thresholds, given by (7) arising from Lemma 1.

A.2.7 Proof of Theorem 4 (See Page 21)
Proof Let Θ 0 ,Θ 0 , R A , and R Y be defined as in the statement of the theorem. Additionally, denote the following quantities r 0 6 := r 6 Θ 0 , C 0 1 := C 1 Θ 0 , C 0 2 := C 2 Θ 0 and where the terms r 6 and C i for i = 0, 1, 2 are defined as in Sect. 3.2. By the monotonicity of the square foot function, it follows that

A.2.8 Proof of Theorem 5 (See Page 25)
Proof If we consider the dynamical systemẋ = g(x, ω), where ω is a bifurcation parameter and the vector field g is C 2 in both x and ω, then the disease-free equilibrium can be viewed as the manifold (x d f e ; ω) where the local stability of x d f e changes at the point (x d f e ; ω ). Now we shall investigate the existence and stability of non-trivial equilibrium states in a neighborhood of the bifurcation point. We focus on the diseasefree equilibrium x d f e and study the occurrence of a transcritical bifurcation at R 0 = 1.
Since R 0 consists of the square root of a complicated combination of parameters, it is not practical to use as a bifurcation parameter. However, observe that R 0 = 1 if and only if In lieu of Lemma 1, it follows that x d f e is locally asymptotically stable when β M < β M and unstable if β M > β M . Thus, the combination of parameters β M is a bifurcation value. To simplify the notation, we rewrite system (4) asẋ = g(x, β M ) where x = (S, E, Y , A, R, M S , M E , M I ) T so that x i is the ith component of x and g = (g 1 , g 2 , g 3 , g 4 , g 5 , g 6 , g 7 , g 8 ) T , or more explicitly Denote J (x d f e , β M ) to be the Jacobian of g evaluated at the DFE x d f e and threshold β M . Thus, J (x d f e , β M ) is given by If we invoke the following positive change of variables then β H can be written as .
As a result, J (x d f e , β M ) can be written as Due to the equivalency of the two conditions R 0 = 1 and β M = β M , it follows that J (x d f e , β M ) contains information of the linearized system evaluated at the disease-free equilibrium and threshold value. Utilizing the machinery covered in "Summary of Stability Theorems" of appendix below, the Jacobian evaluated at the threshold value, i.e., J (x d f e , β M ), has a zero simple eigenvalue with all others having negative real parts. Therefore, the hypothesis of Theorem 7 is satisfied. We proceed by calculating the a and b terms (18) and (19) appearing in Theorem 7. In observance of the conclusions in the theorem, it follows that the SEYAR model (4) will undergo a super-critical bifurcation if a > 0 and b > 0 and a sub-critical bifurcation if a < 0 and b > 0. The main ingredients in calculating a and b are the generalized right and left eigenvectors of the matrix J (x d f e , β M ) and their corresponding nonzero Hessian entries, evaluated at the DFE x d f e . In this fashion, we let w = (w 1 , w 2 , w 3 , w 4 , w 5 , w 6 , w 7 , w 8 ) and For system (14), the nonzero partial derivatives of g evaluated at x d f e are Since D i ∈ R >0 for i = 1, 2, it follows that b is always positive; hence, the bifurcation behavior of system (14) is completely determined by the sign of a. To simplify the pending calculation, we invoke the following change of variables: so that the generalized eigenvectors can be written as: In view of the relatively strong regularity assumptions mentioned in the beginning of Sect. 3.1 and natural symmetry of the system, it follows that and v k = 0 for k = 1, 5, 6. As a result, the terms that contribute to the sum correspond to k = 2, 3, 4, 7. Thus, a can be written as

A.3 Summary of Stability Theorems
The goal of this section of the appendix is to provide the reader with a brief collection and overview of the theorems that are widely used in determining the stability of the equilibrium points for nonlinear dynamical systems. The first theorem will be needed in order to determine the local asymptotic stability of the DFE corresponding to the SEYAR model, while the second is used to investigate the existence of non-trivial subthreshold equilibrium states of the model. In the setting of dynamical systems, one cannot usually pinpoint a solution exactly, but only approximately. As a result, an equilibrium point must be stable to be physically meaningful. A stable equilibrium point Provided an epidemiological model can be grouped into n homogeneous compartments, the local asymptotic stability of the equilibrium states can be established by utilizing the next-generation method, appearing in Van den Driessche and Watmough (2002). By making use of the center manifold theorem (John and Philip 1997), Van den Driessche and Watmough provided a simple prescription for determining the local asymptotic stability of DFE points of a given system. This criterion is given in terms of the reproductive number R 0 of the system which acts as a threshold value. This effectively relates R 0 to the DFE of the system. As a result, this has proven to be very useful in disease control. To cast the above discussion into a mathematical framework, we let x = (x 1 , . . . , x k ) T ∈ R k + and define the space of disease-free states for the compartmental model to be X := x ∈ R k + : x i = 0 for i = 1, · · · , m, m < n . Then, for Φ ∈ C 2 (R k ), we form the following dynamical system: Theorem 6 (Van den Driessche and Watmough 2002) Define R 0 = ρ F V −1 and consider the disease transmission model given by (17) such that Φ satisfies the follow criteria: i If x ≥ 0, then so are F, V + , and V − ,  (17), then x is locally asymptotically stable provided R 0 < 1 and unstable if R 0 > 1.
The above theorem is proved by making use of the following lemma.
Lemma 2 If x is a DFE of system (17) and Φ satisfies assumptions (i)-(v), then the derivatives DF(x ) and DV(x ) can be partitioned as follows: where F and V are defined as: In the above matrix partitioning, F is non-negative, V is an invertible M-matrix, and J 3 , J 4 are sub-matrices of the Jacobian associated with various transmission terms. This theorem provides a convenient epidemiological interpretation of the reproductive threshold R 0 corresponding to a given dynamical system in the SIR family. Additionally, due to the partitioning of the Jacobian mentioned above, the stability of the system is determined by det F V −1 − λI , where I is the identity matrix. If the matrix F containing transmission probabilities and contact rates is set to zero, then all eigenvalues of −V have negative real part. As a result, the stability, or lack of, experienced by the system in question depends on the entries of F. Due to this reason, the bifurcation parameters are chosen from the entries of F. In the case that the transmission probabilities are relatively high, an endemic could occur and additional non-trivial equilibrium can arise. Regarding the SEYAR model, besides the human biting rate σ , the only possible parameters are the transmission probabilities β Y , β A and β M . From a epidemiological perspective, we have more control over the mosquito-to-human transmission probability β M . Due to this reason, β M is chosen as the parameter for the bifurcation analysis. One could perform a similar analysis for β Y and β A ; however, all of the transmission probabilities involved are related through the reproductive number R 0 of the model. As a result, the analysis would be similar.
Let s(A) and ρ(A) stand for the spectral abscissa and radius, respectively. The proof of Theorem 6 hinges on M-matrix theory. A matrix B is said to have the Z -sign pattern, provided all of its off diagonal entries are non-positive. If B = sI − P, where P ≥ 0. If s > ρ(P), then B is a non-singular M-matrix; if s = ρ(P), then B is a singular M-matrix. This observation is then combined with a linear algebra lemma, which we restate below for convenience of the reader. To this end, they make use of the following argument: define J 1 := F − V , then V is a non-singular M-matrix and F is non-negative. It follows that −J 1 := V − F has the Z -sign pattern. By the nonnegativity of F V −1 , it is a direct consequence that −J 1 V −1 := I − F V −1 also has the Z -sign pattern. Now we apply Lemma 3 with H = V and −J 1 := V − F, to conclude that −J 1 is a non-singular M-matrix if and only if I − F V −1 . Also, since all of its eigenvalues have magnitude bounded above by ρ(F V −1 ), we have that I − F V −1 is a non-singular M-matrix if and only if ρ(F V −1 ) < 1. Therefore, s(J 1 ) < 0 if and only if R 0 < 1. Through making use of a similar argument in combination with Lemma 6 of Appendix A, found in Van den Driessche and Watmough (2002), they arrive at the following trichotomy, establishing R 0 as a threshold parameter In Sect. 5 use is made of the following variant of the center manifold theorem specifically adapted to the case of bifurcation analysis for nonlinear systems. The utility of this theorem resides in the classification of the bifurcation due to the sign of the parameters a and b, defined below.
Theorem 7 (Castillo-Chavez and Song 2004) Consider the following dynamical system with real parameter ω: Without loss of generality, we assume that x = 0 is an equilibrium point for the above system for all ω, i.e., f (0, ω) ≡ 0 for all ω. Provided that the following assumptions are satisfied: Then, the local dynamics around the equilibrium point 0 are completely determined by the signs of a and b. More precisely, i a > 0, b > 0. When ω < 0 such that |ω| 1, 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0 < ω 1, 0 is unstable and there exists a negative locally asymptotically stable equilibrium. ii a < 0, b < 0. When ω < 0 such that |ω| 1, 0 is unstable; when 0 < ω 1, 0 is locally asymptotically stable and there exists a positive unstable equilibrium.
iii a > 0, b < 0. When ω < 0 such that |ω| 1, 0 is unstable and there exists a locally asymptotically stable negative equilibrium; when 0 < ω 1, 0 is stable and a positive unstable equilibrium exists. iv a < 0, b > 0. When ω changes sign, 0 changes its stability from stable to unstable.
As a result, a negative unstable equilibrium becomes positive and locally asymptotically stable. In particular, if a > 0 and b > 0, then a backward bifurcation occurs at ω = 0.
The a and b terms appearing in the above theorem depend on generalized eigenvectors, i.e., zero entries are allowed. In the proof of Theorem 7, the a and b terms arise from a differential equation obtained from a parameterization of a one-dimensional center manifold c(t) given byċ = a 2 c 2 + bωc.
Observe how a transcritical bifurcation occurs in the above equation at ω = 0 and can be classified according to the signs of the a and b terms, defined above. These terms depend on the Kronecker product of generalized eigenvectors and entries of the Hessian evaluated at the DFE. As pointed out by Castillo-Chavez and Song (2004), negative components of the generalized eigenvectors are permitted through a modification of Theorem 7. The essence of the argument hinges on the fact that the theorem can still be applied; one only has to compare the negative entries of the eigenvectors with their corresponding entries of the nonnegative equilibrium of interest. Therefore, one has to consider the original parameterization of the center manifold, prior to the coordinate change where the DFE is assumed to be zero. Notice how the negative entries in the generalized eigenvectors calculated in Sect. 5 correspond to the positive entries of the DFE of interest.

A.4 Parameter Values
Presented below are tables of numerical rates corresponding to the following three high transmission sites: Kaduna in Nigeria, Namawala in Tanzania, and Butelgut in Papua New Guinea. The data have been collected from multiple sources, all of which are noted in the footnotes. Due to the additional prevalence of P. vivax in Butelgut, the corresponding data represent combined estimates of both species (Killeen et al. 2000). The data for the other sites exclusively correspond to the P. falciparum species. There are a variety of vector species that inhabit these sites. The data listed here correspond to the dominant vector species of the area being considered. The dominant vector species of each site is listed in the footnotes. Furthermore, as mentioned in Sect. 2 and reported by Laishram et al. (2012a), the parasites carried by asymptomatic human hosts can be more infectious than those of symptomatic. For the parameter values displayed in "Parameter Values" of appendix, we invoke the assumption that the asymptomatic carriers corresponding to each site transmit at a lower rate than that of the symptomatic and that asymptomatic individuals recover faster than symptomatic, i.e., β A < β Y and λ AR > λ Y R . Table 4 contains demographic data for the countries that each site is contained in. These numerical values are used to calculate the recruitment rates of the human populations Ω H corresponding to each region under consideration (Tables 5, 6, 7).  (2015) The human population data displayed in Table 4 are from the Central Intelligence Agency (CIA). Ω LE , Ω BR , and Ω MR denote the life expectancy, birth rate, and migration rate of the population under consideration. Life expectancy is measured in years. The birth rate entries appearing in column 3 of the above table are crude birth rates. The crude birth rate measures the average quantity of live births during a year, per 1000 people, and is given in units of total births per 1000 people per year. The migration rates are net migration rates, which measure the difference of immigrants and emigrants in a given population over the span of a year, and are given in units of humans per year, per 1000 people   (Service 1965). Let Ω L stand for the average An. gambiae life expectancy, which is dependent upon the region under consideration. The data provided by Molineaux et al. (1979) and appearing in row nine column one of Table A.3 in Chitnis et al. (2008, p. 19) correspond to An. gambiae activity in Nigeria. Due to this, we select the entry contained in row one column nine, so that Ω L = 9. Therefore, the average daily An. gambiae mosquito mortality rate is ξ M = 1/Ω L = 1/9 ≈ 0.11 e The assumption is made that the rate of transmission from asymptomatic humans to susceptible mosquitoes is one-tenth the transmission probability corresponding to symptomatic humans, i.e., β A = 0.048 f We adopt the convention utilized by Chitnis et al. (2008), where an estimate of 0.48 will be used for high transmission areas and an estimate of 0.24 will be used for low. Kaduna is a high transmission area; thus, we let β Y = 0.48. An alternate choice would be to choose the average value of the parameters corresponding to the dominant species of parasite. The average value of the parameters appearing in rows one through five in column one of Table A .6 in Chitnis et al. (2008, p. 21) is approximately 0.36. The origin of the parameters used in this calculation is Thomson et al. (1957), Boyd (1949) and Draper (1953) g The effective daily vector-to-human contact rate is obtained by dividing the value in row two column four of Table 3 in Killeen et al. (2000, p. 541) by Ω L . Thus, β M = 0.29/Ω L = 0.29/9 ≈ 0.032 h The human latent periodγ , measure in days, corresponding to P. falciparum infection is taken from row three column one in Table A .7 in Chitnis et al. (2008, p. 21). The range 9-10 is given. The average value is then chosen, i.e.,γ ≈ 9.5. It follows that the average duration of the intermediate host latent period is γ = 1/γ = 1/9.5 ≈ 0.11. The parameter source is Molineaux and Gramiccia (1980) i Letτ denote the Plasmodium incubation period, i.e., the number of days required for parasite development. Thus, by using the entry in row nine column three of Table 2 in Killeen et al. (2000, p. 539), it follows that the average duration of the definitive host latent period is τ = 1/τ = 1/11.6 ≈ 0.09. Technically, this parameter was calculated from the mean and median temperatures listed in the original source (Craig et al. 1999) j The malaria death rates are taken from World Life Expectancy (2016) and are given in units of per 100, 000 people per year. As in the case of the human demographic data listed in Table 4, these rates correspond to the overall country that the region is contained in. Using the data provided by World Life Expectancy (2016), it follows that δ = 60.46/365.25/100, 000 ≈ 1.66 × 10 −6 k The average human biting rate is estimated by dividing the entry in row one column four of Table 3 in Killeen et al. (2000, p. 541) by Ω L . More precisely, we have that σ = 3.8/9 ≈ 0.42 l,m For the purpose of numerical simulations, we assume that asymptomatic carriers transmit malaria to a lesser extent than symptomatic carriers and recover faster. Due to this, we assume that the recovery rate λ Y R for symptomatic individuals is one-tenth the rate λ AR of asymptomatic individuals. Let the quantity 1/λ denote the average duration of the human infectious period, provided that the individual has had no treatment. If Ω I stands for the average duration of the P. falciparum infectious period in humans, we select the average of the entries appearing in row four column one of Table A .9 in Chitnis et al. (2008, p. 21), so that Ω I = 18. Therefore, λ Y R ≈ λ = 1/Ω I = 1/18 ≈ 0.06, and it follows that λ AR = 0.6. The original source of the parameter value is Bloland and Williams (2002) n The temporary immunity loss rate λ RS is such that 1/λ RS is equal to the average duration of the human immune period. As in Chitnis et al. (2008), we assume that the immune period lasts for an average of 5 years in areas of high transmission and 1 year in areas of low transmission. For Kaduna, it follows that λ RS = 1/(5)365.25 ≈ 5.48 × 10 −4 o This baseline assumption is due to the fact that the population under consideration is A-dominant  Table 4, it follows that Ω H = (36.39 − 0.54)/365.25/1000 ≈ 9.82 × 10 −5 b By utilizing the data provided in Table 3 in Killeen et al. (2000, p. 541), it follows that Ω M = (1.8 × 10 6 )/365.25 ≈ 4928.13 c The average natural human mortality rate is calculated in a similar fashion, as in the case of Kaduna, by utilizing the data provided in Table 4 d The dominant vector species of Namawala at the time the field measurements were taken was An. gambiae (Smith et al. 1993). The data provided by Gillies and Wilkes (1965) and appearing in row six column one of Table A.3 in Chitnis et al. (2008, p. 19) correspond to An. gambiae activity in Tanzania. Using these data, it follows that Ω L = 11.26. Therefore, the average daily An. gambiae mosquito mortality rate is ξ M = 1/11.26 ≈ 0.09 e As previously mentioned, the asymptomatic human-to-susceptible mosquito transmission probability is assumed to be β A = 0.048 f Since Namawala is a high transmission area, we let β Y = 0.48 g By dividing the value in row two column five of Table 3 in Killeen et al. (2000, p. 541) by Ω L , it follows that β M = 0.017/11.26 ≈ 0.002 h The average duration of the intermediate host latent period is taken to be γ = 1/γ = 1/9.5 ≈ 0.11. The parameter valueγ is taken from row three column one in Table A .7 in Chitnis et al. (2008, p. 21). The average value is then chosen, i.e., γ ≈ 9.5. The original parameter source is Molineaux and Gramiccia (1980) i By using the entry in row eight column four of Table 2 in Killeen et al. (2000, p. 539), it follows that the average duration of the definitive host latent period is τ = 1/11.6 ≈ 0.09 j By the data provided by World Life Expectancy (2016), it follows that δ = 42.42/365.25/100, 000 ≈ 1.16 × 10 −6 k By using the data in row one column five of Table 3 in Killeen et al. (2000, p. 541), we have that σ = 1.5/11.26 ≈ 0.13 l As in the previous table, the recovery rate pertaining to asymptomatic individuals is assumed to be λ AR = 0.6 m The average of the entries appearing in row four column one of Table A.9 in Chitnis et al. (2008, p. 21) is selected, so that Ω I = 18. Therefore, λ Y R ≈ λ = 1/Ω I = 1/18 ≈ 0.06. The original source of the parameter value is Bloland and Williams (2002) n As Namawala is a high transmission area, the temporary immunity loss rate is taken to be λ RS = 1/(5)365.25 ≈ 5.48 × 10 −4 o This baseline assumption is due to the fact that the population under consideration is A-dominant  Table 4, it follows that Ω H = (24.38 − 0.00)/365.25/1000 ≈ 6.67 × 10 −5 b By utilizing the data provided in Table 3 in Killeen et al. (2000, p. 541), it follows that Ω M = (1.5 × 10 6 )/365.25 ≈ 4106.78 c Similarly, the average natural human mortality rate is calculated by utilizing the data provided in Table 4 d The dominant vector species of Butelgut at the time the field measurements were taken was An. punctulatus (Burkot et al. 1988). The data provided by Peters and Standfast (1960) and appearing in row thirteen column one of Table A.3 in Chitnis et al. (2008, p. 19) correspond to An. punctulatus activity in Papua New Guinea. Using these data, it follows that Ω L = 7.1. Therefore, the average daily An. punctulatus mosquito mortality rate is ξ M = 1/7.1 ≈ 0.14 e As previously demonstrated, the asymptomatic human-to-susceptible mosquito transmission probability is assumed to be β A = 0.048 f Since Butelgut is a high transmission area, we let β Y = 0.48 g By dividing the value in row two column six of Table 3 in Killeen et al. (2000, p. 541) by Ω L , it follows that β M = 0.042/7.1 ≈ 0.006 h The average duration of the intermediate host latent period is taken to be γ = 1/γ = 1/9.5 ≈ 0.11. The parameter valueγ is taken from row three column one in Table A .7 in Chitnis et al. (2008, p. 21). The average value is then chosen, i.e., γ ≈ 9.5. The original parameter source is Molineaux and Gramiccia (1980) i By using the entry in row eight column five of Table 2 in Killeen et al. (2000, p. 539), it follows that the average duration of the definitive host latent period is τ = 1/9.6 ≈ 0.1 j By the data provided by World Life Expectancy (2016), it follows that δ = 37.57/365.25/100,000 ≈ 1.03 × 10 −6 k By using the data in row one column six of Table 3 in Killeen et al. (2000, p. 541), we have that σ = 0.99/Ω L = 0.99/7.1 ≈ 0.14 l Due to similar reasoning as above, the recovery rate pertaining to asymptomatic individuals is assumed to be λ AR = 0.6 m The average of the entries appearing in row four column one of Table A.9 in Chitnis et al. (2008, p. 21) is selected, so that Ω I = 18. Therefore, λ Y R ≈ λ = 1/Ω I = 1/18 ≈ 0.06. The original source of the parameter value is Bloland and Williams (2002) n As Butelgut is a high transmission area, the temporary immunity loss rate is taken to be λ RS = 1/(5)365.25 ≈ 5.48 × 10 −4 o This baseline assumption is due to the fact that the population under consideration is A-dominant