Modelling Typhoid Fever Transmission in the Far North of Cameroon: Simulation-Based Risk Analysis and Environmental Impact

Kikmo Wilba Christophe^1*, Gnassiri Simon¹, Sone Enone Bertin², Batambock Samuel¹, Abanda Andre¹

¹National Advanced School of Engineering of Douala, University of Douala, Douala, Cameroon.
²Faculty of Medicine and Pharmaceutical Sciences, University of Douala, Douala, Cameroon.

DOI: 10.4236/jamp.2025.139171                

Abstract

A comprehensive mathematical framework modelling transmission dynamics of typhoid fever exists for Far North Cameroon where unsanitary conditions significantly exacerbate Salmonella Typhi spread rapidly. Analysis incorporates a deterministic model rooted in ordinary differential equations and a stochastic methodology factoring in uncertainties somewhat randomly. Dual modelling strategy highlights dominant role of water-related factors and climatic variables in shaping epidemic trajectory quite significantly over time. Seasonal disease pattern exhibits two pronounced incidence peaks in April-May and July-August corresponding respectively to drinking water scarcity periods and increased surface runoff facilitating pathogen dissemination. Advanced Bayesian techniques particularly Markov Chain Monte Carlo algorithm and variational inference enable estimation of key epidemiological parameters accurately with Markov processes. Analysis reveals that the basic reproduction number exceeds epidemic threshold during critical periods remarkably often under certain conditions. Simulations of multiple scenarios pretty effectively demonstrate efficacy of targeted control measures like vaccination programs and public awareness crusades nationwide. Such interventions drastically curtail transmission rates and stabilise epidemic trends somewhat effectively meanwhile. Findings contribute valuable insights into epidemiological dynamics of typhoid fever amidst climate variability and offer a robust foundation for public health risk management strategies. Strategic integration of real-time epidemiological data and water-quality surveillance systems holds great promise for enhancing sustainable control of this nasty waterborne disease.

Keywords

Typhoid Fever, Compartmental Model, Water Transmission, Control Strategies

1. Introduction

Typhoid fever remains a pressing public health issue in low-income countries, particularly where sanitation and hygiene conditions are poor. In the far north of Cameroon, the situation is especially concerning as the disease has reached alarming levels . This region, characterized by harsh Sahelian climate conditions, faces a complex mix of challenges: water scarcity, seasonal flooding, and inadequate sanitation infrastructure. These factors, combined with the lack of water treatment systems, increase health risks and lead to frequent contamination of water sources, which fuels community transmission of the disease. Departments such as Mayo-Sava, Logone-et-Chari, and Mayo-Tsanaga are repeatedly hit by typhoid outbreaks. The persistence of these epidemics is largely due to poor wastewater management, communities living in flood-prone areas, limited access to clean drinking water, and a general lack of public awareness regarding prevention [3]-. Together, these factors create a serious health crisis where human behavior and environmental conditions interact in ways that are difficult to manage and understand. One powerful tool for analyzing, predicting, and controlling infectious diseases is mathematical modeling. Compartmental models such as the classic SIR (Susceptible-Infectious-Recovered) and their extensions (SEIR/SEIRS) have been widely used to simulate the spread of diseases. These models rely on key epidemiological parameters like transmission rates, [7] incubation periods, recovery rates, and mortality rates, and are typically constructed using Ordinary Differential Equations (ODEs). Studying these equations allows researchers to evaluate the stability of an outbreak and identify critical thresholds, like the basic reproduction number, $R_0$ . Given the uncertainties associated with natural environments, stochastic models are a useful complement to deterministic ones. These models account for random events such as unpredictable human behavior, weather patterns, and other variables that influence disease spread [8]-[10]. Furthermore, hybrid models that combine ODEs with partial differential equations have been developed to explore how pathogens move through aquatic environments. These models factor in topography, the presence of contaminated sources, and environmental changes over time. Much of the existing research on typhoid fever risk has focused on urban environments, particularly in countries with strong health surveillance systems. Studies by [12] highlight the value of including rainfall and temperature data in models of waterborne disease transmission. Similarly, [13] [14] showed that typhoid epidemics in tropical regions often follow seasonal rainfall patterns. Although SEIR models have been used to study typhoid in several African countries, most research has concentrated on urban or semi-urban areas where data is more readily available. There is still a gap when it comes to integrating climate data—like rainfall, temperature, and flooding—into models designed for semi-arid rural settings, where data is often sparse or unreliable [15]. Moreover, many models are too generic and fail to capture the specific realities of the Sahel region. The lack of real-time health data in rural northern Cameroon further limits the accuracy of these models. Additionally, little attention has been given to models that address both extreme climatic factors and the role of contaminated water in disease transmission in high-risk areas. Few studies have tried to measure the combined impact of seasonal flooding, water stress, and extreme temperatures on the survival of Salmonella Typhi in water systems. And while both deterministic and stochastic modeling approaches exist, they are rarely integrated into a single framework, despite the potential benefits for predicting outbreaks and guiding public health responses. This study proposes a comprehensive and context-sensitive modeling approach tailored to the unique epidemiological, climatic, and public health conditions of northern Cameroon. The model introduced in this research stands out for several reasons:

• It is deterministic and uses ordinary differential equations to describe how the disease evolves over time.

• It also includes a stochastic component to capture the uncertainties in both epidemiological and climate-related data.

• It explicitly incorporates environmental factors like extreme rainfall, seasonal flooding, temperature, and water contamination.

• It models both waterborne and human-to-human transmission, unlike conventional models that often ignore environmental pathways.

The goals of this research are twofold: first, to analyze and quantify how typhoid fever spreads in a high-risk Sahelian setting, and second, to develop a decision-making tool that helps local authorities improve their health responses and build epidemic resilience. By using a hybrid modeling approach, the study aims to generate actionable recommendations for enhancing healthcare coverage and guiding public health policies rooted in scientific evidence. Ultimately, this work addresses a major gap in existing research by providing a reliable, localized model that can adapt to the realities of northern Cameroon. It contributes to a broader effort to strengthen universal health coverage in the country, aligning with both national priorities and global health goals.

2. Deterministic Model Formulation

Typhoid fever also known as abdominal typhus is a bacterial disease caused by infection that Pierre Bretonneau described in eighteen eighteen. Bacteria from genus Salmonella specifically serotypes Typhi and Paratyphi A B and C belonging to Enterobacteriaceae family cause it. Bacterial proliferation occurs pretty rapidly in digestive systems mainly affecting stomach and intestines with rather chaotic consequences unfolding subsequently [14] [16] [17]. Typhoid fever primarily spreads through faecal-oral route involving ingestion of water or food sullied by faeces urine or vomit from bacterial carriers. Transmission happens directly through bodily contact with infected people or indirectly via food tainted with viral particles somehow. Typhoid fever spreads rapidly across Mayo-Sava Mayo-Tsanaga and Logone-et-Chari departments in Cameroon’s Far North region due largely to water scarcity [18] [19]. Underserved populations exposed to tainted water sources are particularly vulnerable thereby amplifying risk of epidemics and muddying implementation of control measures greatly.

A mathematical epidemiological model that accounts for disparities in drinking water access and bacterium transmission dynamics within diverse populations is hereby proposed. A population gets categorised into five distinct compartments under this model.

• S (Susceptible): Individuals who have not contracted the infection and are exposed to the risk of contamination.

• I (Symptomatically infected): Individuals who are contagious and exhibit clinical signs.

• C (Chronic carriers): Individuals who are infected but do not exhibit symptoms and excrete bacteria.

• R (Withdrawn): Individuals who have been immunised for life following infection.

• W (Contaminated environment): A concentration of pathogenic bacteria in water.

There are two main routes of transmission:

• Direct transmission (human-to-human contact): Susceptible individuals can become infected through contact with symptomatic (I) or chronic (C) people, depending on a transmission rate ${\beta }_p$ .

• Environmental transmission (water contamination): A fraction of the population is exposed to indirect transmission through contaminated water at a rate ${\beta }_{w}W$ .

Thus, the susceptible population is divided into two subgroups:

• $S_1$ : The population with access to quality drinking water, with a risk of infection limited to ${\beta }_p\left(I+\omega C\right)$ .

• $S_2$ : The population exposed to contaminated water, with an increased risk of infection ${\beta }_p\left(I+\omega C\right)+{\beta }_{w}W$ .

The following dynamics of infection are postulated:

1) Infected individuals remain contagious for an average duration of $1/\delta$ .

2) A fraction $\theta$ of the infected become chronic (lifetime) carriers, contributing to transmission at a relative rate $\omega <1$ .

3) Pathogenic bacteria excreted by the infected (I) and chronic carriers (C) contaminate the water, but their viability declines at a rate $\xi$ .

The population is regenerated by births at a rate $b$ , with probability $p$ of belonging to $S_1$ and $1-p$ of belonging to $S_2$ . Natural mortality is denoted $\mu$ , and a disease-induced mortality is considered with a rate $\alpha$ . Let $N=S_1+S_2+I+C+R$ denote the size of the total population. Based on the above assumptions, the following system of differential equations is formulated:

