MODELING CO-INFECTION OF IXODES TICK-BORNE PATHOGENS

Ticks, including the Ixodes ricinus and Ixodes scapularis hard tick species, are regarded as the most common arthropod vectors of both human and animal diseases in Europe and the United States capable of transmitting a large number of bacteria, viruses and parasites. Since ticks in larval and nymphal stages share the same host community which can harbor multiple pathogens, they may be co-infected with two or more pathogens, with a subsequent high likelihood of co-transmission to humans or animals. This paper is devoted to the modeling of co-infection of tick-borne pathogens, with special focus on the co-infection of Borrelia burgdorferi (agent of Lyme disease) and Babesia microti (agent of human babesiosis). Considering the effect of coinfection, we illustrate that co-infection with B. burgdorferi increases the likelihood of B. microti transmission, by increasing the basic reproduction number of B. microti below the threshold smaller than one to be possibly above the threshold for persistence. The study confirms a mechanism of the ecological fitness paradox, the establishment of B. microti which has weak fitness (basic reproduction number less than one). Furthermore, co-infection could facilitate range expansion of both pathogens.

1. Introduction.The health burden of tick-borne diseases, such as Lyme disease, human babesiosis and human granulocytic, is increasing partially due to a greater exposure of humans to the infected vectors [12,19] and tick expansion [22].Lyme disease, one of the most common tick-borne diseases worldwide, is caused by the spirochete Borrelia burgdorferi, claiming around 25,000 cases annually during 2008-2011 in the United States.Human babesiosis, caused by the blood protozoan parasite Babesia microti, has also been increasing in the northeastern United States.The cases of human granulocytic anaplasmosis, also known as human granulocytic ehrlichiosis and caused by the bacterium Anaplasma phagocytophilum, has increased significantly over the last decade [15].
All these three diseases, Lyme disease, human babesiosis and human granulocytic ehrlichiosis, are transmitted by Ixodes ticks.The blacklegged tick, Ixodes scapularis, can be infected with a large and diverse array of human pathogens, such as A. phagocytophilum, B. microti, and B. burgdorferi, or all of them simultaneously [13,15].Ticks co-infected with multiple microorganisms have been documented in previous studies [19,28].For example, 92 questing Ixodes ricinus ticks were collected in northern France to test for three micro-organisms, Bartonella sp., Borrelia burgdorferi sensu lato and Babesia sp., which are known as suspected tick-borne pathogens.Seven among 92 samples (7.6%) were positive for at least two of the pathogens and one sample was positive for all three pathogens [13].In another study, 394 questing adult blacklegged ticks, Ixodes scapularis Say (Acari: Ixodidae), collected at four sites were analyzed for five microbial species: Anaplasma phagocytophilum, Babesia microti, Babesia odocoilei, Borrelia burgdorferi, and the rickettsial I. scapularis endosymbiont.In 55% of infected ticks (193/351), a single agent was detected.In 45% (158/351), two or more agents were detected, among which 37% harbored two agents and 8% harbored three agents [26].The reference [15] produces a wonderful graph illustrating the mean level of co-infection prevalence of Anaplasma phagocytophilum, Babesia microti, and Borrelia burgdorferi in questing Ixodes scapularis nymphs in various sites.
Since several of these pathogens can be harbored by the same tick species and carried by the same rodent reservoir hosts, they can be co-transmitted to humans in concert.In theory, when an individual is bitten by an infected tick with multiple pathogens, or sequentially bitten by multiple ticks each transmitting a different pathogen, the individual can get the co-infection.Several studies in various regions have documented the co-infection of multiple tick-borne pathogens [15].For example, it was reported that 3 out of 19 patients with Human Granulocytic Ehrlichia (15.8%) showed evidence of co-infection: 1 (5.3%) with B. burgdorferi, 1 (5.3%) with B. microti, and 1 (5.3%) with both microorganisms.One patient diagnosed with babesiosis was also seropositive for ehrlichiosis [19].Among patients with a confirmed tick-borne infection, the percentage of co-infection rates can reach as high as 39%.The most commonly recognized co-infection in the eastern United States is Lyme disease and human babesiosis, accounting for around 80% of co-infection while Lyme disease and human granulocytic ehrlichiosis co-infection is less common, occurring in 3%-15% of patients in Connecticut and Wisconsin [4].
Co-infection of multiple tick-borne diseases becomes an emerging problem since the interaction of these pathogens may affect the severity and duration of symptoms in humans, making diagnosis and treatment more challenging [12,15,27].For example, simultaneous Lyme disease and human babesiosis infection is shown to be correlated with a more severe clinical progression and persistent symptoms than either condition alone [1,6,15].
Co-infections from tick-borne diseases pose a threat to human health in the northeastern and midwestern United States, but the risk of acquiring a co-infection is not fully understood [15].Field studies in ticks, reservoir hosts, and humans indicate that co-infection with B. burgdorferi and B. microti is common, which could promote transmission and emergence of B. microti in the enzootic cycle.Coinfection by Borrelia burgdorferi, the primary agent of Lyme disease, and Babesia microti, the primary agent of babesiosis, may serve as a diagram for the study of other vector-borne pathogen interactions [6].There is an increasing body of literatures, such as [3,11,21] among others, investigating the co-infection dynamics of various diseases, which greatly improve our understanding of pathogen evolution.To the best of our knowledge, very few theoretical studies, except [24,29] have been performed to address the co-infection of pathogens for vector-borne diseases, especially for tick-borne pathogens.The possible co-infection of Zika virus and several other mosquito-borne viruses of clinical importance, for example, the dengue virus and West Nile virus, and yellow fever virus, poses a public health emergency [23].Identifying ecological drivers of pathogen emergence and host factors that fuel disease severity in co-infected individuals will help to design effective preventive and therapeutic strategies [6].
Although B. microti was identified earlier than B. burgdorferi and the geographic range of both pathogens has expanded in the United States, the spread of babesiosis has lagged behind that of Lyme disease.It is also illustrated that the geographic expansion of B. microti has been restricted to those areas where Lyme disease is already endemic [6].Furthermore, it is generally regarded that the ecological fitness of B. microti is weak and it is not easy to establish in a habitat in the sense that its basic reproduction number is below the threshold for persistence.Therefore, the persistence and geographic expansion pose an ecological paradox [6] in disease transmission: why can the pathogen with reproduction number smaller than one establish?The investigation of this paradox requires the incorporation of pathogen interactions within realistic epidemiological and ecological contexts.Measuring and identifying the effect of co-infection on Babesis emergence should help to predict the spatial and temporal gains of this epidemic and design interventions to reduce disease risk.In this study, we are going to test the hypothesis whether B. burgdorferi can promote the fitness of B. microti.To do so, we develop a mathematical model to assess the effect of co-infection on the disease dynamics of both pathogens.
2. Model formulation and dynamics.Generally, Ixodes scapulars ticks undergo a four-stage life cycle: egg, larva, nymph and adult.When larvae take blood meals from infected reservoirs, the pathogen fuses into and develops in the body.After larvae molting to nymphs and feeding during late spring or early summer in the following year, reservoir hosts may become infected.Humans are accidental hosts and nymph is the primary vector for the transmission of Babesia microti and Borrelia burgdorferi [6].To develop a co-infection model, we start from describing a single pathogen transmission with stage-structured tick population growth.

