Standard routine techniques of modeling of tick-borne encephalitis

1 Introduction

Vector-borne disease is one of the most discussed viral diseases in the contemporary epoch. When vector such as mosquito, tick and bug bites a susceptible host, the infection paves the way for the severe disease. Tick-borne disease is the worst disease among vector-borne diseases. Recently, a mind-boggling disease tick-borne encephalitis (TBE) has been found as the widespread one in European states as well as in some Asian countries [38]. This is an emerging arboviral morbidity characterized by severe acute and chronic neurological infections in humans. The disease is caused by a flavivirus (TBE virus [TBEV]), which is transmitted mainly by Ixodes persulcatus ticks from a wild vertebrate host to humans, domestic animals or through consumption of unpasteurized milk products [5,7].

TBEV is endemic in central, northern and eastern Europe, Russia and the Far East, including Mongolia, the northern part of China and Japan from Asia [10]. Globally, the number of clinical cases of TBE is estimated to be amidst 10,000 and 12,000 a year and in Europe more than 3,000 human TBE cases are hospitalized each year [13].

Various parts of Europe, Central Asia and East Asia are the endemic countries of TBE. The people who travel, hike and do camping in forested areas are at high risk. Specifically, the span of March to November is crucial for the spread of this disease [22]. Outbreak of TBE in Germany by the pathogen Ixodes ricinus has been analyzed in a recent article [17]. This disease has gained momentum and attracted biologists as well as mathematicians.

The clinical symptoms of the infection are biphasic [8,10,33]. The incubation time is between 2 and 28 days, with an approximate average of 7 days. The first phase known as viremic phase lasts between 2 and 8 days, and symptoms of this phase are slight fever, flu, headache etc. The second phase is critical in which the virus spreads to the central nervous system and causes anorexia, fever, headache, vomiting, photophobia and often sensory changes such as visual disturbances, paresis, paralysis or even state of coma. Death might occur as soon as 1 week after the onset of clinical morbidity.

World Health Organization recommends immunization for all people living in the areas where the disease is common or supposed to be outbroke [11]. There are two types of immunity one can develop after vaccination [22]. Partial immunity refers to half or less than the required vaccination rate, and partially immune people live under the risk of getting TBE. Complete immunity refers to full immunity after getting the required vaccination rate [8,10,33]. However, the dose of vaccination is different in different countries. For instance, in Austria, the recommended basic vaccination schedule consists of two vaccinations ≈4 weeks apart followed by a third vaccination after 5–12 months and a fourth vaccination after >3 years. For persons <60 years of age, additional booster immunizations are recommended every 5 years; this interval is reduced to 3 years for persons >60 years of age [32].

The beauty of mathematics is well explored when it is applied on real life phenomena. It has been said that both Mathematics and Biology are opposite terms out-and-out with respect to their applications in real life. However, mathematical modeling on different contagious diseases has helped scientific research to get the solution more accurately and commendably, and for more details see refs. [3,4,14,15,24,40]. In the field of biosciences, mathematical modeling on different diseases got momentum in the twentieth century. Models on diseases such as Zika virus [23] and mathematical model on gonorrhea provide ample assimilation regarding disease as well as the way to control it [2,6,9,20,21,25,28,29,30,35].

Mathematical models are of great importance in natural sciences, particularly in physics. Physical theories are almost invariably expressed using mathematical models. In the early twentieth century, William Hamer and Ronald Ross applied the law of mass action to explain epidemic behavior of diseases. Two aspects are involved in the initial formulation of the law of mass action: (1) the equilibrium aspect concerning the composition of a reaction mixture at equilibrium and (2) the kinetic aspect concerning the rate equations for elementary reactions. Hence, physics plays an elementary role to determine the epidemic behavior of infectious diseases such as chicken pox, measles and TBE.

In a previous study, a mathematical model on TBE was constructed [17]. Question arises if there was already a model constructed on the said disease, what was the need to construct another one or specifically what makes this model unique and different to the cited one? The answer lies in the fact that the kind of extension of susceptible, infected, recovered (SIR) used in both papers is different. The assumptions made to construct a model were different than the assumptions made by us, which is totally based on the facts taken from authentic platform such as World Health Organization. TBE along-with this new proposed model, we have given a try to derive formula for basic reproduction number or more specifically Ro-like threshold which explains the required threshold yielding to get disease outbroken in the area premises.

In Section 2, methodology is discussed, which is opted while modeling TBE. Next-generation method is used to find out the general formula for basic reproductive number. Transmission matrix and transition matrix are used to calculate the spectral radius, which further helps us in deriving the formula for basic reproductive number.

