Non-linear Dynamics of Two-Patch Model Incorporating Secondary Dengue Infection.

In this paper, the impact of human migration on the dynamics of dengue epidemic has been discussed. The vector-host model considers two patches with different dengue serotype in each patch. The model considers the constant rate of migration in susceptible and recovered class from one patch to other. Recovered migrants from prior infection are exposed to secondary infection in the patch where different serotype is present. The basic reproduction number is computed and analyzed in terms of migration parameters. The model is analyzed for the existence and local stability of various equilibrium states in terms of migration parameters. The numerical simulations for the choice of relevant data from literature have been performed to verify analytical results and to further explore the dynamics of the system. The sensitivity analysis of basic reproduction number with respect to migration parameters is carried out. It is found that immigration in a patch increases the basic reproduction in respective patch and vice-versa. The basic reproduction number has been estimated for the two states of Brazil which verifies the occurrence of severe epidemic in one of the states of Brazil. © Springer (India) Private Ltd., part of Springer Nature 2017.

Introduction

Dengue was confined to only nine countries before 1970. However, it has now been spread to more than 100 countries in the tropical, subtropical and temperate areas of North America, South America, Africa and Southeast Asia. According to one estimate, it affects 50–100 million people every year [1]. In fact, World Health Organization has declared dengue as the fastest communicable mosquito-borne disease in the world in terms of human morbidity and mortality [2, 3]. The causes of emergence and resurgence of dengue disease are attributable to multiple factors including urbanization, fast transportation, economic development and changes in human behavior. The fast increase in modes of transportation, their reach, speed and efficiency led to the global movement of human population [2, 4]. This may be responsible for the spread of infectious diseases from isolated locations to the new regions. It has been reported that another 100 million people travel to dengue affected areas where they not only increase the risk of contracting the disease but also spreading it further afield.

Dengue fever exhibits up to four closely related and distinct dengue viruses or serotypes. The origin of all four dengue serotypes is found to have in Asian forests and subsequently they have spread worldwide because of the migration of human and commerce [5]. The mosquitoes Aedes aegypti and Aedes albopictus are the main carriers of these viruses. It is transmitted through infected mosquito bites. Once a person becomes infected by a serotype, he/she will never be reinfected by that serotype due to induced lifelong immunity against it. However, only temporal cross immunity is developed against the other serotypes [6]. Therefore, the recovered individuals from a serotype may get secondary infection when exposed to a different serotype. Primary infections are generally asymptomatic and nonfatal, while secondary infections sometimes lead to the fatal forms of dengue infection, namely, dengue hemorrhagic fever (DHF) or dengue shock syndrome (DSS) [1, 7]. Mathematical modeling has been applied to understand and control the many vector-borne diseases like dengue, malaria, leishmaniasis, chagas diseases etc [8-13]. Some mathematical models on dengue dynamics incorporating primary and secondary infection have been proposed and analyzed by various investigators [11-16].

Due to short flying range of mosquitoes, the role of human migration becomes crucial in dissemination of dengue serotypes from one part to other parts of the region [17]. Some studies on dengue have also confirmed that human migration is the main cause of its spatial spread [18, 19]. In particular, human movement was identified as an important factor for the spread of dengue serotype-1 from the southern region of Vietnam to its northern and central regions [20]. Also, human migration from one patch (area) to another increases the chances of secondary infection [6, 21]. During 1977–1978, the dengue epidemic in Cuba occurred due to DEN-1 serotype with not a single case of DHF [22]. However, human migration from endemic countries (with DEN-2 serotype) to municipalities in Cuba has been reported to be responsible for DHF/DSS cases of subsequent epidemics in 1981 and 1997 [22, 23]. These facts exemplify the two largest outbreaks of secondary infection cases by DEN-2 in Cuba in 1981 and 1997 [22-24]. Further, human migration may be one of the reasons for the spread of secondary infection.