Stage-structured tick population model.
There have been several stagestructured models proposed, such as those in [9,32], which present good descriptions on tick population growth.Here, to keep the model simple while incorporating the key features of both population growth and disease transmission simultaneously, we derive the following model, which is similar to an autonomous version of the models in [14,17], by ignoring effects of seasonal weather conditions on tick development and activities since the main focus of the current research is on the co-infection dynamics of multiple tick-borne pathogens.Although the analytical approach in [14,17] remains applicable to the model here, for reader's convenience, we perform the model analysis tailored to the autonomous model system.The tick population is stratified into the following stages: eggs (E), questing larvae (LQ), feeding larvae (LF ), questing nymphs (NQ), feeding nymphs (NF ) and adults (A), and the size of each stage satisfies the following equations: It is worth to mention that nonlinear terms D L LF 2 (t) and D N NF 2 (t) in the above system are used to measure the additional deaths of feeding individuals due to the grooming effects of immature tick hosts (normally small animals).All variables and parameters are detailed in Tables 1 and 2, respectively.It is easy to see that the model system (1) fits into the general population growth model in [8], which analyzes the stability and persistence of a general ODE model for population growth.Therefore, the results of Fan et al. [8] can be readily applied here.Specifically, the vector field Using the argument in [8], the vector reproduction number for the tick population can be written as the following form (it can also be computed by using the next generation matrix approach [5,7]) where represents the survival probability of ticks from eggs to the adult stage while b E is the number of eggs produced by one adult per unit time and 1 µ A is the egg reproducing duration.It then follows from [8, Theorem 4.2] that the origin z = 0 is locally asymptotically stable if R T < 1 and unstable if R T > 1.Furthermore, Theorems 4.5 and 4.6 of [8] imply that 0 is globally asymptotically stable when R T < 1.Moreover, the system is dissipative and uniformly persistent according to Theorems 7.3 and 7.4 in [8] when R T > 1.Furthermore, there exists a positive equilibrium in this scenario [8,Theorem 6.2].In addition, the vector field g possesses the following properties: (1) g is cooperative on R 6 + and Dg(z) = ( ∂gi ∂zj ) 1≤i,j≤6 is irreducible for every z ∈ R 6 + ; (2) g(0) = 0 and for all i ∈ {1, 2, 3, 4, 5, 6}, g i (z) ≥ 0 for all z ∈ R 6 + with z i = 0; (3) g is strictly sublinear on R 6 + , i.e., g(αz) > αg(z) for any α ∈ (0, 1) and any z ∈ Int(R 6 + ).Therefore, the positive equilibrium is globally asymptotically stable for system (1) by Corallary 3.2 of [34].Summarizing the above arguments, we claim that the vector reproduction number R T can completely determine the global dynamics of the tick growth.
Theorem 2.1.The following statements are valid: and it is globally asymptotically stable with respect to all nontrivial solutions.
The above theorem indicates that the tick population will die out if R T ≤ 1 while it will eventually stabilize at a positive equilibrium 2.2.Single pathogen infection dynamics.Now, we move to the pathogen transmission dynamics by starting from the appearance of only one type of pathogen.We assume that the host for adult ticks, mainly deer, is an incompetent reservoir for pathogen transmission, that is, the adult host can neither harbor nor transmit the pathogen to ticks.Furthermore, no vertical transmission from adult ticks is observed for B. Burgdoferi and B. microti.Therefore, the pathogen transmission cycle is maintained between immature ticks and the host species.Here we assume that there is only one host species-mice for immature ticks.

