Dynamic Model for the Epidemiology of Diarrhea and Simulation Considering Multiple Disease Carriers

Abstract

Diarrhea is responsible for killing around 525,000 children every year, even though it is preventable and treatable. This research focuses on both houseflies’ roles and humans’ roles in carrying pathogens causing diarrhea as multiple disease carriers. Both human and fly compartmental models are simulated with five diseases control strategies in order to identify the epidemic dynamics. The framework considers the life cycle of flies modeled into eggs, larvae, pupae, susceptible flies, and carrier flies, while the human system follows a compartment model as susceptible, infected, recovered, and back to susceptible again (SIRS). The relationships are modeled into an ordinary differential equation-based compartmental system. Then, the control parameters of the compartmental framework are analyzed. In order to propose effective control methods, five control strategies are considered: (1) elimination of flies’ breeding site, (2) sanitation, (3) installation of UV light trap, (4) good personal and food hygiene, and (5) water purification. Then, overall, ten control scenarios using the five control strategies are analyzed. Among them, effective control solutions considering various dynamic epidemiology are provided with the simulations and analyses. The proposed framework contributes to an effective control strategy in reducing the number of both flies and infected humans, since it minimizes the spread of the disease and considers cost-effectiveness.

1. Introduction

Infectious diarrhea, a disease causing fluid loss and dehydration, is the eighth leading cause of death, responsible for around 525,000 children deaths globally every year [1], mostly children in developing countries. Diarrhea is caused by infectious organisms, including viral and bacterial pathogens [2]. These are typically transmitted from the stool of one individual. This means the pathogen spreads through contaminated food and water or from an infectious person to a healthy person as a result of poor hygiene. Although diarrhea is both preventable and treatable, it is still fatal. This is due to several reasons, including a lack of awareness and shortages of existing lifesaving interventions.

Several existing studies showed that bacteria causing diarrhea are carried by houseflies [3,4] mainly. The house fly, Musca domestica L., is known to carry pathogens including bacteria, viruses, fungi, and parasites which cause life threatening diseases in humans and animals [5,6]. Houseflies breed in human feces [7] and the existing studies [8,9,10,11,12] have shown increased incidences of diarrhea during the periods of high fly density. Most importantly, several existing studies [8,13,14,15] have shown that a reduction in housefly density may affect the reduction in diarrhea incidence as well.

Mathematical models [16] are widely used to convert the real case into mathematical epidemiologic representations and predict the dynamics of infectious disease transmission so that they have vital roles in developing public health strategies for disease control and prevention. The formulation process of diarrhea considers multiple carriers (e.g., houseflies and infected humans). Then, the mathematical model is used for relevant disease control. Effective disease control has become an important part of computational epidemiology [17] that can provide useful guidelines for designing effective disease intervention strategies while balancing the costs of the control measures.

The objective of this paper is to understand and to formulate diarrhea’s dynamic epidemiology through a set of differential equation-based mathematical models and to predict the possible future for the effectiveness of disease control strategies not only for reducing the infected population, but also decreasing or even eliminating fly population as principal carriers of the disease as well.

The following section presents the relevant epidemiologic network of diarrhea and its mathematical models along with multiple controls and their theoretical solutions. Then, the effectiveness of the proposed framework and simulation results are analyzed in Section 3 and Section 4. In order to show the influence of each intervention towards the spreading of the disease, numerical simulations under several control scenarios are provided in Section 5.

2. Background and Epidemiologic Network Model

Figure 1 shows a flow diagram of diarrhea transmission through multiple carriers, including humans, flies and other environmental carriers. The system contains humans’ and flies’ epidemiologic systems, which are transformed into both (human and flies) compartmental models. As shown in Figure 1, the flies system shows the lifecycle of flies from egg stage $({\mathrm{E}}_{\mathrm{f}})$, larva $({\mathrm{L}}_{\mathrm{f}})$, pupa $({\mathrm{U}}_{\mathrm{f}})$, until adult flies, which divides into susceptible $({\mathrm{S}}_{\mathrm{f}})$, and carrier flies $({\mathrm{C}}_{\mathrm{f}})$. A carrier means an agent that carries a pathogen causing the disease but it does not show any symptoms of illness.

The human system follows susceptible human $({\mathrm{S}}_{\mathrm{h}})$, infectious $({\mathrm{I}}_{\mathrm{h}})$, recovered $({\mathrm{R}}_{\mathrm{h}})$, and back to susceptible again, since no one has immunity against diarrhea. Figure 2 shows an epidemiologic network where human can be infected by the disease by consuming contaminated food or water due to carrier flies laying pathogens on it.

As shown in Figure 2, a set of mathematical models are derived from the fly-human compartmental model. The model parameters in Figure 2 are provided in

Table 1. Notations of variables and parameters.

