Cross-immunity in a Multi-strain Cholera Model

By Leah LeJeune

In mid-1800s London, physician John Snow developed one of the first epidemiological models in order to combat a cholera outbreak. By mapping the locations of reported cases, he observed that the majority of infected people were located around a specific community water pump. Once the pump handle was removed and the water source was inaccessible, the disease swiftly died out. Since this initial breakthrough, the world has continued to make much progress in the understanding and control of cholera. However, the disease is still present in many countries—recent outbreaks occurred in Haiti, Zimbabwe, and Bangladesh—and requires further study to limit its spread. 

Snow’s historic discovery highlights the benefits of data analysis and disease modeling and identifies one of cholera’s most interesting characteristics: its dependency on a water source. The bacteria Vibrio cholerae can survive, compete, and proliferate in an aquatic environment, meaning that models must consider water sources to accurately capture transmission dynamics. Cholera can also be transmitted from person to person (although the prevalence of this scenario relative to environment-based transmission seems quite low). Because the disease causes severe dehydration via diarrhea and vomiting, contact with waste from an infectious individual transmits the bacteria. Places with poor sanitation and infrastructure—like some developing nations or areas that are damaged by natural disasters—thus tend to experience cholera outbreaks most frequently.

Before creating a model, we must consider another aspect of cholera biology: the presence of distinct serotypes. Serotypes are variants of bacteria that interact with the immune system in different ways. During recent cholera outbreaks, researchers have observed the competition of serotypes Inaba and Ogawa within the infectious population. Both serotypes induce a certain amount of self-immunity within a host, but they also induce a degree of cross-immunity — that is, immunity to each other. In a stroke of generosity for all mathematical modelers, serotypes have identical life histories and thus allow the assumption of identical parameters (e.g., transmission and recovery rates) between the two strains. 

Figure 1. Compartmental diagram of a two-strain cholera model with an aquatic reservoir. Figure courtesy of Leah LeJeune.

Snow created an incredibly effective plan for control measures by simply mapping the locations of infectious individuals. Yet given our present knowledge of disease biology—including multiple modes of transmission and the interplay of Inaba and Ogawa—current cholera models require significantly more consideration, especially when the environmental reservoir is something like a river (rather than a pump with a removable handle). 

We use ordinary differential equations (ODEs) to construct our model. ODE disease models often take the form of susceptible-infected-recovered (SIR) models, wherein individuals move between \(S\), \(I\), and \(R\) compartments at given rates. We adapt this model structure for our purposes by making two major changes. First, we consider an additional compartment that represents the pathogen \((P)\) population in the environment, which acts on the other SIR compartments. And second, we account for the influence of two strains of the bacteria in our system, rather than just one. Some existing studies have explored cholera with a single-strain SIRP model [2], while others have utilized an SIR model for multiple strains of the disease. Our model looks at both components simultaneously (see Figure 1).

For \(j,k = 1,2\), the compartments represent the proportions of susceptible hosts \((S)\), strain \(j\) infectious individuals \((I_j)\), strain \(j\) recovered individuals \((R_j)\), strain \(j\) infectious hosts with a previous strain \(k\) infection \((I_{kj})\), and the quantity of the strain \(j\) pathogen in the environment \((P_j)\). The parameters \(\beta\) and \(\beta_2\) respectively represent host-to-host transmission after a first and secondary infection. Likewise, \(\alpha\) and \(\alpha_2\) are the shedding rates—the rates at which bacteria enter the aquatic reservoir from an infectious individual—for the respective infections. Finally, \(\delta\) is the environment-to-host transmission coefficient, \(\mu\) is the population turnover rate (i.e., it represents both the birth and death rate), \(\nu\) is the rate of recovery, \(r\) represents pathogen decay in the environment, and \(\gamma\) represents cross-immunity. We exclude a growth component for the bacteria in the environment because this process often occurs very slowly relative to the infection timescale. We also normalize the model, so the maximum host population size is fixed at \(1\). The double lines between strains in Figure 1 designate the multiple modes of transmission, the blue compartments represent the pathogen populations, and the dashed lines denote the populations’ effects on the model. The blue and red crossing lines respectively distinguish secondary infections from strain 1 and strain 2.

We make one final simplification to our initial model in Figure 1 by assuming that transmission rates and shedding rates from first and secondary infections are the same. This assumption allows us to define a new compartment, \(y_j = I_j + I_{kj}\), to represent the proportion of strain \(j\) infectious individuals. Doing so consolidates the model (which would be nine equations) into a friendlier system:

\[\frac{dS}{dt} = -\beta(y_1 + y_2)S - \delta(P_1 + P_2)S + \mu(1-S)\]

\[\frac{dy_j}{dt} = (\beta y_j + \delta P_j)S + (1 - \gamma)(\beta y_j + \delta P_j)R_k - (\nu + \mu)y_j\]

\[\frac{dR_j}{dt} = \nu y_j - (1 - \gamma)(\beta y_k + \delta P_k)R_j - \mu R_j\]

\[\frac{dP_j}{dt} = \alpha y_j - r P_j.\]

Note that we can likely extend the results of this simplified system to a system with more complexities (a common and efficient strategy). 

Figure 2. Uniform persistence of the disease in the environment. The function \(\rho\) represents the maximum of infectious compartments \(y_1\) and \(y_2\). The axes represent the values of \(y_1\) and \(y_2\). The origin point is the disease-free equilibrium, where neither strain is present. The \(y_j\) axis represents where the strain \(j\) infectious population is nonzero and the strain \(k\) infectious population is nonexistent. Uniform persistence means that \(\rho\) will stay near these boundaries (i.e., within \(\epsilon\) of the axes) for only a finite period of time before it is drawn to the interior set, where both strains coexist. Figure courtesy of [5].

