First, we model the growth of domestic vectors by the using a delayed logistic term:

where the delay τ is the gestation time of the vectors, K_Vj is the carrying capacity of the vectors in village j, and d_hj(t−τ) is the egg hatching rate in village j at time t−τ. Since adult female vectors lay eggs after having a complete blood-meal, the egg hatching rate is dependent on the biting rate B_j(t−τ), which is the average number of bites each vector makes per day. In addition to the biting rate, the egg hatching rate is naturally assumed to be dependent on the number of eggs a female vector lays after a blood-meal ϕ_l, the proportion of adult females in the population v, and the proportion of eggs that hatch ϕ_h. Thus, the egg hatching rate is given by

Since a larger blood supply results in more biting (up to a maximum per day), the biting rate B_j(t) is dependent on the domestic blood supply b_supj(t) of village j at the time when the female lays eggs, and it is also dependent on the season [29]. In order to accurately capture the seasonal dependence, we construct a seasonal biting rate b(t) (defined in Table 3) based on data from [28]. However, this function b(t) is obtained from data in a setting with a particular known blood supply b_supknown(t) [defined in equation (43)] and does not capture the blood meal dependence, so we also use a Holling Type II response to obtain the final biting rate

where b_supj(t) is the domestic blood supply in village j, β is the maximum possible daily feedings of a vector, and b_max is the maximum value of b(t). A_b is chosen such that

This is done so that if the blood supply in village j matches the conditions in [28] (that is, b_supj = b_supknown), then B(t) = b(t) and the biting rate agrees with the empirical results.

We now consider the domestic blood supply b_supj(t), which is composed of total humans, dogs, and chickens. Since the vectors prefer certain blood sources over others, each one is translated into a number of humans, so that the unit of measurement is human factors [28, 30]. More specifically, the human factors for a dog and a chicken are represented by d_f and c_f, respectively. Also, the blood supply is dependent on the availability of humans, dogs, and chickens in the domestic region in village j at time t, given by a_Nj(t), a_Dhj(t), and a_Chj(t), respectively. These time-dependent values represent the proportion of the respective populations available for biting in the domestic region, which may vary due to various factors, such as dogs sleeping outside or people being unavailable for biting because of their protection (at night) when they are using bednets [28]. Therefore, the total blood supply is given by

which can be thought of as the equivalent number of human factors available for biting.

We now consider the deaths of vectors. Natural death is modeled by

where d_m is the natural death rate. We also assume that the vectors die due to overpopulation, which is modeled by

where d_k is the death rate due to the population being over the carrying capacity. In addition, death can be caused by insecticide spraying. We assume there is a sub-population of vectors V_resj that does not die from spraying, being protected in cracks in the walls or by insecticide resistance. Recent field surveys confirm that domestic vectors that survive spraying continue to be a source of reinfestation [18]. For simplicity, we assume that V_resj is constant for a given village. Therefore, the term that models death due to spraying is given by

where r_j(t) is the mortality rate due to spraying of non-resistant vectors.

As was done in [26], we assume the net movement into (or out of) the domestic region from (or to) the peridomestic region of a village is given by

Here, ρ is a constant parameter with units vectors per day, and the term inside the parentheses captures the density dependence by comparing the ratios of the vector population to the carrying capacity in the respective regions. We note that the movement of vectors between the domestic and peridomestic region is still poorly understood, and a more thorough analysis of this term can be found in [26].

A novel aspect of our model is the migration of vectors between villages. Vectors move between villages primarily in two ways: passively via human movement and actively via flight [15]. However, even with villages within approximately 500 m of each other there is sometimes no flight of vectors between them [15], so we consider active transport negligible in our village system. For passive migration, which occurs when vectors and their eggs are carried by traveling individuals, we assume that the rate is dependent on the number of vectors and humans present in the domestic region of each village in the following way. We let α_{Nj, k} be the daily rate of people traveling from village j to village k, and use η to denote the average ratio of vectors that live in the luggage or other travel items in the domicile. Then the number of people traveling daily from village j to village k is given by α_Nj,kN_j(t), and on average, each person is carrying ηV_j(t)/N_j(t) vectors. Thus, the number of vectors per day passively transporting into (and out of) village j, and from (and into) village k, is

where, p_Vj,k = ηα_Nj,k.

Finally, the complete equation used to model domestic vectors is

A domestic vector becomes infected by biting infected humans and infected dogs [31]. Thus, the term that models the growth of the infected domestic vector population is