Symbol	Description	Initial Values
Variables
Ef	The number of eggs of flies	100
Lf	The number of larvae of flies	100
Uf	The number of pupae of flies	10
Sf	The number of susceptible flies	10
Cf	The number of carrier fly	10
Sh	The number of susceptible humans	1000
Ih	The number of infected humans	1
Rh	The number of recovered humans	0
Parameters
$\mathsf{\lambda }$	An influx rate of susceptible flies	0.1
$\mathsf{\gamma }$	A rate of susceptible flies to become carrier	0.2
$\mathrm{p}$	Probability of female fly	0.5
${\mathsf{\psi }}_{\mathrm{e}}$	Average maturation rate from egg to larva	0.4
${\mathsf{\psi }}_{\mathrm{l}}$	Average maturation rate from larva to pupa	0.6
${\mathsf{\psi }}_{\mathrm{u}}$	Average maturation rate from pupa to adult fly	0.7
${\mathsf{\beta }}_{\mathrm{ef}}$	Natural death rate of eggs	0.1
${\mathsf{\beta }}_{\mathrm{lf}}$	Natural death rate of larvae	0.1
${\mathsf{\beta }}_{\mathrm{uf}}$	Natural death rate of pupae	0.1
${\mathsf{\beta }}_{\mathrm{sf}}$	Natural death rate of susceptible flies	0.1
${\mathsf{\beta }}_{\mathrm{cf}}$	Natural death rate of carrier flies	0.1
$\mathsf{\delta }$	An oviposit rate of adult female flies	0.3
$\mathrm{D}{1}_{\mathrm{sf}}$	Diffusion parameter among susceptible flies	0.001
$\mathrm{D}{1}_{\mathrm{cf}}$	Diffusion parameter among carrier flies	0.001
$\mathsf{\xi }$	Carrier fly’s laying rate of pathogen on water or food	0.6
$\mathsf{\theta }$	Influx rate of susceptible humans	0.1
$\mathsf{\mu }$	Rate from “susceptible” status to “infected” status in humans	0.3
$\mathsf{\epsilon }$	Rate from “infected” status to “recovered” status in humans	0.0008
$\mathsf{\sigma }$	Rate from “recovered” status to “susceptible” status in humans	0.001
${\mathsf{\beta }}_{\mathrm{sh}}$	Natural death rate of susceptible humans	0.0008
${\mathsf{\beta }}_{\mathrm{ih}}.$	Natural death rate of infected humans	0.0008
${\mathsf{\beta }}_{\mathrm{rh}}$	Natural death rate of recovered humans	0.0008
$\mathrm{D}{2}_{\mathrm{sh}}$	Diffusion parameter among susceptible humans	0.1
$\mathrm{D}{2}_{\mathrm{ih}}$	Diffusion parameter among infected humans	0.3
$\mathrm{D}{2}_{\mathrm{rh}}$	Diffusion parameter among recovered humans	0.1
$\mathsf{\eta }$	Rate of contaminated water or food to be consumed by susceptible humans	0.5
${\mathsf{\omega }}_{\mathrm{ih}}$	Disease-induced death rate of infected humans	0.3

.

Table 2. Ordinary differential equations (ODE) considering both the fly and the human system.

Compartment	ODE
${\mathrm{E}}_{\mathrm{f}}$	$\frac{{\mathrm{dE}}_{\mathrm{f}}}{\mathrm{dt}}= \mathrm{p}.{\mathrm{S}}_{\mathrm{f}}.\mathsf{\delta }+\mathrm{p}.{\mathrm{C}}_{\mathrm{f}}.\mathsf{\delta }-{\mathsf{\beta }}_{\mathrm{ef}}.{\mathrm{E}}_{\mathrm{f}}-{\mathsf{\psi }}_{\mathrm{e}}.{\mathrm{E}}_{\mathrm{f}}$
${\mathrm{L}}_{\mathrm{f}}$	$\frac{{\mathrm{dL}}_{\mathrm{f}}}{\mathrm{dt}}={ \mathsf{\psi }}_{\mathrm{e}}.{\mathrm{E}}_{\mathrm{f}}-{\mathsf{\beta }}_{\mathrm{lf}}.{\mathrm{L}}_{\mathrm{f}}-{\mathsf{\psi }}_{\mathrm{l}}.{\mathrm{L}}_{\mathrm{f}}$
${\mathrm{U}}_{\mathrm{f}}$	$\frac{{\mathrm{dU}}_{\mathrm{f}}}{\mathrm{dt}}={\mathsf{\psi }}_{\mathrm{l}}.{\mathrm{L}}_{\mathrm{f}}-{\mathsf{\beta }}_{\mathrm{uf}}.{\mathrm{U}}_{\mathrm{f}}-{\mathsf{\psi }}_{\mathrm{u}}.{\mathrm{U}}_{\mathrm{f}}$
${\mathrm{S}}_{\mathrm{f}}$	$\frac{{\mathrm{dS}}_{\mathrm{f}}}{\mathrm{dt}}=\mathsf{\lambda }.{\mathrm{S}}_{\mathrm{f}}+{\mathsf{\psi }}_{\mathrm{u}}.{\mathrm{U}}_{\mathrm{f}}-{\mathsf{\beta }}_{\mathrm{sf}}.{\mathrm{S}}_{\mathrm{f}}-\mathsf{\gamma }.{\mathrm{S}}_{\mathrm{f}}+\mathrm{D}{1}_{\mathrm{sf}}.{\mathrm{S}}_{\mathrm{f}}$
${\mathrm{C}}_{\mathrm{f}}$	$\frac{{\mathrm{dC}}_{\mathrm{f}}}{\mathrm{dt}}=\mathsf{\gamma }.{\mathrm{S}}_{\mathrm{f}}-{\mathsf{\beta }}_{\mathrm{cf}}.{\mathrm{C}}_{\mathrm{f}}+\mathrm{D}{1}_{\mathrm{cf}}.{\mathrm{C}}_{\mathrm{f}}$
${\mathrm{S}}_{\mathrm{h}}$	$\frac{{\mathrm{dS}}_{\mathrm{h}}}{\mathrm{dt}}= \mathsf{\theta }.{\mathrm{S}}_{\mathrm{h}}-{\mathsf{\beta }}_{\mathrm{sh}}.{\mathrm{S}}_{\mathrm{h}}-\mathsf{\mu }.{\mathrm{S}}_{\mathrm{h}}+\mathsf{\sigma }.{\mathrm{R}}_{\mathrm{h}}+\mathrm{D}{2}_{\mathrm{sh}}.{\mathrm{S}}_{\mathrm{h}}$
${\mathrm{I}}_{\mathrm{h}}$	$\frac{{\mathrm{dI}}_{\mathrm{h}}}{\mathrm{dt}}= \mathsf{\mu }.{\mathrm{S}}_{\mathrm{h}}-{\mathsf{\beta }}_{\mathrm{ih}}.{\mathrm{I}}_{\mathrm{h}}-{\mathsf{\omega }}_{\mathrm{ih}}{\mathrm{I}}_{\mathrm{h}}-\mathsf{\epsilon }.{\mathrm{I}}_{\mathrm{h}}+\mathrm{D}{2}_{\mathrm{ih}}.{\mathrm{I}}_{\mathrm{h}}$
${\mathrm{R}}_{\mathrm{h}}$	$\frac{{\mathrm{dR}}_{\mathrm{h}}}{\mathrm{dt}}= \mathsf{\epsilon }.{\mathrm{I}}_{\mathrm{h}}-{\mathsf{\beta }}_{\mathrm{rh}}.{\mathrm{R}}_{\mathrm{h}}-\mathsf{\sigma }.{\mathrm{R}}_{\mathrm{h}}+\mathrm{D}{2}_{\mathrm{rh}}.{\mathrm{R}}_{\mathrm{h}}$