In Section 3, final results are described, which comprise the final formula for R ₀-like threshold. This formula is amply capable to decide whether the disease will be outbroken in the regional premises or not on the basis of parameters collected from the area concerned. In Section 4, equilibrium points and their local stability are discussed. Numerical computations and graphical analysis are presented in Section 5. In Section 6, the concern and importance of this model are discussed. Some future perspectives are described in Section 7.

2 Description of the basic model of TBE

SIR model is commonly used for mathematical modeling of such type of problems [27]. This model is a modified and extended version of an SIR model in which there are two categories of the population, human population and tick population, where three compartments of human population and two of ticks persist. For the sake of simplicity, different parameters are used to represent different characterizations for human as well as tick population. The parameters are described as follows:

2.2 Compartments of tick population:

1. susceptible ticks; S _t,

2. infectious ticks; I _t.

Both the species are restricted to their own compartments and are strictly prohibited to move from one compartment to the compartment of other specie, which means that neither humans nor ticks would leave their respective specie compartment. It is pertinent to mention that immunity from the disease after the second phase of infection is not possible. If an individual has gone through the required dose of vaccination, then he is said to be attained complete immunity and is free from all sorts of the risk emerging from TBEV. However, a person who is partially immune and has not attained the required dose of vaccination is at risk and there are certain chances for him to get affected. There is a chance for an exposed person to recover and return to the susceptible compartment after recovery before becoming infectious. The susceptible human “S _h” becomes exposed when an infectious Ixodes tick “I _t” interacts with him and transfers the TBEV to him at the rate of “β ₁.” After this, susceptible human “S _h” becomes exposed “E _h” and moves to the exposed compartment. A person who has not been given proper treatment before the virus finds the cell where it would get multiplied becomes infectious “I _h” and moves to the infectious compartment at the rate of “α _h.” On the other hand, if an exposed person gets timely treatment, he would be able to return to the susceptible compartment after getting recovered at the rate of “R _h.” The infectious individuals after getting recovered will then return to the susceptible compartment but they would not get immunity. Some of the infectious ones can die of TBE at the rate of “e _h.” However, at any stage natural death may occur at the rate of “d _h.” (Here natural death means any kind of death other than death due to TBE.) Ixodes susceptible tick “S _t” gets infected when he bites an infectious human “I _h” at the rate of infection “β ₂.” At the rate of “α _t,” tick moves to the infectious compartment and then dies out. The recovery of ticks is not possible albeit. The flow chart given in Figure 1 has been created to get the hang of our model, and arrows are used to show the interactions. Based on the aforementioned assumptions we have the following systems of equations for human and tick compartments:

(1) d S h d t = μ h N h − β 1 S h I t N h − v S h − d h S h + R h E h + γ h I h , d E h d t = β 1 S h I t N h − R h E h − d h E h − α h E h , d I h d t = α h E h − γ h I h − d h I h − e h I h , d S t d t = μ t N t − β 2 S t I h N t − d t S t , d I t d t = β 2 S t I h N t − d t I t .

Figure 1

Description of the basic model of TBE.

The system (1) of differential equations is named the TBE system of ordinary differential equations (ODEs), provided that total human population N _h(t) and N _t(t) satisfy the following condition:

S h ( t ) + E h ( t ) + I h ( t ) = N h ( t ) and S t ( t ) + I t ( t ) = N t ( t ) ,

where the parameters used in system of differential equations (1) and (5) are given in Table 1.

Table 1

Parameters

S h	Susceptible humans
E h	Exposed humans
I h	Infectious human
S t	Susceptible ticks
I t	Infectious ticks
N h	Total human population
N t	Total tick population
μ h	Human birth rate
β 1	Rate of human infection
β 2	Rate of tick infection
V	Human vaccination rate
d h	Human natural death rate
e h	Human infected death rate
R h	Rate at which exposed humans get recovered and become susceptible again
γ h	Rate at which human after getting recovered become susceptible again
α h	Rate at which exposed humans become infectious
d t	Ticks natural death rate

A threshold quantity which determines the possibility of a disease to occur or dies out is known as R 0 -like threshold. R 0 is the number of secondary infections caused by a single infective introduced to a susceptible population [19]. In epidemiological models, if the value of R 0 > 1 , then the disease will outbreak and the value of R 0 < 1 implies that disease will die out. The value of R 0 -like threshold for different diseases is different. For instance, R 0 for HIV disease is 3 and R 0 for SARS is 5. Both the values were calculated by using the same method, otherwise we would not have been able to conclude that SARS is worse than HIV. From the docket of few important methods for calculating R ₀-like threshold such as Anderson and May method, Jacobian method and next-generation matrix method, we have used the next-generation matrix method [12].

