In a first approach, it is assumed that each house has the same contact rate with “fertilized females”, i.e., Γ_i = Γ for all i. Following the analysis of Bottomley et al. [52], the random variable F (size of the population of “fertilized females”) is replaced by its mean m_F for simplification and it is therefore possible to derive differential equations for the mean number of “fertilized females” m_F(t), the mean number of colonies m_C(t) and its variance σ²_C(t). It is reasonable to consider that the life duration of the “fertilized female” stage (which corresponds to the house-seeking phase) is shorter than the duration of its corresponding colony. Therefore dm_F/dt can be set to zero and the model simplifies to the two following equations for the mean number of colonies and its variance (see Bottomley et al. [52] for details):

In a mathematical sense, m_C(t) and σ²_C(t) are limiting values of the mean number of colonies and variance when the number of houses tends to infinity. However, these terms also well describe large but finite populations (i.e., > 100 houses), as observed in field situations.

Epidemiologists define the basic reproduction number (R₀) of a disease as the expected number of secondary infections arising from a single individual during his entire infectious period, in a population of susceptible individuals where there is no immunity and in the absence of interventions to control the infection [60]. The R₀ concept has also been widely used in ecology where it measures individual reproductive success under ideal conditions (i.e., the mean number of offspring produced over the lifetime of an individual) [61]. From these definitions, it emerges that when R₀ < 1, each individual (or infected individual if a disease is concerned) produces, on average, less than one new (infected) individual, and therefore the population (or the disease) will not grow. If R₀ > 1, the population will grow (or, from the disease's point of view, the “pathogen” is able to invade the susceptible population of hosts). This threshold notion is one of the most important and useful characteristic of the R₀ concept. This concept can be applied to the development of triatomine colonies (= the “pathogen”) in the population of houses (= the “hosts”). In the present triatomine model, the basic reproduction number for the mean number of colonies can be computed from Eq. (1) as:

When R₀ < 1, the triatomine species is unable to establish itself in the houses and therefore, m_C = 0. This value is a stable point. When R₀ > 1, the point m_C = 0 is unstable and the number of colonies can increase, up to an equilibrium point where m_C = (1/η) (1 – R₀⁻¹). The mean number of triatomine colonies is therefore dependent on the basic reproduction number R₀ and on the strength of density-dependence η (i.e., the “per-colony reduction in probability of starting a new colony when other colonies are already established”). At this equilibrium, assuming that the rate of disappearance of houses μ_H is much smaller than the other rates, the variance to mean ratio becomes R₀⁻¹. Indeed, when equating Eqs. (1) and (2) to zero at equilibrium, and equating R₀ = δ′/(μ_C + μ_H) ≈ δ′/μ_C, straightforward calculations give σ²_C/m_C = R₀⁻¹. Then, if R₀ > 1 (which enables triatomine colonies to grow in number), then σ²_C/m_C < 1, which implies that the distribution of the number of triatomine colonies across the houses is distributed more evenly than at random. However, field observations prove that triatomine bugs are generally aggregated amongst houses, and therefore some houses harbour more colonies than others [62].

The model should therefore be modified to take into account the aggregated distribution of triatomines and therefore to enable σ²_C/m_C > 1. Several mechanisms can generate such heterogeneity; among them, ecological factors that make some habitats less suitable than others. For example, triatomine establishment can be limited if triatomine refuge quantity and quality diminish [62] or when house walls are improved (this can be assimilated to a kind of “immunity” of the houses). Another likely mechanism is the heterogeneous exposure of houses to colonization which can also mimic the above mechanisms. Among various factors that can influence heterogeneous exposure, one can think for example to close favourable infested peridomestic habitats that can generate higher domestic exposure in the vicinity [18,47]. In the model, heterogeneity can be achieved by enabling a random contact rate of houses with “fertilized females”, i.e. replacing in the model the contact rate Γ, which was fixed, by independent and identically distributed random variables Γ_i (i = 1 to n) with mean m_Γ and variance σ²_Γ (and therefore a coefficient of variation of the exposure rate of v_Γ = σ_Γ²/m_Γ that measures the degree of heterogeneity). In that case, the model's development gives the following equation for the mean number of colonies at stable equilibrium:

Again, when R₀ < 1, the triatomine species cannot invade the environment and when R₀ > 1, colonies tend to a non-zero equilibrium given by Eqs. (4a) and (4b).

These equations indicate that at equilibrium, the mean number of triatomine colonies is, as before, dependent on:

• • the basic reproduction number R₀ of the colonies;
• • the strength of density-dependence η, and now also on;
• • the degree of heterogeneity v_Γ, through which the aggregative distribution of triatomine colonies within houses can be mimicked.

More precisely, Eq. (4a) indicates that the mean number of colonies increases as the degree of heterogeneity σ_Γ decreases. In fact, no house can harbour more than 1/η colonies, and even houses with very high rates of exposure cannot have a corresponding high number of triatomine colonies. On the contrary, houses with rates of exposure close to zero will harbour a very low number of colonies. Thus, the mean number of colonies decreases as heterogeneity in exposure increases. The maximum number of colonies m_Cmax = (1/η) (1 – R₀⁻¹) is reached when σ_Γ = 0, and m_C ≈ 0 when σ_Γ → ∞.

Following the same approach as for the one species model, and first considering that exposure to each species is the same for all houses (i.e., Γ_1i = Γ₁ and Γ_2i = Γ₂ for i = 1 to n) and setting dm_F1/dt = dm_F2/dt = 0 (i.e., life expectancy of “fertilized females” is much smaller than that of a colony), the mean number of colonies for species 1 and species 2 can be modelled as:

Theses equations are symmetrical in regard to species 1 or 2, and are generalizations of Eq. (1).

Once again, the behaviour of the system is dependent on the basic reproduction numbers of the two species colonies and the inter- and intra-specific interactions parameters (i.e., the η_ij). For colonies of species i, the basic reproduction number is R_0i = δ’_i/(μ_Ci + μ_H). As in the “one-species model”, the mean number of colonies of species i will grow if the reproductive number R_0i > 1. If reproductive numbers for both species are greater than 1, then two scenarios are possible:

• • the exclusion of one species by the other and therefore the system approaches a single-species equilibrium, or;
• • the coexistence of both species in a mixed equilibrium.

As a first approach, the probability of establishment of both species can be considered as dependent on the availability of resources in the house (blood access and space). In other words, the effect of species i (i = 1, 2) on the establishment of species i or species j is the same. Therefore in the model, η₁₁ = η₁₂ and η₂₂ = η₂₁. With this assumption, competitive exclusion and coexistence can be understood by analysis of Eqs. (5) and (6) in the plane (m_C1, m_C2) (Fig. 1). In this plane, points P₁, P₂, P₃ and P₄ are computed as: P₁ = (1/η₂₁)(1 – R₀₁⁻¹); P₂ = (1/η₂₂)(1 – R₀₂⁻¹); P₃ = (1/η₁₁)(1 – R₀₁⁻¹); and P₄ = (1/η₁₂)(1 – R₀₂⁻¹).

P₁ corresponds to the mean number of colonies of species 2 (i.e., m_C2) when the isocline dm_C1/dt = 0 crosses the y-axis (m_C2). P₂ corresponds to the mean number of colonies of species 2 (i.e., m_C2) when the isocline dm_C2/dt = 0 crosses the y-axis (m_C2). P₃ corresponds to the mean number of colonies of species 1 (i.e., m_C1) when the isocline dm_C1/dt  =  0 crosses the x-axis (m_C1). P₄ corresponds to the mean number of colonies of species 1 (i.e., m_C1) when the isocline dm_C2/dt = 0 crosses the x-axis (m_C1).

Fig. 1a describes competition of the two species under the assumption η₁₁ = η₁₂ and η₂₂ = η₂₁ (equality of effects of species 1 (2) on the establishment of its own species or of species 2 (1))_. The two isoclines dm_C2/dt = 0 and dm_C1/dt = 0 are parallel and therefore a mixed equilibrium cannot exist. The species with the highest basic reproduction number will exclude the other. In the figure, species 2 will exclude species 1 if R₀₂ > R₀₁, since this implies that P₂ > P₁. This dynamics of competitive exclusion, based on the relative importance of the R₀, has been well described in other biological systems [65,66, for example].

Fig. 1b describes the second scenario, i.e., coexistence. This scenario is possible if species i affects the establishment of its own species more than the other species j. In the model, this implies that η₁₁ > η₁₂ and η₂₂ > η₂₁. In the phase plane analysis, coexistence is achieved when P₁ > P₂ and P₄ > P₃.

However, as for the one species model, heterogeneity in house exposure to infestation should be taken into account. Following Bottomley et al. [52], the model can be slightly modified and therefore, the variation with time of the mean number of colonies m_ci (i = 1,2) is:

As in the single species model, the Eq. (7) describing the colony burden of species 1 in one house (and similarly for species 2) now also depends on the covariance between exposure to house infestation and the burden of colonies in houses (σ_ΓiCj).

Coexistence appears when mutual invadability is possible, i.e., when each species can invade an equilibrium where only the other species is present [67]. Conditions of coexistence can therefore be derived from the model as follows: considering equilibrium e₁ and e₂ of each species, mutual invadability implies that e₁ and e₂ must be unstable to permit invasion of one species by the other. At e₂, if few “infective females” of species 1 start (few) colonies and as such disturb the equilibrium, invasion will be successful if after some time the mean number of colonies of species 1 (m_C1) is growing. In a mathematical sense, dm_C1/dt (Eq. (7)) > 0. At the beginning of the process, the number of colonies of the invading species 1 is small as compared to species 2 and therefore, from Eq. (7), the rate of increase of species 1 is:

When r₁ > 0, equilibrium (2) is unstable, enabling the establishment of species 1. In the same way, r₂ determines the stability of equilibrium 1. The effective reproductive numbers of species 1 and 2 at equilibrium e₂ and e₁ respectively are: R₁ = R₀₁ r₁ and R₂ = R₀₂ r₂, where R_0i = (δ′_i/μ′_c) (i = 1,2). Therefore, the condition for coexistence is R₁ > 1 and R₂ > 1.

The two effective reproductive numbers R₁ and R₂ can be reformulated taking into account that at e_i (i = 1, 2), the maximum possible value of the mean number of colonies for species i is attained when v_Гi (the coefficient of variation for the rate of exposure of species i) is zero. Then, this maximum value is mc_imax = (1/η_ii) (1 – R_0i⁻¹). Now, introducing ρ = σ_{Γ1 Γ2/}σ_Γ1 σ_Γ2 (the correlation between house susceptibility to species 1 and house susceptibility to species 2) and developing r₁ and r₂, the effective reproductive numbers R₁ and R₂ at equilibrium e₂ and e₁, respectively, can be expressed as:

These effective reproductive numbers appear to be the product of the basic reproductive number (R₀) and of a term representing the probability of establishment. This probability depends on the burden of resident triatomines in a house (m_cmax – m_c) and also on the covariance between the burden of the resident species and the susceptibility of the house to the invading species.

These equations indicate that at each respective equilibrium, the effective reproductive numbers are a decreasing function of ρ (the correlation between house susceptibility to species 1 and house susceptibility to species 2). Therefore, intuitively, decreasing the correlation between exposure rates will encourage coexistence.

When heterogeneity is the same for both species (i.e., ν_Г1 = ν_Г2 = ν_Г), then R_is are increasing functions of v_Г for ρ < 1 and, therefore, increasing the heterogeneity in house exposure to both species of triatome will facilitate coexistence. From these equations, it may also be the case that if ν_Г1 >> ν_Г2 (i.e., heterogeneity in house exposure to species 1 is high compared to heterogeneity in house exposure to the other species), then R₁ will be small because ν_Г1/ν_Г2 >> 1, and it will be very difficult for species 1 to establish itself. Therefore, coexistence is unlikely to succeed when the degree of heterogeneity in house exposure greatly differs between the two species of triatomine.

The model captures the observation that the mean number of triatomine colonies in a house is dependent on the basic reproduction number R₀ of colonies, the strength of density-dependence processes (which govern, for example, access to blood sources and/or shelter for triatomes), and the degree of heterogeneity (i.e., the degree of exposure of houses to infestation by triatomines), through which the aggregative distribution of triatomine colonies within houses can be mimicked. Introducing mechanisms that can generate aggregation is of crucial importance in T. infestans dynamics, because field observations indicate clearly the existence of such a phenomenon [2,5,54,61,68].

The model reveals that if density-dependence processes are strong, it will be difficult for a new colony to start. Triatomine dynamics are dependent on the R₀ value, which is:

• • directly proportionate to the rate at which a colony produces new “fertilized females” and the probability of survival of “fertilized females” of starting a colony, and;
• • inversely proportionate to the per capita rate of disappearance of a colony.

Thus, demographic parameters are of importance, in particular, as one might expect intuitively, fecundity and survival [69,70]. However, these parameters are offset by density-dependence processes. In the model, the mean number of colonies decreases as heterogeneity in exposure increases. Heterogeneity of exposure can be understood in terms of an infested peridomicile or an infested nearby house or structure which can be a source of infestation for a house. Indeed, as heterogeneity increases, there are more houses with very high and very low rates of exposure. However, since the number of colonies in a house is limited (as assumed in the model), houses with very high rates of exposure cannot have a corresponding high number of colonies. Thus, in total, the mean number of colonies decreases as heterogeneity increases. In the field, T. infestans exhibit such an aggregative distribution with highly infested houses and low infested ones, but always with a limited maximum number of colonies per house. Moreover, when the peridomestic environment is infested (i.e., high exposure rate for the house), and if the “immunity” of the house is low (i.e. if triatomines can find shelter and accessible blood sources inside) [62], the house also hosts domestic populations of triatomines.

A partial concluding remark would be that for a species to invade and colonize an ecotope, the model points out the predominant role of intuitive parameters such as the survival of “infective females”, the survival of a colony, a good contact rate (exposure of that ecotope to triatomines) and a good fecundity of colonies producing “infective females”.

In the model, competitive exclusion or coexistence can be described. They depend on the R₀ values of both species colonies and the values of the γ_ii and γ_ij, i.e., the “effects” of each species on the establishment of other individuals of its own species and the establishment of individuals of the other species. Generally, the species with the highest basic reproduction number (R₀ for colonies) will exclude the other, but coexistence may also occur in particular conditions (see below). Therefore, because R₀s are directly dependent on fecundity and mortality rates, it would be interesting to compare these life traits between the two species.

Coexistence may occur if one species affects the establishment of its own species more than the other species. Indeed, the phase plane analysis of Fig. 1b demonstrated that the two species will coexist if P₁ > P₂ and P₄ > P₃ [i.e., γ₁₁ (1 – R₀₂⁻¹) > γ₁₂ (1 – R₀₁⁻¹) and γ₂₂ (1 – R₀₁⁻¹) > γ₂₁ (1 – R₀₂⁻¹)]. The developed formulas indicate that for R₀₁ >> 1 and R₀₂ >> 1, coexistence will exist when interspecific effects are weaker than intra-specific ones. Interspecific effects might be weaker than intra-specific ones if, for example, there is site segregation between the two species within human dwellings. Such situations have been observed for the peridomestic distribution of Triatoma garciabesi and Triatoma guasayana in northwest Argentina [21]. The two species have slightly different ecotopes in the peridomicile and, therefore, coexistence situations have been described. At a larger geographic scale, it has also been described for Rhodnius neglectus and R. nasutus [71] and for Triatoma sordida and T. garciabesi [72]. As for T. infestans however, it does not seem that R₀ values are >> 1, and there is no evidence that one group (domestic for example) may affect its proper establishment more than that of the other group if this latter enters the same site (house, in this case). Even if coexistence is possible, a species may still exclude another one if the R₀ of the competitively inferior species is sufficiently close to 1. In that case, for various situations where the R₀s of the triatomine species are low and species are in low densities, the competitive advantage of one species on the other will vary as R₀ values change between situations. Then, it might appear that in one site, one competitor will be present exclusively, while in another site, it will be the other one.

The model clearly shows that the colonization of houses by wild T. infestans populations is unlikely. Personal field observations by the mark–release–recapture technique indicate that T. infestans individuals do not move easily between wild and domestic ecotopes, even if they are close (50 m for example). Therefore, mutual exposure (i.e., exposure of houses to wild triatomines and exposure of sylvatic environment to domestic triatomines) seems to be close to zero while exposure of each type of ecotope by each respective “species” is high and homogenous. If so, the model indicates that (see Eqs. (11) and (12)) it is hard for a wild population of T. infestans for which host heterogeneity in exposure is high (i.e., a low value of exposure for houses and a high value for a wild environment) to invade the already established domestic population where host heterogeneity is relatively low. For example, for the wild species invading a house, Eq. (11) indicates that its reproductive number will be small because ν_Γ1/ν_Γ2 >> 1. Invasion will therefore be very hard. Coexistence of wild and domestic populations of T. infestans is possible, but very unlikely to occur. Indeed, the model demonstrates that coexistence may occur if one species affects the establishment of its proper species more than the other species when both share the same site. For wild/domestic populations, this is unlikely as each group seems to have slightly different ecological needs and thus occupies different ecological niches. Genetics studies tend to distinguish two groups which do not mix [40]. However, a recent study indicated that gene flows can exist between the two groups [41], but, as stated, may represent invasion of sylvatic habitats by domestic populations finding here refuges after insecticide sprayings of houses. In that case, there is no competition, but invasion of an empty niche followed by intra-specific movements of dispersion as one might expect. Conversely, if domestic individuals have survived vector control actions and still occupy their houses, the model tends to minimize the possible role of wild population in invading such a “not totally empty” niche.

Wild populations of T. infestans are mostly found close to the domestic/peridomestic environment (Brenière, comm. pers.) and the definition of “wild” is indeed difficult. Blood meal sources may help in defining a “wild” population which then feed on wild animals. However, at present, the wild foci of T. infestans discovered are mostly under human influence, close to human settlings, except maybe in the Chaco region of Bolivia. For most of these foci, then, it seems that the domestic–peridomestic populations of T. infestans have invaded the wild environment and our model can account for such behaviour. If so, a new category of triatomine environment could be proposed: the “PARA-domestic” environment, in addition to the “domestic”, “peridomestic” and “sylvatic” terms yet in use. The para-domestic environment is then under the influence of the domestic–peridomestic environments, with some different characteristics (such as blood sources for triatomines that may be mainly of wild origin) and may permit population to undergo some level of differentiation.

When two groups are considered for coexistence or competition, our model assumption is that they do not interbreed (i.e., are real species). Below this working hypothesis, our model indicates that coexistence of wild and domestic populations of T. infestans is still possible, but very hard to achieve. If so, the colonization of domestic habitats by wild populations is unlikely. However, if further studies indicated that gene flows amongst wild and domestic T. infestans are significant (i.e., that they interbreed and mix), other models should be used to describe their interactions. Models taking into account source-sink dynamics [73] or better, ecological traps [74,75] should then be of great value and would likely explain the present heterogeneities observed amongst wild and close domestic T. infestans populations, such as genes and phenotypic differences or insecticide resistance. If gene flows exist amongst the two groups, then the population is single and the main question would be to identify likely source-sink dynamics where wild groups cannot persist without close domestic ones (and then without threats for vector control) or to identify mechanisms that enable the persistence of wild populations jeopardizing vector control. If so, distant wild populations without any anthropogenic influence should also be found.