$\mathsf{\mu }=\mathsf{\xi }·\mathsf{\eta }$.

shows the differential equation-based transition models with respect to the lifecycle of the fly and the epidemiologic processes of a human. Humans are classified with the Susceptible-Infected-Recovered-Susceptible (SIRS) reaction-diffusion compartmental model. The ordinary differential equations (ODE) are generated based on each of states in the network. The time unit is “day” in this research.

The solutions of the compartment models, as shown in Table 2, have to be non-negative (positiveness) and exist (boundedness). The below conditions show the positiveness and boundedness of the compartment model.

Let the initial dataset be $S_h\left(0\right)>0, I_h\left(0\right)>0, S_f\left(0\right)>0, C_f\left(0\right)>0, E_f\left(0\right)>0$ and $({\mathrm{L}}_{\mathrm{f}}\left(0\right),{ \mathrm{U}}_{\mathrm{f}}\left(0\right),{ \mathrm{R}}_{\mathrm{h}}\left(0\right))\in \coprod$. $\coprod$ is positive and bounded interval for all time $t>0$.

Consider the Inequality (1) at time $t$.

$$\frac{{\mathrm{dE}}_{\mathrm{f}}}{\mathrm{dt}}\ge -\left({\mathsf{\beta }}_{\mathrm{ef}}+{\mathsf{\psi }}_{\mathrm{e}}\right){\mathrm{E}}_{\mathrm{f}}$$

(1)

$$\int \frac{{\mathrm{dE}}_{\mathrm{f}}}{{\mathrm{E}}_{\mathrm{f}}}\ge -\int ({\mathsf{\beta }}_{\mathrm{ef}}+{\mathsf{\psi }}_{\mathrm{e}})\mathrm{d}\left(\mathrm{t}\right)$$

(2)

$${\mathrm{E}}_{\mathrm{f}}\left(\mathrm{t}\right)\ge {\mathrm{E}}_{\mathrm{f}}\left(0\right){\mathrm{e}}^{-{\displaystyle \int }({\mathsf{\beta }}_{\mathrm{ef}}+{\mathsf{\psi }}_{\mathrm{e}})\mathrm{d}\left(\mathrm{t}\right)}\ge 0$$

(3)

Inequality (1) proves that the solution set is positive for all time $t>0$. The same rule can be applied for Inequality (2) until Inequality (6) at time $t$.

$$\frac{{\mathrm{dL}}_{\mathrm{f}}}{\mathrm{dt}}\ge -\left({\mathsf{\beta }}_{\mathrm{lf}}+{\mathsf{\psi }}_{\mathrm{l}}\right){\mathrm{L}}_{\mathrm{f}}$$

(4)

$${\displaystyle \int }\frac{{\mathrm{dL}}_{\mathrm{f}}}{{\mathrm{L}}_{\mathrm{f}}}\ge -{\displaystyle \int }({\mathsf{\beta }}_{\mathrm{lf}}+{\mathsf{\psi }}_{\mathrm{l}})\mathrm{d}\left(\mathrm{t}\right)$$

(5)

$${\mathrm{L}}_{\mathrm{f}}\left(\mathrm{t}\right)\ge {\mathrm{L}}_{\mathrm{f}}\left(0\right){\mathrm{e}}^{-{\displaystyle \int }({\mathsf{\beta }}_{\mathrm{lf}}+{\mathsf{\psi }}_{\mathrm{l}})\mathrm{d}\left(\mathrm{t}\right)}\ge 0$$

(6)

3. Epidemic Model Dynamics

The dynamics of the models provided in Table 2 depend on the basic reproduction number, which is defined as the average number of secondary infections of an infectious human [18,19].

