Dynamics of Two Pathogens in a Single Tick Population Dynamics of Two Pathogens in a Single Tick Population

A mathematical model for a two-pathogen, one-tick, one-host sys-temispresentedandexplored.Thegoalofthismodelistodetermine howlonganinvadingpathogenpersistswithinatickpopulationinwhicharesidentpathogenisalreadyestablished.Thenumericalsim-ulationsofthemodeldemonstratetheparameterrangesthatallowforcoexistenceofthetwopathogens.Sensitivityanalysishighlights theimportanceofvector-borne,tick-to-host,transmissionratesontheinvasionreproductivenumberandpersistenceofthepathogens overtime.ThemodelisthenappliedtoacasestudybasedonareclaimedswamplandfieldsiteinsoutheasternVirginiausingfield andlaboratorydata.Theresultspinpointthethresholdsrequiredforpersistenceofbothpathogensinthelocaltickpopulation.However,theinvadingpathogenisnotpredictedtopersistbeyondthreeyears.Understandingthepersistenceandcoexistenceoftick-borne pathogenswillallowpublichealthofficialsincreasedinsightintotick-bornediseasedynamics.


Introduction
Tick-borne pathogens are an increasing threat to human health worldwide. In the United States, human diseases caused by tick-transmitted Rickettsia species have been on the rise (Dahlgren, Paddock, Springer, Eisen, & Behravesh, 2016). Rickettsia parkeri, the causative agent of Rickettsia parkeri Rickettsiosis (Paddock et al., 2008(Paddock et al., , 2004, is increasingly reported within Amblyomma maculatum, Gulf Coast tick, populations throughout the southeastern United States, with recently established populations reported from Virginia (Fornadel et al., 2011;Wright et al., 2011), Maryland and Delaware (Florin, Jiang, Robbins, & Richards, 2013). Although the primary vector of R. parkeri is A. maculatum, another potential vector is Amblyomma americanum, the lone star tick (Goddard, 2003), which is widespread within the southeastern U.S. and commonly reported parasitising humans (Stromdahl & Hickling, 2012). Rickettsia parkeri has been detected at low levels ( < 5%) within A. americanum populations (Cohen et al., 2009;Gaines et al., 2014), suggesting the potential for pathogen spillover from A. maculatum populations. In addition to seemingly low levels of R. parkeri spillover into A. americanum populations, A. americanum across the southeast are known to harbour high rates ( > 55%) of another Rickettsia species, Rickettsia amblyommatis, which is generally believed to be nonpathogenic (Nadolny, Wright, Sonenshine, Hynes, & Gaff, 2014). Recent surveys of A. americanum populations in Virginia have indicated R. parkeri spillover at multiple locations, including locations where A. maculatum is thought to be absent (Gaines et al., 2014), suggesting that A. americanum may be able to maintain R. parkeri independently.
Amblyomma americanum ticks have punctuated life histories with four life stages: egg, larva, nymph, and adult. Each non-egg life stage will take a blood meal from a host and have a chance to acquire a pathogen. If a pathogen is acquired, it has the ability to stay with the tick between life stages, transstadially. Unlike other tick species that utilize small hosts for the first bloodmeal, medium hosts for the second, and large hosts for the third, A. americanum primarily feed on large mammals such as Odocoileus virginianus, whitetailed deer, in all life stages within the summer months (Childs & Paddock, 2003). Because of this our model has been simplified to represent all life stages as one population of ticks.
A number of mathematical models have been created to examine the dynamics of ticks and tick-borne diseases building from the standard SIR framework of Anderson and May (Anderson, 1981). The model we use is developed from differential equationbased model for A. americanum and Ehrlichia chaffeensis (Gaff & Gross, 2007;Gaff, Gross, & Schaefer, 2009;Gaff & Schaefer, 2010). In addition, many mathematical models have been created to examine the dynamics of competing pathogens specifically with interest in establishing parameters that would allow for invasions of new pathogens and coexistence of two different pathogens. The majority of these studies have focused on a single host with multiple strains of a single pathogen (Adams & Sasaki, 2007;Feng & Velasco-Hernandez, 1997;Kirupaharan & Allen, 2004;Vasco, Wearing, & Rohani, 2007).
There are a few studies that have looked at multiple pathogens within a vector-borne disease model (Ackleh & Allen, 2003), but this will be one of only a few models to look at potentially competing pathogens for a tick-borne disease (Dunn et al., 2014).
The purpose of this effort is to investigate the dynamics between a single tick species and two related tick-transmitted pathogens in order to identify the parameters necessary for pathogen maintenance and coexistence. Although this effort focuses broadly on a onetick, two-pathogen system, it is motivated by the case study presented, which is based on field and laboratory data from southeastern Virginia (Nadolny et al., 2014;. This paper is organized as follows. In Section 2, we introduce the two pathogens model. In Section 3, we calculate the basic reproductive number and invasion reproductive number. In Section 4, we consider the sensitivity of the model simulations as parameters vary within a biologically feasible parameter space. Then in Section 5, we apply the model to a specific case study based upon a field site in southeastern Virginia we have been surveying since 2010 to explore the role of site variation in pathogen competition. Finally, we present overall conclusions and next steps in Section 6.

Model
In Equations (1)-(8), we present a deterministic model that allows us to quantify the conditions under which one or two pathogens may experience long-term survival within single host and single vector populations. The pathogenicity of R. amblyommatis has not been established in humans, but will be referred to as a pathogen in this model as it infects both the ticks and the hosts. Although the formulation and discussion is intentionally general, the model is constructed to consider the transmission cycle of R. parkeri and R. amblyommatis in A. americanum, with Odocoileus virginianus, white-tailed deer, as the host. Because A. americanum are more frequently associated with a single type of host and are found questing at all life stages during the same time of year, we used a single host population in modelling its population dynamics. This model is an extension of the model in Gaff and Gross (2007) and includes deterministic differential equations that track host and vector population dynamics while also tracking the pathogen infections within each population. Note that variable and parameter definitions and baseline values are summarized in Table 1. Baseline values were chosen based on previous models (Gaff & Gross, 2007;Gaff & Schaefer, 2010;Gaff et al., 2009) and published data (Nadolny et al., 2014;.

Population dynamics
Host population dynamics are taken from Gaff and Gross (2007), with the total number of hosts is denoted by N, and the change in the host population is based on a logistic growth model with an established carrying capacity K and growth rate β. This logistic growth model incorporates density-dependent mortality, and additional density-independent mortality b from hunting or other external factors is considered as well.
All Hosts: The vector population V is based on a logistic growth model, with birth rateβ and densitydependent death, but instead of employing a standard population carrying capacity, as seen in the host population dynamics, tick populations are limited by the maximum amount of ticks per host M times the number of hosts N. Background density-independent mortalitŷ b is included within the population as well. All Vectors:

Transmission
The model allows three forms of pathogen transmission. With vector-borne transmission, the pathogen is spread from infected host to susceptible tick (Figure 1(a)) or from infected tick to susceptible host during bloodmeals (Figure 1(b)). Transovarial and transstadial transmission refers to the transfer of the pathogen from an adult female tick to her offspring and the maintenance of a pathogen as a tick moults from one stage to the next, respectively (Figure 1(c)). Based on related laboratory experiments, our model assumes that ® Host background density-independent mortality rate 0.0 b Tick background density-independent mortality rate 0.001 Notes: Baseline values are estimated based on unpublished field and lab research except where noted by † (Gaff & Gross, 2007). All rates are per month except as noted by *.
transovarial transmission will not result in coinfected offspring . Lastly cofeeding transmission can be described as the transmission of a pathogen or pathogens from one tick to another by feeding upon the same host at the same time in close proximity (Figure 1(d)). For cofeeding, the pathogen is transmitted through the host tissues from an infected tick to a susceptible tick feeding simultaneously next to each other on the same host without infecting the host. In the following, we denote the hosts infected by pathogen 1, 2, or both by Y 1 , Y 2 , and Y 12 , respectively. Likewise, we denote the vectors infected by pathogen 1, 2, or both by X 1 , X 2 , and X 12 , respectively. The equations describing host infections are given in Equations (3)-(5). Tick-to-host vector-borne transmission for pathogen i, i = 1, 2, is modelled by The rate of transmission from tick-to-host is A i , which includes biting rate, probability of transmission, and proportion of hosts to ticks. This rate is multiplied by the proportion of the susceptible host population to the total population and by the number of vectors that carry pathogen i, which includes coinfected vectors. Similar dynamics are considered for hosts who become coinfected with both pathogens, with the exception that the proportion of susceptible hosts who can move to the coinfected class is simply Y 1 /N for transmission of pathogen 2 and Y 2 /N for transmission of pathogen 1. There is no evidence for cross immunity of these pathogens and thus susceptibility does not change after contracting a single pathogen. Host recovery, ν i , from either pathogen is included within the model, allowing recovery from a single pathogen to the susceptible class. In the case of coinfection, ν 12 is the rate of host recovery from coinfection into the pathogen 1 infected class, and similarly ν 21 represents the rate of recovery from coinfection into the pathogen 2 infected class.
Hosts infected with only pathogen 1: Hosts infected with only pathogen 2: Coinfected hosts: (5) While there is only one mode of pathogen transmission for hosts, all three modes of pathogen transmission are possible with tick populations. The differential equations modelling the infected tick populations are given in Equations (6), (7), and (8) and explained as follows. Similar to transmission to the host population, ticks also acquire pathogens through vector-borne host-to-tick transmission. Host-to-tick vector-borne transmission of pathogen i, i = 1, 2, is modelled bŷ where the host-to-tick vector-borne transmission rate of pathogen i isÂ i , and is multiplied by the number in the vector population that are susceptible and by the proportion of hosts infected with pathogen i. Similarly, the termsÂ 12 ((Y 2 + Y 12 )/N)X 1 and A 21 ((Y 1 + Y 12 )/N)X 2 describe the host-to-tick transmission from either pathogen class into the coinfected class.
Additionally, ticks acquire pathogens through transovarial or transstadial transmission. Transovarial and transstadial transmission is described for pathogen i, i = 1, 2 bŷ where the proportion of infected females' (X i or X ij ) offspring with transovarial transmission is γ i (γ ij for coinfected females) and is regarded the same as ticks maintaining the pathogen between life stages. As explained at the start of the transmission subsection, we choose not to include transovarial cotransmission in this model, and observe γ 12 + γ 21 ≤ 1 to avoid double transmission by coinfected ticks.
Lastly, we consider cofeeding transmission of the pathogens between ticks. Cofeeding transmission of pathogen i, i = 1, 2, is modelled by where the cofeeding transmission rate of pathogen i, μ i , is multiplied by the size of the susceptible population of ticks and the proportion of ticks infected with the pathogen i. Cofeeding transmission resulting in coinfection is likewise modelled by (μ 12 + μ 21 )(X 1 X 2 /V) where the cofeeding coinfection transmission rate for a tick with pathogen 1 then becoming coinfected is μ 12 , and a corresponding definition is made for μ 21 . Ticks infected with pathogen 1: ® 56 A. WHITE ET AL.

Model analysis
The infection dynamics of pathogens 1 and 2 can be described by well-defined epidemiological threshold quantities, the basic reproductive numbers R 1 , R 2 , and the invasion reproductive numbersR 1 andR 2 . The basic reproductive number (BRN) for each pathogen describes the ability of a single pathogen to persist in a completely susceptible population, i.e. the second pathogen is not present. The invasion reproductive number (IRN) describes the ability for a second pathogen to invade when the first pathogen is already present. A pathogen is able to invade and persist in a susceptible population if the basic reproductive number exceeds 1. If both pathogens' basic reproductive numbers exceed 1 then only the pathogen(s) with an IRN exceeding 1 will be able to invade a population with an already present pathogen.
The pathogens within our model are treated as having identical mechanistic structure, although their parameter values and initial conditions may differ. The basic reproductive numbers are derived using next-generation operator approaches (Diekmann, Heesterbeek, & Metz, 1990;Van den Driessche & Watmough, 2002). The basic reproductive numbers for each pathogen depend on direct vector-to-vector transmission as well as on host-vector transmission and i = 1, 2 provides the corresponding terms for each pathogen. Then ® The expressions for R i have the same form seen in some other studies of vector-borne infections with vertical transmission (Kribs-Zaleta & Mubayi, 2012;Pelosse et al., 2013) and can be seen to have a less-than-additive effect on the two components, with max(R Vi , R Ni ) < R i < R Vi + R Ni . Finally, We now consider the invasion reproductive numbersR 1 andR 2 for pathogens 1 and 2, respectively. Again, the IRN gives us a measure of the ability of one pathogen to invade when the other pathogen is already present. The IRN values are also calculated via nextgeneration operator approaches, but now only the invading pathogen is considered to be an infection, and we require the endemic equilibrium for the resident pathogen rather than the disease-free equilibrium. The IRN will depend on this equilibrium, and we present the values following. For convenience, we represent the following notation for susceptible hosts and vectors, respectively: We note that at equilibrium the total numbers of hosts and vectors are given by whereM is given in Equation (9). We then find that the equilibrium values for pathogen i-infected hosts and vectors (in the absence of pathogen j) are given by, respectively (i = 1, 2): and where Some algebra using standard methods shows that there is a unique equilibrium in which only pathogen i is present, if and only if R i > 1 (i = 1, 2). The next-generation matrix forR i (i = 1, 2) yields a cubic characteristic equation whose largest root is the IRN. Although the expression for this root is too complicated to interpret term by term, it is straightforward to find numerically, and a special case offers some insight into how each IRN compares to the corresponding BRN. Under the simplifying assumptions that coinfected vectors are equally likely to transmit the pathogen as singly infected vectors, and coinfected hosts are equally likely to recover from the pathogen as singly infected hosts, i.e. γ ij = γ i and ν ji = ν i (i = 1, 2, j = i), we have that where, as before, the components express the efficiency of vector-to-vector transmissioñ again using the notation (forR i ) Note that since the expressions in parentheses inR Vi andR Ni are weighted averages of transmission rates,R Vi < R Vi if (and only if) μ ji < μ i , andR Ni < R Ni if A ji < A i andÂ ji <Â i . ThereforeR i < R i if coinfection is (uniformly) disadvantageous relative to primary infection, andR i > R i if it is instead enhanced. In the more general case without the simplifying assumptions, the IRNs also contain terms which are weighted averages of the γ 's and of the ν's. If we do not make the simplifying assumptions above (note that we do not meet these criteria in our parameter assumptions in Table 1), then the IRN can be calculated numerically. The IRNR 2 is the largest root of the following equation: where (under the realistic assumption that The equation forR 1 simply interchanges 1 and 2 in the subscripts. We calculated the BRN and IRN numerically for multiple parameter sets of interest to our field work, and in these cases we found that the values were equivalent. If we assume as a baseline that pathogens 1 and 2 have equal transmission parameters, as in Table 1, then following the derivations above we find R i = 1.2360 which tells us that each pathogen would survive individually. Additionally, in this case we see thatR i = 1.23596 which means that each pathogen is able to invade and coexist in the presence of the other pathogen. Although slight, the difference in the IRN and BRN can be accounted for biologically by coinfected tick reproduction only producing offspring that have a single infection, in proportions γ 12 and γ 21 , thereby reducing marginally each pathogen's replication in the presence of the other. On the other hand, if we assume that pathogen 2 has a host-to-tick transmission rate that is decreased four-fold -soÂ 2 =Â 1 /4 -then the BRN remains unchanged for pathogen 1, but both BRN and IRN decrease to 0.7516 for pathogen 2 (making pathogen 1's IRN undefined altogether, since pathogen 2 cannot establish itself).

Baseline simulations and results
In this study, we work to quantify the relative ability of the transmission pathways to allow long-term persistence of either or both pathogens, as well as the sensitivity of the IRN to changes within these pathways. First, we considered the effect on the model using baseline parameter values except for varying pathogen 1 values for tick-to-host transmission (A 1 ), host-to-tick transmission (Â 1 ), and cofeeding transmission rate (μ 1 ). Coinfection transmission rates (A 12 ,Â 12 and μ 12 ) were simulated at half the respective value of a single-infection transmission rate. See Table 1 for the baseline values, which are based on the work by Gaff and Gross as well as field and lab derived data (Gaff & Gross, 2007;Nadolny et al., 2014;. To look at the dynamics of these parameters, we conducted a series of simulations comparing the predicted outcome of pathogen competition. For each set of simulations, we randomly sampled 5000 values from a uniform distribution allowing variation of one or more of the following: • tick-to-host transmission for pathogen 1 between 0 and 0.1 hosts/tick/month, • host-to-tick transmission for pathogen 1 between 0 and 0.05/month, and • cofeeding transmission between 0 and 0.1/month. I I ® Figure 2. These graphs show the importance of both transmission rates in the short and to a greater degree long term. Red circles represent outcomes with coexistence of pathogens. Green circles are when pathogen one only exists, and blue is when pathogen two only exists. Magenta circles represent when neither pathogen exists. While there is coexistence of both pathogens in the first year (a), pathogen 1 begins to disappear in the second year for the lowest transmission rates (b). Only the higher transmission rates are able to have coexistence of the pathogens in the longer times (c) and (d). Standard baseline values are noted by the + symbol.
Ranges for transmission rate parameters were calculated based on reasonable limits determined from experimental and field observations (Nadolny et al., 2014;. Our model is deterministic, thus a pathogen will never truly go extinct in real time. Therefore, extinction in our model was assumed when the prevalence rate goes below 0.05% of ticks infected with a given pathogen. This threshold was chosen based on the likelihood of collecting an infected tick through active surveillance combined with the current limits of detection thresholds for molecular pathogen identification techniques. The scatter plots in Figure 2 show the relationship between tick-to-host transmission rate (A 1 ), host-to-tick transmission rate (Â 1 ), and pathogen persistence. Here we see that variations in both host-to-tick transmission and tick-to-host transmission rates are important to the long-term survival of the pathogen.  short-term coexistence of the two pathogens before the loss of the pathogen with the lowered transmission parameters. Lastly, a two-way sensitivity was performed between tick-to-host (A 2 ) and cofeeding transmission (μ 2 ) on the IRN of pathogen 2. All other parameters remained constant at baseline levels. If cofeeding transmission can control the ability of pathogen 2 to invade with low levels of tick-to-host vector transmission this could highlight the importance of cofeeding in the persistence and prevalence of pathogen 2. However, our results in Figure 3 support the findings of Gaff and Gross (2007) that this system is more sensitive to tick-to-host vector transmission than cofeeding. For this set of simulations values altered were: • tick-to-host transmission for pathogen 2 between 0.05 and 0.55 hosts/tick/month and • cofeeding transmission for pathogen 2 between 0.002 and 0.02/month.

Case study simulations and results
While the model is intentionally described in general terms, as we suggested previously we are immediately interested in exploring the transmission cycle of R. parkeri, represented as pathogen 1, and R. amblyommatis, represented as pathogen 2, infections in the A. americanum populations found in the southeastern region of Virginia. Here we applied the model to a specific case study. This case study takes place at a reclaimed swampland adjacent to the Great Dismal Swamp National Wildlife Refuge. The site is located in Chesapeake, Virginia, and is managed by the Nature Conservancy. The plants and habitat  of this site have changed greatly over the past decade from farm fields to natural swampland. This site was one of the first locations in Virginia to have an established A. maculatum population, which was found to have a high prevalence of R. parkeri. Although A. maculatum ticks population disappeared in 2014, there are A. americanum ticks, which can transmit the residual R. parkeri and be a concern for public health officials. By applying the model to this case we hope to find how long R. parkeri is expected to persist within the population of ticks and be a risk for humans. To simulate this case study using our model we included specific pathogen dynamics for transmission of R. parkeri and R. amblyommatis by the single tick species A. americanum. Amblyomma americanum is not the primary vector of R. parkeri, and thus in our system this bacteria is considered to be less likely to survive long term in A. americanum, when compared to R. amblyommatis. We altered parameters in reference to data collected from previous work as well as laboratory findings based on known transmission dynamics of the two bacteria (Nadolny et al., 2014;. Parameter changes and reproductive numbers can be seen in Table 2. Rickettsia amblyommatis is prevalent and well maintained within ecosystems in southeastern Virginia with A. americanum ticks present. Therefore, while the exact values for the hostto-tick and tick-to-host transmission for R. amblyommatis are not known, we assume that they are much greater than those for R. parkeri. Furthermore, transovarial and transstadial transmission of R. amblyommatis are also more likely in A. americanum and were increased from the baseline value accordingly. Coinfection transmission values are half of the respective single pathogen transmission, as well as transovarial and transstadial transmission for coinfection are half of the respective single pathogen rate. All other values remain unchanged including initial values for the state variables. ® Figure 4. These scatter plots show the results of the scenario assessing how long a less competitive pathogen, R. parkeri could persist. Red circles represent outcomes with coexistence of pathogens. Green circles are when pathogen one only exists, and blue is when pathogen two only exists. Magenta circles represent when neither pathogen exists. While all parameter sets allow for coexistence in the first year (a), the less competitive pathogen is completely eliminated by the fifth year (c). The more competitive pathogen is able to persist for nearly all parameter sets.
Using the applied case study parameters, we conducted a series of simulations comparing the predicted outcome of pathogen competition. We followed the design of the simulations as described in Section 4, and we repeat the design here for clarity. For each set of simulations, we randomly sampled 5000 values from a uniform distribution allowing variation of one or more of the following: • tick-to-host transmission for R. parkeri, between 0 and 0.1 hosts/tick/month, • host-to-tick transmission for R. parkeri, between 0 and 0.05/month, and • cofeeding transmission for R. parkeri, between 0 and 0.1/month.
The results of the case study simulation are shown in Figure 4 and demonstrate the BRN/IRN results presented in Section 3 that predicted no long-term persistence of R. parkeri, with only the host-to-tick transmission value decreased. The scatter plots also demonstrate that R. parkeri is maintained for less than five years in the model, while  R. amblyommatis is consistently maintained for most parameter choices. Amblyomma americanum as a vector alone does not allow for persistence of R. parkeri, which is not unexpected based on our observational and experimental data, but the one to three-year time frame presented does provide new insight. The continuing persistence of R. amblyommatis is also consistent with what we expect. Based on our results, R. amblyommatis will remain in the A. americanum population, but R. parkeri will present in these ticks for less than 5 years. Therefore, without any other vectors in this system, it is unlikely that R. parkeri would be a long-term concern for public health officials. However, if other factors such as host recovery or competent reservoirs for pathogens were introduced the disease would be likely to persist and could be analysed in the future.

Conclusions
A mathematical model for two-pathogens, one tick, and one host was presented and explored numerically. Not surprisingly, the model simulation results suggest that tick-tohost transmission is important for persistence of pathogens on any time scale. Cofeeding as well as host-to-tick transmission plays less of a role in pathogen survival. To add to these results, the model was applied to a case study where available data suggest approximate initial conditions and parameter estimates. The BRN and IRN for R. parkeri indicate that even if it is introduced into the system it will only be a public health issue for a few years, regardless of how infectious R. amblyommatis is. The results from this site demonstrate that coexistence of both pathogens long term is unlikely. Also the long-term persistence of R. amblyommatis is found to be quite likely, which is not a concern to public health officials as it is nonpathogenic to humans. Future field study data and laboratory results will allow us to pinpoint parameter estimates and give further understanding into prevalence of these pathogens at specific sites potentially with other vectors species, and thus allow additional understanding of potential risk of tick-borne diseases in those areas.