$\{\begin{array}{l}{\dot{S}}_1=pbN-\mu S_1-{\beta }_p\left(I+\omega C\right)S_1\hfill \\ {\dot{S}}_2=\left(1-p\right)bN-\mu S_2-\left({\beta }_p\left(I+\omega C\right)+{\beta }_{w}W\right)S_2\hfill \\ \dot{I}=-\left(\delta +\mu \right)I+{\beta }_p\left(I+\omega C\right)\left(S_1+S_2\right)+{\beta }_{w}W{S}_2\hfill \\ \dot{C}=-\mu C-\delta \theta I\hfill \\ \dot{R}=-\mu R+\delta \left(1-\alpha -\theta \right)I\hfill \\ \dot{W}=-\xi W-\gamma \left(I+\omega C\right)\hfill \end{array}$

The variables and parameters are summarised in Table 1 and Table 2.

Table 1. Describing the state variables of the epidemic model.

State value	Describe
$S_1$	People expected to be able to access safe water
$S_2$	People not expected to be able to access safe water
$I$	Number of infected individuals with clinical signs
$C$	Infected persons without clinical signs (chronic carriers)
$R$	Number of people with durable immunity
$N$	Total population
$W$	CFU/ml of pathogenic bacteria in water

Table 2. Describing the model parameters.

Parameters	Describe
$b$	Fertility rate
$p$	Probability of neonates in $S_1$
$\mu$	Death rates
$1/\delta$	How long it lasts
${\beta }_p$	Direct transfer rate (human to human)
${\beta }_w$	Indirect transmission (waterborne)
$\alpha$	Proportion of infected persons with disease-related mortality
$\theta$	Percentage of infected people who become chronic carriers
$\omega$	Relative infectiousness of chronic patients (contribution of chronic carriers to disease spread)
$\gamma$	Rate of leaching to watercourses
$\xi$	Infectious particle degradation in water

2.1. Parameter Calibration

Parameter calibration is critical to ensure that model simulations accurately reflect the dynamics of typhoid fever in Far North Cameroon. The arid climate and landlocked terrain severely compromise health infrastructure in departments such as Mayo-Tsanaga, where access to drinking water remains limited [20]-[22]. Demographic parameters such as the birth rate $b$ and natural mortality rate $\mu$ are well-documented in regional reports. In the Far North region, the birth rate is approximately 37.5 per 1000 inhabitants, while the natural mortality rate is around 11.8 per 1000 inhabitants, reflecting inadequate sanitary conditions.

The duration of Salmonella Typhi infectiousness depends on patient age, bacterial strain, and access to antibiotics [23]-[25]. In rural areas of the Far North, the period of contagiousness typically ranges from two to five weeks. Transmission rates, ${\beta }_p$ (direct person-to-person) and ${\beta }_w$ (environmental via water), are context-dependent and exhibit substantial variability. Estimates for ${\beta }_p$ in urban areas with relatively safe water are below 1.2 per person per day, whereas in rural settings with poor water access, it can exceed 2.0 per person per day. The indirect transmission rate via contaminated water, ${\beta }_w$ , typically ranges from 0.3 to 0.7 per person per day in rural areas and increases sharply during the rainy season due to runoff and environmental contamination.

Bacterial excretion in water ($\gamma$ ) and environmental degradation ($\xi$ ) are strongly influenced by temperature, sunlight exposure, and the presence of other microorganisms. Salmonella Typhi demonstrates high survival rates in shallow, stagnant water commonly found in dry tropical climates [26]-[28]. Field studies in sub-Saharan Africa report excretion rates ranging from 0.05 to 1.2 per day and degradation rates from 0.02 to 1.0 per day, depending on local water physico-chemical characteristics. To improve rigor, specific citations from controlled studies in comparable Sahelian contexts should be included here—or example, references providing measured bacterial persistence in rivers, wells, or storage tanks [29].

Typhoid-induced mortality varies with the speed of diagnosis and the availability of treatment. In Cameroon’s Far North, mortality rates can reach 5% under limited healthcare access, whereas in better-served regions, rates remain below 1%. The chronicity rate $\theta$ , representing asymptomatic carriers, ranges from 1% to 4% in existing literature. Chronic carriers serve as a significant reservoir, contributing substantially to disease persistence. The relative infectivity of chronic carriers, $\omega$ , is complex to quantify; literature suggests a range of 0.3% - 25% depending on hygiene and sanitation practices. Accurate local estimates require targeted field studies.

Simulations will explore diverse scenarios to optimize parameter calibration and enhance understanding of transmission dynamics, particularly in vulnerable departments. Incorporating local epidemiological and environmental data is crucial to improve model fidelity and provide actionable insights for public health interventions.

2.2. Calculation of $R_0$

Proposition 1: The basic reproduction number $R_0$ is determined by the spectrum of the Jacobian matrix evaluated at DFE = $\left(S_1^{\text{*}},S_2^{\text{*}},0,0,0,0\right)$ . It is

$R_0=\frac{{\beta }_p\left(S_1^{\text{*}}+S_2^{\text{*}}\right)}{\delta +\mu }+\frac{{\beta }_p\omega \left(S_1^{\text{*}}+S_2^{\text{*}}\right)\delta \theta }{\mu \left(\delta +\mu \right)}+\frac{{\beta }_w{S}_2^{\text{*}}}{\xi }.$

Proof. Let $F$ be the new infection rate matrix and $V$ the transition rate matrix between infected compartments given by:

$F=\left[\begin{array}{c}{\beta }_p\left(I+\omega C\right)\left(S_1+S_2\right)+{\beta }_{w}W{S}_2\\ 0\\ 0\end{array}\right]\text{and}V=\left[\begin{array}{c}\left(\delta +\mu \right)I\\ \mu C+\delta \theta I\\ \xi W+\gamma \left(I+\omega C\right)\end{array}\right]$

The Jacobians of $F$ and $V$ evaluated in the state free of infection, i.e., when $I=C=W=0$ (Figure 1), to analyze the local stability of the system.

Figure 1. Showing the evolution of $R_0$ as a function of the waterborne transmission rate (${\beta }_w$ ) and the relative infectivity of chronically ill patients ($\omega$ ).

$F=\frac{\partial F}{\partial \left(I,C,W\right)}=\left[\begin{array}{ccc}{\beta }_p\left(S_1^{\text{*}}+S_2^{\text{*}}\right)& {\beta }_p\omega \left(S_1^{\text{*}}+S_2^{\text{*}}\right)& {\beta }_w{S}_2^{\text{*}}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$

$V=\frac{\partial V}{\partial \left(I,C,W\right)}=\left[\begin{array}{ccc}\delta +\mu & 0& 0\\ -\delta \theta & \mu & 0\\ -\gamma & -\gamma \omega & \xi \end{array}\right]$

$R_0$ denoted by a rather obscure symbol can be expressed thusly $R_0$ equals $\rho$ of some matrix $A$ where $\rho \left(A\right)$ signifies its rather dominant eigenvalue. Latter gets defined on condition that $S_P\left(A\right)$ represents $A$ ’s spectrum quietly. Eigenvalues of $-F{V}^{-1}$ emerge as solutions of characteristic equation $\text{det}\left(-F{V}^{-1}-\lambda I\right)=0$ rather quietly.

$V^{-1}=\left[\begin{array}{ccc}\frac{1}{\delta +\mu }& 0& 0\\ \frac{-\delta \theta }{\mu \left(\delta +\mu \right)}& \frac{1}{\mu }& 0\\ \frac{\gamma }{\left(\delta +\mu \right)\xi -\gamma \omega \delta \theta }& \frac{\gamma \omega }{\left(\delta +\mu \right)\xi -\gamma \omega \delta \theta }& \frac{1}{\xi }\end{array}\right]$

$R_0=\frac{{\beta }_p\left(S_1^{\text{*}}+S_2^{\text{*}}\right)}{\delta +\mu }+\frac{{\beta }_p\omega \left(S_1^{\text{*}}+S_2^{\text{*}}\right)\delta \theta }{\mu \left(\delta +\mu \right)}+\frac{{\beta }_w{S}_2^{\text{*}}}{\xi }.$

Direct transmission between healthy and infectious individuals is signified by first term quite explicitly; second term reckons with chronic cases $C$ effect thoroughly and third term corresponds vaguely to environmental transmission $W$ via various routes.

The stochastic model presented in the subsequent section can be regarded as an extension of the deterministic framework, designed to incorporate the intrinsic randomness associated with typhoid transmission and environmental contamination. While both models share the same fundamental mechanisms, the stochastic parameters can be directly mapped onto their deterministic counterparts as follows:

In this manner, the stochastic model preserves the core transmission dynamics of the deterministic system while explicitly accounting for random fluctuations that may arise in small populations or under conditions of data uncertainty. The rates $\lambda$ and $\nu$ can be interpreted as individual-level realizations of the deterministic infection forces ${\beta }_p\left(I+\omega C\right)$ and ${\beta }_{w}W$ , scaled according to the effective population size. (Table 3)

Table 3. Mapping between stochastic and deterministic model parameters.

Stochastic Parameter	Deterministic Parameter	Interpretation
$\lambda$	${\beta }_p$	Direct human-to-human transmission rate. In the stochastic model, $\lambda$ corresponds to ${\beta }_p$ in the deterministic model, representing the risk of infection through contact with symptomatic individuals (I) or chronic carriers (C).
$\nu$	${\beta }_{w}W$	Environmental transmission rate. The parameter $\nu$ captures random exposure to contaminated water, analogous to ${\beta }_{w}W$ in the deterministic model, which quantifies the impact of bacterial concentration $W$ combined with the environmental transmission coefficient ${\beta }_w$ .
$\alpha$	$\alpha$	Incubation rate. In both models, $\alpha$ represents the rate at which exposed individuals progress to the infectious stage.
$\mu$	$\delta$	Removal or isolation rate. The stochastic parameter $\mu$ is approximately equivalent to $\delta$ in the deterministic model, reflecting the average duration an individual remains infectious prior to recovery or isolation.

This explicit mapping not only elucidates the conceptual rationale behind the stochastic extension but also establishes a direct connection between deterministic analyses and probabilistic simulations, thereby offering a coherent and robust framework for interpreting epidemic dynamics.

3. Stochastic Model of Typhoid Transmission

To better capture the dynamics of Salmonella Typhi transmission, we extend the deterministic model by introducing an additional compartment, $E$ , representing individuals who have been exposed to the pathogen but are not yet infectious nor symptomatic. The progression from exposure to infectiousness is governed by the incubation rate $\alpha$ , corresponding to an average latent period of $1/\alpha$ days.

The stochastic model has three primary objectives: 1) to analyse the temporal evolution of exposed and infectious populations, 2) to investigate theoretical properties such as stability and extinction probabilities, and 3) to enable probabilistic simulations based on a continuous-time Markov process.