The basic reproduction number is denoted by ${\mathrm{R}}_0$ as the number of secondary infections caused by an infected individual. If the value of ${\mathrm{R}}_0<1$, then the disease dies out. While the value of ${\mathrm{R}}_0>1$, then the number of infectious individuals increases and the disease invades the population. Let $\mathrm{x}={\left({\mathrm{x}}_1,{\mathrm{x}}_2,\dots ,{\mathrm{x}}_{\mathrm{n}}\right)}^{\mathrm{T}}$ be the number of individuals in each compartment where the first $m<n$ compartments contain infected individuals. The disease-free equilibrium (DFE) is given by (S,E,I,R) = ($S_0,0,0,0$). The terms S, E, I and R mean Susceptible, Exposed, Infected, and Recovered, respectively. In this model dynamics, it is assumed that the DFE exists and is stable in the absence of disease. Consider those equations written in the form $\frac{{\mathrm{dx}}_{\mathrm{i}}}{\mathrm{dt}}={\mathrm{F}}_{\mathrm{i}\left(\mathrm{x}\right)}–{\mathrm{V}}_{\mathrm{i}\left(\mathrm{x}\right)}$ for $\mathrm{i}=1,2,\dots ,\mathrm{m}$ where, ${\mathrm{F}}_{\mathrm{i}\left(\mathrm{x}\right)}$ is the rate of appearance of new infections in compartment i and ${\mathrm{V}}_{\mathrm{i}\left(\mathrm{x}\right)}$ is the rate of other transitions between compartment i and other infected compartments. F ($\mathrm{F}=\left[\frac{\partial {\mathrm{F}}_{\mathrm{i}\left({\mathrm{x}}_0\right)}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right]$) is entry wise non-negative and V ($\mathrm{V}=\left[\frac{\partial \mathrm{Vi}\left({\mathrm{x}}_0\right)}{\partial {\mathrm{x}}_{\mathrm{j}}}\right]$) is a non-singular matrix for $1\le \mathrm{i},\mathrm{j}\le \mathrm{m}$. $\mathrm{F}\left(\mathrm{X},\mathrm{Y}\right)$ denotes a vector of new infection rates (flows from X to Y) and $\mathrm{V}\left(\mathrm{X},\mathrm{Y}\right)$ is a vector of all other rates. For each compartment, an inflow in V is negative and an outflow in V is positive. It is assumed that $\mathrm{F}\left(0,\mathrm{Y}\right)=0$ and $\mathrm{V}\left(0,\mathrm{Y}\right)=0$. ${\mathrm{FV}}^{-1}$ is called the next generation matrix where the spectral radius of it is equal to ${\mathrm{R}}_0$, which is the largest eigenvalue of ${\mathrm{FV}}^{-1}$. ${\mathrm{FV}}^{-1}$ is derived using Equations (7)–(10).

$$\mathrm{F}=\left[\frac{\partial {\mathsf{\mu }\mathrm{S}}_{\mathrm{h}}\mathrm{I}}{\partial \mathrm{I}}\right]=\left[{\mathsf{\mu }\mathrm{S}}_{\mathrm{h}}\right]$$

(7)

$$\mathrm{V}=\left[\frac{\partial \left({\mathsf{\beta }}_{\mathrm{ih}}+{\mathsf{\omega }}_{\mathrm{ih}}+\mathsf{\epsilon }+\mathrm{D}{2}_{\mathrm{ih}}\right){\mathrm{I}}_{\mathrm{h}}-\partial {\mathsf{\mu }\mathrm{S}}_{\mathrm{h}}}{\partial {\mathrm{I}}_{\mathrm{h}}}\right]$$

(8)

$${\mathrm{V}}^{-1}=\frac{1}{{\mathsf{\beta }}_{\mathrm{ih}}+{\mathsf{\omega }}_{\mathrm{ih}}+\mathsf{\epsilon }+\mathrm{D}{2}_{\mathrm{ih}}}$$

(9)

$${\mathrm{FV}}^{-1}=\left[\frac{{\mathsf{\mu }\mathrm{S}}_{\mathrm{h}}}{{\mathsf{\beta }}_{\mathrm{ih}}+{\mathsf{\omega }}_{\mathrm{ih}}+\mathsf{\epsilon }+\mathrm{D}{2}_{\mathrm{ih}}}\right]$$

(10)

Then, ${\mathrm{R}}_0$ is determined as the basic reproduction number using ${\mathrm{FV}}^{-1}$.

$${\mathrm{R}}_0=\frac{{\mathsf{\mu }\mathrm{S}}_{\mathrm{h}}}{{\mathsf{\beta }}_{\mathrm{ih}}+{\mathsf{\omega }}_{\mathrm{ih}}+\mathsf{\epsilon }+\mathrm{D}{2}_{\mathrm{ih}}}$$

(11)

After modeling the infection processes of diarrhea provided in Section 2 and Section 3, the derived differential equations are solved with the initial conditions of ${\mathrm{E}}_{\mathrm{f}}=10$; ${\mathrm{L}}_{\mathrm{f}}=10$; ${\mathrm{U}}_{\mathrm{f}}=10$; ${\mathrm{S}}_{\mathrm{f}}=10$; ${\mathrm{C}}_{\mathrm{f}}=10$; ${\mathrm{S}}_{\mathrm{h}}=1000$; ${\mathrm{I}}_{\mathrm{h}}=1$; ${\mathrm{R}}_{\mathrm{h}}=0$. The initial parameters are supposed with real-world examples. Figure 3 shows the disease’s infection simulations without any control method, in terms of humans’ and flies’ status. The simulation is performed using Matlab©. According to Figure 3, the disease starts with an outbreak at t = 0 and, the number of infected people is around 2700 at t = 20. The time unit is day in this research.

The provided ODE-based compartment model is used to check the epidemiologic processes of diarrhea and the effectiveness of control methods, which are provided in the following section.

4. Effective Control Framework Considering Multiple Disease Carriers

In this section, the previous fundamental infection model is extended with several control strategies on the spread of the disease. The term “effectiveness” is evaluated with the blocks of additional infections and a cost-effective concept. This consideration may help to block the spread of SARS-2 or COIVID-19. The cost-effective criteria are determined using Pontryagin’s [20] maximum principle on five variations of control methods. Pontryagin [20,21] introduced the idea of adjoint functions which has a similar purpose as Lagrange multipliers, to append the differential equation to the objective function.

This section focuses on effective control frameworks to prevent additional infections of humans. In order to control additional spreads of diarrhea, five control methods are considered: (1) elimination of fly’s breeding site, (2) sanitation-related investment, (3) installation of UV light traps for killing flies, (4) good personal and food hygiene, and (5) water purification. In general, these methods intend to prevent additional infections among humans or to remove one of primary disease carriers—flies.

As the first control method, “Elimination of flies’ breeding site” is considered. As flies feed on garbage, food waste, and animal feces, therefore one way to prevent additional infections is to keep the places clean, and to spray pesticide on prone areas. This intervention could reduce or even eliminate the maturation rate of fly eggs, larva, and pupa, which eventually impacts the number of adult flies. Let $\mathsf{\alpha }$ denote the level of a breeding site elimination strategy ($0\le \mathsf{\alpha }\left(\mathrm{t}\right)\le 1$). The effect of breeding site elimination will decrease the maturation rate of eggs, larvae, and pupae, which is modelled as $\mathsf{\alpha }·{\mathrm{E}}_{\mathrm{f}}, \mathsf{\alpha }·{\mathrm{L}}_{\mathrm{f}}$and $\mathsf{\alpha }·{\mathrm{U}}_{\mathrm{f}}$, respectively.

The second control option is sanitation efforts [22], which include the sterilization of cooking utensils, washing, and drinking water by boiling it properly. The parameter $\mathsf{\pi }\left(\mathrm{t}\right)$ ($0\le \mathsf{\pi }\left(\mathrm{t}\right)\le 1$) denotes the level of sanitation strategy. The effect of sanitation is to decrease the number of pathogens, which is modelled as a reduction in the rate of the carrier fly population by the term $\mathsf{\pi }·{\mathrm{C}}_{\mathrm{f}}$.

The third control method is to install light traps. When flies see ultra-violet (UV) light, they are naturally lured in the direction of the source of the UV rays. This strategy is more desirable than spray for indoor fly control as it keeps the surroundings clean. $\mathrm{\vartheta } \left(0\le \mathsf{\vartheta }\left(\mathrm{t}\right)\le 1\right)$ denotes the effectiveness of the installation of UV light trap strategy. The effect of this effort is to decrease the number of susceptible and carrier flies by the terms $\mathsf{\vartheta }·{\mathrm{S}}_{\mathrm{f}}$ and $\mathsf{\vartheta }·{\mathrm{C}}_{\mathrm{f}}$.

The next consideration is to keep good personal and food hygiene. Good personal and food hygiene can be defined as handling, preparing, and storing food in a way that reduces the risk of becoming contaminated. Let $\mathsf{\rho } \left(0\le \mathsf{\rho }\left(\mathrm{t}\right)\le 1\right)$ be the good personal and food hygiene effort. The effect of this effort is to suppress the number of infected people by the term $\mathsf{\rho }·{\mathrm{I}}_{\mathrm{h}}$.

The final option is to purify water. Chlorine is often a choice for water purification, since it effectively inactivates the bacteria causing diarrhea, leaves residual protection, has low cost and is easy to transport and use. There are two chlorine-based options [23] used in diarrhea outbreaks, such as tablets and liquid. $\mathsf{\tau } \left(0\le \mathsf{\tau }\left(\mathrm{t}\right)\le 1\right)$ denotes the level of water purification strategy. The effect of water purification is to inactivate pathogens, which is modelled as a reduction in the disease transmission rate by the term $\mathsf{\tau }·{\mathrm{I}}_{\mathrm{h}}.$

Table 3. Control parameters.

Symbol	Description	Initial Values
Control Parameters
$\mathsf{\alpha }$	Effective control using eliminations of fly’s breeding site	0.03
$\mathsf{\pi }$	Effective rate using sanitation methods	0.1
$\mathsf{\vartheta }$	Effective rate using installation of UV light trap	0.04
$\mathsf{\rho }$	Effective rate using good personal and food hygiene	0.01
$\mathsf{\tau }$	Effective rate using water purification	0.02

explains each control parameter for each control option.

Based on the control parameters shown in Table 3, the initial transmission models are modified with additional terms to Equations (12)–(17). The additional term is underlined in each equation.

$$\frac{{\mathrm{dE}}_{\mathrm{f}}}{\mathrm{dt}}=\mathrm{p}.{\mathrm{S}}_{\mathrm{f}}.\mathsf{\delta }+\mathrm{p}.{\mathrm{C}}_{\mathrm{f}}.\mathsf{\delta }-{\mathsf{\beta }}_{\mathrm{ef}}.{\mathrm{E}}_{\mathrm{f}}-{\mathsf{\psi }}_{\mathrm{e}}.{\mathrm{E}}_{\mathrm{f}}-\underset{\_}{\mathsf{\alpha }.{\mathrm{E}}_{\mathrm{f}}}$$

(12)

$$\frac{{\mathrm{dL}}_{\mathrm{f}}}{\mathrm{dt}}={\mathsf{\psi }}_{\mathrm{e}}.{\mathrm{E}}_{\mathrm{f}}-{\mathsf{\beta }}_{\mathrm{lf}}.{\mathrm{L}}_{\mathrm{f}}-{\mathsf{\psi }}_{\mathrm{l}}.{\mathrm{L}}_{\mathrm{f}}-\underset{\_}{\mathsf{\alpha }.{\mathrm{L}}_{\mathrm{f}}}$$

(13)

$$\frac{{\mathrm{dU}}_{\mathrm{f}}}{\mathrm{dt}}={\mathsf{\psi }}_{\mathrm{l}}.{\mathrm{L}}_{\mathrm{f}}-{\mathsf{\beta }}_{\mathrm{uf}}.{\mathrm{U}}_{\mathrm{f}}-{\mathsf{\psi }}_{\mathrm{u}}.{\mathrm{U}}_{\mathrm{f}}-\underset{\_}{\mathsf{\alpha }.{\mathrm{U}}_{\mathrm{f}}}$$

(14)

$$\frac{{\mathrm{dS}}_{\mathrm{f}}}{\mathrm{dt}}=\mathsf{\lambda }.{\mathrm{S}}_{\mathrm{f}}+{\mathsf{\psi }}_{\mathrm{u}}.{\mathrm{U}}_{\mathrm{f}}-{\mathsf{\beta }}_{\mathrm{sf}}.{\mathrm{S}}_{\mathrm{f}}-\mathsf{\gamma }.{\mathrm{S}}_{\mathrm{f}}+\mathrm{D}{1}_{\mathrm{sf}}.{\mathrm{S}}_{\mathrm{f}}\underset{\_}{-\mathsf{\pi }.{\mathrm{S}}_{\mathrm{f}}-\mathsf{\vartheta }.{\mathrm{S}}_{\mathrm{f}}}$$

(15)

$$\frac{{\mathrm{dC}}_{\mathrm{f}}}{\mathrm{dt}}=\mathsf{\gamma }.{\mathrm{S}}_{\mathrm{f}}-{\mathsf{\beta }}_{\mathrm{cf}}.{\mathrm{C}}_{\mathrm{f}}+\mathrm{D}{1}_{\mathrm{cf}}.{\mathrm{C}}_{\mathrm{f}}\underset{\_}{- \mathsf{\pi }.{\mathrm{C}}_{\mathrm{f}}-\mathsf{\vartheta }.{\mathrm{C}}_{\mathrm{f}}}$$

(16)

$$\frac{{\mathrm{dI}}_{\mathrm{h}}}{\mathrm{dt}}=\mathsf{\mu }.{\mathrm{S}}_{\mathrm{h}}-{\mathsf{\beta }}_{\mathrm{ih}}.{\mathrm{I}}_{\mathrm{h}}-{\mathsf{\omega }}_{\mathrm{ih}}{\mathrm{I}}_{\mathrm{h}}-\mathsf{\epsilon }.{\mathrm{I}}_{\mathrm{h}}+\mathrm{D}{2}_{\mathrm{ih}}.{\mathrm{I}}_{\mathrm{h}}\underset{\_}{-\mathsf{\rho }.{\mathrm{I}}_{\mathrm{h}}-\mathsf{\tau }.{\mathrm{I}}_{\mathrm{h}}}$$

(17)

In this research, an effective control model is considered as a cost-effective control method within a controllable infection size. The objective function that aims to minimize cost in the control strategy $\mathrm{J}=\left\{\mathsf{\alpha }, \mathsf{\pi }, \mathsf{\vartheta }, \mathsf{\rho }, \mathsf{\tau }\right\}$ is obtained by Equation (18). In Equation (18), the parameter ${\mathrm{W}}_1$ denotes the cost for elimination of breeding site, ${\mathrm{W}}_2$ is the cost for sanitation-based works, ${\mathrm{W}}_3$ is the cost for installing UV light trap, ${\mathrm{W}}_4$ is the cost for isolation of infected person, and ${\mathrm{W}}_5$ is the cost of water purification. Each ${\mathrm{C}}_{\mathrm{i},\mathrm{i}\in \left[1,5\right]}$ denotes each control parameter-based cost.

$$\begin{array}{ll}{\mathrm{C}}_{\mathrm{stf}}\left(\mathrm{T}\right)=\min {{\displaystyle \int }}_0^{\mathrm{T}}& [{\mathrm{W}}_1\cdot \mathsf{\alpha }\left(\mathrm{t}\right)\cdot K_1+{\mathrm{W}}_2\cdot \mathsf{\pi }\left(\mathrm{t}\right)\cdot {\mathrm{K}}_2+{\mathrm{W}}_3\cdot \mathsf{\vartheta }\left(\mathrm{t}\right)\cdot {\mathrm{K}}_2+{\mathrm{W}}_4\cdot \mathsf{\rho }\left(\mathrm{t}\right)\cdot {\mathrm{I}}_{\mathrm{h}}\left(\mathrm{t}\right)\\ & +{\mathrm{W}}_4\cdot \mathsf{\rho }\left(\mathrm{t}\right)\cdot {\mathrm{I}}_{\mathrm{h}}\left(\mathrm{t}\right)+{\mathrm{W}}_5\cdot \mathsf{\tau }\left(\mathrm{t}\right)\cdot {\mathrm{I}}_{\mathrm{h}}\left(\mathrm{t}\right)+\frac{1}{2}\cdot {\mathrm{K}}_3]\mathrm{dt}\end{array}$$

(18)

where ${\mathrm{K}}_1={\mathrm{E}}_{\mathrm{f}}\left(\mathrm{t}\right)+{\mathrm{L}}_{\mathrm{f}}\left(\mathrm{t}\right)+{\mathrm{U}}_{\mathrm{f}}\left(\mathrm{t}\right)$, ${\mathrm{K}}_2={\mathrm{S}}_{\mathrm{f}}\left(\mathrm{t}\right)+{\mathrm{C}}_{\mathrm{f}}\left(\mathrm{t}\right)$, and ${\mathrm{K}}_3={\mathrm{C}}_1{\mathsf{\alpha }}^2+{\mathrm{C}}_2{\mathsf{\pi }}^2+{\mathrm{C}}_3{\mathsf{\vartheta }}^2+{\mathrm{C}}_4{\mathsf{\rho }}^2+{\mathrm{C}}_5{ \mathsf{\tau }}^2$.

Due to the fact that the cost function is nonlinear with the infection trends, time and other conditions, the mathematical programming belongs to the differential equation-based nonlinear mathematical programming. When denoting the prices associated with their respective classes by ${\mathsf{\varphi }}_{{\mathrm{E}}_{\mathrm{f}}}, {\mathsf{\varphi }}_{{\mathrm{L}}_{\mathrm{f}}}, {\mathsf{\varphi }}_{{\mathrm{U}}_{\mathrm{f}}}, {\mathsf{\varphi }}_{{\mathrm{S}}_{\mathrm{f}}}, {\mathsf{\varphi }}_{{\mathrm{C}}_{\mathrm{f}}}, {\mathsf{\varphi }}_{{\mathrm{S}}_{\mathrm{h}}}, {\mathsf{\varphi }}_{{\mathrm{I}}_{\mathrm{h}}} and {\mathsf{\varphi }}_{{\mathrm{R}}_{\mathrm{h}}}$, the nonlinear mathematical programming is solved using Pontryagin’s maximum principle and the derivation of necessary conditions. As explained in [24,25], the existence of an optimal control is a sequence of the convexity of the integrand $\mathrm{J}$ with respect to $\mathsf{\alpha }, \mathsf{\pi }, \mathsf{\vartheta }, \mathsf{\rho }, \mathsf{\tau }$; a priori boundedness of the state variables, and the Lipschitz property of the state system with respect to the state variables. The differential equations are obtained by differentiating the Hamiltonian function, as shown in Appendix A. Finally, optimal control values are obtained as shown in Appendix B.

$${\mathsf{\alpha }}^{\mathrm{c}}=\left[\left({\mathrm{E}}_{\mathrm{f}}^{*}{\mathsf{\varphi }}_{{\mathrm{E}}_{\mathrm{f}}}+{\mathrm{L}}_{\mathrm{f}}^{*}{\mathsf{\varphi }}_{{\mathrm{L}}_{\mathrm{f}}}+{\mathrm{U}}_{\mathrm{f}}^{*}{\mathsf{\varphi }}_{{\mathrm{U}}_{\mathrm{f}}}\right)-{\mathrm{W}}_1\left({\mathrm{E}}_{\mathrm{f}}^{*}+{\mathrm{L}}_{\mathrm{f}}^{*}+{\mathrm{U}}_{\mathrm{f}}^{*}\right)\right]/{\mathrm{C}}_1$$

(19)

$${\mathsf{\pi }}^{\mathrm{c}}=\left[\left({\mathrm{S}}_{\mathrm{f}}^{*}{\mathsf{\varphi }}_{{\mathrm{S}}_{\mathrm{f}}}+{\mathrm{C}}_{\mathrm{f}}^{*}{\mathsf{\varphi }}_{{\mathrm{C}}_{\mathrm{f}}}\right)-{\mathrm{W}}_2\left({\mathrm{S}}_{\mathrm{f}}^{*}+{\mathrm{C}}_{\mathrm{f}}^{*}\right)\right]/{\mathrm{C}}_2$$

(20)

$${\mathsf{\vartheta }}^{\mathrm{c}}=\left[\left({\mathrm{S}}_{\mathrm{f}}^{*}{\mathsf{\varphi }}_{{\mathrm{S}}_{\mathrm{f}}}+{\mathrm{C}}_{\mathrm{f}}^{*}{\mathsf{\varphi }}_{{\mathrm{C}}_{\mathrm{f}}}\right)-{\mathrm{W}}_3\left({\mathrm{S}}_{\mathrm{f}}^{*}+{\mathrm{C}}_{\mathrm{f}}^{*}\right)\right]/{\mathrm{C}}_3$$

(21)

$${\mathsf{\rho }}^{\mathrm{c}}=\left[{\mathrm{I}}_{\mathrm{h}}^{*}{\mathsf{\varphi }}_{{\mathrm{I}}_{\mathrm{h}}}-{\mathrm{W}}_4{\mathrm{I}}_{\mathrm{h}}^{*}\right]/{\mathrm{C}}_4$$

(22)

$${\mathsf{\tau }}^{\mathrm{c}}=\left[{\mathrm{I}}_{\mathrm{h}}^{*}{\mathsf{\varphi }}_{{\mathrm{I}}_{\mathrm{h}}}-{\mathrm{W}}_5{\mathrm{I}}_{\mathrm{h}}^{*}\right]/{\mathrm{C}}_5$$

(23)

Therefore, it can be concluded by the standard control arguments involving the bounds on the controls in

Table 4. The optimal control parameter of each strategy.

Control Parameter	The Optimal Value of Control Parameter
${\mathsf{\alpha }}^{*}$	$\{\begin{array}{c}0\\ {\mathsf{\alpha }}^{\mathrm{c}}\\ 1\end{array}$$\begin{array}{c}\mathrm{if}\\ \mathrm{if}\\ \mathrm{if}\end{array}$$\begin{array}{c}{\mathsf{\alpha }}^{\mathrm{c}}\le 0\\ 0<{\mathsf{\alpha }}^{\mathrm{c}}<1\\ {\mathsf{\alpha }}^{\mathrm{c}}\ge 1\end{array}$
${\mathsf{\pi }}^{*}$	$\{\begin{array}{c}0\\ {\mathsf{\pi }}^{\mathrm{c}}\\ 1\end{array}$$\begin{array}{c}\mathrm{if}\\ \mathrm{if}\\ \mathrm{if}\end{array}$$\begin{array}{c}{\mathsf{\pi }}^{\mathrm{c}}\le 0\\ 0<{\mathsf{\pi }}^{\mathrm{c}}<1\\ {\mathsf{\pi }}^{\mathrm{c}}\ge 1\end{array}$
${\mathsf{\vartheta }}^{*}$	$\{\begin{array}{c}0\\ {\mathsf{\vartheta }}^{\mathrm{c}}\\ 1\end{array}$$\begin{array}{c}\mathrm{if}\\ \mathrm{if}\\ \mathrm{if}\end{array}$$\begin{array}{c}{\mathsf{\vartheta }}^{\mathrm{c}}\le 0\\ 0<{\mathsf{\vartheta }}^{\mathrm{c}}<1\\ {\mathsf{\vartheta }}^{\mathrm{c}}\ge 1\end{array}$
${\mathsf{\rho }}^{*}$	$\{\begin{array}{c}0\\ {\mathsf{\rho }}^{\mathrm{c}}\\ 1\end{array}$$\begin{array}{c}\mathrm{if}\\ \mathrm{if}\\ \mathrm{if}\end{array}$$\begin{array}{c}{\mathsf{\rho }}^{\mathrm{c}}\le 0\\ 0<{\mathsf{\rho }}^{\mathrm{c}}<1\\ {\mathsf{\rho }}^{\mathrm{c}}\ge 1\end{array}$
${\mathsf{\tau }}^{*}$	$\{\begin{array}{c}0\\ {\mathsf{\tau }}^{\mathrm{c}}\\ 1\end{array}$$\begin{array}{c}\mathrm{if}\\ \mathrm{if}\\ \mathrm{if}\end{array}$$\begin{array}{c}{\mathsf{\tau }}^{\mathrm{c}}\le 0\\ 0<{\mathsf{\tau }}^{\mathrm{c}}<1\\ {\mathsf{\tau }}^{\mathrm{c}}\ge 1\end{array}$

A *: optimal value of parameter of A.

as follows:

5. Simulation and Analysis of Control Model Considering Multiple Disease Carriers

With the provided optimal control frameworks, numerical simulations are performed using the parameter values given in Table 1 and Table 3. Several control scenarios are considered with combinations of intervention strategies for non-adult flies, adult flies, and infected humans, as summarized in

Table 5. Scenarios with combinations of control strategies.

Scena-Rio	Strategy\Control Parameters	Egg Flies	Larva	Pupa	Susceptible Flies	Carrier Flies	Infected Human
-	Initial condition (without controls)	126	70	51	163	273	2795
I	Elimination of breeding site $\left(\alpha =0.03\right)$	5	0	0	133	223	1991
II	Sanitation $\left(\pi =0.1\right)$	103	57	42	27	4	1837
III	Installation of UV light trap ($\vartheta =0.04$)	69	38	28	4	33	1739
IV	Good personal and food hygiene ($\rho =0.01$)	103	57	42	133	223	1599
V	Water purification ($\tau =0.02$)	103	57	42	133	223	1641
VI	Combination of I, II and IV ($\alpha =0.03,\pi =0.1,\rho =0.01$)	50	13	10	30	50	1251
VII	Combination of I, II and V ($\alpha =0.03,\pi =0.1,\tau$=0.02)	23	13	10	30	50	1480
VIII	Combination of I, III and IV ($\alpha =0.03,\vartheta =0.04,\rho =0.01$)	11	7	5	15	21	1182
IX	Combination of I, III and V ($\alpha =0.03,\vartheta =0.04,\tau =0.02$)	17	10	7	22	36	961
X	Combination of I-V ($\alpha =0.03,\pi =0.1,\vartheta =0.04,\rho =0.01,\tau =0.02)$	0	0	0	0	0	877

. Table 5 shows the highest number of flies and humans using 50 days’ of simulation.

Each control is applied and is combined with more than one control to examine the impact on the human population. The simulation results are shown in Figure 4.

Table 6. Cost values for control simulations.

Parameter	Unit Cost ($)
Control costs
Eliminations of fly’s breeding site	100
Sanitation methods	60
Installation of UV light trap	240
Good personal and food hygiene	1138
Water purification	0.46
Medical treatment cost
Hospitalization	207.7

shows the assumed control costs and medical treatment costs. The costs are assigned with the consideration of existing relevant literatures [26,27,28,29,30].

After comparing all of those ten scenarios, the most effective strategy is considered Scenario X, which combines all control method under the assumptions.

Table 7. Result of control scenario.

Control Scenario	Relevant Figure	Cost ($)	Results
I	Figure 4a	413,632	- Maturation rates of eggs, larva, and pupa are decreased since they died before even matured. - The number of eggs of fly, larva, and pupa vanish at t = 42 while the number of infected humans reaches its peak at t = 26 with 1991 people and then started to decrease compared with the initial simulation. - An outbreak at t = 23 along with the increasing number of flies.
II	Figure 4b	381,605	- Sanitation aim to increase the death rate of adult flies, both the susceptible and carrier flies. - Applying this intervention impacts to the number of susceptible and carrier flies significantly decrease to 27 and 4 respectively while for initial simulation, there are 163 susceptible and 273 carrier flies. - At t = 24, the number of infected human starts decreasing after reaching its maximum with 1837 people get the disease.
III	Figure 4c	361,430	- This intervention will trap and catch the adult flies which means the death rate of adult flies is increased. - Thus, the population number of susceptible and carrier flies reaches its lowest at t = 50 with 4 and 33 flies. - The highest number of infected humans is 1739 people which is lower that Strategy II.
IV	Figure 4d	333,250	- By keeping good personal and food hygiene, it does not give opportunity for the pathogen which cause the disease to spread among humans. - Applying this control leads to almost half of reduction number of infected humans at t = 13 with 1373 people compared to the initial simulation
V	Figure 4e	340,836	- Ensure every healthy people to only consume hygiene water is like cutting the chain of the disease to be spreading in the environment. - It will affect the reduction rate of contaminated water to be consumed by human. Therefore, the number of infected people will decrease at t = 14 after experiencing its peak with 1641 people.
VI	Figure 4f	261,131	- Combination of these strategies will tackle each of non-adult flies (eggs, larva, pupa), adult flies, and infected human population. - At time t = 70, the number of infected people is almost zero meanwhile the disease infected the highest number 1251 people at t = 10 before the transmission becomes slower after that.
VII	Figure 4g	307,556	- After performing this scenario, the number of non-adult flies has much reduction at t = 48 which is good for indicating the decreasing of the adult flies’ population as well. - The largest population of infected people is 1.480 at t = 12 then the infection subsided.
VIII	Figure 4h	246,979	- Combination of these strategies could lead to the reduction number of all individual, where at t = 36, the population of flies is almost zero and the highest number of infected humans is 1182 at t = 9.
IX	Figure 4i	199,940	- These combinations provide a good result in decreasing the maturation rate of non-adult flies and increasing the death rate of adult flies which affects to the reduction rate of susceptible status becomes infected status in human. At t = 7, the number of infected humans is 961 people but the population becomes disease-free at t = 50 since no one is infected.
X	Figure 4j	183,691	- The population of non-adult and adult flies are zero at t = 27 meanwhile the highest number of infected people is 877 before reaching the diseases free at t = 43. - Combination of all strategies surely give the best result in terms of number of each individual. - In addition, the scenario achieves the lowest cost among all scenarios.

shows the simulation result for each control scenario I to scenario X.

In terms of cost and the number of infected persons, Scenario X shows the best performance. This indicates that the overall control methods are required for effective disease control and their portions influence its performance. It also denotes the significant reduction in number of infected people from 2795 to 877 as a result of the decreasing number of disease carriers. It means there is a 65.7% reduction in the number of infected people compared to the simulation without applying any control.

The provided framework can be used to check the effectiveness with different assumptions and numerical simulations considering multiple disease carriers.

6. Conclusions

This paper presents the ordinary differential equation-based epidemiologic models to understand diarrhea’s infection dynamics and with the consideration of multiple carriers. Numerical simulation and mathematical analyses were performed to identify the relationship and the status of flies and humans when the disease outbreaks.

Then, each control strategy is applied, as shown in Scenarios I to V, and is combined in some combinations of these strategies to obtain the optimal condition, as shown in Scenarios VI to X. Under a certain simulation scenario, Scenario X, which combines all control methods, is proven to be the most effective strategy in reducing both of the number of flies and the infected population, as it minimizes the spread of the disease.

In the future, further study about broader conditions such as the effect of temperature on the transmission of a disease is required, as the population of housefly reaches its peak when the temperature is warmer.

The proposed framework is considered an effective disease control considering multiple carriers with epidemiologic dynamics and the relevant parameters.