Next-generation matrix is a method involving transmission and transition matrices, and spectral radius of the negative product between transmission matrix and inverse of transition matrix gives the value of R ₀. For further elaboration of transmission matrix and transition matrix see ref. [12].

3 Methodology and results

In this section, the reproductive number R ₀ for TBEV is calculated. Next-generation matrix is the widely used method to find R ₀-like threshold [41]. It precisely deals with the transmission and transition rates inside the infected subsystem. For the stability analysis of mathematical model of TBE, next-generation matrix is used. In the method, two kinds of matrices are used to calculate the value of R ₀-like threshold. One is the transmission matrix and the other is the transition matrix. The transmission matrix, T, denotes the pathogen passing between all stages of the infection subsystem. It is calculated by taking Jacobian of all the transmission rates:

The transition matrix, Σ , denotes all other transitions to and from the infection subsystem. It is calculated by taking Jacobean of all the transition rates.

where − T Σ − 1 is given by the following matrix:

R 0 is the spectral radius of − T Σ − 1 (matrices derived above). Hence, we get the following:

The formula computed in the equation is the required formula to calculate the basic reproductive value of TBEV. From the derived formula of R ₀, it can be seen that:

1. R 0 2 is directly proportional to β 1 β 2 , which shows that with the increase in the rate of transmission between humans and ticks the value of R 0 increased with squared escalation.

2. R 0 2 is inversely proportional to the product of two factors ( R h + d h + α h ) and ( γ h + d h + e h ) which is of exposed population and of infected population, respectively, reflecting the fact that exposed are becoming either susceptible or infective and infective are becoming either susceptible or going to die. So, there is de-escalation of population in both compartments, i.e., exposed and infective. The higher the speed of population jumping into other compartment, the lower the value of R 0 2 .

3. R 0 2 is inversely proportional to the death rate of tick, which shows that the increase in death rate of ticks will result in decrease in R 0 2 .

4. Due to the unavailability of precise data of the disease, the accurate value of R 0 2 could not be calculated; one of the major reasons is that the disease is so rare and data are not provided in the literature so far.

4 Equilibrium points and their local stability

To discuss the equilibrium points of the mathematical models, first we do some modifications of system (1) by considering the assumption of fractions as: x 1 = S h N h , x 2 = I h N h , x 3 = E h N h , x 4 = I t N t , x 5 = S h N t , which further simplify system (1) as:

For the ease of handling calculations, the following notation is assumed for the modified system (3): μ h = α 1 , υ = α 2 , d h = α 3 , R h = α 4 , I h = α 5 , α h = α 6 , r h = α 7 , e h = α 8 , α = α 9 , d t = α 10 , μ t = α 11 , h t = α 12 and the following system is obtained:

The biologically meaningful equilibria are disease-free equilibrium (DFE) and endemic equilibrium, depending on S h , I h , E h and S t and I t . To find DFE, all time derivatives of x 1 , x 2 , x 3 , x 4 and x 5 are set equal to zero, whose biological meaning is “without infection;” then, we obtain:

The Jacobian technique at DFE is used for system (4). By adopting this technique, the local stability of DFE is found.

The characteristic equation implies that the first eigenvalue λ = − ( μ h + v ) is negative. To further evaluate, the submatrix of matrix J is taken as:

where β 1 = − α 3 − α 4 − α 6 , β 2 = α , β 3 = − α 7 − α 3 + α 8 , β 4 = α 10 , β 5 = α 9 and β 6 = − α 9 − α 10 . Characteristic equation corresponding to J 1 is as follows: λ 4 + 3 β 4 + β 6 λ 3 + 3 β 1 β 3 + β 4 β 6 λ 2 + ( β 1 β 3 β 4 + β 3 β 1 β 6 ) λ + 3 β 1 β 3 β 4 β 6 = 0 . For the sake of ease, equation (5) can be rewritten as:

where

To prove the local asymptotic stability through the Routh–Hurwitz criterion [34,36], we prove that all coefficients in equation (5) are positive and a 1 a 2 a 3 > a 4 are greater than zero for R 0 < 1 . It is quite clear that a i ' s are positive; hence, the first condition is fulfilled. For the second condition that a 1 a 2 a 3 > a 4 are greater than zero for R 0 < 1 , of the Routh–Hurwitz criterion, we proceed as follows:

It is effortless to prove the second condition. Hence, both the conditions are fulfilled for R ₀ less than one and the infection-free equilibrium point is LAS. Hence, it proves the theorem.

where A 1 = ( a 2 2 − 3 a 1 a 3 + 12 a 4 ) and A 2 = 2 a 2 3 − 9 a 1 a 2 a 3 + 27 a 3 2 + 27 a 1 2 a 4 − 72 a 2 a 4 .

By determining the sign of λ , λ 1 , λ 2 , λ 3 and λ 4 conclusion about stability can be made. If all eigenvalues are negative or purely imaginary or real part of complex eigenvalues are negative, then the system is asymptotically stable otherwise unstable for all possibilities of sign of eigenvalues.

5 Numerical computation and graphical analysis

To get the deeper insight of the role of R ₀-like threshold, we solved the coupled system of ODEs given by (1). Runge-Kutta (RK) method of order 4 with initial conditions for the tick compartment is as follows:

and initial conditions for the human compartment are as follows:

The approximated data collected from Austria’s 26 years TBE vaccination program are used [31]. RK method of order 4 is implemented for numerical solution of the coupled system of ODEs in MATLAB environment. MATLAB 2016 with double precision floating point operations is used to carry out all computations. TBE is solved. In order to avoid any error encountering, the values of the parameters are approximated and are not exact as given in the program. In order to keep the content precise and lucid, a population of anonymous state suffering from TBE has been taken into consideration. Figure 2 shows that the susceptibility to disease in humans decreases with time. Figure 2 shows that it is highest at initial time. But as the time passes by the human would leave the compartment either becoming exposed or immune.

Figure 2

Graph of susceptible vs infected human.

Figure 3 represents the human population of the infected, and as per the aforementioned assumptions prior modeling, initially infected humans are zero but as the infection travels from exposed to infected, population of infected goes on increasing.

Figure 3

Graph of susceptible & infected vs exposed human.

Figure 4 shows that the exposed humans are the same as the infected but the difference between both is that in the exposed compartment the growth of population is a bit faster than that of infected, where the growth of population would be steady.

Figure 4

Graph of infected & exposed vs susceptible human.

It is evident from Figures 2–5 that the solution to the coupled system of ODEs exists and is inline with numerical results obtained. Clearly, the susceptibility decreases as the infection increases to maximum. Likewise, the exposed ones are those who lie between susceptible and infected humans. The increase and decrease are witnessed from the graphs that our model confirms the growth and decay of population in different compartments.

Figure 5

Graph of susceptible & exposed vs infected human.

6 Discussion

Pakistan is an agrarian state and history of few decades bears witness to the fact that Pakistan has been suffering from vector-borne disease more than other viral diseases since the twentieth century [1,37]. Mosquitoes have been the primary vectors serving for devastating vector-borne diseases such as dengue, malaria and Japanese encephalitis. TBE has been a very rare disease in Pakistan premises [39]. However, TBEV can follow the pattern of Crimean–Congo hemorrhage fever (CCHF) spread in Pakistan. The outbreak of this disease occurred due to a couple of reasons given as follows:

1. Unavailability of developed health-care system and insufficient equipment in the institutes to control CCHF [16].

2. Unawareness of the spread of this disease among general public because of the vastness of agricultural sector of Pakistan.

3. The migration of numerous refugees from endemic areas to Pakistan.

This disease started spreading in Pakistan via an infected domestic animal which later started affecting humans and vice versa.

Likewise, Pakistan is under the risk of TBE, which is already very common in its neighborhood countries i.e. China, Russia, Kazakhstan etc. and in European countries. Pakistan is not an isolated state, it is a trading partner of China and many other European countries. The import of animals such as pets, sheep and camels is very common in Pakistan. Risks are high to a considerable extent. Even unpasteurized milk can pave the way for this virus between an infected sheep and susceptible human and Ixodes ticks [18]. Some other factors such as climate change, global warming, increasing interdependence and urbanization can also become the reason for the outbreak of this disease. Pakistan has been found in the docket of top ten countries vulnerable to climate change [26], which is really dangerous and point of concern for Pakistan.

7 Future perspective

The future perspective of this work resides as follows:

1. After finding the precise data one can be able to predict accurately whether the disease is going to be endemic in the area premises or not.

2. More work on vector-borne disease and mathematical modeling will be promoted. Biological mathematical modeling will be promoted and enhanced.

3. Before getting extremely endemic, disease will be controlled.