Model equations
Our host-parasite mathematical model involves the following host population components: ‘susceptible’ hosts denoted by (S), and hosts infected by k distinct types of parasites ((k=1,2,…,n)), the corresponding population numbers of infected hosts are denoted by (I_{i_1,i_2,…,i_k}), where each index (i_j) can take a value from 1, …, n (to avoid repeated counting of the same infection configuration, we require throughout the paper that (i_1<i_2<i_3<…<i_k)). Altogether our model contains (2^n) equations for the densities of S and (I_{i_1,i_2,…,i_k}). Importantly, we consider that parasites are distinct species, or are very different strains of the same species, which can not mutate into each other, so we do not explore the effects of kin selection. Note that our model is an extension of previous host-parasite models with multiple infections, such as the model by Choisy and de Roode25,38, or, more generally, well-known co-infection models in epidemiology31,39. The flowcharts of the model for (n=2) are shown in the supplementary material (SM2).
The model equations for S and (I_{i_1,i_2,…,i_k}) (with (k=1,…,n)) read as follows:
$$begin{aligned}{}&begin{aligned} frac{dS}{dt}=F(N)S-Sleft( sum _{j=1}^{n}phi _{i_j}+mu right) , end{aligned} end{aligned}$$
(1)
$$begin{aligned}{}&begin{aligned} frac{dI_{i_k}}{dt}=Sphi _{i_k}-I_{i_k}left( sum _{j=1,i_jne i_k}^{n}phi _{i_j}+mu +alpha _{i_k}right) , end{aligned} end{aligned}$$
(2)
$$begin{aligned}{}&begin{aligned} frac{dI_{i_{k_1},i_{k_2},…,i_{k_{{bar{n}}}}}}{dt}=sum _{i_j in {i_{k_1},i_{k_2},…,i_{k_{{bar{n}}}}}} I_{{i_{k_1},i_{k_2},…,i_{k_{{bar{n}}}}}backslash {i_j}}phi _{i_j}-I_{i_{k_1},i_{k_2},…,i_{k_{{bar{n}}}}}left( sum _{i_j notin {i_{k_1},i_{k_2},…,i_{k_{{bar{n}}}}}}phi _{i_j}+mu +alpha _{{i_{k_1},i_{k_2},…,i_{k_{{bar{n}}}}}}right) , end{aligned} end{aligned}$$
(3)
$$begin{aligned}{}&begin{aligned} frac{dI_{1,2,…,n}}{dt}=sum _{j in {1,2,…,n}} I_{{1,2,…,n}backslash {j}}phi _{j}-left( mu +alpha _{{1,2,…,n}}right) I_{1,2,…,n}{,} end{aligned} end{aligned}$$
(4)
where
$$begin{aligned} begin{aligned}{}&phi _{i_j}=I_{i_j}beta ^{{i_j}}_{i_j}+sum _{j_1=1,j_1ne j}^{n}Biggl (beta ^{{i_j}}_{i_j,i_{j_1}}I_{i_j,i_{j_1}}\&quad +sum _{j_2=j_1+1,j_2ne j}^{n}left( beta ^{{i_j}}_{i_j,i_{j_1},i_{j_2}}I_{i_j,i_{j_1},i_{j_2}}+…+sum _{j_{n-1}=j_{n-2}+1,j_{n-1}ne j}^{n}beta ^{{i_j}}_{{i_j,i_{j_1},i_{j_2},…,i_{j_{n-1}}}}I_{i_j,i_{j_1},i_{j_2},…,i_{j_{n-1}}}right) Biggr ), end{aligned} end{aligned}$$
(5)
is the infection rate of strain (i_j) infecting a host currently not infected by (i_j).
Our model assumes the mass action (bi-linear) mechanism of infection40. Here we suggest that individuals containing co-infections (i_1,i_2,…,i_k) can infect healthy or infected hosts in a way that the acquisition of only one new type of parasite is possible at a time, i.e. we neglect simultaneous double, triple, quadruple, etc. infections. In the model equations, (beta _i) is the transmission coefficient of the infection by host individuals having only a single type of parasite ((I_i)), whereas (beta _{i_1,i_2,…,i_k}^{i_j}), (i_jin (i_1,i_2,…,i_k)) denote the transmission coefficient of the parasite of type (i_j) from an infected host (I_{i_1,i_2,…,i_k}) to a healthy host or to a host which does not contain the parasite of type (i_j). The parameters (alpha _{i_1,i_2,..,i_{k}}) denote the virulence (i.e. an extra mortality of the host) due to the presence of parasites of types (i_1,i_2,..,i_{k}). Following previous studies14,25, we assume a trade-off between the virulence and the transmission rate (details provided below). Note that in the above model, we consider parasites to be obligate ones; in most cases, we assume that parasites are pathogens, however, the generic nature of the model allows its application to non-pathogenic types of parasites.
In the equation for S, F(N) is the host’s per capita reproduction rate (N denotes the total number of hosts). Here assume that F(N) is a decreasing linear function. (mu) is the background (unrelated to parasitism) mortality of the host. For simplicity, we disregard the possible recovery of infected hosts, the corresponding extension of the model can be easily done (e.g. see25,38). A summary of the dynamical variables, functions and model parameters is in Table 1. For the particular cases with double, triple and quadruple infections, i.e. for (n=1,2,3,4), the corresponding explicit equations are provided, for simplicity, in the supplementary material.
Parameterisation of model terms
We implement the standard assumption about the existence of a trade-off between virulence and transmission of pathogens25,38. In particular, we will use the following well-known hyperbolic parameterisation to describe the trade-off between transmission and virulence for a single infection of any type i
$$begin{aligned} begin{aligned} beta ^{i_1}_{i_1}=frac{beta _{0i}alpha _{i_1}}{K_{0i}+alpha _{i_1}}, end{aligned} end{aligned}$$
(6)
where (beta _{0i}) and (K_{0i}) are constants (for simplicity, in most cases we assume them to be independent of i). Parameterisation of the trade-off function in the case of co-infection is a more complicated matter since this should include mutual interactions between different types of pathogens competing for the same host. To the best of our knowledge, there is no universally accepted way of parametrising (beta ^{i_j}_{i_1,i_2,…,i_k}). Here we will use the following expression using a hyperbolic parameterisation
$$begin{aligned} begin{aligned} beta ^{i_j}_{i_1,i_2,…,i_k}=frac{beta _{0i}alpha _{j_i}}{K_{0i}+A^{i_j}_{i_1,i_2,…,i_k}}, end{aligned} end{aligned}$$
(7)
where (A^{i_j}_{i_1,i_2,…,i_k}) has the meaning of an effective virulence, i.e. a certain function of single-pathogen virulence (alpha _i), which describes the presence of all co-infecting pathogens in (I_{i_1,i_2,…,i_k}) in the ability to transmit pathogen (i_j) (we assume that (i_j in ( i_1,i_2,…,i_k))). For simplicity, we consider the following linear combination of (alpha _i)
$$begin{aligned} begin{aligned} A^{i_j}_{i_1,i_2,…,i_k}=sum _{kin (i_1,i_2,…,i_k)}C_kalpha _k{.} end{aligned} end{aligned}$$
(8)
Note that in the particular case with (C_k=0) (and (C_{i_j}=1)) we obtain (beta ^{i_j}_{i_1,i_2,…,i_k}=frac{beta _{0i}alpha _{j_i}}{K_{0}+alpha _{j_i}}), which coincide with the same expression as in25,38 which was suggested to take place in the absence of competition between the pathogens inside the host and without phenotypic plasticity. However, we should stress that even without a direct competition of pathogens for resources inside the host, co-infections can still affect the transmission of pathogens, for example, by reducing the contact rate of hosts. For the mentioned reason, we assume that generally all (C_k>1), even in the absence of competition for host resources: in most of our simulations, we considered (C_k=1). However, for comparison purposes, we also explored the scenario with (C_k=0), (C_{i_j}=1).
The parameterisation of the virulence for multiple infections is given by
$$begin{aligned} begin{aligned} alpha _{i_1,i_2,…,i_k}=sum _{kin (i_1,i_2,…,i_k)}alpha _k. end{aligned} end{aligned}$$
(9)
Note that the above expression describes the scenario of the absence of competition for host resources25: for a more general case, one needs to introduce weights in the above expression. In this paper, for the sake of simplicity, we do not include explicit competition for host resources which should be done elsewhere.
In this study, we explore long-term evolutionary dynamics using the Adaptive Dynamics framework, which considers ecological (epidemiological) time scales to be much faster than the slow evolutionary dynamics13,33,41. The essence of the method employs game theory by the introduction of a rare mutant with slightly different traits into the resident population at ecological equilibrium, or generally, an ecological attractor. This process occurs iteratively, with successive successful mutant invasions excluding the resident41,42,43,44. Following numerous such invasions and substitutions, the species evolve towards an evolutionary (and convergent) stable singular point called an Evolutionarily Stable Strategy (ESS). Similarly, simultaneous co-evolution of life traits of several pathogens results in a co-Evolutionarily Stable Strategy (co-ESS). The invasion fitness characterises the initial exponential growth (or decay) of a rare mutant and can be utilised in the analytical derivation of an ESS strategy: the evolutionarily singular point signifies the vanishing of the gradient of invasion fitness. We should stress, however, that for the number of co-infections (n>2), analytical expressions for invasion fitness become intractable, thus we use direct numerical simulations, where we successively introduce rare mutants after the system reaches the close vicinity of its ecological equilibrium. We apply a superinfection framework for mutants and the resident strains of the same type of pathogen38, i.e. different strains of the same pathogen cannot coexist in the same host. The equations defining invasion fitness for (n=1,2) and the corresponding flowcharts are provided in the supplementary material.
Measuring success of biological control
We need a measure of the success of biological control when using co-infections: we recall that the original goal is to suppress, as much as possible, the target host (pest) species by introducing co-infecting pathogens. The simplest measure of biological control by pathogens is the total number of individuals N of the host, which is expected to decrease as a result of control. However, this does not account for the effects of debilitation of the host due to infections. Indeed, a heavily infected pest generally produces considerably less damage to the ecosystem (e.g. less consuming of other species) compared to a healthy pest individual. Reduction in the fitness of the host caused by parasites should arguably be a certain function of infection load related to the virulence, which generally increases in the case of co-infections. As such, the damage of infected individuals (I_{i_1,i_2,…,i_n}) of the considered pest species on the environment should include some weighting. By following the logic above, we can use the efficient population size of the pest (N_e) defined as follows as a proxy for the success of biological control
$$begin{aligned} begin{aligned}{}&N_e=S+sum _{m=1}^{n}I_m w(alpha _m)+sum _{m_1=1}^{n}sum _{m_2=m_1+1}^{n}I_{m_1,m_2} w(alpha _{m_1,m_2})\&quad +…+sum _{m_1=1}^{n}sum _{m_2=m_1+1}^{n}…sum _{m_k=m_{k-1}+1}^{n}I_{m_1,m_2,…,m_k} w(alpha _{m_1,m_2,…,m_k})+…+ I_{1,2,…,n}w(alpha _{1,2,…,n}), end{aligned} end{aligned}$$
(10)
where (w(alpha )) are weighting functions. These functions should be non-linear and have the following properties: (w(0)=1) (in the absence of infection, the damage is maximal and given by 1), and for (alpha gg 1), we have (w(alpha )approx 0) (i.e. a heavily infected pest produces almost no damage to the environment). Various non-linear functions (w(alpha )) satisfy the required above properties. In this study, we consider the following parameterisations of the weighting functions (w(alpha ))
$$begin{aligned}{}&begin{aligned} w_1(alpha )=frac{a}{1+exp (b(alpha -alpha _0))}, end{aligned} end{aligned}$$
(11)
$$begin{aligned}{}&begin{aligned} w_2(alpha )= exp (-calpha ). end{aligned} end{aligned}$$
(12)
The above functions are the sigmoidal and the exponential dependencies, respectively; (a,b,c,alpha _0) are positive parameters (note that the parameter a is not independent, it can be found from the condition (w_1(0))=1). According to the scenario described by the sigmoidal function, the impact of the infected individuals on the environment is close to that of the healthy ones until the virulence attains some threshold value determined by (alpha _0): with the virulence higher than (alpha _0), the infected organism becomes too debilitated and its negative impact on the other species and the environment is close to zero. For the scenario modelled by (w_2(alpha )), the negative impact on the environment caused by an infected pest exponentially gradually drops with an increase of virulence starting already from low values of (alpha). Note that both expressions can describe the particular case, where only healthy individuals of the host would produce damage to the environment, in this case, we constrain the coefficients such that (b,cgg 1) and (alpha _0ll 1).
In this study, we also perform a simple cost-effectiveness analysis of biological control with co-infections. Here we apply two well-known cost-effectiveness metrics45, the incremental and average cost-effectiveness ratios, denoted by ICER and ACER, respectively. The mathematical expressions for the ICER and the ACER are as follows
$$begin{aligned}{}&begin{aligned} ICER=frac{Delta C}{Delta B}, end{aligned} end{aligned}$$
(13)
$$begin{aligned}{}&begin{aligned} ACER=frac{C}{B}, end{aligned} end{aligned}$$
(14)
where (Delta C) is the difference between two control strategies to be compared, (Delta B) is the difference between the corresponding effectiveness, C is the cost of a single control strategy, and B is the effectiveness of the considered strategy. For simplicity, we consider the scenario where the cost of biological control of introducing a natural enemy of each type is the same (given by (C_0), thus adding a new type of enemy to the existing ones results in an extra cost (Delta C=C_0)). The cost is additive, i.e. implementation of n types of enemies requires the cost of (C=nC_0), we can always assume that (C_0=1). We must stress that there can be other scenarios, with a non-additive cost of natural enemies, for example, in the situation where the cost of fieldwork to release natural enemies does not largely depend on the total number of enemies used, and the cost of breeding of natural enemies is low. Such cases, however, need to be assessed in some separate study.
To measure the effectiveness of control B, we use the normalised efficient population size (N_e(n)) for n co-infections, namely, (B=1-frac{N_e(n)}{N_e(0)}). Here (N_e(0)) denotes the pest density without biological control ((n=0)). The rationale behind this formulation is that the reduction in the negative impact of pests after applying the biological control is arguably proportional to the difference between the numbers (N_e(0)) and (N_e(n)), in particular, in the case (N_e(n) approx 0) we have (B approx 1), i.e. the maximal effectiveness. For the incremental effectiveness -while adding a new parasite to the existing (n-1) types- we have (Delta B =frac{N_e(n-1)-N_e(n)}{N_e(0)}).