where P_NV and P_MV are the proportion of vectors that become infected after taking a blood meal from infected humans and mammals (including dogs). Note that B_j(t)(V_j(t)−V_ij(t)) is the daily number of bites by non-infected vectors and that the latter half of (12) is the fraction of bites that cause infection. The vectors do not pass T. cruzi to their young, so this is the only term for the growth of the infected vector population [32].

The death rates for infected vectors are assumed to be the same as those for all vectors with no preference for the infected or susceptible populations, so the terms accounting for mortality are the same as those for total domestic vectors but multiplied appropriately by V_ij/V_j. The transport rates between populations of vectors are subject to the same assumption, and they are correspondingly multiplied by the appropriate ratios of infected vectors.

Thus, the complete equation is

Naturally, the dynamics of the peridomestic vectors are similar to those of the domestic vectors; however, there are some minor differences due to the different setting. First, the egg hatching rate is updated to depend on the peridomestic blood supply, b~supj(t), which is given by

Here, a_Dpj(t) is the peridomestic availability of dogs, a_Cpj(t) is the peridomestic availability of chickens, a_Mj(t) is the availability of mammals, and m_f is the mammal factor analagous to d_f and c_f [28]. Thus, the peridomestic biting rate B~j(t) is given by

and the peridomestic egg hatching rate is given by

Second, we assume that luggage containing items such as clothing is what harbors vectors for passive transport [15], so as there is no luggage in the peridomestic region there is correspondingly no term for passive transport (We relax this assumption in Simulation 4 and consider a case with transport of peridomestic vectors). Therefore, we have

The dynamics of infected peridomestic vectors are derived from those of the total peridomestic vectors just as for the domestic vectors, so the equation is given by

We now consider susceptible humans. We note that all children of susceptible humans are born susceptible, and a certain fraction P_NN of children born to infected mothers are also born infected due to congenital transmission with the rest being born susceptible [33–35]. Also, due to the relatively long time it takes for Chagas disease to cause death, the birth rate of infected mothers is assumed to be the same as that of susceptible ones. Thus, using a logistic model, the growth term is given by

where G_N is the daily growth rate of humans with unlimited resources. For the death term, we use

where γ_Ns is the per day death rate of susceptible humans.

Along with death, we consider the loss of susceptible humans due to infection. The rate of infection is taken to be the product of the number of bites per day by infected vectors and the fraction of those bites that cause human infection, which yields the following term:

Finally, as we consider multiple villages, we include terms for the movement of people between them. We assume this movement depends on a per-person travel rate, so that the rate of people entering and leaving village j (from and to village k) is given by

where α_Nj,k is the daily rate of movement of people from village j to k.

Thus, the complete equation is

The model for infected humans is similar to that of susceptible ones. However, the birth term reflects that only infected mothers can give birth to infected children and that only a certain percentage of their children are born infected. Also, the susceptible humans that become infected join the infected population, and the death rate of infected humans is assumed to be higher. However, the travel rates are assumed to be the same. Thus, the equation is

We first consider the birth rate of infected dogs. To do this, we consider the total population of dogs in village j, D_j(t), which we assume is known. Then, assuming a growth rate of α_j(t) for the total dogs and that susceptible and infected dogs die at the same rate of γ_D, we have Dj′(t)=αj(t)Dj(t)-γDDj(t). Finally, since D_j(t) is known, and hence, Dj′(t) is also known, we also have α_j(t) in terms of known quantities:

Therefore, as infected dogs are only born to infected mothers, the birth term for infected dogs is given by

where P_MM is the proportion of dogs infected by vertical transmission.

In addition, the infected dog population could grow due to infected vector bites on susceptible dogs. Therefore, the term we use is similar to that used to model infection of humans. However, since dogs inhabit both domestic and peridomestic regions, there are two separate terms to model infection using appropriate respective biting terms, blood supplies, populations, and availability values. Thus, the terms used are

and

for the domestic and peridomestic regions, respectively.

Then, including the death term, the complete equation is given by

Finally, we consider infected mammals. This equation is similar to that for dogs, but it differs slightly in that mammals do not enter the domestic region and thus there are no terms for infection from domestic vectors. Therefore, we have

Thus, the full model for village j is given by the following eight differential equations and the initial conditions found in Table 2:

We now establish the baseline parameter values. We separate these values into two categories: those considered to be biological in nature and generally consistent in all villages and those that depend on the structure of the particular collection of villages being modeled. For example, we consider the vector gestation time and the proportion of vector bites that result in a human being affected to be the same in any village, whereas values such as the number of houses in a village and the availability of dogs in the domestic region clearly depend on the particular village in question. The parameter values used are based upon existing data of Chagas disease and the Gran Chaco region when available, but we note that these values have a high degree of aleatory variability and epistemic uncertainty. The constant biological parameters are found in Table 1. We note that the human growth rate G_N was chosen so that the total human population grows at an approximate annual rate of 1.0% over the 30 years of simulation. This is consistent with recent population data for Bolivia, Paraguay, and Argentina, the three countries that contain the vast majority of the Gran Chaco region [46–48]. In addition, the infection proportion P_VN is unknown, so we choose a value that agrees with the T. cruzi seroprevalence in humans of 51.7% found in surveys in the Bolivian Chaco [9]. The village dependent values are found in Table 2 with citations where the value was chosen based on empirical data.

We now consider the time-dependent parameters. Based on data from [28], the number of goats in village j is given by

where H_j is the number of houses in village j, which yields a (smooth) function with a period of half a year. This function is extended periodically. Similarly, the number of pigs is also derived from [28] and it is taken to be yearly periodic. The (smooth) function for the first year is given by

Combining these two functions, we have that

For availability, we first define the total chicken availability a_Cj(t), goat availability a_Gj(t), and pig availability a_Pj(t), all of which are obtained from data in [28]. We use continuous, yearly-periodic, piecewise-linear functions defined by the values listed in Table 3 and shown in Figure 1. We assume that there are an equal number of chickens in the peridomestic region and the domestic region, so we take a_Cpj(t) = 0.5*a_Cj(t) for the peridomestic chicken availability and a_Chj(t) = 0.5*a_Cj(t) for the domestic availability. We then take a_Mj(t) to be the weighted average of a_Pj(t) and a_Gj(t), yielding

For the biting term, we again use a yearly-periodic piecewise-linear function as defined in Table 3. For the known blood supply, we use

where a_C(t) and a_M are defined the same as a_Cj(t) and a_Mj, and C_known(t) and M_known(t) are given by H_known = 74. Both of these functions are also derived from [28]. Finally, we define a yearly periodic active spraying mortality function, r(t), with the first year given by

with t₁ = 212.5, t₂ = 303.75, and

as defined in [25]. This corresponds to spraying at t₁ with residual effects that diminish until disappearing at t₂. We take r_j(t) = r(t) in spraying years for village j and r_j(t) = 0 otherwise. The time-dependent parameters, with the exception of the active spraying mortality function, are shown over the course of 1 year in Figure 1.

For the baseline case with resistant vectors, the simulation is run for 30 years, during which insecticide spraying in village 1 occurs from year 11 to year 20. This is done so that the populations can reach steady oscillations before spraying begins and so that the behavior of the populations before, during, and after insecticide use can be analyzed. The results of this baseline simulation may be seen in Figure 2.

Before spraying begins, all the non-human populations demonstrate almost steady oscillations. These oscillations continue for the full 30 years in village 2. In village 1, after spraying starts, the vector populations quickly reach new smaller steady oscillations while the infected populations decrease. Once insecticide use stops, the non-human populations rapidly return to the same oscillating patterns that they demonstrated before it began, and the number of infected humans begins to increase again. This is largely in line with the qualitative results found in the baseline simulations in [26], although there is a more noticeable decrease in the infected populations in the current work as a result of the updated parameters.

In this section, we simulate the populations of two villages that are connected by human travel, including the passive transport of vectors. We use the baseline parameters in both villages and include insecticide spraying from year 11 to year 20 only in village 1. This scenario investigates how isolated rural villages with different control measures might affect each other, and in particular, the role village 2 might play as a reinfestation source for village 1. Here we simulate the case where travel between the villages is common, with the average person traveling to the other village once per week (α_N1,2 = α_N2,1 = 1/7). Additionally, we specify that 10% of the vectors live in luggage or other travel items in the homes. The results are shown in Figure 3.

We see that during the spraying years, the domestic vector populations in village 1 and village 2 are higher and lower, respectively, as compared to the baseline case. These changes reflect the mixing of the vector populations in the two villages due to the passive transport of vectors. More strikingly, there is a complete mixing of the human populations between the villages so that they have the same number of infected humans despite the difference in control measures. While the travel between the villages does increase the number of human infections in village 1, it is important to note that the total number of human infections in the two villages after 30 years is identical to the baseline case. That is, the additional infections in village 1 are offset by fewer infections in village 2. Thus, the movement of infected humans and vectors from village 2 to village 1 spreads the disease burden between the villages but does not increase the total number of infected humans.