Key parameters include:

• $\lambda$ : human-to-human transmission rate,

• $\nu$ : environmental transmission rate through contaminated water,

• $\alpha$ : incubation rate (progression from exposed to infectious),

• $\mu$ : removal or isolation rate of infectious individuals.

Several simplifying assumptions are made to render the model analytically and numerically tractable. Population dynamics, including births and natural deaths, are considered negligible over the short time frame of the epidemic, implying that the total population remains approximately constant. Chronic carriers are not explicitly included in this stochastic model, as waterborne transmission is assumed to be the dominant infection pathway. Individuals withdrawn following hospitalization or treatment are incorporated implicitly through the removal rate $\mu$ .

The stochastic model is based on a continuous-time bivariate Markov process ${\left(E_t,I_t\right)}_{t\ge 0}$ , taking its values in the set $S=\mathbb{N}\times \mathbb{N}$ , where $E_t$ and $I_t$ represent, respectively, the number of exposed and infectious individuals at time $t$ within the population.

At a given state $\left(e,i\right)=\left(s,t\right)$ in $S$ , the process evolves according to transitions determined by the disease transmission and recovery mechanisms. The possible transitions and their associated rates are defined as follows:

• $\left(e,i\right)\to \left(e+1,i\right)$ with a rate of $\lambda i+\nu$ .

• $\left(e,i\right)\to \left(e-1,i+1\right)$ with a rate of $\alpha e$ .

• $\left(e,i\right)\to \left(e,i-1\right)$ is possible with a rate of $\mu i$ .

The infinitesimal generator $Q$ of the process ${\left(E_t,I_t\right)}_{t\ge 0}$ is thus expressed as follows:

$\{\begin{array}{ll}{Q}_{\left(e,i\right),\left(e+1,i\right)}=\lambda i+\nu ,\hfill & \text{for}e,i\ge 0,\hfill \\ Q_{\left(e,i\right),\left(e-1,i+1\right)}=\alpha e,\hfill & \text{for}e\ge 1,i\ge 0,\hfill \\ Q_{\left(e,i\right),\left(e,i-1\right)}=\mu i,\hfill & \text{for}e\ge 0,i\ge 1,\hfill \\ Q_{\left(e,i\right),\left(e,i\right)}=-\left(\left(\lambda +\mu \right)i+\alpha e+\nu \right),\hfill & \text{for}e,i\ge 0,\hfill \end{array}$

where the parameters $\lambda ,\nu ,\alpha$ and $\mu$ are positive.

It is assumed that the initial state of the process is given by $\left(E_0,I_0\right)=\left(e,i\right)\in S$ . The following notation is employed:

$p_{\left(e,i\right),\left(e^{\prime },i^{\prime }\right)}\left(t\right)=\mathbb{P}\left(\left(E_t,I_t\right)=\left(e^{\prime },i^{\prime }\right)|\left(E_0,I_0\right)=\left(e,i\right)\right),\left(e^{\prime },i^{\prime }\right)\in S,t\ge 0.$

The transition probabilities of the stochastic process ${\left(E_t,I_t\right)}_{t\ge 0}$ satisfy the following progressive Kolmogorov equation:

$\begin{array}{c}\frac{\text{d}}{\text{d}t}{p}_{\left(e,i\right),\left(e^{\prime },i^{\prime }\right)}\left(t\right)=\left(\lambda i^{\prime }+\nu \right)p_{\left(e,i\right),\left(e^{\prime }-1,i^{\prime }\right)}\left(t\right)\\ +\alpha \left(e^{\prime }+1\right)p_{\left(e,i\right),\left(e^{\prime }+1,i^{\prime }-1\right)}\left(t\right)\\ +\mu \left(i^{\prime }+1\right)p_{\left(e,i\right),\left(e^{\prime },i^{\prime }+1\right)}\left(t\right)\\ -\left(\left(\lambda +\mu \right)i^{\prime }+\alpha e^{\prime }+\nu \right)p_{\left(e,i\right),\left(e^{\prime },i^{\prime }\right)}\left(t\right)\end{array}$

with the initial condition

$p_{\left(e,i\right),\left(e,i\right)}\left(0\right)=1\text{and}{p}_{\left(e,i\right),\left(e^{\prime },i^{\prime }\right)}\left(0\right)=0\text{for}\left(e,i\right)\ne \left(e^{\prime },i^{\prime }\right).$

The model establishes a robust and systematic framework for investigating the stochastic dynamics of typhoid fever, encompassing both direct and environmental transmission pathways while accounting for the intrinsic randomness of infection events. This formulation enables simulation-based analysis of epidemic scenarios and supports the derivation of essential epidemiological metrics, including the probability of disease extinction and the expected outbreak size.

3.1. Mathematical Properties of the Model

Proposition 2: Mean process behaviour ${\left(E_t,I_t\right)}_{t\ge 0}$ . Let $e_0,i_0\ge 0$ . If $\lambda =\mu$ , then

$\{\begin{array}{l}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]=\frac{\lambda \nu }{\lambda +\alpha }t-k\left({\text{e}}^{-\left(\lambda +\alpha \right)t}+\frac{\lambda }{\alpha }\right)+i_0\frac{\lambda }{\alpha }+\frac{\nu }{\lambda +\alpha }\\ {\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]=\frac{\lambda \nu }{\lambda +\alpha }t-k\left({\text{e}}^{-\left(\lambda +\alpha \right)t}-1\right)+i_0\end{array}$

and if $\lambda \ne \mu$ , then

$\{\begin{array}{l}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]=k_1{\text{e}}^{\frac{-\left(\mu +\alpha \right)+\sqrt{{\left(\mu +\alpha \right)}^2+4\alpha \lambda }}{2}t}+k_2{\text{e}}^{\frac{-\left(\mu -\alpha \right)+\sqrt{{\left(\mu +\alpha \right)}^2+4\alpha \lambda }}{2}t}+\frac{\mu \nu }{\alpha \left(\mu -\lambda \right)}\\ {\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]=k_3{\text{e}}^{\frac{-\left(\mu +\alpha \right)+\sqrt{{\left(\mu +\alpha \right)}^2+4\alpha \lambda }}{2}t}+k_4{\text{e}}^{\frac{-\left(\mu -\alpha \right)+\sqrt{{\left(\mu +\alpha \right)}^2+4\alpha \lambda }}{2}t}+\frac{\nu }{\mu -\lambda }\end{array}$

with

$k=\left(e_0+i_0-\frac{\nu }{\mu -\lambda }\right)\frac{\alpha }{\alpha +\lambda },$

$k_1=\frac{1}{a}\left(i_0+\frac{\mu -\alpha +a}{2}{e}_0-\frac{\nu }{\mu -\lambda }\left(1+\frac{\mu \left(\mu -\alpha +a\right)}{2\alpha }\right)\right),$

$k_2=e_0-k_1-\frac{\mu \nu }{\alpha \left(\mu -\lambda \right)},$

$k_3=\frac{\mu -\alpha +a}{2a}\left(i_0+\frac{\mu -\alpha +a}{2}{e}_0-\frac{\nu }{\mu -\lambda }\left(1+\frac{\mu \left(\mu -\alpha +a\right)}{2\alpha }\right)\right),$

$k_4=-\frac{\mu -\alpha +a}{2}{k}_2,a=\sqrt{{\left(\mu +\alpha \right)}^2+4\alpha \lambda }.$

Proof: Let $\left(e_0,i_0\right)\in S$ , and let $f$ be a positive or bounded function on $S$ . For all $e,i\ge 0$ , we have:

${\mathbb{E}}_{\left(e_0,i_0\right)}\left[f\left(E_t,I_t\right)\right]={\displaystyle \sum }_{\left(e,i\right)\in S}{p}_{\left(e_0,i_0\right),\left(e,i\right)}\left(t\right)f\left(e,i\right).$

The transition probabilities of the stochastic process ${\left(E_t,I_t\right)}_{t\ge 0}$ satisfy the following progressive Kolmogorov equation:

$\begin{array}{c}\frac{\text{d}}{\text{d}t}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[f\left(E_t,I_t\right)\right]={\displaystyle \sum }_{\left(e,i\right)\in S}{p^{\prime }}_{\left(e_0,i_0\right),\left(e,i\right)}\left(t\right)f\left(e,i\right)\\ ={\displaystyle \sum }_{\left(e,i\right)\in S}\left({\displaystyle \sum }_{\left(e^{\prime },i^{\prime }\right)\in S}{p}_{\left(e_0,i_0\right),\left(e^{\prime },i^{\prime }\right)}\left(t\right)Q_{\left(e^{\prime },i^{\prime }\right),\left(e,i\right)}\right)f\left(e,i\right)\\ ={\displaystyle \sum }_{\left(e^{\prime },i^{\prime }\right)\in S}{p}_{\left(e_0,i_0\right),\left(e^{\prime },i^{\prime }\right)}\left(t\right)\left({\displaystyle \sum }_{\left(e,i\right)\in S}{Q}_{\left(e^{\prime },i^{\prime }\right),\left(e,i\right)}f\left(e,i\right)\right)\\ ={\displaystyle \sum }_{\left(e^{\prime },i^{\prime }\right)\in S}{p}_{\left(e_0,i_0\right),\left(e^{\prime },i^{\prime }\right)}\left(t\right)(\left(\lambda i^{\prime }+\nu \right)\left(f\left(e^{\prime }+1,i^{\prime }\right)-f\left(e^{\prime },i^{\prime }\right)\right)\\ +\alpha e^{\prime }\left(f\left(e^{\prime }-1,i^{\prime }+1\right)-f\left(e^{\prime },i^{\prime }\right)\right)\\ +\mu i^{\prime }\left(f\left(e^{\prime },i^{\prime }-1\right)-f\left(e^{\prime },i^{\prime }\right)\right)).\end{array}$

Thus,

$\begin{array}{c}\frac{\text{d}}{\text{d}t}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[f\left(E_t,I_t\right)\right]={\mathbb{E}}_{\left(e_0,i_0\right)}[\left(\lambda I_t+\nu \right)\left(f\left(E_t+1,I_t\right)-f\left(E_t,I_t\right)\right)\\ +\alpha E_t\left(f\left(E_t-1,I_t+1\right)-f\left(E_t,I_t\right)\right)\\ +\mu I_t\left(f\left(E_t,I_t-1\right)-f\left(E_t,I_t\right)\right)].\end{array}$

For $f\left(E_t,I_t\right)=E_t$ , we have:

$\frac{\text{d}}{\text{d}t}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]={\mathbb{E}}_{\left(e_0,i_0\right)}\left[\lambda I_t+\nu -\alpha E_t\right]=-\alpha {\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]+\lambda {\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]+\nu .$