A very limited but relevant metapopulation epidemic models involving movement of individuals between discrete spatial patches have been proposed by investigators [25-31]. In particular, such compartment models for multi-patches have been discussed for measles, influenza and SARS [25, 28, 29, 31]. Prosper et al. [32] have discussed Ross-Macdonald type two patch SI metapopulation model for malaria by considering different degrees of transmission in patches. Gao and Ruan [27] have derived a multi-patch model to study the impact of population dispersal on spatial spread of malaria between patches. A limited work on metapopulation model incorporating host-vector dynamics has been carried out. A two patch vector-host model incorporating secondary infection of dengue disease has been discussed by Gakkhar and Mishra [33]. They have assumed no vector migration and human migration is allowed only in one direction.

In this paper, a two-patch vector-host model for the dengue epidemic incorporating secondary infection has been proposed and analyzed. The model considers human migration between two discrete patches having different dengue serotypes. Due to imperfect cross immunity towards heterologous serotype, the recovered migrants of patch-1 may get secondary infection in patch-2. The impact of human migration on persistence or extinction of disease in patches has been discussed by considering all possible combinations of migration parameters. Stability analysis and numerical simulations have been carried out for various equilibrium states. Further, sensitivity indices of basic reproduction numbers for migration parameters have been computed. In the case study, the estimation of basic reproduction has also been performed for the two states of Brazil. The conclusion is given in the last section.

Formulation of Model