The dynamics of Borrelia infection only.
When there is only Borrelia in the environment, the pathogen transmission can be described as where In this system, we use subscript 0 to represent susceptible class while 1 to denote Borrelia infected class.Then the total population size for mice follows the logistic growth which admits a globally asymptotically stable positive equilibrium To be population persistent, we assume b M − µ M > 0.
A standard argument can show that every solution of (3) with a nonnegative initial value remains nonnegative (see for example [25, Theorem 5.2.1]).Moreover, the tick population sizes of eggs (E), questing larvae (LQ), feeding larvae (LF =LF 0 + LF 1 ), questing nymphs (NQ=NQ 0 + NQ 1 ), feeding nymphs (NF =NF 0 + NF 1 ) and adults (A=A 0 + A 1 ) satisfy the tick population growth system (1).When R T ≤ 1, we have and similarly, we have The equation for M 1 gives the following asymptotically autonomous system By the theory of asymptotically autonomous semiflows (see [30,Corollary 4.3]), we have lim Then we have the following asymptotically autonomous system to describe the Borrelia transmission dynamics The basic reproduction number of the Borrelia can be derived through the next generation matrix method proposed by [5,7] as Using a similar argument as in Theorem 2.1, we have the following result indicating that R 1 determines the disease dynamics.Lemma 2.2.Assume that R T > 1, then the following statements are valid: then Borrelia-free solution is globally asymptotically stable for system (4) in R 3 + ; (ii) If R 1 > 1, then system (4) admits a unique positive endemic equilibrium and it is globally asymptotically stable for all nontrivial solutions.
A rigorous proof similar to that in [18, Theorem 2.2] via the theory of internally chain transitive sets [16,33] shows that the dynamics of the model system (3) are determined by the two threshold parameters: the reproduction number of ticks (R T ) and the basic reproduction number of Borrelia (R 1 ).The theory of internally chain transitive sets has been widely applied to lift the dynamics of the limiting system to the whole model system, and readers may refer to [10,14,18] for some wonderful applications.

