The Seasonal Reproduction Number of p.vivax Malaria Dynamics in Korea

Understanding the dynamics of Malaria can help in reducing the impact of the disease. Previous research proved that including animals in the human transmission model, or ’zooprophylaxis’, is effective in reducing transmission of malaria in the human population. This model studies plasmodium vivax malaria and has variables for animal population and mosquito attraction to animals. The existing time-independent Malaria population ODE model is extended to time-dependent model with the differences explored. We introduce the seasonal mosquito population, a Gaussian proﬁle based on data, as a variant for the previous models. The seasonal reproduction number is found using the next generation matrix, endemic and stability analysis is carried out using dynamical systems theory. The model includes short and long term human incubation periods and sensitivity analysis on parameters and all simulations are over three year period. Simulations show for each year larger peaks in the infected populations and seasonal reproduction number during the summer months and we analyze which parameters have more sensitivity in the model and in the seasonal reproduction number. Analysis provides conditions for disease free equilibrium (DFE) and the system is found to be locally asymptotically stable around the DFE when R 0 < 1, furthermore we ﬁnd the uniqueness of the endemic equilibrium point. The sensitivity analysis for the parameters shows that the model was not sensitive to the exact values of the long or short term periods as it was to the average number of contacts between host and mosquito or rate of disease progression for mosquitoes. This model shows that inclusion of variable mosquito population informs how domestic animals in the human population can be more effectively used as a method of reducing the transmission of malaria. The most relevant contribution of this work is including the time evolution of mosquito population and simulations show how this feature af-fects human infection dynamics. An analytical expression for the endemic equilibrium point is provided for future work to establish conditions over which an epidemic may be prevented.


Introduction
In the study of the characteristics of the disease of Malaria effort has been put into modeling it's transmission. This research seeks to expand on previous efforts and motivation of this work based on two such efforts, one which demonstrated how domestic animals attract mosquitoes over humans and the other which considers the incubation period as a variable categorizing short versus long periods. The dilution effects of the domestic animal population on the transmission of plasmodium vivax malaria has been explored by several authors [1], [5], [8], [10], [12], [13], [14], [16], [20], [21].
The fact that domestic animals seem to attract this variety of mosquito with subsequent decrease in mosquito bite rate on humans was expanded to show how much the number of domestic animals influences the rate of malaria transmission. The dynamics of malaria transmission without morbidity was modeled with a SEIS epidemiological system. This study assumed a constant mosquito population. However it did use data to fit time 200 The Seasonal Reproduction Number of p.vivax Malaria Dynamics in Korea dependent mosquito population with a resulting distribution of a Gaussian-like curve where the mosquito population reaches the peak during July and August, demonstrating seasonality as seen in [2], [11], [18].
Including variable incubation periods shows to be necessary and a justified complexity, see [9] for more details, the long and short incubation periods of p.vivax Malaria in Korea are considered where the incubation period was obtained based on the difference between date of exposure to onset. Using the estimated mean incubation period, cases were classified as short incubation period or long incubation period where symptom onset before 25 weeks post-infection was defined as a short incubation period and symptom onset after 33 weeks post-infection was defined as a long incubation period as in [15], [17].
In particular, the incubation period included in [17] for p.vivax malaria transmission model, called p(t) is represented as a weighted sum of an exponential function and a step function, which means short-term incubation periods is modelled by an exponential distribution with long term incubation period assumed to have fixed lengths as per step function. As following : P (u) = pP s (u) + (1 − p)P l (u) where P s (u) = e −uτ , P l (u) = u − u i .
A secondary transmission is addressed, a discrete version of renewal process for p.vivax malaria modeling is applied and it is then fitted to an extrapolated parametric function of the yearly component of the reproduction numbers to statistically detect the trend of secondary transmission. Analysis led to the conclusion that the best fitted model appeared to be the model with abrupt decline in secondary transmission by 34% from 1998 and this coincides with the expansion of chemoprophlaxi.
Monthly incidence of malaria has the peak in July and August and the highest incidence in the area near the DMZ(demilitarized zone) in the northern part of Korea. The incidence was the highest in 20-25 years old group and the mean of age was 31.0 in males but 45.3 in females. Reports adaptation of HGLM (hierarchical generalized linear model) to estimate the annual pattern of malaria incidence by age and sex in Korea. And this model considered spatial and temporal correlations for the best estimation, [18]. There is no preceding research for basic reproduction number of the ODE model with seasonal characteristics to include spread in Summer and Autumn. The seasonal reproduction number is obtained and the basic reproduction number is calculated including an exponential functional to model the incubation period of malaria. In changing to existing time-independent ODE model to time-dependent model, differences are explored. Also, using bi-modal malaria incubation, changes to the infectious population in contrast to constant incubation period of the ODE model are discussed. Stability results for the disease free equilibrium were obtained. An explicit expression for the endemic equilibrium point is reached and stability results are obtain under certain constraints on the parameters by using similar methods such as those seen in [3], [4], [6], [7], [19], [22]. Furthermore, simulations are performed, confirming the seasonality of the disease and impact of bi-modal incubation on malaria transmission.