Keeping the above facts into consideration, let there be two distinct serotypes of dengue namely serotype-1 (DEN-1/DEN-3/DEN-4) and serotype-2 (DEN-2) prevalent in two discrete patches: patch-1 and patch-2 respectively. Human (host) population is compartmentalized into and for susceptible, infected and recovered individuals respectively in the ith patch, . It is assumed that the infected individuals are unable to migrate due to disease-induced weakness. The susceptible individuals of patch-1 (patch-2) may become primary infected with serotype-2 (serotype-1) when they migrate to patch-2 (patch-1). The recovered migrants from patch-1 to patch-2 may become secondary infected when exposed to serotype-2 and form a new compartment of secondary infection in patch-2. As, the serotype DEN-2 is more associated with secondary infection [7, 23, 24, 34, 35], the almost perfect cross immunity has been assumed for the recovered migrants from patch-2 to patch-1. The probability of getting secondary infection from serotype-1 is assumed to be negligibly small in patch-1. Accordingly, the recovered migrants from patch-2 to patch-1 will not get secondary infection. The secondary infected population is not being considered ( in patch-1. Further, this simplifying assumption reduces the dimension of the model.

Let be the constant rate of recruitment in susceptible class and be the natural death rate of host in both the patches. The disease-induced death rate e is considered only in the secondary infection compartment of patch-2. Let and ; be the rates of transmission of infection and recovery of the host population respectively in the ith patch. The proportion of recovered migrants from serotype-1 become secondary infected by serotype-2 at the rate in patch-2. Let and be the rates of migration from patch-1 to patch-2 and vice-versa respectively.

Let and be the susceptible and infected vector population in the ith patch. Due to short life span of vector, once infected they never recover. Accordingly, no recovery class has been considered for them. Further, let be the recruitment rate of mosquitoes in the absence of vertical transmission. Let and be the constant natural death rates of mosquitoes in patch-1 and patch-2 respectively. The transmission rate of infection to mosquitoes by infected host is assumed to be . The mosquitoes are not migrating and their dynamics remains the same irrespective of the patch. Considering all parameters to be positive, the following model has been formulated: The model is associated with following non-negative initial conditions:The schematic diagram has been drawn in Fig. 1 showing the transmission dynamics of disease in two patches. Note that for simplicity, the mortality rates of mosquitoes are considered to be equal in both the patches i.e. . Consider the following non-dimensional parameters:The dimensionless model is given as: The model is associated with following non-negative initial conditions:

Fig. 1

Schematic diagram representing the transmission dynamics of disease in two patches

Analysis of the Model

Consider a system of differential equations in as

Lemma 1

The closed setis positive invariant for the non-linear system (12)–(22).

Proof

Let N(T) be the total host population with where = and = be the total host population at time T in patch-1 and patch-2 respectively. Further, let M(T)= be the total vector population where, ; .

By adding the host dynamics from the model (12)–(22) gives,Using standard comparison theorem [36],Again, by adding vector dynamics in from the model (12)–(22) gives,Using standard comparison theorem [36],Also, and for and respectively. This shows that the solutions of the system (12)–(22) converge towards the set .

The Lemma 1 shows that the all solutions of the non-linear system (12)–(22) are non-negative and bounded. Therefore, the model is mathematically as well as biologically well behaved. In next subsection, the existence of equilibria for the model (12)–(22) has been discussed.

Equilibrium Points

The non-linear system of Eqs. (12)–(22) have four equilibrium states :

The disease-free state () = ,

The disease in patch-1 state () = exists for The expressions for non-zero state variables at equilibrium point are given as follows:

The disease in patch-2 state () = exists for The expressions for non-zero state variables at equilibrium point are obtained as follows: where, and

The endemic state ()= exists under the following conditions: where, The expressions for the state variables at equilibrium level are omitted as they are lengthy and complex.

Basic Reproduction Number

The basic reproduction number is computed by next generation approach [37]. The details of the new infections and transfer matrices are given in “Appendix”. The basic reproduction number, (say) is given asorwhere, and are the basic reproduction numbers of patch-1 and patch-2 respectively.

Observe that there will be no effect of migration on the basic reproduction numbers and when .

Stability of Equilibrium States

This section analyzes the stability of the non-linear system (12)–(22) for various equilibrium states by computing eigenvalues of the Jacobian matrix J[E]. The general matrix J[E] of the system is given in “Appendix”.

Theorem 1

The disease-free state is locally asymptotically stable for

About disease-free state, seven of the eigenvalues of the Jacobian matrix () are negative and computed as (multiplicity 4), (multiplicity 2) and -1 - . The remaining four eigenvalues are given as:Further simplifications give all eigenvalues with negative real part under condition (27). Therefore, the disease-free state is locally asymptotically stable under condition (27).

Further, for global stability of disease-free state, the following theorem is concluded:

Theorem 2

The locally asymptotically stable state is also globally stable for

For arbitrarily chosen positive constants A, B, C and D, consider the positive definite function as:Taking derivative of with respect to t and its simplifications yield,For to be negative, the conditions areLet us choose and . Their substitution in inequalities and further simplifications giveLet and for positive and .

then the inequalities can be combined to giveThe positive arbitrary constants C and D can now be chosen to satisfy the above inequality ensuring that is negative. Accordingly, the function is a Lyapunov function for the condition (28). At , , and , becomes zero. If , , and , then is the only largest invariant set that contains a subset in which all these variables are zero. By applying LaSalle’s invariance principle [38], all trajectories in the closed set approach the equilibrium point . Hence, the locally asymptotically stable disease-free state is also globally asymptotically stable. Hence, the result is proved.

It is observed that may still be globally stable even when one or both of the conditions in (28) are not satisfied.

Further, when the state is unstable, existence of some other states may become possible. These possibilities are explored next.

It may be noted that the existence of requires [see condition (23)] which implies the instability of disease-free state .

Lemma 2

When then

Assuming leads to which is the contradiction. Therefore, whenever .

Theorem 3

For given in (29), the state will be locally asymptotically stable when

For the local stability of , the seven of the eigenvalues of the Jacobian matrix are computed as:

(multiplicity 3), , , .

Note that all of these seven eigenvalues are negative provided condition (30) is satisfied. The other four eigenvalues are the roots of the polynomialwhere,Now, substituting the values of state variables at and simplifying, it has been found that all the four roots of the polynomial are having negative real part for , if the Routh–Hurwitz conditions satisfy. Hence, the state will be locally asymptotically stable if the Routh–Hurwitz conditions hold.

Remark 1

Accordingly, state is stable when . However, when , the state may be stable/unstable subject to the condition (30). These are further explored in terms of migration:

Observe that in absence of migration (). Therefore, in absence of migration the state will be locally asymptotically stable when .

When migration is allowed only in patch-1 i.e. but , it is observed that and the state is locally asymptotically stable for .

When and then by simplifying the expression for , it may be noted that still remains smaller than 1. Therefore, the state is again locally stable for .

Remark 2

Stability of the state is possible even though since the stability condition may still be satisfied for sufficiently small in presence of migration.

Local stability of (disease in patch-2 state) state is discussed below:

Lemma 3

When then

Assuming leads to which is the contradiction. Therefore, whenever .

Theorem 4

For given in (31), the state will be locally asymptotically stable when

For the local stability of =, the seven eigenvalues of the Jacobian matrix are given as follows:

(multiplicity two), , and .

Note that all of these seven eigenvalues are negative for condition (32). The other eigenvalues are the roots of following polynomial of degree fourwhere,Now putting the values of variables at state , all the four roots of the above polynomial are having negative real part for provided the Routh–Hurwitz conditions satisfy. Hence, the state will be locally asymptotically stable if the Routh–Hurwitz conditions hold.

Remark 3

Accordingly, state is stable when . However, when , the state may be stable/unstable. These are further explored in terms of migration:

When both the patches are isolated i.e. and Assuming , the state will be locally asymptotically stable when .

When migration is allowed only in patch-2 i.e. =0 but , then it is observed that gives the local stability of for .

When and then by simplification yields . Therefore, the state is again locally stable for .

Remark 4

Stability of the state is possible even though since the stability condition may still be satisfied for sufficiently small in presence of migration.

Remark 5

When and , the states and both exist. when or , the states or respectively are locally stable while the state does not exist. The state exists for and .

For the local stability of the endemic state , the three of the eigenvalues of Jacobian matrix are (multiplicity two) and . Two of the eigenvalues are the negative roots of second degree polynomialFurther, the remaining eigenvalues are the roots of the polynomial of degree 6:The expressions for the coefficients of above six degree polynomial omitted from the text as they are lengthy and complex. The numerical simulations have been performed for the stability of the endemic state in the next section.

Numerical Simulation

Numerical simulation of the system (1)–(11) are carried out for the data given in Table 1. For different combination of migration parameters and , the thresholds responsible for the existence and stability of states , and have been computed in Table 2.

Table 1

Parameters values

Parameters	Values
β_1	(0, 0.05) [39]
β_2	(0, 0.05) [39]
γ_1	0.3428 [40]
γ_2	0.3667 [40]
μ	0.0000457 [41]
μ_1	0.0714 [39]
ω	0.001
ω_1	400 [41]
ρ	(0, 1)
σ	(0, 0.05) [39]
e	0.001 [41]

Table 2

Computation of and for different combination of migration parameters

	m_1	m_2	R_10	R_01	R012ξ	R102η
Case 1:	0.001	0.005	9.99268	4.17131	0.636259	8.20645
Case 2:	0.005	0.001	4.50934	9.25027	6.24823	0.775578
Case 3:	0.004	0.005	8.16928	6.76781	1.1458	2.1211

For case 1, the existence condition (25) of the state is violated while both the states and exist and are given below:According to stability conditions (30) and (32) of the states and respectively, is found to be locally stable while gets unstable.

A projection of phase plot in hyperplane has been drawn in Fig. 2 for the set of initial conditions , , and chosen in the neighborhood of . The solution trajectories are found to be converging to the state in hyperplane .

Fig. 2

Starting from initial conditions and chosen in the neighborhood of converges to the state

Again, for case 2 in Table 2, the existence condition (25) of the state is violated while both the states and exist.The Fig. 3 shows a projection of phase plot in hyperplane for the initial conditions , , and chosen in the neighborhood of . It is observed that the all the trajectories starting from the neighborhood of converge to the state .

Fig. 3

The initial conditions , and chosen in the neighborhood of showing the stability of the state

For case 3, the states , and exist:The initial conditions and have been considered in the neighborhood of and the and have been considered in the neighborhood of . It is observed from a 3D phase plot drawn in Fig. 4 that all the solution trajectories converge to in hyperplane . This verifies the instability of the states () and (). Further, this confirms the stability of .

Fig. 4

Starting with the initial conditions , in the neighborhood of and , in the neighborhood of , the solution trajectories converge to the state

Sensitivity Analysis

In disease modeling, sensitivity analysis is important to know the change in basic reproduction number () in response to the parameters it involves. The control measures can then be targeted on the basis of the sign and magnitude of the sensitivity indices with respect to various parameters. To compute the sensitivity of parameters of different scale, the normalized forward sensitivity index has been used [42]. The normalized forward sensitivity index of for the parameter p (say) is defined asConsider basic reproduction numbers and in terms of dimensional parameters asThe sensitivity indices with respect to the migration parameters () for the two patches are obtained as follows:Due to emigration from patch-1 is found to be negative. The same is true for . However, and are observed to be positive. This is due to immigration in respective patches.

From the above expressions, it can be observed that sensitivity indices for the migration parameters depend on host death rate also. For , the sensitivity indices with respect to migration parameters are given in Table 3. It can be concluded that when rate of immigration is higher than emigration in a patch, the basic reproduction number of respective patch increases. Particularly, from given Table 3, if increases by 10%, the basic reproduction number for patch-1 decreases by 1.238%. However, there will be increase of 3.664% in patch-2.

Table 3

Sensitivity indices of and for and

Parameters	Sensitivity index for R102	Sensitivity index for R012
m_1=0.001, m_2=0.003	Υm_1R_10=-0.1238, Υm_2R_10=0.1254	Υm_1R_01=0.3664, Υm_2R_01=-0.3713
m_1=0.003, m_2=0.001	Υm_1R_10=-0.3718, Υm_2R_10=0.3670	Υm_1R_01=0.1254, Υm_2R_01=-0.1237

Case Study: Estimation of Basic Reproduction Number

In this paper, a two patch model incorporating human migration is applied for 2003 dengue outbreak in the two states of Brazil namely, Ceara and Rio de Janeiro. The first case of DEN-2 serotype in Rio de Janeiro was reported in 1990 then spread over to other states of the country namely Ceara, Bahia, Rio Grande do Norte, Alagoas, So Paulo, Mato Grosso do Sul, Mato Grosso, with the highest incidence rate in northeast region [43]. During summer 2002 a large outbreak of dengue due to DEN-3 serotype occurred in Rio de Janeiro [44]. In 2003, the human migration rate in Ceara was quite high as compared to Rio de Janeiro [45]. The sequential infection by DEN-2 serotype led to severe secondary infection cases in the form of DHF in the state Ceara [46, 47]. Keeping these facts in mind, let us consider the state Rio de Janeiro as patch-1 with DEN-3 serotype prevalent at time t and the state Ceara as patch-2 where DEN-2 serotype is present. The estimation of has been performed for the two patches from the initial growth phase of the epidemics [48]. Let us consider that at the beginning of the epidemic, the cumulative number of cases, K(t), varies exponentially aswhere, is the force of infection and is constant. Accordingly, the infected host and vector population of both the patches are assumed asAgain, the number of non-susceptible hosts and vectors can be assumed negligible,Now using Eq. (34) for infected host and vector population in the model (1)–(11) gives, Multiplying Eqs. (36) and (37) together for patch-1 and Eqs. (38) and (39) for patch-2 respectively give the following: The and are estimated reproduction numbers for the patch-1 and patch-2 respectively.

Now, from the Eq. (33), it can be easily seen that, at the beginning of the epidemic number of new cases in a month in the patch-1 and patch-2 (say and respectively) would be proportional to the cumulative number of cases, i.e. and . Again by plotting the number of new cases per month against the cumulative number of cases K(t) for the two patches separately, the force of infection which would be the slope of respective curves, can be obtained [49].

Using the monthly data of the two states of Brazil namely, Rio de Janeiro (say patch-1) and Ceara (say patch-2) of 2003 from “WHO Dengue NET”, the commutative number of cases have been plotted against new number of cases in Figs. 5 and 6 respectively. The force of infection of patch-1 and patch-2 have been given in Table 4 by computing the slope of respective curves. The time series for the data of respective states has been given in Tables 5 and 6 in “Appendix”. The host parameters (host death rate and recovery rate) have been taken from literature [40]. The estimation of mosquito mortality rates are described in next subsection.

Fig. 5

Force of infection for the state Rio de Janeiro

Fig. 6

Force of infection for the state Ceara

Table 4

Parameter estimation of mosquito mortality rates and basic reproduction number for two states of Brazil

Parameters	Rio de Janeiro	Ceara	Source
Force of infection	0.0072/month	0.2937/month	Estimated
Mosquito mortality rate	0.0319026(=μ_1)per day	0.03075193(=μ_2)per day	Estimated
Basic reproduction number	1.009178(=R10E)	1.423395(=R01E)	Estimated

Table 5

Monthly dengue cases in Rio de Janeiro in 2003 from WHO DENGUE NET

Months	Number of infectives
1	1765
2	2311
3	1742
4	972
5	523
6	245
7	119
8	77
9	69

Table 6

Monthly dengue cases in Ceara in 2003 from WHO DENGUE NET

Months	Number of infectives
1	1195
2	2335
3	4143
4	4264
5	4850
6	2844
7	1159
8	548
9	336

Estimation of Mosquito Mortality Rate

The mortality rate of mosquitoes is temperature dependent. By using the average temperature records of nine months of the states Rio de Janeiro and Ceara, the estimation of mortality rate of the mosquito Aedes aegypti has been explored by the enzyme experiment [50, 51]. The following formula is obtained to calculate the mortality rate of Aedes aegypti mosquito:where, T is the temperature and is the mortality rate of mosquito.

By using the monthly temperature records of the states Ceara [52] and Rio de Janeiro [53] and substituting in Eq. (42), the range of mortality rate for the two respective states has been calculated. The mortality rate of Aedes aegypti mosquito is found to be in the range of 0.02987767–0.0363767 for the state Rio de Janeiro. However, for the state Ceara, it is coming in the range of 0.02714614–0.03403341. Using the average value of the mortality rates, the estimated value of basic reproduction number of the respective states are given in Table 4.

It is concluded from the Table 4 that the basic reproduction number for the Ceara state of Brazil is greater than that of the Rio de Janeiro state. This verifies the fact that severe epidemic was occurred in Ceara in 2003. It supports that human migration is one of the reasons for the severe outbreak of DHF.

Conclusion

In this paper, a non-linear two patch dynamic model has been proposed to study the dynamics of dengue transmission. It is assumed that the two distinct serotypes are predominant in respective patches. The main emphasis is given on inter-patch migration of human population and its impact on spatial spread of primary as well as secondary infection. The existence conditions for equilibrium states are obtained in terms of . The disease-free state is found to be globally stable provided the basic reproduction numbers () of both the patches are less than one. It is found from analysis that for the equilibrium state (), where disease persists in patch-1 (patch-2), it will persist there irrespective of migration if the basic reproduction number of patch-2 (patch-1) is less than 1. Numerical implementations have been carried out for the relevant data from the literature to explore the stability of endemic state. For the set of data where exists, it is found that starting from the initial conditions in the neighborhood of the states and , the solution trajectories converge to the state . This shows the stability of the state .

The magnitude of migration parameters and their mutual relationships are very crucial to predict the stability of various states. When equal rates of migration are considered in both the patches, the status of disease will not be affected due to migration. In other words, the disease-free state will not become endemic or vice-versa in any patch due to migration. However, the level of primary/secondary infection changes without affecting the disease-free/endemic status. On the other hand, when rate of migration is different in the patches, the persistence/extinction of infection will depend on the magnitude of the basic reproduction number of patches as well as the migration parameters.

Sensitivity indices of the basic reproduction number for the migration parameters are computed to analyze the effect of emigration and immigration in respective patches. Accordingly, it is found that emigration from a patch decreases the basic reproduction number of respective patch which consequently increases the chances to establish the infection-free state. However, higher rates of immigration in a patch may increase the infection level.

The case study for the two states of Brazil namely, Rio de Janeiro and Ceara is carried out to validate the model for the sequential infection by DEN-2 and DEN-3 serotypes occurred during 2003. It has been found from literature that Ceara had more immigration as compared to Rio de Janeiro in 2003 because of this Ceara had more severe dengue epidemic in 2003. By estimation of basic reproduction numbers for these two states, the basic reproduction number for Ceara state is found to be higher than the Rio de Janeiro state. This is in line with the fact that severe epidemic was occurred in Ceara in 2003.