For $f\left(E_t,I_t\right)=I_t$ , we have:

$\begin{array}{c}\frac{\text{d}}{\text{d}t}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]={\mathbb{E}}_{\left(e_0,i_0\right)}\left[\lambda I_t+\nu -\alpha E_t\right]\\ =-\alpha {\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]+\lambda {\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]+\nu .\end{array}$

The following system of differential equations is obtained by combining the two equations:

$\{\begin{array}{l}\frac{\text{d}}{\text{d}t}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]=-\alpha {\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]+\lambda {\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]+\nu \\ \frac{\text{d}}{\text{d}t}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]=\alpha {\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]-\mu {\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]\end{array}$

Upon solving the system, the following expectation expressions are obtained: In the event that $\lambda =\mu$ , then

$\{\begin{array}{l}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]=\frac{\lambda \nu }{\lambda +\alpha }t-k\left({\text{e}}^{-\left(\lambda +\alpha \right)t}+\frac{\lambda }{\alpha }\right)+i_0\frac{\lambda }{\alpha }+\frac{\nu }{\lambda +\alpha }\\ {\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]=\frac{\lambda \nu }{\lambda +\alpha }t-k\left({\text{e}}^{-\left(\lambda +\alpha \right)t}-1\right)+i_0\end{array}$

and if $\lambda \ne \mu$ , then

$\{\begin{array}{l}{\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right]=k_1{\text{e}}^{\frac{-\left(\mu +\alpha \right)+\sqrt{{\left(\mu +\alpha \right)}^2+4\alpha \lambda }}{2}t}+k_2{\text{e}}^{\frac{-\left(\mu -\alpha \right)+\sqrt{{\left(\mu +\alpha \right)}^2+4\alpha \lambda }}{2}t}+\frac{\mu \nu }{\alpha \left(\mu -\lambda \right)}\\ {\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]=k_3{\text{e}}^{\frac{-\left(\mu +\alpha \right)+\sqrt{{\left(\mu +\alpha \right)}^2+4\alpha \lambda }}{2}t}+k_4{\text{e}}^{\frac{-\left(\mu -\alpha \right)+\sqrt{{\left(\mu +\alpha \right)}^2+4\alpha \lambda }}{2}t}+\frac{\nu }{\mu -\lambda }\end{array}$

The constants $k,k_i$ for $i=1,2,3,4$ are determined from the initial condition $\left(E_0,I_0\right)=\left(e_0,i_0\right)$ .

Corollary. Asymptotically, under the initial condition $\left(E_0,I_0\right)=\left(e_0,i_0\right)$ , the average behaviour of the process ${\left(E_t,I_t\right)}_{t\ge 0}$ depends on the ratio between the parameters $\lambda$ and $\mu$ :

• If $\lambda >\mu$ : The average numbers of exposed and infected individuals grow exponentially. This phenomenon can be interpreted as an uncontrolled spread of the epidemic, characterised by a rapid increase in infections over time.

• If $\lambda =\mu$ : The expected numbers of exposed and infected individuals increase linearly over time, indicating a proportional progression of the epidemic.

• If $\lambda <\mu$ : The average process stabilises, leading to convergence towards a constant level of the average number of exposed and infected individuals, defined by:

$\underset{t\to \infty }{\text{lim}}\left({\mathbb{E}}_{\left(e_0,i_0\right)}\left[E_t\right],{\mathbb{E}}_{\left(e_0,i_0\right)}\left[I_t\right]\right)=\left(E^{\text{*}},I^{\text{*}}\right)=\left(\frac{\mu \nu }{\alpha \left(\mu -\lambda \right)},\frac{\nu }{\mu -\lambda }\right).$

The equilibrium values $E^{\text{*}}$ and $I^{\text{*}}$ are attained in instances wherein the transmission rate, denoted by $\lambda$ , is less than the recovery rate, symbolised by $\mu$ . (Figure 2)

Proposition 3: (Asymptotic behaviour of ${\left(E_t,I_t\right)}_{t\ge 0}$ )

If $\lambda <\alpha$ and $\frac{\lambda +\alpha }{2}<\mu$ , then ${\left(E_t,I_t\right)}_{t\ge 0}$ is positive recurrent and admits a unique invariant law.

Proof: The process ${\left(E_t,I_t\right)}_{t\ge 0}$ is clearly irreducible, its generator $L$ is given by: for any function $f:{\mathbb{N}}^2\to ℝ$ ,

Figure 2. Simulation and analysis of the process $\left(E_t,I_t\right)$ and its mean for parameters $\lambda =0.23$ , $\mu =0.2$ , $\alpha =0.22$ and $\nu =0.02$ (for comparison, $\lambda =\mu =\alpha =0.2$ and $\nu =0.02$ ).

$\begin{array}{c}Lf\left(e,i\right)=\left(\lambda i+\nu \right)\left(f\left(e+1,i\right)-f\left(e,i\right)\right)\\ +\alpha e\left(f\left(e-1,i+1\right)-f\left(e,i\right)\right)\\ +\mu i\left(f\left(e,i-1\right)-f\left(e,i\right)\right).\end{array}$

Applying this generator $L$ to the function $f\left(e,i\right)=e^2+i^2$ , we obtain:

$\begin{array}{c}Lf\left(e,i\right)=\left(\lambda i+\nu \right)\left({\left(e+1\right)}^2-e^2\right)+\alpha e\left({\left(e-1\right)}^2-e^2\right)\\ +\alpha e\left({\left(i+1\right)}^2-i^2\right)+\mu i\left({\left(i-1\right)}^2-i^2\right)\\ \le 2\lambda ei-2\alpha e^2+2\alpha ei-2\mu i^2\\ \le \left(\lambda -\alpha \right)e^2+\left(\lambda +\alpha -2\mu \right)i^2\\ \le \text{max}\left\{\lambda -\alpha ,\lambda +\alpha -2\mu \right\}\left(e^2+i^2\right),\end{array}$

where $\text{max}\left\{\lambda -\alpha ,\lambda +\alpha -2\mu \right\}<0$ because $\lambda <\alpha$ and $\frac{\lambda +\alpha }{2}<\mu$ . More precisely,

$Lf\left(e,i\right)\le \text{max}\left\{\lambda -\alpha ,\lambda +\alpha -2\mu \right\}\left(e^2+i^2\right)+C{1}_{\left(e,i\right)\ge \left(e_0,i_0\right)}\text{with}C>0.$

Consequently, the process ${\left(E_t,I_t\right)}_{t\ge 0}$ satisfies the Foster-Lyapunov condition, and thus the Markov chain $\left(E_t,I_t\right)$ admits a unique invariant law and converges in $L^2$ when $t$ tends to infinity.

Estimating parameters of the exposed-infected model plays a crucial role in quantifying transmission dynamics and evolution of epidemics fairly accurately nowadays. Available data is often shrouded in uncertainty necessitating rather robust estimation methods amidst considerable ambiguity and limited information. A combination of Markov Chain Monte Carlo (MCMC) and variational inference underpins this study’s approach enabling estimation of probability distributions on key model parameters such as transmission, incubation, and cure rates. This strategy has been shown to take better account of epidemiological uncertainty and provide a more reliable estimate beneath noisy data.

Numbers of exposed individuals and infected ones aren’t directly observable at any given time $t$ . Daily new case counts accumulate sharply and jumps occur suddenly with amplitude −1 in $I_t$ process. The expression $O_{\left]s,t\right]}$ denotes the number of jumps of amplitude −1 in the process $I_t$ over the interval $\left]s,t\right]$ .