\[R_0 = \frac{\beta}{\mu + \nu} + \frac{\delta}{(\mu + \nu)}\left(\frac{\alpha}{r}\right).\]

When broken down into these two terms, \(R_0\)'s dependence on multiple modes of transmission is clear. The first term indicates the contribution of host-to-host transmission, and the second indicates environment-to-host transmission with a magnitude that is regulated by shedding rate \(\alpha\) and decay rate \(r\). Using standard analysis, we can demonstrate that the disease will completely die out after a period of time if \(R_0 < 1\), regardless of the number of initial infections. 

But when \(R_0 > 1\), the dynamics are less clear. We can certainly establish the existence of a unique coexistence equilibrium (where both strains achieve a positive steady state). Numerically, many simulations exhibit convergence to identical steady states for the two infectious compartments \(y_1\) and \(y_2\). And analytically, global results from persistence theory show that the disease will move away from the boundary—defined when one or both serotypes are not present—and towards an attracting set in the interior, which could either be an equilibrium or a periodic solution (see Figure 2). 

One of the most interesting biological options for coexistence behavior corresponds to a periodic solution for our system. For instance, a 2006 study explored cholera serotype dynamics over a period of 40 years in Matlab, Bangladesh [4]. The two serotypes alternated between dominance and subservience, but neither strain reached extinction or a steady state (see Figure 3). This begs the following two-part question: What drives the dominance switch, and can we replicate this behavior in our model? Prediction of when and why the switch occurs would allow decision-makers to take efficient control measures, such as vaccinating at the optimal time and against the correct serotype.

Figure 3. Proportion of symptomatic cases of Inaba (blue) and Ogawa (red) serotypes from January 1983 to October 2005. Figure courtesy of [4].

We do have some idea as to what may drive the switch. A previous study used a model that was similar to ours and contained a cross-immunity component but excluded a pathogen compartment [1, 3, 4]. This model exhibited periodic solutions with anti-phase (i.e., dominance cycling) oscillations between the two strains that were driven by the cross-immunity component (\(\gamma\) in our model). Figure 4 depicts changes in our model’s oscillatory behavior for both the infectious populations and pathogen populations when \(\gamma\) varies from \(0\) (full cross-immunity) to \(1\) (no cross-immunity). As \(\gamma\) approaches \(1\), oscillations are barely present — both infectious populations quickly converge to the same equilibrium value. With no cross-immunity and identical rates of infection and recovery, the two strains behave identically in the system after a sufficient period of time to account for the difference in initial infections. This trend is common in our simulations, although we have not yet characterized the specifics of how—and by how much—\(\gamma\) controls this process.

Figure 4. Here, \(\beta = 312/365.2422\), \(\delta = 0.1\), \(\mu = 1/(60*365.2422)\), \(\nu = 52/365.2422\), \(\gamma = 0.2\), \(\alpha = 1/5\), and \(r = 0.9\). As \(\gamma\) ranges from \(0\) to \(1\), anti-phase oscillatory behavior is most present around \(\gamma = 0.5\) (i.e., when the protection to infection from the other strain is about 50 percent) and least present as \(\gamma\) approaches \(1\). Notice that our equilibria values for the pathogen compartments can be written as \(P_j = (\alpha/r)y_j\). Since the pathogen population is always proportional to the infectious host population of the same strain, the graphs will exhibit the same general behavior after enough time. Figure courtesy of Leah LeJeune.

In each instance of Figure 4, every population eventually converges to an equilibrium value; no simulations exhibit the undamped, anti-phase oscillatory behavior seen in Figure 3. To address this disparity, we adjust the environment-to-host transmission coefficient \(\delta\) from a constant parameter to a periodic function of time \(\delta_t\), since cholera—a water-borne bacteria—spreads rampantly during rainy seasons. We define this periodic function as

\[\delta_t = \delta(1 + \delta_0 \sin(t 2\pi/365.2422)).\]

Though this modification further complicates the model, it introduces the effect of seasonality in the context of disease spread. Simulations of this forced, nonautonomous model do indeed exhibit oscillatory behavior that is similar to that in Bangladesh [4]. 

In the future, we will continue to analyze this forced system and attempt to understand cross-immunity’s effects on solution behavior. Incorporating all of these considerations will advance our knowledge of cholera’s spread.

Acknowledgments: The author would like to thank her advisor, Cameron Browne, for his support and guidance of this research. This material is based upon work that is supported by the National Science Foundation under grant DMS-1815095.

References
[1] Adams, B., & Boots, M. (2007). The influence of immune cross-reaction on phase structure in resonant solutions of a multi-strain seasonal SIR model. J. Theor. Biol., 248(1), 202-211. 
[2] Bani-Yaghoub, M., Gautam, R., Shuai, Z., van den Driessche, P., & Ivanek, R. (2012). Reproduction numbers for infections with free-living pathogens growing in the environment. J. Biol. Dyn., 6(2), 923-940. 
[3] Kamo, M., & Sasaki, A. (2002). The effect of cross-immunity and seasonal forcing in a multi-strain epidemic model. Physica D, 165(3-4), 228-241. 
[4] Koelle, K., Pascual, M., & Yunnus, M. (2006). Serotype cycles in cholera dynamics. Proc. R. Soc. B, 273(1603), 2879-2886. 
[5] Smith, H. (2015). Persistence theory in population dynamics. Retrieved from https://math.la.asu.edu/~halsmith/PersistenceTutorial.pdf.

Leah LeJeune is a Ph.D. candidate at the University of Louisiana at Lafayette. Her research interests include mathematical modeling with applications to epidemiology and ecology.