The dynamics of disease mediated invasions by hosts with immune reproductive tradeoff

The modern world involves both increasingly frequent introduction of novel invasive animals into new habitat ranges and novel epidemic-causing pathogens into new host populations. Both of these phenomena have been well studied. Less well explored, however, is how the success of species invasions may themselves be affected by the pathogens they bring with them. In this paper, we construct a simple, modified Susceptible-Infected-Recovered model for a vector-borne pathogen affecting two annually reproducing hosts. We consider an invasion scenario in which a susceptible native host species is invaded by a disease-resistant species carrying a vector-borne infection. We assume the presence of abundant, but previously disease-free, competent vectors. We find that the success of invasion is critically sensitive to the infectivity of the pathogen. The more the pathogen is able to spread, the more fit the invasive host is in competition with the more vulnerable native species; the pathogen acts as a ‘wingman pathogen,’ enhancing the probability of invader establishment. While not surprising, we provide a quantitative predictive framework for the long-term outcomes from these important coupled dynamics in a world in which compound invasions of hosts and pathogens are increasingly likely.

Biological invasions play a critical role in shaping the ongoing community dynamics of established ecosystems 1,2 . While invasions are a natural process of population growth, they are also facilitated by increasingly common influences of land use change, climate change, and anthropogenic transport (whether purposeful or accidental) of novel species to non-native habitats [3][4][5][6] . The ability of native populations to repel invasion attempts, or to survive successful invasions, depends on a multitude of factors 7 and are, of course, actively influenced by conservation/ restoration efforts 8 .
Invasions can be even more complicated when the invasive species is itself a pathogen that can use members of native species as hosts. Recent work has explored the fascinating dynamics of parasitic invasions 9,10 . While clearly a special case of a well-studied scenario in which a novel predator (the infectious agent, whether parasite or pathogen) arrives to potentially decimate a native prey (the host), this has been much more thoroughly explored within the scope of epidemiological dynamics of the introduction of novel pathogens, rather than through the lens of invasion ecology. As a result, most of the focus has been on establishing criteria for whether a novel pathogen will either force the native hosts extinct (and thus likely die out in the new environment itself), establish itself in the native host population as a new endemic disease, or sweep through the native population in one or multiple outbreaks only to result in herd immunity that eradicates the disease from a remaining host population (which could then experience another outbreak after sufficient migration or the birth of new hosts) [11][12][13][14] . Some work has gone so far as to explore the evolutionary implications. Research has suggested that introduction of a pathogen into a novel host should select for the evolution of changes in virulence over time 15 . When endemic competition is high, selection should favor increased virulence, which however also increases the probability of either host extinction or self-limiting outbreaks [16][17][18] . On the side of the host, introduction of novel pathogens has been proposed as a sufficient sudden selective pressure that it could lead to evolutionary rescue effects [19][20][21] . While it is unusual to discuss host evolution without also considering pathogen evolution due to the relative mutability and generation times of most host/pathogen pairs, there has now been empirical evidence of such dynamics in wildlife populations; e.g., 22 .
Returning to an ecological perspective, however, highlights a thus far less well studied aspect of the likely scenarios around the introduction of novel pathogens: they are likely to arrive as passengers within hosts who are themselves potential invaders, and subsequently infect the native hosts as well. Although some excellent