The joint process ${\left(E_t,I_t,O_{\left]s,t\right]}\right)}_{t\ge 0}$ is characterised by the following transition probabilities:

For $t>0$ ,

$p_{\left(e_0,i_0\right),\left(e^{\prime },i^{\prime },o\right)}\left(t\right)=\mathbb{P}\left(\left(E_t,I_t\right)=\left(e^{\prime },i^{\prime }\right),O_{\left]s,t\right]}=o|\left(E_0,I_0\right)=\left(e_0,i_0\right)\right).$

These probabilities are in accordance with the Kolmogorov equation and the lemma below.

Theorem 1: When $t>0$ and $e,i,e^{\prime },i^{\prime },o\in \mathbb{N}$ , we have

$\begin{array}{c}\frac{\text{d}}{\text{d}t}{p}_{\left(e,i\right),\left(e^{\prime },i^{\prime },o\right)}\left(t\right)=\left(\lambda i+\nu \right)p_{\left(e+1,i\right),\left(e^{\prime },i^{\prime },o\right)}\left(t\right)+\alpha e{p}_{\left(e-1,i+1\right),\left(e^{\prime },i^{\prime },o\right)}\left(t\right)\\ +\mu i{p}_{\left(e,i-1\right),\left(e^{\prime },i^{\prime },o-1\right)}\left(t\right)-\left(\left(\lambda +\mu \right)i+\alpha e+\nu \right)p_{\left(e,i\right),\left(e^{\prime },i^{\prime },o\right)}\left(t\right),\end{array}$

for $i\le i^{\prime }+o$ , and $p_{\left(e,i\right),\left(e^{\prime },i^{\prime },o\right)}\left(t\right)=0$ when $i>i^{\prime }+o$ .

The initial condition is as follows:

$p_{\left(e,i\right),\left(e\text{'},i\text{'},o\right)}\left(0\right)={\delta }_{i=i^{\prime }}{\delta }_{o=0}.$

Proof: We use the infinitesimal characterization of the distribution of the joint process. For $h>0$ and $i\le i^{\prime }+o$ , then we have:

The process ${\left(E_t,I_t\right)}_{t\ge 0}$ is a homogeneous Markov process, which implies that:

$\begin{array}{c}{p}_{\left(e,i\right),\left(e^{\prime },i^{\prime },o\right)}\left(t+h\right)=\left(\lambda i+\nu \right)p_{\left(e+1,i\right),\left(e^{\prime },i^{\prime },o\right)}\left(t\right)+\alpha eh{p}_{\left(e-1,i+1\right),\left(e^{\prime },i^{\prime },o\right)}\left(t\right)\\ +\mu ih{p}_{\left(e,i-1\right),\left(e^{\prime },i^{\prime },o-1\right)}\left(t\right)-\left(\left(\lambda +\mu \right)i+\alpha e+\nu \right)p_{\left(e,i\right),\left(e^{\prime },i^{\prime },o\right)}\left(t\right)+o\left(h\right).\end{array}$

Estimation problem arises here from hidden information model in this particular study. Cumulative new cases reported daily are utilized in lieu of direct data on exposed and infected individuals corresponding to discrete variations in process $I_t$ . Variations represented by amplitude jumps of −1 are observed roughly over time interval $\left]s,t\right]$ . Process $N_t$ equals summation of $O_k$ cumulatively counting new cases up to time $t$ essentially. Primary objective of this study involves estimating quadruplet parameters $\left(\lambda ,\mu ,\alpha ,\nu \right)$ fairly accurately using various novel methods. Initially methodology proposed entails consideration of a scenario where variables $E_t$ and $I_t$ get observed at each discrete time instant. Average behaviour of joint process $\left(E_t,I_t,N_t\right)$ gets analysed and its moments are expressed explicitly as some function of parameters. Stability condition defined previously guarantees invertibility of obtained system of expectations enabling explicit expression of parameters in terms of process moments $\left(E_t,I_t,N_t\right)$ . A plug-in method employing moment estimator facilitates estimation of these parameters thereafter. Estimation’s second phase incorporates missing data via hidden Markov chain framework where triplet $\left(E_t,I_t,N_t\right)$ figures as Hidden Markov Multichain Model.

Theorem 2: If $\lambda <\mu$ , then the parameters $\lambda$ , $\mu$ , $\alpha$ et $\nu$ are given by the following formulas:

$\lambda =\frac{N^{\star }\left(\frac{R^{\star }}{E^{\star }{I}^{\star }}-1\right)\left(E^{\star }+I^{\star }\right)}{I^{\star }\left(\frac{R^{\star }}{E^{\star }{I}^{\star }}-1\right)\left(E^{\star }+I^{\star }\right)},\mu =\frac{N^{\star }}{I^{\star }},\alpha =\frac{N^{\star }}{E^{\star }},\nu =I^{\star }\left(\mu -\lambda \right),$

where:

$\begin{array}{l}{N}^{\star }=\underset{t\to +\infty }{\text{lim}}\frac{{\mathbb{E}}_{\left(e_0,i_0,n_0\right)}\left[N_t\right]}{t},\\ E^{\star }=\underset{t\to +\infty }{\text{lim}}{\mathbb{E}}_{\left(e_0,i_0,n_0\right)}\left[E_t\right],\\ I^{\star }=\underset{t\to +\infty }{\text{lim}}{\mathbb{E}}_{\left(e_0,i_0,n_0\right)}\left[I_t\right],\\ R^{\star }=\underset{t\to +\infty }{\text{lim}}{\mathbb{E}}_{\left(e_0,i_0,n_0\right)}\left[E_t{I}_t\right].\end{array}$

with $\left(e_0,i_0,n_0\right)\in \mathbb{N}$ .

Proof: Let ${\left(E_t,I_t,N_t\right)}_{t\ge 0}$ be a Markov process with values in $\tilde{S}=\mathbb{N}\times \mathbb{N}\times \mathbb{N}$ whose generator $\tilde{Q}$ is given by: for all $e,i,n\ge 0$ ,