It is not clear from this simulation to what extent the change in the infected human populations is a result of the humans traveling vs. the vectors being transported. To isolate the effects, we ran a similar simulation with the same level of human travel but with no vector transport. The results for the human populations were identical to those seen in Figure 3 while the rest of the populations had dynamics similar to the baseline case. We see that the transport of vectors has little effect on the number of infected humans when the human populations are mixing at such a high rate. These results were further confirmed in scenarios with much lower rates of human travel. For example, when humans travel once per year between the villages, it takes about 22 years for the two human populations to fully mix, but the overall number of infections remains the same. And, removing vector transport in this case also has no effect on the total number of human infections in either village.

Overall, we see that vector migration is not a significant source of reinfestation in this scenario and that travel redistributes the number of human infections, but does not increase them. The human travel in this scenario is similar to migration as people travel to and stay in the other village. As such, the eventual complete mixing of the populations is likely unrealistic, potentially obscures the differences between the two villages, and does not increase the total number of human infections. Thus, to further investigate the potential effects of vector migration, we will consider in Simulations 2–7 human travel that passively transports the vectors while the humans themselves do not move to the new village. These simulations represent scenarios where humans transport objects containing vectors but do not stay in the other village.

In this section, we simulate the populations of two villages with the same parameters and insecticide spraying as in Simulation 1, but with the travel terms removed from the human equations. That is, we consider the baseline simulation with the addition of passive vector migration but not human travel. The results are shown in Figure 4 where we see that, except for infected humans, all the population dynamics are similar to those in Simulation 1 (Figure 3).

Because we have removed human population mixing, the differences in the number of human infections are attributable solely to the passive transport of vectors. As in Simulation 1, we see more vectors in village 1 during spraying years as compared to the baseline case. These additional vectors in village 1 that come from village 2 do result in five more infected humans in village 1 after 30 years as compared to the baseline. However, there are six fewer infected humans in village 2 in this scenario. So while the disease burden has slightly shifted from village 2 to village 1, the overall effect is small and the number of human infections across both villages remains nearly constant. We see that the transport of vectors does not have a significant impact on human infections in this scenario.

We further investigate the effect of the travel parameters in this section by considering the number of infected humans in both villages after 30 years as a function of P_V1,2 = P_V2,1 = ηα_N1,2. We recall that α_Nj,k is the daily travel rate of humans from village j to village k and η is the average ratio of vectors that live in the luggage or other travel items. As in Simulation 2, to isolate the effect of vector migration, the humans transport the vectors between villages but do not stay in the other village. All other parameters are held at the baseline values in both villages, with insecticide spraying in village 1, but not in village 2. The results are shown in Figure 5.

Overall, we see that vector migration has a relatively small effect on the infected human populations across a broad range of travel scenarios. Small changes in the travel parameters have negligible effects on the infected human populations and even a change from no travel to high levels of travel and vector transport results in only a handful of additional infections in village 1 after 30 years. Furthermore, any increase in human infections in village 1 is offset by a similar decrease in infections in village 2.

Next we consider a simulation where vectors may be passively transported from the peridomestic regions of a village in addition to being transported from the homes. In this scenario, vectors may be hiding in materials (e.g., straw, hay, or animal feed) that are transported from the peridomestic region in one village to the peridomestic region in another during agricultural activities. To accomplish this in the model, we modify the equations for the peridomestic vector equations to include transport terms that are analogous to the transport terms in the domestic vector equations. The travel parameters are set to the same level as those in Simulation 2, while all other parameters are at their baseline values and insecticide spraying occurs only in village 1. As in Simulations 2–3, the humans transfer terms are removed from the equations, though the humans still transport the vectors. The results are shown in Figure 6.

As expected, we now see a mixing of the peridomestic vectors, just as we saw with the domestic vectors in Figures 3, 4. Additionally, we see differences in the infected dog and infected mammals populations as compared to those in Figures 3, 4, where there was no peridomestic vector transport. The decline and rebound of the infected human population is qualitatively similar to those in Figure 4, though the rebound after spraying ceases is more pronounced. After 30 years, the number of people infected in village 1 is about nine higher than the baseline case and four higher than Simulation 2, the case with only domestic vector transport. Notably, the number of infected humans in village 2 does not decline to offset the increase in village 1, but rather stays about the same as in Simulation 2. Nevertheless, the total infected humans across both villages is only four higher than the baseline case.