The model
Following the work in 36 , we construct an epidemiological model which tracks the disease dynamics and population of two species of hosts following the introduction of a pathogen. The native host (hereafter simply referred to as "type 1") is vulnerable to the disease, but due to being well adapted to the native habitat has high fecundity when uninfected. The invasive host (hereafter referred to as "type 2"), has coevolved defenses to the pathogen that increase both its tolerance of and resistance to the disease, but is not inherently as well-adapted to the habitat in the absence of infection (i.e., its intrinsic rate of growth in the new habitat is lower than that of the native).
Our initial conditions correspond to a population of uninfected type 1 hosts with a small number of both uninfected and infected type 2 hosts, representing an invasion by a novel competitor carrying a novel pathogen into the type 1 population. We consider a vector-borne pathogen, and make the simplifying assumption that there is an already abundant competent vector species in the habitat. (For this initial formulation, we considered a scenario of mosquito-borne infections in birds, such as avian malaria 37 or West Nile virus 38 , to motivate concrete choices.) The model couples two biological dynamics: the daily vector-borne spread of the disease among hosts, and a yearly host breeding cycle. We simulate in discrete time-steps that represent days using an SIR model taking into account the interactions between the disease, the two species of host, and the vectors. The model also includes a passive death rate for hosts of vectors, which increases for hosts while infected. While the vectors are assumed to breed daily, the hosts reproduce as part of an assumed annual breeding season, every t c time-steps (typically equal to 365). These dynamics were informed by considering an annually breeding bird population in a tropical environment, however, they are not meant to reflect the realism of any one biological system. They are chosen here merely to allow a clean interpretation of modeled scenarios. Future models should explore the impact of greater variety in the dynamics of possible vector and host reproductive patterns.
Epidemiological model. The model tracks eight variables corresponding to combinations of host species and vectors with their infection status. Hosts may be of type 1 or 2, and are either susceptible to the disease ( S 1 , S 2 ), currently infected ( I 1 , I 2 ), or recovered ( R 1 , R 2 ). We assume that recovery is complete and recovered www.nature.com/scientificreports/ individuals suffer no residual effects from their infection aside from a lifelong immunity to becoming reinfected.
(We later set the recovery rate for host type 1 to 0, so R 1 = 0 at all times, but leave it defined for the sake of generality.) For simplicity, we model using only one stage of infection in which individuals are both infectious and symptomatic. The model also tracks the status of the vector population, which may either be susceptible ( S v ) or infected ( I v ). We assume that vectors do not recover from the disease, but also suffer no negative effects from being infected, acting only as carriers.
For convenience of notation, we denote the total number of hosts and the relative frequencies of infection within their respective population which allows some equations to be written more compactly. Table 1 shows a summary of these variables. The model also has several constant parameters that affect the dynamics. β j determines the probability that hosts of type j become infected when bitten by a single infected vector. We typically set β 1 > β 2 , making type 2 hosts less likely to become infected.
Likewise, δ j determines the probability that a vector becomes infected when biting an infected host of type j. b j determines the bite rate for vectors on host type j. We assume that each vector bites the same number of hosts per day, so each vector's probability of becoming infected depends only on the frequency of infection among hosts, while each host will be bitten more if there are more vectors.
γ j determines the proportion of infected hosts of type j that recover from the disease each day. We typically set γ 1 = 0 < γ 2 , meaning infected hosts of type 1 do not recover, while infected type 2 recover after an average of 1/γ 2 days. µ j− determines the daily death rate for uninfected hosts of type j and µ j+ determines the death rate for infected host of type j. We typically set µ 1− = µ 2− < µ 2+ < µ 1+ , meaning uninfected hosts have the same death rate regardless of type, infected type 2 have a higher death rate than uninfected hosts, and infected type 1 have the highest. (Both susceptible and recovered hosts are considered to be uninfected.) Table 2 shows a summary of parameters related to the SIR dynamics. Equation 1 shows continuous ordinary differential equations approximating the dynamics. Note that the actual model instantiates these in discrete time-steps using the forward Euler method with h = 1.  www.nature.com/scientificreports/ Following a standard SIR model, susceptible hosts can become infected, and infected hosts become recovered, but each equation also contains a negative term corresponding to deaths. Thus, the total population of hosts is strictly decreasing in this time-frame. We assume that the vectors breed on a much shorter timescale than hosts, so we include a term for their births here, while host births are implemented by a yearly breeding event. We assume no vertical disease transmission, so all new vectors begin in the susceptible category. We assume that the daily birthrate for each vector increases with access to hosts, and decreases with competition among other vectors for hosts and breeding sites, so we set it equal to α v H S v +I v , where α v is a constant scaling factor. Since the birthrate for each vector contains the total number of vectors in its denominator, the total number of vector births in the population will simply be α v H.
A population with a larger number of hosts will be able to sustain a larger number of vectors. For a population with a constant number of hosts, the equilibrium vector population will be proportional to the number hosts: aH where a = α v µ v is the equilibrium vector density (number of vectors per host). For instance if a = 2 , then in equilibrium there will be twice as many vectors as hosts. Given a fixed number of hosts, the population of vectors will asymptotically approach the equilibrium value. In practice the total number of hosts is constantly changing, so the population of vectors will chase after this moving equilibrium, though for our standard parameters α v and µ v are sufficiently large such that this will occur on a short timescale, and the population of vectors remains close to the current equilibrium value.
Breeding event. Table 3 shows a summary of parameters related to the breeding event. Every t c days (typically 365), a breeding event occurs according to the following process. Let be the number of new host offspring of each type born this generation. In order to maintain consistency of temporal units among the parameters, each birthrate parameter is multiplied by t c . Let H be the current total number of hosts. Let be the proportion of offspring that survive to adulthood. (None, if the population is already above carrying capacity. All, if the difference between the reproducing population size and the carrying capacity exceeds the new births. If the population is approaching carrying capacity, juvenile mortality scales proportionally so that the population will hit carrying capacity but not exceed it.) Then www.nature.com/scientificreports/ We assume there is no vertical disease transmission, so all new hosts begin in the susceptible category. We assume that the host population is iteroparous, such that the new offspring and the existing adult population both carry over to the next generation. If the new population would exceed the carrying capacity, we assume the limited space or supplies reduces the number of successful offspring so that the population exactly reaches the carry capacity by reduction in juvenile survival rather than population-wide competition that could also reduce the adult population. The carrying capacity is therefore what drives the interspecific host competition. Because births of both species are summed and then normalized by the total number of births, the higher the birthrate of one host, the larger a fraction of the available space it will capture during the breeding event. Similarly, the lower the death-rate of a host, the less space it frees up for the next breeding event. Even if one host species would be able to sustain a stable population on its own, the presence of a more fit competitor can lead to the extinction of the less fit type by driving its effective birth rate down.
Immune-reproductive trade-offs and boundary conditions. We assume that host type 1 is evolutionarily stable in the absence of the disease; an uninfected monoculture population below the carrying capacity will have at least as many births as deaths each cycle. In a continuous version of this model where births and deaths happened simultaneously, this might be defined by α 1− ≥ µ 1− . However in our model, the population spends many days decreasing due to deaths before the next breeding event occurs. The population exponentially decays throughout the cycle, and then jumps up during the breeding event. The number of new host births is proportional to the number of hosts at the start of the breeding event, which will be the lowest value of any other time during the cycle. Thus, the birth rate needs to be high enough that the surviving hosts can compensate despite their diminished numbers. Taking this into account, we get the condition Which is a higher bound on α 1− than the simpler one above, but will be close to it if µ 1− and t c are small.
To implement the scenario in which type 2 has increased resistance and tolerance to the disease at the expense of overall fecundity, we implement the following boundary conditions: Type 2 hosts are less likely to contract the disease, and are able to recover from it, while type 1 lack the immunological strength to eradicate it completely. Additionally, while both types of host are weakened by the disease, type 2 suffer fewer negative effects. However, this stronger immune response comes at the cost of reducing their birth rate when compared to healthy type 1 hosts.
Due to the heterogeneous population, there is ambiguity in defining R 0 for the disease. The two types of host have different transmission rates and durations of infection, and will therefore be responsible for different amounts of disease spread. To resolve this, we define several related values. Let R j 0 be the R 0 of the disease in a homogeneous population of type j hosts: the average number of hosts infected (indirectly, through vectors) from a single infected host in a population consisting entirely of type j hosts.
We simplify the equation for R 1 0 since γ 1 = 0 . We define w to be the frequency of host type 1: w := (S 1 + I 1 )/H . Then R 0 for the vectors is which will also be the effective R 0 of the disease for the hosts in the mixed population.
For simplicity of results, we restrict to the case where type 1 is more infectious overall than type 2, in particular R 1 0 > R 2 0 . This allows us to avoid edge cases in simulation outcomes which are beyond the scope of this paper. We intend to lift this restriction and study these outcomes in future work.
Note. Although usual epidemiological model formulations can rely on the value 1 as the boundary condition for R 0 to determine the epidemic potential of an outbreak, in this case we are calculating effective R 0 in a dynamic host population, such that the decrease in disease spread due to saturation from recovered hosts and already infected hosts increases the actual thresholds. More accurate criteria require a technical and somewhat cumbersome analysis, which we leave for a future paper.

