The four-patch cofactor model presented here builds upon and inherits all of the parameter names (*μ, α, β, v, ϵ* Table 1) from the two-patch SI model of Anderson et. al. (1988) [18] describing the impact of HIV on overall host population growth in developing countries. That model tracks the susceptible (*X*(*t*)) and infected (*I*(*t*)) subpopulations of the host population (*N*(*t*)) through time (Figure 1a), and is completely tractable from an analytic viewpoint. Incorporation of the cofactor into this basic model framework requires two new population patches (Figure 1b). The susceptible subpopulation is now divided into *S*(*t*) - susceptible to the primary disease and to the cofactor, and *C*(*t*) - susceptible to the primary disease with the cofactor already present. The infected subpopulation is now divided into *Y*(*t*) - infected by the primary disease but not by the cofactor, and *Z*(*t*) - infected by the primary disease and possessing the cofactor. The expanded model includes parameters controlling differences in the infectivity of and susceptibility to the primary infection caused by the cofactor (non-negative parameters *δ*
_{2} and *δ*
_{1}, respectively, Table 1), as well as parameters for the rates of transmission and loss of the cofactor (*ζ* and *γ*, respectively, Table 1).

Whereas the two-patch model has a single transition (Figure 1a), the four-patch cofactor model includes eight new transitions among subpopulations (Figure 1b). These new transitions arise because the cofactor can be cured, introducing bidirectionality; because the two diseases can be transmitted either singly or together; and because of the greater number of players in the system. Therefore, while the two-patch cofactor-free model admits an exact solution, the additional non-linearity in the four-patch cofactor model makes mathematical analysis of the latter much more difficult. Complete results pertaining to the circumstances under which the primary disease can successfully invade a population in the four-patch cofactor model have been obtained, and are described in the next subsection. We then proceed to apply the four-patch cofactor model through simulations of a specific disease/cofactor system.

### Primary disease invasion in the analytical model

#### Invading the disease-free/cofactor-free state

Suppose that the primary disease agent attempts to invade a population which is entirely free of both the primary disease and the cofactor. In the two-patch model, this invasion occurs through the introduction of one individual who is infected by the primary disease (from population *I*, Figure 1a). In that model, the primary disease successfully invades when the number of new *I*s (susceptible individuals infected) produced by an invading *I* individual (i.e. the basic reproductive number) is greater than one. In the four-patch cofactor model, the primary disease may invade through the introduction of an individual who is infected either by the primary disease alone (from subpopulation *Y* , Figure 1b) or by both the primary disease and the cofactor (subpopulation *Z*, Figure 1b). Analysis of the Jacobian matrix of the linearization of the subsystem for *Y* and *Z* [25] shows that because an invading *Y* cannot produce *Z*s in the absence of *C*s, the primary disease successfully invades when either the number of new *Y*s produced from an invading *Y* or the number of new *Z*s produced from an invading *Z* is greater than one (Appendix, first subsection). This leads to an *R*
_{0}(0) (the basic reproductive number of the population which starts with no individuals possessing the cofactor) value of max{*β*, λ_{2}ζ - γ}/(*μ+α*).

Since patches *Z* and *C* in the four-patch cofactor model are initially empty, for an invasion into a disease-free and cofactor-free population the basic reproductive number of the two-patch model is the same as the number of new *Y*s produced from an invading *Y* in the four-patch cofactor model. Therefore, it is never harder for the primary disease to invade the disease-free/cofactor-free population in the four-patch cofactor model when compared with the two-patch model. In other words, our *R*
_{0}(0) is at least as large as the basic reproductive number of the two-patch Anderson et. al. (1988) model. The ability of the primary disease to invade the disease-free/cofactor-free population is facilitated by the cofactor when the basic reproductive number of the two-patch model is less than one, and the number of *Z*s produced from an invading *Z* is greater than one. This occurs when (λ_{2}
*ζ* - γ)/(*μ* + *α*) > 1. Keeping in mind that *λ*
_{2}
*ζ* is the rate at which *S*s convert to *Z*s, we see that except in the event that the disease and the cofactor are transmitted together substantially more often than singly during a single interaction (or λ_{2}
*ζ > β*), the cofactor will have no effect on the ability of the infection to invade the disease-free/cofactor-free population.