Model formulation
In Korea p.vivax Malaria is endemic disease. p.vivax Malaria ended at 1970's in Korea. But it re-occurred in 1993 at DMZ and then, the military took the precaution of Chemo-prophylaxis. Malaria transmits only between human and mosquito, that is, it does not transfer from human to human or mosquito to mosquito. When an infectious mosquito (I M ) bites a susceptible human (S H ), then susceptible human is able to become exposed human having either short (E s H ) or long (E l H ) term incubation. After the incubation, exposed human can become an infectious human(I H ). Similarly, when a susceptible mosquito (S M ) bites an infectious human, then the susceptible mosquito becomes an exposed mosquito(E M ) and then an infected mosquito. After an infected human or mosquito recovers, the model assumes return to the susceptible state. As mentioned before, we classified human (H) to four classes (susceptible, exposed long, exposed short, infected) and mosquito (M ) to three classes (susceptible, exposed, infected). Research has shown that mosquitoes have higher probability of being attracted to animals for a blood meal than to humans [10], therefore, the mosquito contact rate can be assumed as σ(1 − a)H/((1 − a)H + aA)(0 ≤ 1). Thus, the proportion of bites on human host (σ) is modified to be (1−a)H 2 ((1−a)H 2 +aĀ 2 ) . Since p.vivax malaria shows seasonal infectious dynamics, we considered dividing the mosquito population (M ) into two categories. In previous research it is assumed that the malaria population is constant, taken as the average number,M . Then, since each population remains constant, this is time-independent malaria dynamics. This is extended by assuming the malaria population has a Gaussian function. The time-dependent malaria dynamics are characterized by a sharp increase in Summer and Autumn followed by a decrease each year. The seasonal reproduction number is presented to find the number of secondary transmissions of time-dependent malaria. Accordingly, comparison of the differences between time-independent and time-dependent malaria in the number of secondary infectious human and mosquitoes and investigations with parameter changes are presented. The four parameters b m , r, σ, T m in table 1, are varied and as these parameters change and with those changes the impact on the model and the number of infectious humans and mosquitoes informs sensitivity of the model to certain parameters. The basic reproduction number (R 0 ) and seasonal reproduction number (R s ) with the difference of R 0 and R s for three years were calculated. In addition, consideration of the bi-modal incubation periods of malaria provides a more robust model as the incubation of exposed human has short and long period and rather than using the average value and several values are considered. Using this, comparisons of how the incidence of malaria changes. See Figure 1. for a flow diagram representing the interaction and transmission between the human and mosquitoes populations. , susceptible human (S H ), exposed human with short (E s H ) or long (E l H ) term incubation, exposed becomes an infectious human(I H ), susceptible mosquito (S M ) becomes an exposed mosquito(E M ) and then an infected mosquito; classified human (H) to four classes (susceptible, exposed long, exposed short, infected) and mosquito (M ) to three classes (susceptible, exposed, infected).

Model System of Ordinary Differential Equations
A system of ordinary differential equations is used to model the dynamics of malaria:

Model Invariance
In this model, rather than specific populations, the proportion of human, mosquito and animal population is used. That is, the following variables are dimensionless : (1), The rate of change for the total human population is given by similarly,the rate of changes for the total mosquito population is given by It is observed that all solution of System (1) starting in remain in D for all t ≥ 0, then D is a positively invariant and it is enough to consider solutions in D.

The threshold dynamics
In this section, the analysis of System (1) is developed, regular theory of dynamical systems is used in order to analyze the system locally around the equilibrium points. The disease free equilibrium E 0 and its stability is included. The endemic equilibrium of System (1) is provided . The basic reproduction number R 0 , defined as the number of secondary infections generated by a typical infectious mosquito in a fully susceptible population, is evaluated and the next generation matrix method [22] is used to establish the basic reproduction number R 0 , furthermore a seasonal reproduction number R s which varies with respect to time is established.  Per capita rate of progression of humans from long term of exposed state to the infectious state(days −1 )