Results
The long-term behavior of the model is sensitive to parameter values, but does not depend on the initial conditions, provided the starting size for each population is nonzero. Thus, we focus on presenting analysis of the parameter space in the competition between hosts, rather than sensitivity to initial conditions. We classify outcomes for the system into one of four categories: 1. Failure to establish The invading host 2 population asymptotically goes to zero, while the host 1 population remains near the carrying capacity. 2. Coexistence Both host types survive at a stable level without going extinct. 3. Competitive exclusion The host 1 population decreases asymptotically to zero and is replaced completely by type 2 hosts. 4. Extinction Introduction of infection alters the system such that both host populations asymptotically go to zero.
We define a set of parameters that lead to coexistence, which we refer to as the 'default parameters' , shown on Table 4. All figures and numerical results are made using the default values for each parameter except when otherwise specified.
Additionally, we set the carrying capacity κ = 15000 , days per year t c = 365 , bite rate b j = 1 , and as initial conditions set S 1 = 14000, S 2 = 1300, I 2 = 200, S v = 14000 , and all other initial populations to 0. Although in general the vector transmission rate from the host types, δ 1 and δ 2 , need not be equal, for simplicity here we set them both equal to 0.05.
As intended, our default parameters yield Coexistence between the two host types (Fig. 1).
To observe the longer-term trends, we use the same default parameters and sample data points once each year immediately after the breeding event (Fig. 2), thereby smoothing out the yearly cycles in the population. Under www.nature.com/scientificreports/ this default scenario, the initial infection grows into an epidemic which reduces the host 1 population, which then causes the outbreak to recede. This in turn allows the host 1 population to recover until it triggers another smaller epidemic, again reducing their population. These oscillations gradually decrease in magnitude and the population approaches a stable equilibrium (we leave analytic characterization of these dynamics to future work). Because the total host population reaches the carrying capacity after each breeding event, the host 2 population varies inversely with the host 1 population. Similar behavior is observed over a wide range of parameters, with the equilibrium frequency of host 1 depending primarily on parameters that influence the spread of infection. Figure 3 shows the 200 year projected results for simulations using default parameters for every parameter except α v , which we set equal to 0.2a, where a is the desired vector density.
We observe that the host outcome is strongly dependent on vector density. Low vector density leads to the Failure to Establish outcome. As vector density increases, Coexistence occurs, with the frequency of each host population changing continuously with vector density. For high density, we observe the Competitive Exclusion outcome.  www.nature.com/scientificreports/ These outcomes are rooted in the infection dynamics. Figure 4 shows infection frequencies as a function of vector density. When the vector density is low, the pathogen is unsuccessful at spreading. At the same threshold observed in Fig. 3, there is a discontinuous jump in infection success. Afterwards, contrary to expectation, the infection rates for vectors and all hosts actually decreases as vector density increases. This decrease can be attributed to something akin to Simpson's paradox: each species of host maintains a constant level of infection in this region ( Fig. 4; we leave discussion of the technical details causing this constant level for a future work). Because type 1 hosts are more likely to be infected than type 2, the average infection level of the whole population increases or decreases alongside the type 1 frequency. Since vectors interact with hosts proportionally to their frequency in the host population, this in turn causes the vector infection rate to decrease as well. Once the type 1 hosts go extinct, the rate of infection in the host 2 population increases with vector density, as we would normally expect.
Extinction never occurs under the range of parameters shown in these figures, since even if every host becomes infected, the host 2 population can replace itself faster than it dies. Extinction can occur under different birth or death rates that do not guarantee demographic replacement for host 2. When populations are below the carrying capacity, the birth and death equations are proportional to the current population size, so the host populations will grow approximately exponentially, assuming the infection level is stable. If the exponent is positive, the population will increase until it reaches the carrying capacity. If the exponent is negative, the population will asymptotically approach 0. An approximation for this exponent would be the average birth rate of the host type, given the average frequency of infected and uninfected individuals minus the average death rate. (Note: While this is a reasonable approximation for most parameters, it is not quite accurate since breeding happens after a year of cumulative host death, therefore the size of the population that reproduces is smaller than its average size throughout the year.) In order to better show how the epidemic interacts with the host frequency equilibrium and extinction, we allow the birth rates of both host types to vary. In particular, we multiply the default values for α 1 and α 2 by the same fixed value , and construct a plot that shows the population state for a simulation after 200 years, as and vector density vary (Fig. 5). For generality, we show outcomes where varies from 0 to 1, although only values above 0.61 satisfy the boundary condition for uninfected host 1 being evolutionarily stable.
All four possible outcomes are achievable through different combinations of the and vector density (Fig. 5). When birth rates are low, Extinction can occur, regardless of vector density, though the threshold for survival is seen to be vector density dependent. When birth rates are high and vector density is low, no epidemic occurs and we see Failure to Establish. Under intermediate values of vector density, we see Coexistence between both host types. Additionally, we observe a continuous gradient of host frequencies, with more type 2 hosts as vector density increases. When density is high, we see Competitive Exclusion, with only type 2 hosts surviving.
We also observe a phase transition between extinction and reaching carrying capacity as the birth rate varies, with a small transition region between them. This occurs since populations near the phase transition will exponentially grow or decay with an exponent very close to 0, so the time required to reach equilibrium will exceed the 200 year horizon presented in Fig. 5. As the time horizon increases, the boundary between the extinction (white) and non-extinction (colored) regions in the figure becomes increasingly sharp (not shown).