It is also possible for the primary disease to invade an exponentially growing (*ν* > *μ*) disease-free (*Y* + *Z* = 0) population even though the proportion of infected individuals approaches zero. By analyzing the Jacobian matrix of the proportions *Y*/*N* and *Z*/*N* of individuals with primary infection only and of individuals with primary infection and the cofactor near the disease and cofactor free state, we find that the proportion of individuals with the primary disease always approaches zero unless *R*
_{0}(0) > (*ν+α*)/(*μ+α*) >1. This means that the percentage of the population with the primary infection will only become bounded away from zero in the long run when the disease grows fast enough to outpace not only the general death rate, *μ*, but also the population growth rate, *ν*. [18, 26]

#### Invading a cofactor endemic population

The primary disease could be introduced to a population in which a cofactor is already established. This would occur if an emerging disease arises in a population in which an established behavior, condition or pathogen exists. For HIV, this scenario best represents the conditions under which the primary disease was introduced into most human populations. In this case it is still true that if either the number of new *Y*s produced from an invading *Y* or the number of new *Z*s produced from an invading *Z* is greater than one, then the primary disease will successfully invade the population (Appendix, second subsection). Since, in the presence of *C*s, the *Z*s can produce *Y*s and the *Y*s can produce *Z*s, the off-diagonal terms in the Jacobian matrix of the linearization of the *Y*-*Z* subsystem are now both non-zero. This produces a scenario in which invasion may occur even when both of these numbers (*Y*s from a *Y* and *Z*s from a *Z*) are smaller than one. The condition governing this type of invasion is complex and not easily interpreted in biological terms (Appendix, second subsection).

Despite this complication, we analytically verified (Appendix, third subsection) that whenever *R*
_{0}(0) > 1, then the basic reproductive number of the population which starts out with the stable proportion *c** of individuals possessing the cofactor, *R*
_{0}(*c**), is also greater than one, meaning that the disease-free/cofactor-*endemic* population can also be invaded by the primary disease. We show that because *R*
_{0}(*c**) ≥ *R*
_{0}(0), it is never the case that the disease-free/cofactor-free population can be invaded, while the disease-free/cofactor-endemic population cannot. This supports the idea that the cofactor cannot exhibit a preventative effect on disease invasion in this model. The demonstration of this inequality is somewhat subtle, since the starting equilibrium state affects the number of new *Y*s that can be produced by an invading *Y*. During an invasion of the disease-free/cofactor-endemic population, this number is smaller than in an invasion of the disease-free/cofactor-free population due to the smaller available proportion of cofactor-free susceptibles (*S*s) in the population.

When the impact of the cofactor (as measured by *δ*
_{1} and *δ*
_{2}) is small, whenever the disease-free/cofactor-free population resists invasion the disease-free/cofactor-endemic population will also resist invasion by the primary infection. This means that when the basic reproductive number of the two-patch model is smaller than one, there will be a domain of values of *δ*
_{1} and *δ*
_{2} for which *R*
_{0}(*c**) < 1 as well. However, if the impact of the cofactor is sufficiently large, *R*
_{0}(*c**) will be larger than one and the primary disease will successfully invade the disease-free/cofactor-endemic population even when it would fail to invade the disease-free/cofactor-free population. Specifically, infection of cofactor-possessing individuals (C to Z transitions) can be quite common in the disease-free/cofactor-endemic case. This group has the potential to have the largest transmission rate in the system, because both the risk of infection in the *C* individual and the infectivity of the *Z* individual have been increased by the cofactor condition. As a result each invading *Z* may produce more new *Z*s than an invading *Y* can produce *Y*s, and will in fact produce more than one new *Z* when the values of *δ*
_{1} and/or *δ*
_{2} are sufficiently large. The precise criteria for this scenario appears in the second subsection of the Appendix.

In typical ecological models, having once come to understand the boundary equilibria (in which one subpopulation is absent), the analysis would proceed via persistence theory [27, 28] to consider the existence of a coexistence equilibrium state at which all subpopulations are present. However, the four-patch cofactor model is somewhat unusual among ecological models because it includes so many intrinsic interactions among subpopulations. As a result persistence theory cannot be applied, and it is impossible to explicitly solve the general system or to precisely determine the long term behavior of the solutions. It is even difficult to locate (in terms of the parameters) an equilibrium concentration state in which the cofactor-only (*C*s), the primary disease-only (*Y*s), and the cofactor with primary disease subpopulations (*Z*s), are all present. This situation likely results from the interactive nature of the *Y* and *Z* populations - that is, *Z*s create *Y*s from interacting with *S*s and *Y*s create *Z*s from interacting with *C*s. It is the non-linear terms in the model which result from these interactions that prevent the model from being solved in closed form for the positive equilibrium state. (Note that the model does admit a positive equilibrium, but that it cannot be found analytically.)