2.2.2.
The dynamics of Babesia infection only.In system (3), replacing β 11 with β 22 which is the Babesiosis transmission probability of ticks from infected mice to susceptible questing larvae and replacing γ 11 with γ 22 which is the Babesia transmission probability of mice from infected questing nymphs to susceptible mice, we obtain a model system to describe the Babesia transmission between ticks and reservoirs.Then the basic reproduction number for Babesia R 2 can be figured out by employing the same method, which becomes The similar result as Theorem 2.3 also holds for Babesia transmission when other tick-borne pathogens do not appear.

2.3.
Co-infection dynamics.We consider co-infection of B. microti and B. burgdorferi between reservoir mice and vector ticks.The co-infection disease transmission process between ticks and their hosts can be described by the diagram in Figures 1 and 2. Each post-egg age group is divided into four subgroups based on their infection statuses with various infectious pathogens: susceptibles to both B. microti and B. burgdorferi (subscript 0), infected with B. burgdorferi only (subscript 1), infected with B. microti only (subscript 2), and co-infected with B. microti and B. burgdorferi (subscript 3).The same subscripts are used for mice which are stratified into four compartments: M 0 , M 1 , M 2 , and M 3 .In theory, hosts can get infected with multiple pathogens either through a single bite by a tick harboring both pathogens, or sequential bites of ticks with each transmitting a different pathogen [15].All variable are summarized in Table 1, with 18 variables for various tick stages and disease statuses, and 4 variables for mice with different disease statuses.Therefore, we can write down a system (6) of 22 equations, which is postponed to the appendix.
Theorem 2.4.For any initial value x 0 ∈ R 22 + , system (6) has a unique nonnegative and bounded solution x(t, x 0 ) for all t ≥ 0.
When R T ≤ 1, an argument similar to that used in the previous subsection shows that all variables related to ticks go to zero and the disease dies out.When R T > 1, we will figure out the co-infection basic reproduction number R 0 .According to Theorem 2.1, we have  Linearizing the system (6) at the disease-free equilibrium, we can obtain the following system for the infected compartments: A schematic diagram of co-infection in the tick population.Here E (eggs), LQ (questing larvae), LF (feeding larvae), NQ (questing nymphs), NF (feeding nymphs) and A (adults) represent the stages of tick population with subscripts denoting the infectious status for each pathogen.Subscript 0: no pathogen in ticks; 1: Borrelia only; 2: Babesia only; 3: both pathogens.
A schematic diagram of co-infection in mice M with subscripts denoting the infectious status for each pathogen.
Note that in this system, infected feeding nymphs and adults are not included as they play no role on pathogen transmission between immature ticks and reservoirs.
From the above linearized system, we can figure out the basic reproduction number using the ideas of [5,7].In particular, we set the disease transmission matrix and the transition matrix Then, the next next generation matrix can be written as: Based on [5,7], the spectral radius of F V −1 is the basic reproduction number, which can be computed as where, similar to R 1 and R 2 , We can see that γ 31 , γ 32 , γ23 , γ33 , γ13 , γ23 do not appear in the expression of the basic reproduction number.These parameters measure the probabilities of further infection of infected mice due to sequential tick bites.Hence sequential infection of another pathogen for infected reservoirs does not increase the basic reproduction number.
3. Numerical simulations and discussion.To perform the numerical simulations, it is pivotal to choose biologically acceptable parameter values.The birth and death rates for mice are taken from [2].The density-dependent death rate for mice is chosen such that the mice population size is at the steady state 200.In reality, the tick development and feeding rates are strongly affected by biotic and abiotic factors, in particular, seasonal temperature variations [32].Since the main focus of the current study is to investigate the co-infection dynamics of tick-borne pathogens, we formulate the model in a simpler setting that all related parameters are temperature-independent, and hence time-independent by taking the annual average.For this purpose, the averaged tick development and feeding rates are adopted from [17].Previous studies measure the transmission probabilities for Borrelia and we take reasonable values as in Table 2.The transmission probabilities related to the Babesia transmission are set such that the basic reproduction number for Babesia (R 2 ) is less than one.Since there are few studies on the effect of coinfection on promoting the transmission probabilities, we assume the co-infection can enlarge the transmission probabilities for any pathogen 1.5 times.The model parameters and their values are summarized in Table 2.In all the simulation plots, we focus on nymphs as this stage is responsible for the majority of human infections with tick-borne diseases.Based on the parameter values in Table 2, the vector reproduction number for ticks is R T = 7.2892, and therefore the tick population size will stabilize at a constant level (Figure 3(a)).When co-infection is not considered, the reproduction numbers of Borrelia and Babesia are R 1 = 1.2148 and R 2 = 0.9111, respectively.According to Theorem 2.1, single Babesia transmission can not be established as R 2 < 1 while single Borrelia transmission can be maintained since R 1 > 1.The single pathogen transmission dynamics are simulated in Figures 3(b) and 3(c).
However, when co-infection is considered, both Borrelia and Babesia can establish in the habitat.In this case, the basic reproduction number is R 0 = R 1 = 1.2148.Some ticks are infected with Borrelia only while some infected with Babesia only (Figures 3(d) and (e)).Some other ticks may get infected with both pathogens (Figure 3(f)).The co-infection can increase the infected tick size, in particular, the size of infected ticks with Borrelia is also promoted.Since R 0 ≥ R 1 ≥ R 2 , co-infection promotes the transmission of human babesiosis, which provides one mechanism to resolve the ecological paradox of babesisa transmission: it has weak fitness (R 2 < 1) while gets established in the region where Borrelia remains endemic.
Ticks transmit more pathogens than other arthropods.Previous studies have demonstrated that co-infection occurred in almost half of the infected ticks, and that ticks could be infected with up to five pathogens [20].The current work studies the co-infection dynamics of tick-borne pathogens, which highlights the co-infection phenomenon in ticks.In particular, it illustrates that co-infection can facilitate the transmission of both pathogens and promote the risk of each disease, which emphasizes the need for new diagnostic tests better adapted to tick-borne diseases and novel control measures for co-infection of tick-borne pathogens.
Pathogens interacting within a single host could in theory facilitate, compete or have no effect on others [11,15].In our numerical simulations, we only assume that co-infection of Borrelia and Babesia can boost the transmission probabilities of both pathogens.However, in a long-term field study, evidence for both positive and negative interactions between B. microti and A. phagocytophilum were reported, with the outcome dependent on the duration of A. phagocytophilum infection [15].It would be interesting to study the negative interaction of co-infection in hosts in the future.Further investigations can also be performed to refine the obtained results here with the consideration of following aspects: (a) the host community of immature ticks containing more than one host species.Since rates of transmission from infected ticks to vertebrate hosts can vary and increased co-infection was observed in larval ticks that fed on small mammal hosts, but not on sciurid, or avian hosts [15,17]; (b) spatial spread of ticks due to host movement, such as deer, mice and birds [14,22,31]; (c) Seasonality on tick population growth which is driven by the seasonal temperature variations [17].Moreover, since diverse results are reported about the prevalence of co-infection and infection by a single pathogen in different regions [15], sites with different epidemiological and ecological features should be tested to make the model more realistic and conclusions more reliable.From the mathematical point of view, the interaction of multiple pathogens in a simple model can induce very complicated dynamics, such as the competition exclusion, backward bifurcation, and so on [11] while here, we haven't performed these analyses for such a complicated model.All these topics remain as our future work.
Appendix: The model system for co-infection dynamics.

Figure 3 .
Figure 3. Solution simulations with the model parameters in Table2.Solutions through different initial values converge to the constant level for ticks (a), constant infected ticks for Borrelia infection only (b) and Babesia transmission cycle can not establish without the co-infection (c).However, on the scenario of coinfection, both pathogens can get established ((d), (e) and (f)).More interestingly, some ticks becomes infected with only Babesia or Borrelia while some others get infected with both pathogens.

Table 1 .
The state variables for the co-infection model.Bo and Ba represent Borrelia and Babesia, respectively.
3 number of adults co-infected with Ba and Bo M 0 number of mice susceptible to both Ba and Bo M 1 number of mice infected with Bo only M 2 number of mice infected with Ba only M 3 number of mice co-infected with Ba and Bo

Table 2 .
Definitions and corresponding values of the model parameters with the daily timescale.Abbreviations: Bo: Borrelia; Ba: Babesia; TP: transmission probability; AS: assumed parameter values.