Discussion
As an increasing number of animal (and plant) species move (or are transported) around the world, the dynamics of biological invasions get more complicated 39 . When those invasions are further complicated by involving pathogens or parasites, they can reshape the nature of entire ecosystems 40,41 . Our model has demonstrated that invasive hosts carrying vectorborne pathogenic "wingmen" can drastically alter the dynamics of invasions, meaningfully shifting the likely outcomes among the options: Failure to Establish, Coexistence, Competitive Exclusion, and Extinction. Although failed invasions are difficult to study in detail, and may go completely unnoticed due to their lack of substantial impact, examples of successful invasions assisted by a disease have been observed in several cases such as in squirrels 27,28 , moose and deer 30 , and crayfish 31 , among others 23 . The study of human infectious diseases further demonstrates how the global spread of vectorborne pathogens may easily be driven by the mobility of infected hosts, rather than solely by the expansion of habitat range of vectors (e.g. the global Existing studies have already considered the opposite scenario from the one here presented 43,44 in which native hosts have increased resistance to infection relative to their more susceptible would-be invaders. Invasion success in this case is unlikely, even if the invading hosts have a competitive fitness advantage in the absence of infection 43 . In that scenario, rather than acting as a wingman to the invaders, the pathogen acts as a protective barrier against invasion. Although their model is different from ours in many ways, (e.g. studying a direct transmission, SI disease model using a stochastic cellular automaton), our conclusions are mutually consistent. The success of either the native or invading host are both possible, as is stable coexistence, but that the outcome depends on the relative demographic and etiological factors in each host type, where increasing pathogen transmissibility shifts selective pressures and competitive advantage to favor the disease-resistant host.
Many models of vector-borne disease spread have also considered a dilution effect, where an increase in host diversity decreases the spread of infection borne by generalist vectors, typically by decreasing the density of highly competent hosts for the disease and mixing them with less competent hosts [45][46][47][48] . While carrying a native pathogen into a habitat with an additional novel host does increase available host diversity, dilution would only occur in the case in which the native host is less susceptible to infection than the invaders. Increased native susceptibility would lead the native host to amplify, rather than dilute, the disease risks to both populations (as our results show; Fig. 5). Of course, this can be further complicated by factors such as vector feeding preferences 49,50 . If vectors focus more attention on a single host type, vector bites will be more concentrated on a small group of hosts, increasing the contact rate between infected hosts and uninfected vectors. Thus, the dilution effect of adding more host types to an ecosystem may be overestimated by models that do not consider feeding preferences in generalist vectors.
Our model contains terms for vector bite rate b j , which are analogous to vector feeding preference but fail to conserve total number of bites per vector as the host frequencies change. In this way, our model highlights the need to consider the full ecological, evolutionary, and epidemiological complexity of systems in being able to understand and predict the interactions among, and trajectories of, host populations.
In our model, when disease outbreaks are similarly likely in both host populations (i.e. when large outbreaks occur in both host types or else in neither), only one host type should ultimately persist. Using vector density as the dial by which to tune the relative force of infection, we see that when the introduced pathogen is unlikely to spread in either host type, type 1 hosts will outcompete type 2. Conversely, when the infection is likely to spread in both host types, type 2 hosts dominate. Of course, vector density yields these observed results due to its action on the force of infection in each host population, but other factors in the model similarly impact the force of infection. Therefore, tuning any of these factors ( δ j , β j , b j ) would result in similar system-wide dynamics.
The evolutionary dynamics, in fact, depend very little on the actual disease severity, except in so far as severity affects the force of infection. For example, in any modeled scenario, if we were to multiply both the type 1 death rate attributable to infection, µ 1+ , and the type 1 infection rate, β 1 , by a factor of 100 (thereby keeping R 1 0 constant), there would be no change in the predicted outcome (excepting edge cases). If the disease fails to spread, relative death rate is irrelevant to the evolutionary outcome. If the disease spreads among the host 2 population then the host 1 population will still die out and the higher death rate will simply hasten this inevitable outcome. Even in the Coexistence outcome, the equilibrium frequency of type 1 hosts won't change significantly; a more deadly disease will lower the equilibrium infection level required to keep the host 1 population in check, but not www.nature.com/scientificreports/ the resulting frequency of type 1 hosts. (Future work is underway to explore the analytic boundary conditions of these dynamics.) An important result from our model is that increasing transmissibility of the infection increases the relative fitness of type 2 hosts, and therefore actually increases their equilibrium frequency (assuming extinction does not occur). The competitive evolutionary benefit outweighs the epidemiological cost. This wingman pathogen dynamic can therefore play a pivotal role in determining whether invasion leads to Coexistence or Competitive Exclusion. The more easily the disease spreads, the higher the frequency of the invasive host species in the resulting equilibrium compared to the native host, and this change in frequencies happens in a continuous way. The distinction between survival and extinction, however, depends more on the birth and death parameters, and happens in a discontinuous way: the population as a whole either goes to the carrying capacity, or to extinction, with no equilibrium in between. We attribute this to our choice for the vector/host interaction rate to depend on the ratio of vectors to hosts. If our model had instead assumed that decreasing host population necessarily implies decreasing host density, then resulting in decreasing opportunities for transmission, this would slow the spread of disease and should lead to the persistence of small host populations in cases where our model leads to extinction.

Conclusion
While invasion success is determined by a complicated and diverse set of environmental and ecological factors, pathogens carried by invasive hosts can alter the competitive landscape and significantly alter their probability of establishment. This is especially true in cases where, either due to accident or coevolutionary selective pressures, the pathogens cause only minor fitness costs in the invaders, but cause substantial costs to native hosts. We have shown how some cases of such vectorborne "wingman" pathogens allow for stable Coexistence of both host types where, in their absence, the invading species would have simply failed to establish a persisting community, and can even shift the balance entirely allowing for the displacement of the native entirely (Competitive Exclusion) instead of failing to establish in their new habitat. These results clearly demonstrate the need for more nuanced, community ecology perspectives on the epidemiological-ecological dynamics of invasions in a global world of increasing species movement of hosts, vectors, and pathogens.

Data availibility
All simulations and figures were generated by original code by the authors, available at https:// github. com/ kazar raha/ SIRVe ctorM odel. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.