$\{\begin{array}{l}{\tilde{Q}}_{\left(e,i,n\right),\left(e+1,i,n\right)}=\lambda i+\nu ,\hfill \\ {\tilde{Q}}_{\left(e,i,n\right),\left(e-1,i+1,n\right)}=\alpha e,\hfill \\ {\tilde{Q}}_{\left(e,i,n\right),\left(e,i-1,n+1\right)}=\mu i,\hfill \\ {\tilde{Q}}_{\left(e,i,n\right),\left(e,i,n\right)}=-\left(\left(\lambda +\mu \right)i+\alpha e+\nu \right).\hfill \end{array}$

Let $\left(e_0,i_0,n_0\right)\in \tilde{S}$ . For any positive or bounded function $f$ on $\tilde{S}$ , and for all $e,i,n\ge 0$ , we have:

${\mathbb{E}}_{\left(e_0,i_0,n_0\right)}\left[f\left(E_t,I_t,N_t\right)\right]={\displaystyle \sum }_{\left(e,i,n\right)\in \tilde{S}}{p}_{\left(e_0,i_0,n_0\right),\left(e,i,n\right)}\left(t\right)f\left(e,i,n\right).$

The progressive Kolmogorov equation gives:

$\begin{array}{c}\frac{\text{d}}{\text{d}t}{\mathbb{E}}_{\left(e_0,i_0,n_0\right)}\left[f\left(E_t,I_t,N_t\right)\right]={\displaystyle \sum }_{\left(e,i,n\right)\in \tilde{S}}{p^{\prime }}_{\left(e_0,i_0,n_0\right),\left(e,i,n\right)}\left(t\right)f\left(e,i,n\right)\\ ={\displaystyle \sum }_{\left(e,i,n\right)\in \tilde{S}}\left({\displaystyle \sum }_{\left(e^{\prime },i^{\prime },n^{\prime }\right)\in \tilde{S}}{p}_{\left(e_0,i_0,n_0\right),\left(e^{\prime },i^{\prime },n^{\prime }\right)}\left(t\right){\tilde{Q}}_{\left(e^{\prime },i^{\prime },n^{\prime }\right),\left(e,i,n\right)}\right)\end{array}$

$\begin{array}{c}f\left(e,i,n\right)={\displaystyle \sum }_{\left(e^{\prime },i^{\prime },n^{\prime }\right)\in \tilde{S}}{p}_{\left(e_0,i_0,n_0\right),\left(e^{\prime },i^{\prime },n^{\prime }\right)}\left(t\right)\left({\displaystyle \sum }_{\left(e,i,n\right)\in \tilde{S}}{\tilde{Q}}_{\left(e^{\prime },i^{\prime },n^{\prime }\right),\left(e,i,n\right)}f\left(e,i,n\right)\right)\\ ={\displaystyle \sum }_{\left(e^{\prime },i^{\prime },n^{\prime }\right)\in \tilde{S}}{p}_{\left(e_0,i_0,n_0\right),\left(e^{\prime },i^{\prime },n^{\prime }\right)}\left(t\right)[\left(\lambda i^{\prime }+\nu \right)\left(f\left(e^{\prime }+1,i^{\prime },n^{\prime }\right)-f\left(e^{\prime },i^{\prime },n^{\prime }\right)\right)\\ +\alpha e^{\prime }\left(f\left(e^{\prime }-1,i^{\prime }+1,n^{\prime }\right)-f\left(e^{\prime },i^{\prime },n^{\prime }\right)\right)\\ +\mu i^{\prime }\left(f\left(e^{\prime },i^{\prime }-1,n^{\prime }+1\right)-f\left(e^{\prime },i^{\prime },n^{\prime }\right)\right)]\\ ={\mathbb{E}}_{\left(e_0,i_0,n_0\right)}[\left(\lambda I_t+\nu \right)\left(f\left(E_t+1,I_t,N_t\right)-f\left(E_t,I_t,N_t\right)\right)\\ +\alpha E_t\left(f\left(E_t-1,I_t+1,N_t\right)-f\left(E_t,I_t,N_t\right)\right)\\ +\mu I_t\left(f\left(E_t,I_t-1,N_t+1\right)-f\left(E_t,I_t,N_t\right)\right)].\end{array}$

In the context of functions $f\left(E_t,I_t,N_t\right)$ , such as $E_t,I_t,N_t$ , $E_t^2$ , $E_t{I}_t$ and $I_t^2$ , the system of differential equations obtained admits the following asymptotic solutions:

Consequently, if the parameter, $\lambda <\mu$ then:

$N^{\star }=\underset{t\to +\infty }{\text{lim}}\frac{E_{\left(e_0,i_0,n_0\right)}\left[N_t\right]}{t}=\frac{\mu \nu }{\mu -\lambda },E^{\star }=\underset{t\to +\infty }{\text{lim}}{E}_{\left(e_0,i_0,n_0\right)}\left[E_t\right]=\frac{\mu \nu }{\alpha \left(\mu -\lambda \right)},$

$I^{\star }=\underset{t\to +\infty }{\text{lim}}{\mathbb{E}}_{\left(e_0,i_0,n_0\right)}\left[I_t\right]=\frac{\nu }{\mu -\lambda },R^{\star }=\underset{t\to +\infty }{\text{lim}}{\mathbb{E}}_{\left(e_0,i_0,n_0\right)}\left[E_t{I}_t\right]=\frac{\mu \nu \left(\left(\mu +\alpha \right)\nu +\alpha \lambda \right)}{\alpha {(\mu -\lambda )}^2\left(\mu +\alpha \right)}.$

Consequently, over time, $I_t$ and $E_t$ attain stationary values, while $N_t$ increases linearly. These results delineate the asymptotic evolution of the system and facilitate the expression of the model parameters analytically. Consequently, estimators (of the moment type) of the parameters $\lambda$ , $\mu$ , $\nu$ and $\alpha$ are derived.

$\lambda =\frac{N^{\star }\left(\frac{R^{\star }}{E^{\star }{I}^{\star }}-1\right)\left(E^{\star }+I^{\star }\right)}{I^{\star }\left(\frac{R^{\star }}{E^{\star }{I}^{\star }}-1\right)\left(E^{\star }+I^{\star }\right)},\mu =\frac{N^{\star }}{I^{\star }},\alpha =\frac{N^{\star }}{E^{\star }},\nu =I^{\star }\left(\mu -\lambda \right),$

Process $E_t$ becomes observable within this study’s context while variables ${\left(E_n,I_n\right)}_{n\in \mathbb{N}}$ stay shrouded in obscurity beneath surface level observations. A Hidden Multi-Chain Markov Model incorporating conditional probabilities coupled between several chains extends classical Hidden Markov Model for estimating $E^{\star }$ , $I^{\star }$ and $R^{\star }$ . Estimating transition probabilities of hidden states ${\left(E_n,I_n\right)}_{n\in \mathbb{N}}$ from observations ${\left(O_n\right)}_{n\in \mathbb{N}}$ via adapted HMM algorithms enables subsequent Monte Carlo exploitation quite effectively.

3.2. The Hidden Multi-Chain Markov Model

Exposed-infected models under scrutiny feature chains signifying distinct population states namely exposed and infected and isolated ones occasionally. Coupling between chains reflects dynamic interaction between these groups; hidden state of infected depends heavily on state of exposed at time $n-1$ and vice versa [30]. On observation depends heavily on both $I_n$ and preceding $I_{n-1}$ thus introducing non-standard transmission scheme capturing epidemic propagation quite precisely.

Theorem 3: The process ${\left(Z_n,O_n\right)}_{n\in {\mathbb{N}}^{\text{*}}}$ is a hidden Markov model characterised by the triplet $M=\left(Q,\psi ,\rho \right)$ , where:

1) Transition matrix $Q$ of the hidden chain $\left(Z_n\right)$ :

$Q_{\left(e,i,j\right),\left(e^{\prime },i^{\prime },j^{\prime }\right)}=\frac{p_{\left(e^{\prime },i^{\prime }\right),\left(\cdot ,j^{\prime }\right)}{p}_{\left(e,i\right),\left(e^{\prime },i^{\prime }\right)}}{p_{\left(e^{\prime },i^{\prime }\right),\left(\cdot ,j^{\prime }\right)}}{\delta }_{i^{\prime }=j},\left(e,i,j\right),\left(e^{\prime },i^{\prime },j^{\prime }\right)\in {\mathbb{N}}^3,$

where $p_{\left(e,i\right),\left(e^{\prime },i^{\prime }\right)}=p_{\left(e,i\right),\left(e^{\prime },i^{\prime }\right)}\left(\text{Δ}t\right)$ is the transition matrix of the chain ${\left(E_n,I_n\right)}_{n\in \mathbb{N}}$ and

$p_{\left(e,i\right),\left(\cdot ,j\right)}={\displaystyle \sum }_k{p}_{\left(e,i\right),\left(k,j\right)}.$

2) The probability of emission of process $O$ given process $Z$ is represented by:

${\psi }_{\left(e,i,j\right)}\left(o\right)=\mathbb{P}\left(O_n=o|Z_n=\left(e,i,j\right)\right)=\frac{p_{\left(e,i\right),\left(\cdot ,j,o\right)}}{p_{\left(e,i\right),\left(\cdot ,j\right)}},e,i,o\in \mathbb{N},$

where $p_{\left(e,i\right),\left(\cdot ,j,o\right)}={\displaystyle {\sum }_k{p}_{\left(e,i\right),\left(k,j,o\right)}\left(\text{Δ}t\right)}$ as defined above.

3) Initial law $\rho$ of the steady state $Z$ :

${\rho }_{\left(e,i,j\right)}=\mathbb{P}\left(Z_1=\left(e,i,j\right)\right)=p_{\left(e,i\right),\left(\cdot ,j\right)}{\pi }_{\left(e,i\right)},$

where $\pi$ is the initial chain law of ${\left(E_n,I_n\right)}_{n\in \mathbb{N}}$ .

Proof:

1) The transition matrix $Q$ is defined as the probability of transition between the various hidden states of the process $\left(Z_n\right)$ . Each state $Z_n$ is defined by a triplet ${\left(E_n,I_n,J_n\right)}_{n\in \mathbb{N}}$ , where $E_n$ denotes the state of the exposed population, $I_n$ the state of the infected population, and $J_n$ the state of the isolated population at time $n$ . The transition probability between two states $\left(e,i,j\right)$ and $\left(e^{\prime },i^{\prime },j^{\prime }\right)$ is given by:

$\begin{array}{c}{Q}_{\left(e,i,j\right),\left(e^{\prime },i^{\prime },j^{\prime }\right)}=\mathbb{P}\left(Z_{n+1}=\left(e^{\prime },i^{\prime },j^{\prime }\right)|Z_n=\left(e,i,j\right)\right)\\ =\frac{\mathbb{P}\left(E_n=e^{\prime },I_n=i^{\prime },I_{n+1}=j^{\prime },E_{n-1}=e,I_{n-1}=i,I_n=j\right)}{\mathbb{P}\left(E_{n+1}=e,I_{n+1}=i,I_n=j\right)}.\end{array}$

Using the Markov property of the chain $\left(E_n,I_n\right)$ , we obtain:

$\begin{array}{c}{Q}_{\left(e,i,j\right),\left(e^{\prime },i^{\prime },j^{\prime }\right)}= \frac{\mathbb{P}\left(I_{n+1}=j^{\prime }|E_n=e^{\prime },I_n=i^{\prime }\right)}{\mathbb{P}\left(I_n=j|E_{n-1}=e,I_{n-1}=i\right)}\\ \times \mathbb{P}\left(E_n=e^{\prime },I_n=i^{\prime }|E_{n-1}=e,I_{n-1}=i\right){\delta }_{i^{\prime }=j}\\ =\frac{p_{\left(e^{\prime },i^{\prime }\right),\left(\cdot ,j^{\prime }\right)} p_{\left(e,i\right),\left(e^{\prime },i^{\prime }\right)}}{p_{\left(e^{\prime },i^{\prime }\right),\left(\cdot ,j^{\prime }\right)}}{\delta }_{i\text{'}=j}.\end{array}$

2) The emission probability $\psi$ is defined as the probability of observing $O_n=o$ , conditional on the hidden state $Z_n=\left(e,i,j\right)$ , i.e., the probability of the observed output for a given hidden state. This follows directly from the definition of conditional probability.

3) The initial law $\rho$ determines the distribution of the hidden state at initial time $n=1$ , i.e. the probability of occurrence of each hidden state $Z_1=\left(e,i,j\right)$ :

$\begin{array}{c}{\rho }_{\left(e,i,j\right)}=\mathbb{P}\left(Z_1=\left(e,i,j\right)\right)\\ =\mathbb{P}\left(I_1=j|\left(E_0,I_0\right)=\left(e,i\right)\right)\mathbb{P}\left(\left(E_0,I_0\right)=\left(e,i\right)\right)\\ =p_{\left(e,i\right),\left(\cdot ,j\right)}{\pi }_{\left(e,i\right)}.\end{array}$

Estimating parameters $\theta =\left(Q,\psi ,\rho \right)$ of an epidemic model requires accounting for noisy, incomplete, or uncertain observed data $O=\left\{O_1,O_2,\cdots ,O_n\right\}$ . Transition parameters $Q$ and emission parameters $\psi$ alongside initial law $\rho$ play a pivotal role in dynamics of epidemic risk evaluation. Direct solution of such problems using conventional methods proves rather arduous especially when available data remains severely limited and highly uncertain. A Bayesian approach leveraging Markov Chain Monte Carlo algorithm with Metropolis-Hastings method playing a pivotal role is adopted quite effectively here. Bayesian formulation of estimation problem aims to obtain a posteriori distribution of MMM parameters given by

$P\left(M|O\right)\propto P\left(O|M\right)P\left(M\right)$

quite elaborately. Likelihood function $P\left(O|M\right)$ defines likelihood of observations conditional on model parameters in a probabilistic sense. $P\left(O|M\right)$ is essentially the probability of observing data given some model. Probability of observing data $O_n$ given hidden transitions $Z_n$ and model parameters MMM is determined through intricate calculations. A priori distribution of parameters is given priority $P\left(M\right)$ . This prior can be defined based on historical data or biological knowledge from epidemiological studies. Estimation involves employing a Markov chain sampling posterior distributions of model parameters thereby ensuring robust solutions for parameters $M=\left(Q,\psi ,\rho \right)$ amidst uncertainty in observed data. The fundamental principle behind the algorithm involves proposing candidate values for parameters at each iteration and accepting or rejecting them based on likelihood. (Figure 3)

Figure 3. Parameterising and cross-validating a statistical model on simulated data.

3.3. Effects of Seasonal Rainfall and Flooding on Typhoid Dynamics: Markov Chain Monte Carlo (MCMC) and Variational Inference

To capture the influence of seasonal rainfall and flooding on typhoid fever dynamics, we employ a sequence of probabilistic and computational methods that build upon one another. The stochastic compartmental model described in Section 3 provides a mechanistic framework for disease transmission, but its parameters are not directly observable. To link the latent epidemic states with actual case reports, a Hidden Markov Model (HMM) framework is introduced. The HMM formalism allows us to model the unobserved compartments (e.g., exposed and infectious individuals) as hidden states, while the reported case data serve as observations, thus bridging the gap between the mechanistic model and empirical data.

Once the HMM structure is defined, we require robust methods to estimate the model parameters from noisy and incomplete data. Markov Chain Monte Carlo (MCMC) is employed to generate samples from the posterior distribution of the parameters $M=\left(Q,\Psi ,\rho \right)$ conditional on observed epidemic data. This approach enables the quantification of uncertainty in parameter estimates and in the resulting epidemic trajectories. At each iteration $t$ , a candidate parameter set ${\theta }^{\prime }$ is proposed from a Gaussian distribution centered on the previous value ${\theta }^{\left(t-1\right)}$ , and the acceptance of ${\theta }^{\prime }$ is determined by the likelihood ratio comparing the new and old parameter configurations.

To complement MCMC and accelerate inference in higher-dimensional parameter spaces, Variational Inference (VI) is employed. VI approximates the posterior distribution by an analytically tractable distribution, optimizing the parameters of this approximation to minimize divergence from the true posterior. By combining MCMC and VI, we obtain both precise posterior samples and a fast approximation, ensuring robust estimation of the impact of environmental factors such as rainfall and flooding on typhoid transmission dynamics.

Markov chains get initialised rather haphazardly for each parameter $\theta$ comprising $Q$ , $\psi$ and $\rho$ mostly under specific conditions. A new value of ${\theta }^{\prime }$ emerges at each iteration $t$ via some function of previous value ${\theta }^{\left(t-1\right)}$ according to some transition rule often involving normal distribution around current value. A new candidate ${\theta }^{\prime }$ gets generated at each step for each parameter from a proposal function incorporating Gaussian noise on parameters. Likelihood is furnished by some pertinent data implicitly:

$L\left({\theta }^{\prime }|O\right)=P\left({\theta }^{\prime }|O\right)={\displaystyle \prod }_{n=1}^N{\psi }_{\left(e,i,j\right)}\left(O_n\right)$

In this context, $O=\left(O_1,O_2,\cdots ,O_N\right)$ denotes the observations.

In the event that the likelihood of the new parameter ${\theta }^{\prime }$ is greater than that of the current parameter ${\theta }^{″}$ , we accept ${\theta }^{\prime }$ . Conversely, if the likelihood of the new parameter is less than that of the current parameter, then the former is accepted with a probability given by:

$\alpha =\text{min}\left(1,\frac{L\left({\theta }^{\prime }|O\right)}{L\left({\theta }^{\left(t-1\right)}|O\right)}\right)$

Following a series of iterations, a parameter distribution is obtained, from which the values of the parameters $Q$ , Ψ, and $\rho$ can be estimated.

3.3.1. Estimation of the Transition Matrix Q

From the samples obtained via MCMC, the transition matrix $Q$ is estimated by calculating the relative frequency of transitions observed between different states $\left(e,i,j\right)$ . These transitions are contingent on both the causal relationships and the probabilities observed in the data.

$Q_{\left(e,i,j\right),\left(e^{\prime },i^{\prime },n^{\prime }\right)}=\frac{{{\displaystyle \sum }}_{n=1}^{N}I\left(Z_n=\left(e^{\prime },i^{\prime },n^{\prime }\right)\right)I\left(Z_{n-1}=\left(e,i,j\right)\right)}{{{\displaystyle \sum }}_{n=1}^{N}I\left(Z_{n-1}=\left(e,i,j\right)\right)}$

The expression $I\left(Z_1=\left(e,i,j\right)\right)$ represents the indicator function, which is equal to 1 if the hidden state $Z_1$ at time $n=1$ is equal to $\left(e,i,j\right)$ , and 0 otherwise.

3.3.2. Estimation of the Emission Probability, Denoted by ψ

In a similar manner, the emission probability, denoted by ${\psi }_{\left(e,i,j\right)}\left(o\right)$ , is estimated as a function of the observations $O$ , and the hidden states $Z_n=\left(e,i,j\right)$ . The calculation of the conditional probability of the observations given the hidden states is achieved by means of the MCMC samples.

${\psi }_{\left(e,i,j\right)}\left(o\right)=\frac{{{\displaystyle \sum }}_{n=1}^{N}I\left(Z_n=\left(e,i,j\right)\right)I\left(O_n=o\right)}{{{\displaystyle \sum }}_{n=1}^{N}I\left(Z_n=\left(e,i,j\right)\right)}$

3.3.3. Estimation of the Initial Distribution ρ

The initial distribution $\rho$ is estimated from the first observed states. The probability of an initial hidden state $Z_1=\left(e,i,j\right)$ is calculated from the observed data, taking into account the observed initial transitions:

${\rho }_{\left(e,i,j\right)}=\frac{{{\displaystyle \sum }}_{n=1}^{N}I\left(Z_1=\left(e,i,j\right)\right)}{N}$

The objective of this study is to identify a distribution $q\left(\theta \right)$ that optimises the Kullback-Leibler (KL) divergence from the exact a posteriori distribution $p\left(\theta |O\right)$ .

$D_{KL}\left(q\left(\theta \right)\parallel p\left(\theta |O\right)\right)={\displaystyle {\int }_{\Theta }q\left(\theta \right)\text{log}\left(\frac{q\left(\theta \right)}{p\left(\theta |O\right)}\right)\text{d}\theta }$

where $\theta =\left(Q,\psi ,\rho \right)$ and $\text{Θ}$ is the space of parameters: $\text{Θ}=\left\{Q,\psi ,\rho \right\}$ . (Figure 4)

3.3.4. Analysis of Q, ψ and ρ Matrices