Timeseries simulations of the model over a variety of parameter ranges have provided no hint that the model can exhibit either oscillatory or chaotic behavior in its component subpopulations. It seems that whenever the primary disease can successfully invade the population in the four-patch cofactor model, then the proportions of each of the subpopulations will tend toward a positive equilibrium value. This means that even though it is not possible to precisely predict the long term behavior of the system in general terms, the model will still have great predictive value through simulations with a given set of parameter values.

The comparison of the four-patch cofactor model with the two-patch model highlights the important analytical changes that occur in this system as a result of the addition of a single player: a linear and explicitly solvable system becomes non-linear and unsolvable. This qualitative change in the system as a result of the presence of a cofactor indicates that in populations where a potential cofactor co-occurs with a primary disease of interest, management approaches based on parallel and independent models can be expected to considerably misrepresent potential effects of interventions to the spread of primary infection. In particular, in populations that are completely naive to both the primary disease and the cofactor and in populations in which a cofactor is established at the time of invasion by a primary disease, a simple, cofactor-free model like that of Anderson et. al. [18] can substantially overestimate the invasion force of the primary infection by attributing the entire rate of spread of the infection to the infection itself, rather than to the combined force of the primary infection and its assisting cofactor.

### Simulation of HIV with *Schistosoma haematobium*cofactor

We further investigate the effect of cofactors on the spread and persistence of a primary disease and the potential impact of intervention policies targeting cofactors rather than the primary disease through simulations. In these simulations, estimates of transmission contact rates for the primary disease (*β*) and the cofactor (*ζ*) are based on those of HIV and *Schistosoma haematobium*, the causal agent of genital schistosomiasis. We chose the cofactor *S*. *haematobium* for this study because it is endemic in many areas with high HIV prevalence [9, 29], has been identified as a cofactor for HIV [30–32], and is easily and inexpensively treatable [9, 33].

In Anderson et. al. (1988) [[18], and refs within], many demographic parameters are chosen directly from measurement data of birth and death rates in the sub-Saharan region which comprises the focus of that study. They also estimated what we call the total transmission rate at the start of the HIV epidemic (*β*
_{
T
}= *β*, Figure 1a) in such a way that the output from their model would match HIV prevalence data from those same locales. By contrast, the four-patch cofactor model makes use of a direct transmission rate (*β*
_{D} = *β*, Figure 1b) from which any effect of the cofactor has been removed. For the following simulations we estimated *β*
_{D}, *β*
_{
T
}, and other HIV-specific parameters using epidemiological data and estimates from HIV-specific models [[7, 18, 21, 22] details in Methods]. Since the schistosomiasis incidence and prevalence data [23, 24] provide a broad range of estimates for schistosome transmission rates (*ζ*), we choose *ζ* in the center of the estimated range with the result that the timeseries of the total number of individuals infected by the primary disease under the four-patch cofactor model (with *β*
_{D}) agrees well with the total number of individuals infected by the primary disease under the two-patch model (with *β*
_{
T
}). All of these parameters are fixed for the following simulations (Table 2, values chosen; methods section, motivation for choices).

The quantitative results of the following simulations are based on very basic models of both HIV and schistosomiasis which suppress many significant details of both conditions. Consequently, the focus here is on the qualitative results of the simulations under the parameter choices that have been made, illustrating the possible importance of the disease/cofactor synergy on the primary disease progression within a population. In particular, the primary importance of the two-patch model of Anderson et. al. (1988) is its ability to give an early indication, based on sparse observational data during the early stages of an epidemic, of the potential demographic impact of the epidemic. With the goal of qualitative evaluation of cofactor impact during the same early stages of such an epidemic in mind, it is therefore reasonable to base both the construction of the four-patch cofactor model and the selection of parameter values for use in simulation on the structure and parameter choices made in that two-patch model.

#### Agreement of the four-patch cofactor model with two-patch model results

If the four-patch cofactor model is to be used to measure the decrease in HIV prevalence resulting from a treatment program for a cofactor condition, it should first be verified that it models the total population and the actual number of infected individuals under the conditions where the cofactor is present in the population and no recovery from the cofactor takes place (*γ* = 0) in a way which is comparable to the two-patch model. For this comparison, in the four-patch cofactor model simulation we allow individuals to occupy all four patches (*S*, *C*, *Y* and *Z*), set *β*
_{
D
}= 0.08 per year, and the other parameters as listed in Table 2. In contrast, in the two-patch solution, *β*
_{
T
}is set to 0.2 per year, with the other two-patch parameters the same as the analogous four-patch parameter values. The parameter *β*
_{
T
}is always at least as large as *β*
_{
D
}because it implicitly incorporates transitions involving *Z*s as well as those involving *Y*s. Figure 2 shows the timeseries of the total population and the number of infected individuals for both the four-patch model (in red), and the two-patch model (in black). The similarity of these timeseries indicates that the four-patch model does indeed correctly model the total population and the actual number of infected individuals under our selection of parameters as compared to the two-patch model.