We ran additional simulations to see if other differences between the villages might affect the role that vector migration plays. We first considered the size of village 2, allowing village 1 to be connected to a smaller or larger village by varying the number of houses in village 2 from 40 to 145. We then considered the number of infected humans after 30 years using different rates of domestic vector migration. We did not include peridomestic vector migration or the transfer of humans between the populations. All other parameters are held at the baseline values in both villages, with insecticide spraying only in village 1.

When vector migration is low (α_N1,2 = α_N2,1 = 1/365 and η = 0.01), the number of infected humans in village 1 is essentially independent of the number of houses in village 2, varying by less than a single person over the whole range of house values in village 2. We then increased the travel parameters to their values from Simulation 2, see Figure 7. In this case, the number of humans infected after 30 years increases by about thirty people over the whole range of house values in village 2. This suggests that the size of the neighboring village with no insecticide spraying may be a significant factor in determining the effects of vector migration in spreading the disease.

Finally, we ran a similar simulation where we allowed the use of bed nets in village 2 to vary from 0% to 100% effectiveness. In both high and low vector migration scenarios, the use of bed nets in village 2 had no meaningful effect on the number of infected people in village 1. Overall, we see under a variety of circumstances that vector migration plays only a minor role in the spread of infection, except when there is substantial travel to a nearby larger village.

In the preceding scenarios where infected populations are already established, travel and vector migration has had only a secondary or marginal effect on the dynamics. Accordingly, in this section we consider scenarios where infections in village 1 are very low so that travel can play a potentially larger role in reintroducing infected vectors.

First, we modify Simulation 2 so that initially in village 1, there are no infected animals or vectors, only infected humans. We modify the corresponding susceptible populations so that the total populations are initially the same and leave the initial human populations, infected and susceptible, as in the baseline case. All other parameters are set to their baseline values, the travel parameters are as in Simulation 2, and insecticide spraying occurs only village 1. The results are shown in Figure 8 where they are contrasted with an otherwise identical scenario that involves no travel between the two villages.

We see that in this extreme scenario, the reinfestation of village 1 caused by the transport of vectors from village 2 does increase the number of infected humans above the level of infection that would occur without the travel between the villages. However, we see that in both cases the other infected populations rebound within a few years and the vector migration from village 2 results in fewer than seven additional human infections at the 30 year mark.

Next, we modify the above scenario so that there are no infected animals or humans initially in village 1 and adjust the susceptible populations accordingly to maintain the same total values at time 0. The travel parameters are set much lower to α_N1,2 = α_N2,1 = 1/365 and η = 0.01. Also, notably we do not include insecticide spraying in either village. The results are shown in Figure 9.

As expected, the introduction of infected vectors from village 2 leads to an eventual explosion of infection throughout all populations in village 1. The infected vector populations reach steady levels within about 5 years and the dogs and mammals reach steady levels soon thereafter. Even though the travel between the villages is very low, within 30 years, around 20% of the humans in village 1 are infected. However, we note that if travel with village 2 is entirely removed from this last scenario and instead we introduce a single infected vector in the domicile at time 0, we get nearly identical dynamics and a slightly higher number of infected humans after 30 years. It is clear that once the infection is introduced from any source, it is the local dynamics in the village that drive the further infection in the village and not the ongoing reinfestation of vectors through transport.

Finally, we look at the sensitivity of the baseline simulation to the number of resistant vectors in the village. Figure 10 depicts the number of infected humans after 30 years in a single, unconnected village as a function of the number of resistant vectors V_res = W_res. All other parameters are set to their baseline values and insecticide spraying occurs from year 11 to year 20, as in the baseline simulation.

We see that the first resistant vector has a substantial impact on human infections, leading to 16 more infections after 30 years. And, as expected, more resistant vectors leads to more infected humans, though at a decreasing rate. Figure 11 shows all the populations in the village in the extreme cases of 0, 1, and 40 resistant vectors. When there are no resistant vectors, the total number of vectors approaches zero during spraying years and does not recover for at least 5 years after spraying. However, with only a single resistant vector, the vector populations are able to rebound quickly to pre-spraying levels with 1–2 years. Additional resistant vectors lead to a quicker rebound of vector populations, but not significantly. Indeed, the presence of any resistant vectors allows the vector populations to rebound quickly.

Figure 11 also shows increased human infections with V_res = 40 as compared to V_res = 1, but this is due almost entirely to the sustained higher level of vectors during the spraying years, and not due to a quicker rebound of vector populations.