Analysis of Q, ψ and ρ matrices revealed structurally defined seasonal dynamic in typhoid fever transmission in Far North of Cameroon directly correlated with fluctuations in water levels. Transitional probabilities reveal two major epidemic

Figure 4. Monitoring the course of an epidemic: Number of people exposed, infected and recovered over 365 days.

peaks with first peak observed in April-May coinciding with end of dry season when populations resort to dubious water sources due to scarcity. Bacterial contamination gets exacerbated by run-off water during July-August when rainy season reaches its zenith and second peak occurs. Between October and December a temporary decline manifests coinciding with stabilisation of post-rainy sanitary conditions rather quietly. Remarkable accuracy of analytical results provides robust scientific foundation for optimisation of prevention strategies by intensifying targeted hygiene control measures during critical periods.

4. Conclusion

A sophisticated mathematical framework was developed analysing transmission dynamics of typhoid fever in Far North Cameroon where unsanitary conditions prevail heavily. Coupling deterministic modeling rooted in differential equations with stochastic processes that capture uncertainties and random fluctuations revealed climatic and hydric factors’ pivotal role in driving epidemic dynamics erratically. Analysis of transition matrices Ψ, emission matrices $\rho$ and Markov processes Σ demonstrated existence of rigorous seasonal structuring of transmission with two major epidemic peaks identified April-May and July-August. April-May period corresponding roughly to end of dry season holds particular significance since scarcity of potable water forces populations towards dubious sources. Rainy season typically intensifies run-off and spreads pathogens resulting in infections surging anew during July and August pretty much everywhere. Health conditions stabilise somewhat after rains subside between October and December bringing a fleeting respite under normally dreary circumstances. Application of Markov Chain Monte Carlo algorithm and variational inference enabled accurate estimation of epidemiological parameters thereby affirming basic reproduction number $R_0$ surpasses epidemic threshold during critical seasons quite frequently. Simulations of various scenarios show that implementation of targeted measures like vaccination programmes and awareness campaigns substantially reduces transmission stabilising epidemic eventually nationwide. Findings signify considerable scientific progress in grasping dynamics of typhoid fever outbreaks amidst substantial regional climate fluctuations obviously. Such findings furnish a crucial underpinning pretty effectively for optimisation of health risk prevention strategies and diverse management tactics subsequently. Integration of water and epidemiological surveillance systems in real time alongside adaptation of health interventions according to seasonal rhythms will greatly facilitate effective response augmentation against water-borne disease.

Acknowledgements

Authors express profound gratitude rather quietly to local health authorities and community stakeholders in Far North region of Cameroon whose active collaboration and field observations greatly facilitated collection of contextual data and calibration of models utilised therein. Research did not receive specific funding from public, private or non-profit organisations notably. It forms part of a broader academic endeavour aimed at bolstering local epidemiological modelling capabilities in regions riddled with environmental and health issues. Authors also gratefully acknowledge valuable contributions made by informal networks drawing heavily on hydrology expertise and applied mathematics in public health.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Anderson, R.M. and MAY, R.M. (1992) Infectious Diseases of Humans: Dynamics and Control. Oxford University Press.
[2]	Brassard, P. and Vidal, L. (2004) Modélisation mathématique des maladies infectieuses: Applications à la fièvre typhoïde. Mathematical Modelling and Simulation, 12, 295-312.
[3]	Chan, C.H. and Hyman, J.M. (2008) Mathematical Models in Epidemiology: A Case Study on Typhoid Fever. Journal of Mathematical Biology, 56, 423-440.
[4]	Che, A.M. and Tchuente, M. (2013) Impact des facteurs climatiques sur la transmission des maladies hydriques au cameroun: Une analyse de modèles dynamique. African Journal of Environmental Studies, 21, 141-158.
[5]	Gomes, M.G.M. and Rojas, D.P. (2017) Mathematical Models for Waterborne Diseases: A Review of Typhoid Fever Models. SIAM Journal on Applied Dynamical Systems, 16, 1021-1043.
[6]	Green, H.K. and Ferguson, M.D. (2010) The Role of Waterborne Diseases in Public Health: Modeling and Control Strategies. International Journal of Environmental Health Research, 15, 449-460.
[7]	Kouadio, J.R. and Silué, K.D. (2015) Fièvre typhoïde et climat: Une approche de modélisation stochastique. Mathematical Biosciences and Engineering, 12, 1045-1062.
[8]	Liu, S.H. and Wang, L. (2019) Seasonal Impact of Floods and Precipitation on Typhoid Fever Dynamics in Sub-Saharan Africa. Journal of Infectious Diseases, 227, 1042-1054.
[9]	Meyers, L.A. and Liu, Y. (2009) The Role of Environmental Factors in Epidemic Disease Spread: A Model for Typhoid Fever in Developing Countries. Epidemiology, 20, 649-658.
[10]	Mounier-Jack, S. and Coker, R. (2010) Health Systems and Climate Change: Implications for Typhoid and Other Waterborne Diseases in Africa. Public Health, 124, 755-763.
[11]	Nabarro, D.F. (2004) Typhoid Fever in the Developing World: A Review of Environmental Risk Factors and Control Strategies. Tropical Medicine and International Health, 9, 534-540.
[12]	Ndiaye, S.M. and Barry, M.A. (2011) Modélisation mathématique des maladies infectieuses dans un contexte de changement climatique: Applications à la fièvre typhoïde en afrique. Mathematical and Computer Modelling, 53, 1299-1310.
[13]	Patel, R. and Singh, A.S. (2016) Linking Water, Sanitation, and Hygiene to Typhoid Transmission: A Mathematical Approach. Water Quality Research Journal of Canada, 51, 214-226.
[14]	Rahman, S. and Khan, A.M. (2018) Impact of Water Contamination on Typhoid Fever: A Mathematical Modeling Approach. Environmental Health Perspectives, 126, Article 057001.
[15]	Rojas, D.P. and Gomes, M.G.M. (2015) Stochastic Models for Typhoid Fever: A Case Study of Sub-Saharan Africa. Mathematical Biosciences, 271, 36-48.
[16]	Sampson, G. and Khanna, S. (2019) Impact des précipitations saisonnaires sur la propagation de la fièvre typhoïde: Modélisation et prévisions. Journal of Climate, 32, 1452-1463.
[17]	Singh, N. and Dutta, S. (2021) Integrated Modeling Approaches for Waterborne Disease Control: The Case of Typhoid Fever. Environmental Modelling & Software, 134, Article 104853.[CrossRef]
[18]	Sissoko, S. and Toure, S. (2014) Modélisation des maladies hydriques et leur sensibilité aux facteurs environnementaux: Cas de la fièvre typhoïde. Mathematical Methods in the Applied Sciences, 37, 1412-1425.
[19]	Smith, K.A. and Davis, L.S. (2007) The Role of Climatic and Hydrological Factors in the Transmission of Typhoid Fever. Environmental Research, 105, 172-182.
[20]	Tchuente, M. and Brice, G. (2011) Analysis of Typhoid Fever Transmission in Sub-Saharan Africa Using Compartmental Models. Mathematical Modelling and Computer Simulation, 9, 207-220.
[21]	Tirado, M.L. and Clarke, K. (2008) Waterborne Pathogens and Their Transmission in Arid Regions: A Mathematical Analysis of Typhoid Fever in the Middle East. Journal of Water and Health, 6, 279-288.
[22]	Vlassoff, C. and Tan, K. (2015) Integrating Climate and Health Modeling: A Study of Typhoid Fever Transmission in the Sahel. Climate Change and Health, 15, 44-56.
[23]	Wang, L. and Liu, J. (2019) Integrating Climate Data into Epidemiological Models for Typhoid Fever Prediction in Rural Areas of Africa. Epidemiology and Infection, 147, e25.
[24]	White, L. and Smith, J. (2006) Typhoid Fever: Modeling and Prediction of Transmission Patterns in Africa. Journal of Global Health, 14, 153-160.
[25]	World Health Organization (WHO) (2018) Global Strategy for Waterborne Disease Control: A Framework for Typhoid Fever Prevention. WHO Technical Report Series, 893.
[26]	Zhou, M. and Zhang, Z. (2017) Modeling the Interaction Between Water Contamination and Disease Spread: The Case of Typhoid Fever in Sub-Saharan Africa. Journal of Mathematical Biology, 74, 91-104.
[27]	Ziegler, T. and Brecht, L. (2009) Impact of Climate Change on Waterborne Disease Dynamics: A Model for Typhoid Fever. Ecological Modelling, 221, 1601-1612.
[28]	Sahli, M. and Lahmidi, S. (2012) Transmission Dynamics of Waterborne Diseases: A Stochastic Approach to Typhoid Fever. Mathematical Methods in the Applied Sciences, 34, 1291-1301.
[29]	Ng, C., Le, T.H., Goh, S.G., Liang, L., Kim, Y., Rose, J.B. and Gin, K.Y. (2010) Investigating the Decay Rates of Escherichia coli Relative to Vibrio parahaemolyticus and Salmonella Typhi in Tropical Coastal Waters. Water Research, 44, 4662-4672.
[30]	Yara, S. and Ouédraogo, N. (2014) Mathematical Modelling of Typhoid Fever and the Role of Water Contamination in Its Transmission. Mathematical Medicine and Biology, 31, 67-81.