Model steady states
T m Per capita rate of progression of mosquitoes from the exposed state to the infectious state (days −1 ) 0.08 matrices F and V defined by F = DF(E 0 ) and V = DV(E 0 ) are given in partitioned form as follow Notice that F is a non-negative matrix of rank 2 and V is non-singular, both were evaluated at E 0 , [16], for the derivation of the spectral radius ρ(F V −1 ) = ρ(F 1 V − 1 1), the expression for the basic reproduction number R 0 is given by The following Lemma shows stability result with R 0 for the disease free equilibrium E 0 by using Theorem 2 from van den Driessche and Watmoungh [22].
the eigenvalues for the matrix F 1 − V 1 have negative real parts given by −(T m +2b m )/2, additionally the eigenvalues for matrix −V 2 are all negative and given by −(T s h + b h ), −(T l h + b h ) and −(r + b h ), therefore R 0 < 1 and the disease free equilibrium E 0 , is locally asymptotically stable.  Proof. By setting the left-hand sides equations of System 1 to zero, the following relations are established:

The seasonal reproduction number
In existing ODE model of malaria, the number of mosquito population M , has been set to the constant value. Using the number of mosquito population follows the Gaussian function such that it has the peak in July and August, we calculate the seasonal reproduction number R S differently than existing basic reproduction number. The basic reproduction number is calculated to using the human population H = 10000, and the mosquito population M = 100. Thus, we used M = M/H = 0.01

M (t) = M (t)/H
This function has the characteristics of rapid increase in summer and autumn followed by decrease in spring and winter. As following R s :

Simulations
The following table contains values that are varied in the following simulations.

The seasonal reproduction number sensitivity
We now consider the seasonal reproduction number, R s with changes to the values of parameters b m , r, σ, T m , p, in the ODE model and whether and now the seasonal reproduction number changes. Using this data the sensitivity of R s with respect to the five parameters can be estimated.
In Figure 2, b m , r, T m and σ are separately varied. The plots show changes in the seasonal R s for each t over a three year period and we note is that we have implications of sensitivity with respect to a parameter if for each value we see different plots. The first parameter investigated is the newly emerging mosquito coefficient which was varied among five values, b m ={0.1, 0.5, 0.7949, 1.25, 1.50}. The sensitivity is the greatest as b m increases from 0.1 to 0.5 as R s decreases as b m increases. In Figure 2 observe that while the value of b m is increasing, R s is decreasing. That is, if the number of newly emerging mosquitoes is larger, then the number of secondary infectious human decreases. Next, vary r =  is increasing, then R s is also increasing. That is, as the number of infected mosquitoes is larger, then the number of second infected human is increasing. In Figure 3, where the parameter p,probability of exposed humans going through short-term incubation periods, is varied p = {0.25, 0.50, 0.75}there is almost zero change with implications of little sensitivity of R s as for each value we see almost identical plots.

Infected mosquito/human sensitivity
Now consider the population of infected mosquitoes or humans over time and vary the parameters b m , r in Figure 4. We observe greater sensitivity for the smaller values of the parameters and then not as much change in the infected populations for the rest of the values. It is not surprising to see an increase in the infected populations at the peak summer months of mosquito prevalence. In Figure 5 where T m , σ are varied the corresponding mosquito and human infection populations are also affected. The conclusion is that the model is also sensitive to these parameters. Again note that there is an increase in I m , I h during peak mosquito months as this model considers variable mosquito population M (t). Figure 6 shows the infected mosquito and human populations as p is varied and these show minimal sensitivity to this parameter.
At first, b m , a per capita rate of newly emerging adult mosquitoes, is varied. The default value of b m = 0.7949 is set and the model is evaluated at a total of five values({0.1, 0.5, 0.7949, 1.25, 1.50}).
As we can see in Figure 4, as the value of b m increases, that is as newly emerging mosquitoes increases, infectious humans and mosquitoes appear to decrease. Next, five values of r = {0.01, 0.15, 0.30, 0.45, 0.50} are considered. As r is a per capita rate of progression of humans from the infectious state to the susceptible state, we see that contrary to b m , the graph of infectious mosquitoes and infectious humans both look like a flat lines except at the lowest value of r = 0.01. As the value of r increases, naturally infectious humans decrease and consequently uninfected mosquitoes are less likely to get infected. Third, in Figure 5, σ ={0.25, 0.31, 0.37, 0.44, 0.50} and observe infected mosquitoes and humans. Note that σ is average number of contacts made to host by a single mosquito. Infectious humans shows a marked increase at the peak mosquito population around the 20th week per year. Infectious mosquito population, for all values of σ follows a similar trend. That is, the effect of σ to infectious human and mosquito population is bell shaped at peak mosquito growth which follows the shape of the seasonal mosquito population considered in this model (2). Also in Figure 5, consider the model with five values of T m = {0.05, 0.08, 0.120.160.20}. T m is a per capita rate of progression of mosquitoes from the exposed state to the infectious state. T m , like σ, has an impact on both infectious human and mosquito populations. As T m gets larger the proportion of infectious status mosquitoes and human population shows a greater increase during the months of greatest mosquito prevalence compared to that at lower values of T m . As T m increases then naturally the number of infectious mosquitoes increases and we observe this reflected in both infected groups.
In Figure 6, where the parameter p,probability of exposed humans going through short-term incubation periods, is varied p = {0.25, 0.50, 0.75} there is change, mostly in the first of three years with implication of some sensitivity of I h and I m , the infected populations. Again note the increase in the summer peak mosquito months.

