A general multi-vector/multi-host model. We extend the SEIR epidemiological framework (e.g., 19), into model that incorporates multiple vector types and multiple host types. This multi-vector/multi-host model depicts multiple transmission pathways of WNV. We explicitly consider three types of vector representing three classes of feeding preferences: ornithophilic (mosquitoes that prefer to feed on birds, e.g. pipiens), mammophilic (mosquitoes that prefer to feed on mammals, e.g., molestus) and generalist feeders (mosquitoes with no particular preference, e.g., hybrid mosquitoes), and two types of hosts representing two competent classes: competent hosts (that can get infected and transmit the virus onwards to a susceptible mosquito, e.g., birds), and incompetent hosts (that can become infected but do not develop sufficient viremia to further transmitting the pathogen, e.g., humans) (Fig. 1A). A full description of the models is provided in the Supplementary Material.
We assume four possible epidemiological states: susceptible, exposed (i.e., infected but not infectious), infectious, and recovered (i.e., immune to re-infection); designated as S, E, I, and R, respectively. For the competent host category and for the three vector types we include a latency period of virus incubation, i.e., the period of time between initially becoming infected to the moment of becoming infectious. This window20 can be particularly important to consider in the case of mosquitoes, since the latency period and mosquito’s life expectancy often have similar time scales and may shift in different temperatures implying that some individuals may die in the period after becoming infected but before becoming infectious due to a relatively short life span. The dead-end host category can only get infected, but never become infectious. After a characteristic duration of infection, hosts are assumed to recover and become indefinitely immune to all future infections. Even though in reality immunity may wane, assuming lifelong immunity is reasonable here considering our time scale of interest. For mosquitoes, however, we do not incorporate a recovery compartment due to their relatively short life span, hence once infected they remain infected.
The transmission of the virus is assumed to only take place from vector to host, and from host to vector (i.e., while some cases of horizontal transmission and vertical transmission have identified in WNV 21, these are sporadic events and hence are not included in our framework as routes of transmission). Hence all new individuals added by birth are assumed to be susceptible. WNV transmission is assumed to be frequency-dependent, in that the number of blood meals taken up by a mosquito per unit of time is limited, and hence independent of host population size 22. We assume equal biting rate (number of bites) for the three vectors23, but we consider that they differ in their feeding preferences.
Because we are interested in transmission dynamics of WNV within a single season, we assume constant birth and mortality rates of hosts and vectors. While all vectors and hosts experience a natural mortality rate, hosts may also suffer from disease induced mortality when infected with WNV.
The case-study. The inclusion of three vector and two host types makes it possible to study: 1) how unequal feeding preference of vector types affects the risk of WNV invasion in different habitat types; and within each habitat type, 2) the relative contribution of each vector variant to the risk of WNV invasion; and 3) which attributes of the multi-vector/multi-host/virus system are most influential for WNV invasion risk. This model explicitly accounts for the pathogen transmission pathways arising between the multiple vector and host types (Fig. 1B).
Our mathematical framework was based with field-data using house sparrows (Passer domesticus) as amplifier hosts. This species was chosen as a model bird because it is a competent reservoir (i.e., amplifier host) for WNV 24, common prey of mosquito bites 25 and widely distributed throughout most continents, especially related to urban environments 26. Sparrows thus play a key role in WNV amplification and transmission to humans 27–29. In addition, humans were considered as dead-end species (i.e., they can become infected but do not develop sufficient viremia to further transmitting the pathogen) so that biting vectors are "wasting" a bite leading to a dilution effect.
Habitats characterization. The relative risk of WNV invasion, expressed by R0, was investigated along a natural-rural-urban gradient. For this purpose, three habitat types were considered as in Martínez-de la Puente et al.,10 characterized by different composition of vector and host abundances, and different feeding preference rates of the vectors. Urban habitat were those that contained a higher abundance of molestus and more densely populated by humans; rural habitats displayed a higher abundance of hybrid mosquitoes than pipiens or molestus; and natural habitats showed a higher abundance of both pipiens and birds, and a lower human density, than the other two habitat types. However, it is possible that our data reflect an idealized scenario, and that their validity is contingent upon the specific environmental conditions present in each habitat.
We also considered different scenarios for each habitat type for a total of nine study cases. The scenarios evaluated included a baseline scenario, where vector abundances were set to the estimated raw data for each habitat, as well as increased climate change and urbanization scenarios, where the number of each vector type and human population density were increased, respectively. Scenarios were chosen to reflect possible future changes in natural, rural and urban environments caused, such as by increasing mosquitoes due to higher temperatures, or by an expansion of human population as a result of anthropization of the environment. All vector and host abundances assumed can be found in Table S1.
Calculating R0. The R0 metric was calculated from the Next-generation matrix 30, which encapsulates all the possible transmission pathways of WNV (i.e., all combinations of hosts to vectors and vectors to hosts). When R0 > 1 this implies that the virus has the potential to invade, while R0 < 1 implies that the virus is not expected to successfully invade. To study the relative contribution of each vector type and their respective combinations to the invasion risk of WNV, we disregarded WNV prevalence at equilibrium 38. Here, WNV was assumed non-endemic in the population (as almost no WNV infected blood-fed mosquitoes were found).
For our multi-vector/multi-host model we define a 5 by 5 matrix (see Supplementary Material, A). The value in each cell is the respective transmission rate, βij , of vector i to host j (/or of host j to vector i). Each rate is a function of: 1) the biting rate, i.e., the number of times a host of type i (/or j) is bitten by a mosquito of type j (/or i) in a unit of time (noted as bp, bm, bh), and 2) transmissibility, i.e., the probability that the virus is successfully transmitted from an infected mosquito to a susceptible host (qij) (/or from an infected host to a susceptible mosquito, qji) (Table S2).
Model parametrization. Vector 19, 57 and host 31 data were obtained in the framework of different studies from natural, rural, and urban habitats (see 32 for further details on the environment characterization) in southwestern Spain, a Mediterranean area with active WNV circulation 33,34. All values used for parametrization of the model are shown in Table 2 (see Supplementary Material, B for further details on the numerical values used for parametrization).
Table 2
Parameter used in the numerical examples and respective sources. Values are taken for pipiens (eco)type, molestus (eco)type, and their hybrid mosquitoes as the vectors, and birds (Passer domesticus) and humans as the hosts.
Parameter
|
Description
|
Estimate
|
Source
|
Δp
|
Birth rate of pipiens
|
0.537
|
1
|
Δb
|
Birth rate of molestus (eco)type
|
0.537
|
1
|
Δh
|
Birth rate of hybrid mosquitoes
|
0.537
|
1
|
Δb
|
Birth rate of birds
|
0.023
|
2
|
Δd
|
Birth rate of humans
|
5.5 x10− 5
|
3
|
1/µp
|
Longevity of pipiens (eco)type (days)
|
1/0.029
|
1
|
1/µm
|
Longevity of molestus (eco)type (days)
|
1/0.029
|
1
|
1/µh
|
Longevity of hybrid mosquitoes (days)
|
1/0.029
|
1
|
1/µb
|
Longevity of birds (days)
|
1/0.0015
|
4
|
1/µd
|
Longevity of humans (days)
|
1/3.91x10− 5
|
5
|
bp
|
Biting rate pipiens (eco)type (bites/day)
|
0.333
|
6
|
bm
|
Biting rate of molestus (eco)type (bites/day)
|
0.333
|
6
|
bh
|
Biting rate of hybrid mosquitoes (bites/day)
|
0.333
|
6
|
qpb
|
Transmission probability from infected pipiens (eco)type to birds
|
0.8
|
6
|
qmb
|
Transmission probability from infected molestus (eco)type to birds
|
0.8
|
6
|
qhb
|
Transmission probability from infected hybrid mosquitoes to birds
|
0.8
|
6
|
qpd
|
Transmission probability from infected pipiens (eco)type to humans
|
0.88
|
7
|
qmd
|
Transmission probability from infected molestus (eco)type to humans
|
0.88
|
7
|
qhd
|
Transmission probability from infected hybrid mosquitoes to humans
|
0.88
|
7
|
qbp
|
Transmission probability from infected birds to pipiens (eco)type
|
0.16
|
6
|
qbm
|
Transmission probability from infected birds to molestus (eco)type
|
0.12
|
6
|
qbh
|
Transmission probability from infected birds to hybrid mosquitoes
|
0.09
|
6
|
1/σp
|
Latency period of infected pipiens (eco)type (days)
|
1/0.107
|
8
|
1/σm
|
Latency period of infected molestus (eco)type (days)
|
1/0.107
|
8
|
1/σh
|
Latency period of infected hybrid mosquitoes (days)
|
1/0.107
|
8
|
1/σb
|
Latency period of infected birds (days)
|
1/0.5
|
9
|
1/σd
|
Latency period of infected humans (days)
|
1/0.25
|
9
|
αb
|
WNV-induced additional death rate of birds
|
0.1
|
10
|
αd
|
WNV-induced additional death rate of humans
|
5x10− 7
|
7
|
1/γb
|
Duration of infection in birds (days)
|
0.18
|
8
|
1/γd
|
Duration of infection in humans (days)
|
0.91
|
11
|
1Wonham et al., 2004; 2Moschini et al., 2017; 3Pawelek et al., 2014; 4Marini et al., 2018; 5Blayneh et al., 2010; 6Vogels et al., 2017; 7 Bowman et al., 2005; 8Hartemink et al., 2007; 9 del Amo et al., 2014; 10Komar et al., 2003; 11Chowell-puente et al., 2004. |
Transmission rate (βij) was parameterized according to Vogel et al. 23, defining the probability as the average between the three temperature levels reported in that study (Table S2).
Host feeding preference (selection) was estimated using the Manly resource selection design II index35, a ratio that uses relative density as the measure of host availability (density-based selection ratio; ŵi) and was estimated for each vector type at follows:
where oi is the proportion of utilized host species i, and π̂i is the proportion of available host species. The measurement of the mosquito feeding preferences, following Hamer et al.36 and Simpson et al.37, relied on the estimation of oi as the proportion of blood-fed mosquitoes of each vector type on a given host (e.g., house sparrows) out of all available hosts (see Table 3 for parameters used). The Manly selection ratio equals 1 when mosquito feeding on host i is in equal proportion to estimated availability (no preference, i.e., random biting); is > 1 when a host is overused (i.e., more frequent feeding than expected by chance), and is < 1 when a host is underused (i.e., less frequent feeding than expected by chance). When estimating the Manly host selection ratio, host species that were not observed as blood meal sources but were identified in vertebrate surveys were given a blood meal value of one. Host species observed as blood meal sources but not observed in vertebrate surveys were given a density equal to the lowest observed hist density at each site.
Table 3
Parameters used in the numerical examples, calculations of the vector feeding preference indexes (wi) and abundances (Ni) in each habitat, with the respective sources.
Variable
|
Variable description
|
Natural
|
Rural
|
Urban
|
Total
|
Source
|
Mi
|
Total Culex complex mosquitoes captures
|
10,598
|
4,886
|
3,784
|
19,268
|
2
|
npb
|
pipiens (eco)type bites on Passer domesticus
|
0
|
3
|
3
|
6
|
1
|
nmb
|
molestus (eco)type bites on Passer domesticus
|
0
|
1
|
7
|
8
|
1
|
nhb
|
hybrid mosquitoes bites on Passer domesticus
|
0
|
5
|
4
|
9
|
1
|
cp
|
Total pipiens (eco)type blood-fed
|
16
|
12
|
8
|
36
|
1
|
cm
|
Total molestus (eco)type blood-fed
|
8
|
10
|
27
|
45
|
1
|
ch
|
Total hybrid blood-fed mosquitoes
|
9
|
14
|
16
|
39
|
1
|
Nb
|
Passer domesticus contacts
|
2,556
|
7,075
|
7,396
|
17,027
|
2
|
Nd
|
Human abundance by census
|
82
|
176
|
2,811
|
3,069
|
2
|
Nt
|
Total host contacts (all birds and mammals)
|
18,500
|
14,829
|
18,834
|
52,163
|
2
|
1Martínez-de la Puente et al., 2016; 2Ferraguti et al., 2016; 3 Ferraguti et al., 2018. |
Data for host parameters were based on i) the bird contacts (rather than raw samples) as an estimate for πi, being π̂i equal to the house sparrows contacts / contacts of all other hosts present in the study area; and ii) human census data. In natural areas, low values for the ŵi of the three mosquito types were assumed by default because no bites were recorded on Passer domesticus following Martínez-de la Puente et al.10 thus not allowing the calculation of feeding preference (Table 4). To account for uncertainty in the feeding preference estimate, due to low bite records, 100 simulations were run based on the fixed parameters (Table 3), and with wi randomly sampled from a normal distribution with means given by the preferences estimated in Table 4, i.e., wi ~ N(ŵi, 0.1). Overall, the way the feeding preferences are calculated reflect the different host compositions in each habitat.
Table 4
Mosquito feeding preference indexes, wi, for the three vectors based on the values obtained from Martínez-de la Puente et al., (2016)10, and Ferraguti et al., (2018)31.
Parameter
|
Variable description
|
Natural
|
Rural
|
Urban
|
ŵp
|
pipiens (eco)type preference on Passer domesticus
|
0.01*
|
0.52
|
0.96
|
ŵm
|
molestus (eco)type preference on Passer domesticus
|
0.01*
|
0.21
|
0.66
|
ŵh
|
Hybrids mosquitoes preference on Passer domesticus
|
0.01*
|
0.75
|
0.64
|
* Assumed values. |
Owing to the limitations of not being able to determine exactly how many mosquitoes of each (eco)type were present in each habitat (due to the high costs of conducting a genetic study characterising the total number of Cx. pipiens complex mosquitoes captured in 57), vector abundance was assumed to follow the proportions found in Martínez-de la Puente et al., 10 (see Table 5). In all cases, the risk of invasion of WNV into that population was assessed given by R0.
Table 5
Mosquito estimated abundance (Ni) for the three Cx. pipiens (eco)types in each habitat based on the values obtained from Ferraguti et al., (2016)32.
Parameter
|
Variable description
|
Natural
|
Rural
|
Urban
|
Np
|
pipiens (eco)type estimated abundance
|
48.5%
(5,138)
|
33.3%
(1,629)
|
15.7%
(594)
|
Nh
|
Hybrids mosquitoes estimated abundance
|
27.3%
(2,890)
|
38.9%
(1,900)
|
31.4%
(1,187)
|
Nm
|
molestus (eco)type estimated abundance
|
24.2%
(2,569)
|
27.8%
(1,357)
|
52.9%
(2,003)
|
Nt
|
Total vector estimated abundance
|
100%
(10597)
|
100%
(4886)
|
100%
(3784)
|
Note: The abundance of pipiens in the natural habitat was assumed to be Np = cp / (cp+ cm+ ch) * Mn = 0.485 * 10,597 = 5,138, where ci is the mosquito count of each pipiens (eco)type at the natural habitat and Mn is the total captured mosquitoes at the natural habitat. The same estimation was applied for the other vector types and habitats. |
Elasticity analysis. To determine the contribution of each vector and host related parameters to increasing or decreasing R0, an elasticity analysis was carried out 39. The elasticity of R0 to each parameter a is given by:
(2)
were kij is the element in row i and column j of the Next-generation matrix. Keeping all other parameters fixed, we measured how R0 reacted proportionally to change in parameter a. We systematically repeated this analysis for all parameters in each scenario and habitat combination (Table 2). The elasticity of R0 to the individual parameters were calculated numerically in R (v. 4.2.0), using the popbio package based on the parameter estimates shown in Table 2, and assuming no vector feeding preferences (wp= wm= wh = 1). A positive value for e(a) was interpreted as the parameter a leading to an increase in R0, while a negative value for e(a) was translated as the parameter a leading to a decrease in R0.