#### The HIV prevalence in a cofactor-free population

In order to determine whether or not the four-patch cofactor model correctly measures the decrease in HIV prevalence resulting from a treatment program for a cofactor condition, it is necessary to estimate the expected HIV prevalence in a cofactor-free population under the two-patch model. In a truly cofactor-free population the total transmission rate *β*
_{
T
}would simply be equal to the direct transmission rate *β*
_{
D
}. Therefore, to measure the difference in prevalence of HIV between the observed cofactor-endemic and a truly cofactor-free context, we compare the solutions of the two-patch model using the total transmission rate (with *β* = *β*
_{
T
}= 0.2 as the HIV transmission rate) just as above, and then using our estimated direct transmission rate (*β* = *β*
_{D} = 0.08).

Figure 3 shows timeseries for the total population and the number of infected individuals for two two-patch model solutions. The total transmission curve uses the estimate of the total transmission rate, *β*
_{
T
}(used in the above model), as the estimate of the parameter *β*. The direct transmission curve uses the estimated direct transmission rate, *β*
_{
D
}, as the estimate of *β*. Under our choices of these parameters, the relative proportion of infected individuals is approximately stable (or near the positive equilibrium) after simulating a time period of 70 years. The contrast between the curves in Figure 3 indicates that cofactor conditions can greatly exacerbate the negative impact of the HIV epidemic on a population. The simulated proportion of the population with HIV is near 75% by 70 years under the total transmission rate while under the direct transmission rate the number infected increases at a much slower rate, and the concentration of the infection decreases to near zero by 70 years (Figure 4, black). In addition, while the total transmission rate results in population decline after only 25 years the direct transmission rate allows the population to continue to grow. A program for treating cofactor conditions within a population where the true direct transmission rate is lower than the total transmission rate, as modeled here, might therefore offer the significant benefit of reducing HIV prevalence.

#### The effect of cofactor intervention on HIV prevalence

We now investigate the predicted decrease in HIV prevalence resulting from a treatment program for a cofactor condition in the four-patch cofactor model. We expect to see four-patch timeseries results which fall between the two different two-patch curves in Figure 4 (thick and thin black), with the exact location depending on the magnitude of the treatement effort (as measured by the induced cofactor recovery rate *γ*). To test this, we simulate the four-patch cofactor model using exactly the same collection of parameters as above (Table 2, Figure 2) but now include a positive rate of recovery from the cofactor (*γ* > 0). The total population timeseries from two such simulations (*γ* = 0.025 and *γ* = 0.075) appear in Figure 4 (red), along with the two-patch direct and total transmission rate curves. The simulation with a greater intervention parameter (*γ* = 0.075) shows spread of primary disease more like the two-patch direct transmission rate curve while the simulation with the smaller intervention parameter (*γ* = 0.025) is closer to the two-patch total transmission rate curve. Therefore, the four-patch model does serve to interpolate between the different two-patch curves in a way which depends on the force of the intervention.

We also simulate the four-patch cofactor model many times over a range of positive rates of recovery from the cofactor (*γ* from 0 to 0.2) to explore the effects of intervention over a range of intensities (Figure 5), focusing on the final concentration of infected individuals after 25 years from each of the simulations and comparing it to the results from the two-patch model with both total and direct transmission rates. Under the total rate in the two-patch model and under low rates of intervention in the four-patch model the total populations are in decline (dotted portions of the curves, Figure 5). Higher rates of intervention, like the direct transmission rate from the two-patch model, result in continued population growth (solid portions of the curves, Figure 5). Notably, when the recovery rate from the cofactor is near zero, the infection concentration from the four-patch cofactor model is quite similar to the infection concentration from the two-patch model with total transmission rates, and as the recovery rate increases, the infection concentration moves closer to that observed under the direct transmission rates, just as suggested by the simulation from Figure 4. Figure 5 can be used to directly estimate the difference in the expected HIV prevalence between a two-patch model system and a four-patch system including treatment of the cofactor at a given constant rate during a 25 year time period.