Long and short term incubation period sensitivity
When infected mosquitoes bite human hosts, the exposed human has an incubation period before developing malaria. This incubation period is largely divided into two types, the first one being a short term incubation period. Historical data shows the short term incubation period has a mean equal 25.5 days and follows the gamma distribution [?, ?] [17]. In previous models a constant value was assumed for short term incubation in the ODE model with the value of short term incubation period set to 25.9. This research extends the previous ODE model by using the three values 7, 25.5, 50 for this parameter for comparison in a sensitivity analysis. The second period considered is long term incubation and has the mean of 329.4 days and follows a distribution similar to the normal distribution [16], [17]. Similarly, rather than use a constant value of 329.4, this research will apply three values 250, 329.4, 360. In using three values each of the two incubation periods, the model is more robust as it allows for variability in the incubation periods and correspondingly estimates how infectious human and mosquitoes interact. In Figure 7 the infected mosquito and human populations are shown as parameters T s h and T l h are varied. We note that the infected human population shows sensitivity to these values with the infected mosquito population less so.
In Figure 7 where the parameters T s h (top two plots) this first  parameter investigated is the short incubation period rate which is varied among three values, 7, 25.5, and 50 (weeks) for corresponding rates, T s h ={1/7, 1/25.5, 1/50} and T l h (bottom two plots), long term incubation period with the exposed human incubation periods 250, 329.4 and 360 (weeks) and corresponding values for this rate are T l h = {1/250, 1/329.4, 1/360}. We note the curves for infected humans do vary as T s h and T l h are varied with peaks at summer months per year. The infected mosquito population appears to be unaffected.
In Figure 8, we vary parameters T s h , T l h and the plots show  seasonal R s for each. We observe that while the value of T s h is decreasing, R s stays the same. That is, if the mean short term rate of incubation is increased from 7 to 25.5 and then to 50, then the number of secondary infectious humans does not seem to change with respect to shorter incubation period. Upon close inspection there is a slight change with the shortest incubation period corresponding to higher peaks in , R s (see third plot in Figure 8). Consider the model with T l h decreasing with corresponding incubation periods 250, 329.4, 360 days, R s again appears unchanged. The constant nature implies that as the incubation period for malaria in humans increases, then the number of secondary infectious human does not seem to be impacted. Again upon zoomed in view (fourth row in Figure  8), we see some differences, however with respect to scale they are negligible.

Conclusions
In this paper an epidemiological model for the reemergence of malaria in Korea is presented. We investigate the impact of variable mosquito population as well as add variables long and short incubation periods on the transmission dynamics of malaria. The SEIS model is based on malaria infected mosquitoes which bite humans or animals. This model studies plasmodium vivax malaria and including variables related to animal population and mosquito attraction to them. The standard basic reproduction number of the ODE model is evaluated using regular methods from epidemiology, furthermore we concentrate more on the seasonal reproduction number for simulation purposes since the mosquito population varies throughout the year. The existing timeindependent Malaria population ODE model is extended to a time-dependent model with simulations reflecting three years of variable mosquito population. Also, using bi-modal Malaria incubation, changes to the infectious population from constant incubation period is extended in the ODE model and also includes seasonal mosquito prevalence over a three year period. Our interest is also to find an explicit expression for the endemic equilibrium as presented, and stability analysis for the model is conducted.
Sensitivity analysis of the parameters included studying R s , seasonal reproduction number, and infection populations I h and I m in response to changes in parameters b m , r, T m , σ and p. In conclusion, b m , r is negative correlation to R s while σ, T m is positive correlation to R s . The seasonal reproduction number (R s ) showed sensitivity to all parameters except the probability of exposed humans going through short-term incubation (p). Infected mosquito populations were not sensitive to the variable of human bi-modal incubation, but infected human population showed some sensitivity to these variables.
The results presented here extend ideas from [16], [17] by performing endemic equilibrium analysis and DFE with stabil-ity and inclusion of variable mosquito population, M (t) in the reproduction number analysis and ODE simulations. Also parameters are varied with the seasonality of the mosquito population in the model with simulations presented. Although both the infectious human and mosquito populations or compartments show differences as the four parameter are changed, the impact of the parameter changes is most noticeable during peak mosquito weeks.