Using mathematical modeling to study the dynamics of Legionnaires' disease and consider management options

1.   Introduction

Legionnaires' disease (LD) was first identified in 1977 in the United States of America (USA) after a large pneumonia outbreak in 1976 ^[1]. It primarily appears in regions with poorly maintained artificial sources of water ^[1]. LD was responsible for 786 cases (37%) of all drinking water related illnesses in the USA from 2015 to 2020. Additionally, it accounted for the vast majority of all hospitalizations and disease induced mortality, with 544 (97%) of all hospitalizations and 86 (98%) of all deaths ^[2]. In the USA, there are approximately 6000 cases reported each year, with an additional 10–15 reported cases per million population each year in Australia and Europe ^[1,3]. However, due to the lack of proper surveillance and detection of LD, these numbers most likely underestimate the true number of cases. In the European Union/European Economic Area (EU/EEA) during 2021, the largest recorded rate of LD was 2.4 cases per 100,000 population ^[4]. While the largest recorded rate of LD by year was more recent (2021), the largest LD outbreak to date was in Murcia, Spain, in July 2001, with 449 total confirmed cases and over 800 suspected cases ^[5]. In this outbreak, the case-fatality rate was 1%. In this paper, we focus on a 15-day subset of this data from June 25 to July 9 ^[6]. During this period, there was a significant spike in the number of confirmed cases by the date of onset symptoms, which rapidly declined after July 6. In the discussed LD outbreak, the likely source of infection was from a cooling tower in the vicinity. Outbreaks such as these demonstrate the need to have proper ways to control an LD outbreak if a body of water is contaminated. This signifies the importance of saniizing of large water sources, as the common sources of LD include artificial sources of water such as cooling towers or potable water ^[7].

Legionella is the bacteria responsible for LD, which is a pneumonic disease. It is transmitted through inhaling water droplets containing Legionella pneumophila ^[7]. The most frequent cause of legionellosis, L. pneumophila, is a specific species of Legionella, which causes about 90% of all reported cases of legionellosis in the USA ^[8]. Legionella is also known to cause Pontiac fever, which is a nonpneumonic disease. While its nonpneumonic counterpart Pontiac fever can go away on its own in most cases, LD requires antibiotics for treatment. The primary symptoms of Legionnaires' include fever, a loss of appetite, headache, shortness of breath, and lethargy, while health complications include shock and multi-organ failure, especially kidney and lung failure ^[1,7]. Additionally, the overall disease-induced mortality rate of LD is 4%–18% through progressive pneumonia ^[9].

While LD surveillance has improved over time, it is still underdiagnosed and underreported. LD cannot be clinically distinguished from other pneumonic diseases ^[8]. Due to this difficulty, it usually requires microbiological testing for confirmation. The most common tests are generally given together: taking a culture from the lower respiratory secretion and a urinary antigen test (UAT) ^[10]. However, this UAT is only 70% accurate for the most common serogroup (serogroup 1) of LD ^[11]. It is important to improve detection since it is so difficult to separate this infection from other variations of pneumonia. This is especially necessary for high risk individuals, which include being 50 years or older, smoking or having chronic lung diseases, or having weakened immune systems ^[7]. From previously reported cases, about 75%–80% of cases are from individuals 50 years or older ^[1]. Additionally, LD is extremely costly for patients. The estimated cost per hospital stay ranges from ＄7950 to ＄149,000 USA dollars, let alone the productivity loss for nonfatal cases requiring hospitalization, making it the second most expensive waterborne disease in the USA to treat ^[9].

There have been quantitative studies performed on LD for different outbreaks. For instance, Smith et al. attempted to use data and clustering analyses to classify sources of Legionella which caused the cases in Genesee County, Michigan ^[12]. Similarly, other researchers used different data analysis techniques on LD data and incidences ^[13,14,15]. Specifically, there was an undergraduate thesis by Miramontes that formulated an initial framework of an ordinary differential equation (ODE) model for LD with high and low-risk individuals ^[16]. The dynamics were based on those of a disease model for a different environmentally transmissible disease by Edholm et al. ^[17].

Due to the lack of mathematical models on LD, we utilized work on similar environmentally transmitted pathogens to inform our model formulation. Another relevant environmentally transmitted disease is Buruli Ulcers (BU), which is a chronic debilitating disease that causes skin damage ^[18]. A significant amount of past research has been done on BU that provides an approach to mathematically model an environmentally transmissible disease. Typically models for BU consider contact with the water environment that contains Mycobacterium ulcerans, which causes BU, or with a vector such as an insect or plant. Some of these models have used ODEs, based on the susceptible-infected-recovered (SIR) individuals framework, but incorporating the environmental transmission or other pathways ^{[17,19,20,21,22]}. Additionally, we consider models of cholera, an environmental pathogen spread to humans by consumption or interaction with a contaminated water source ^{[23,24,25,26]}. The study of cholera gives helpful insights into understanding LD, as it is also a waterborne disease that relies on sanitization methods for control ^[23,24,25].

In this paper, building off of similar mathematical models for BU and cholera, we implement our own susceptible-exposed-infected-recovered (SEIR) individuals ODE-based model to study LD. Using the time series of new cases from the peak of the LD outbreak in Murcia, Spain, in 2001 ^[6], we fit our model based on the cumulative number of infections. For the parameter estimation, we do note that the focus is to capture the overall dynamics so that we can understand how model parameters influence the spread of the disease, and produce a baseline fit rather than creating a method to predict future epidemics. Therefore, based on the results from our parameter estimation, we also perform a global sensitivity analysis for the number of infected individuals. Furthermore, to provide management measures for future outbreaks of LD, we formulate an optimal control problem for our model and conduct numerical simulations of our base strategy, along with varying key parameters.

2.   Model and methods

We developed a novel model for an LD outbreak that considers different compartments of individuals and a water environment containing Legionella. Individuals are divided into four compartments related to different stages of infection: susceptible ($S$), exposed ($E$), infected ($I$), and recovered ($R$), as reflected in Table 1. Additionally, we included a compartment for the concentration of the Legionella pathogen in the environment ($B$), since LD is contracted by exposure and interaction with the environment. Our model dynamics are captured in Figure 1. Individuals are born susceptible at the rate $\pi$ and die naturally at the rate $\mu$. Infected individuals can die from LD-related causes at the rate $\delta$, in addition to the natural mortality rate $\mu$. Susceptible individuals ($S$) may become exposed ($E$) by interacting with a pathogen in the environment ($B$) at the rate $\beta$. After spending an average of $1/\gamma$ days in the exposed state, exposed individuals ($E$) become infected ($I$). Infected individuals ($I$) recover at the rate $\lambda$. The system of ODEs associated with our diagram in Figure 1 is given by the following:

Table 1.  Description of state variables.

Notation	Description
$S(t)$	Number of susceptible individuals at time $t$
$E(t)$	Number of exposed individuals at time $t$
$I(t)$	Number of infected individuals at time $t$
$R(t)$	Number of recovered individuals at time $t$
$B(t)$	Concentration of pathogen in the environment at time $t$

---

 | Show Table

DownLoad: CSV

Figure 1.  Flow diagram for Legionnaire's outbreak, including individuals and environment interaction. The circles represent the human compartments in the model, which flow between the classes defined by the state variables in Table 1 and parameters defined in Table 2. The square represents the concentration of the pathogen in the environment, which interacts with the transmission rate, $\beta$. The associated equations for the model are in Eq (2.1).

DownLoad: Full-Size Img PowerPoint

The parameters are summarized in Table 2, which also includes the values we will use later for the numerical simulations. Note, in Eq (2.1), we include the logistic growth of the pathogen in the environment ($B$) along with natural decay to capture the concept of a carrying capacity for the concentration of the pathogen in the environment. Additionally, our choice was based on models for other environmentally transmitted pathogens, since there are no pre-existing mathematical models for Legionella ^[17,27]. Since there is no evidence of individuals shedding Legionella into the environment, compartment $B$ is not affected by the outbreak. For the interaction of susceptible individuals with the environment, we included a half-saturation function for the concentration of the pathogen in the environment ($B$).

Table 2.  Description of parameters, along with the values used in the numerical simulations. Parameters from the literature have their associated citations, whilst parameters we estimated in Section 3.1 are indicated with an asterisk (^*).

Notation	Description	Value	Reference
$\beta$	Transmission rate from $S$ to $E$	$8.00 \times 10^{-5}$	*
$h_B$	Half-saturation constant of the pathogen	$199.64$	*
$1/\gamma$	Average incubation period in humans	$6$	[28,29]
$\lambda$	Recovery rate from infected individuals	$1/14$	[1]
$r_B$	Growth rate of the pathogen in the environment	$2.52 \times 10^{-5}$	*
$K_B$	Carrying capacity of the pathogen in the environment	$500$	*
$d_B$	Decay rate of the pathogen in the environment	$3.31 \times 10^{-5}$	*

---

 | Show Table

DownLoad: CSV

2.1.   Basic reproduction number

For our model in Eq (2.1), we calculated the basic reproduction number, $\mathcal{R}_0$, using the next generation matrix approach ^[30,31]. We can establish our Jacobian matrices, $F$ and $V$, from new infections $\mathcal{F}$ and transitions $\mathcal{V}$ for states $[E, B]^T$ at the disease-free equilibrium $(S^*, 0, 0, 0, 0) = \left(\frac{\pi}{\mu}, 0, 0, 0, 0\right)$:

Therefore, to find the $\mathcal{R}_0$, we take the following:

Thus, $\mathcal{R}_0 = \frac{r_B}{d_B}$, which signifies that an outbreak will occur when the pathogen growth rate is larger than the natural decay rate, as described in Table 2.

2.2.   Optimal control problem

We formulated an optimal control problem for our model to explore how management solutions could be applied to affect an LD outbreak. Due to the ability of Legionella to survive at temperature ranges between 25 and 42 ℃ ^[32], humans are presumed to be the main cause for LD outbreaks. This information suggests that the primary sources for Legionella growth are from artificial sources of water. Thus, we propose filtering the sources of water individuals interact with as a method of control.

For feasibility purposes, and to cut potential government expenses, we propose to clean water at the individual interaction, $S$ with $B$, instead of the source, $B$. Thus, we introduce the time-dependent control function, $u(t)$, as the percentage ($0\leq u(t)\leq1)$ of cleaned water, which leads to the following:

Next, we construct the following objective functional, which represents the cost of implementing our control procedure that we want to minimize:

where $C_S$ is the cost per new case, $C_u$ represents the cost of water sanitization, $\epsilon_u$ is the nonlinear cost of water sanitization, and $T$ denotes the final time. We set $C_S$ to one, $\epsilon_u$ to 0.00001, and vary the sanitization cost, $C_u$, in our numerical simulations, with the baseline value set to 24. Since the controls, state variables, and derivatives of our state variables satisfy the boundedness, then the standard compactness results infer the existence of an optimal control function ^[24,25,33]. Using the existence result, we apply Pontryagin's Maximum Principle ^[34] to obtain the following optimal control characterization:

The Hamiltonian and adjoint equations constructed from Pontryagin's Maximum Principle are in Appendix A.1.

3.   Numerical simulations

3.1.   Parameter estimation

We used the publicly available data of LD to approximate the recovery rate ($\lambda$) and the average incubation period ($1/\gamma$) in humans. LD treatment lasts about 1 to 3 weeks ^[29], while the incubation period is between 2 to 14 days, usually being 5 or 6 days ^[1,28]. Hence, by averaging each time period, we set $\lambda = 1/6$ and $\gamma = 1/14$ using days as the unit of time.

For the peak of the Murcia, Spain outbreak (June 25 – July 9, 2001), there was only one death due to the disease, which proved difficult to capture with the ODE model dynamics, thus resulting in a small disease-induced mortality rate. Hence, we elected to not include this rate for the numerical simulations (i.e., $\delta = 0$). Because the simulation period was very brief (15 days), we set the natural birth and deaths rates to zero ($\pi = \mu = 0$). We used a report on the dataset related to this outbreak in which the data was given as the number of cases of LD on each date where the onset of symptoms occurred ^[6]. Since the data was given in this manner, fitting our parameters to the number of cumulative cases seemed more logical, as we have data on the growing total number of cases, which is updated each day. Based on this and the amount of data, we focused on performing a general fit as a starting point for the model, understanding that we require additional information in terms of data yet to be produced for a more precise fit. To approximate the rest of the parameters, we used MATLAB's $\texttt{fmincon}$ function from the optimization toolbox to minimize the $L^2$ norm of the distance between the simulated results and the data at the corresponding time points, as seen in Eq (3.1) ^[17,35]:

where $z$ is the vector of (unknown) parameters that we are estimating, $MI(z)$ is the vector with the cumulative number of infections in the numerical simulation, and $MI^*$ is the vector with the corresponding data.

In our simulation, we used a total population of $N =$ 360,000, and the initial conditions $E(0) = I(0) = R(0) = 0$, $S(0) = N-I(0)-R(0)-E(0) =$ 360,000, $B(0) = 1000$. The resulting estimates of the parameter values are listed in the third column of Table 2, and the model simulation with these estimates is shown in Figure 2. The resulting normalized difference between the numerical simulation and the data is $0.2071$. From Figure 2, we can see that the general trends of the outbreak are captured, with the individuals moving from the susceptible class, through the exposed and infected, and to recovered. Furthermore, we can see the different dynamics of the exposed and infected classes, which are related to the incubation period and recovery rates. Our environmental class started above the carrying capacity; therefore, it is decaying towards the carrying capacity, during the short window of the outbreak related to the smaller growth and decay values associated with Legionella. Since we are working with a smaller dataset and a short period, our fit does a good job of capturing the general trend of the outbreak. Here, we are working to create a set of baseline parameter values for the model, which we will use to further explore their implications through a sensitivity analysis.

Figure 2.  Time series of number of susceptible individuals ($S$), exposed individuals ($E$), infected individuals ($I$), recovered individuals ($R$), and the concentration of pathogen in the environment ($B$) from model, Eq (2.1), with the estimated parameter values listed in Table 2 (blue curves). The data used to estimate those parameter values are shown as orange dots in the bottom right panel for the cumulative number of infections, $\mathcal{I}(t)$.

DownLoad: Full-Size Img PowerPoint

3.2.   Sensitivity analysis

We are interested in understanding how model factors, $\theta = (\beta, h_B, \gamma, \lambda, r_B, K_B, d_B)$, affect the number of infections, $I$. We performed a global sensitivity analysis of the model in Eq (2.1) using Latin Hypercube Sampling (LHS) to sample the parameter space and the Partial Rank Correlation Coefficient (PRCC) to evaluate the sensitivity of the model output, $I(t)$, to variations in the model factors ^[36]. Each model factor is sampled 2000 times from a uniform distribution from $50\%$ to $200\%$ of its value listed in Table 2. As shown in Figure 3, all the parameters associated with individuals ($\beta, \gamma, \lambda, h_B$) were statistically significant, while the environmental parameters ($r_B, K_B, d_B$) were not; this was likely due to the assumption that the growth and decay of the pathogen in the environment were unaffected by infected individuals. Specifically, the transmission rate ($\beta$) and the exposure rate ($\gamma$) were positively correlated with the infected population $I$, while the recovery rate ($\lambda$) and the half-saturation constant of $B$ ($h_B$) were negatively correlated with the infected population $I$. These correlations make sense, since increasing the transmission or exposure will result in additional infections, while decreasing the and the half-saturation means that we have less infected individuals at a given time due to a quick recovery and the diminished infection from the environment. Two parameters of specific interest to us were $\beta$ and $h_B$, since they are related in the model dynamics for the exposure of individuals to the pathogen, as seen in Eq (2.1). The magnitude of the PRCC values for $\beta$ and $h_B$ indicate the relative importance of $\beta$ in this relationship, and leads to our further exploration of this parameter with the optimal control problem in Section 3.3.

Figure 3.  Global sensitivity analysis results using PRCC and LHS to examine the correlation between the number of infections, $I$, and the model parameters, with PRCC values in the left panel and $p$ values in the right panel. From the p-values, we have that parameters $\beta, \gamma, \lambda, h_B$ were statistically significant, with $\beta$ and $\gamma$ being positively correlated with the number of infections and $\lambda$ and $h_B$ being negatively correlated.

DownLoad: Full-Size Img PowerPoint

3.3.   Management method results

In Section 2.2, we introduced our objective functional in Eq (2.3) and our control vector $u(t)$ for our defined optimal control problem. The simulations were performed over $T = 20$ days to correspond with the data from our given outbreak and we set $\epsilon_u = 0.00001$. We used the forward-backward sweep method outlined in Lenhart and Workman ^[37] to find the optimal control vector to minimize our objective functional through numerical simulations in MATLAB. Initially, we ran the simulation with our established baseline values for the costs and parameters, with the result shown in Figure 4. We observed that the control starts at the peak allotted value for the percentage of cleaned water, assumed at 0.5; then, this transitions to no control, which makes sense with our optimal control problem formulation. Compared with Figure 2 without a control, we can see the effect of applying a control in Figure 4 especially on the exposed and infected compartments in the shape of the curves. We note that there is no change in the concentration of the pathogen in the environment $B$, as the proposed control method of cleaning the water from an individual interaction does not interact with the environment.

Figure 4.  Numerical simulation results for the model in Eq (2.2) with the base value of $C_u = 24$, including the number of individuals that are susceptible ($S$), exposed ($E$), infected ($I$), and recovered ($R$), the concentration of pathogens in the environment ($B$), and the optimal level of water sanitization effort ($u$) over time. In comparison with Figure 2, we can observe the effect of control implementation, specifically on the exposed individuals.

DownLoad: Full-Size Img PowerPoint

We decided to explore how varying the cost of a applying sanitization control would affect the outcomes by changing the values for $C_u$. This would mean that the cost to apply sanitization is either lower or higher than the baseline value that was initially assumed. The results from different variations of $C_u$ are shown in Figure 5. We varied $C_u$ around the baseline value of 24, thereby going down to 23.98 (Variation 1) and up to 24.04 (Variation 7) by even increments. The first variation, Variation 1, results in a constant control ($u(t) \equiv 0.5$), since we have a lower cost, and will implement control throughout the 20-day outbreak. Meanwhile, our final variation, Variation 7, results in no control being applied ($u(t) \equiv 0$) since the cost of control is too high for implementation. Additional simulations confirmed that these behaviors persisted for $C_u$ values below Variation 1 and above Variation 7. The results for Variation 2-6 show a control that is "turned on" for a portion of the outbreak, based on the cost of the control. Our base value of 24 for $C_u$ is reflected in Variation 3, and we can see the cascade effect surrounding this as the increase and decrease in $C_u$ result in different lengths of control implementation before dipping down to zero.

Figure 5.  Numerical simulation results for the model in Eq (2.2) with different variations of $C_u$, the cost of water sanitization, including the number of individuals that are susceptible ($S$), exposed ($E$), infected ($I$), and recovered ($R$), and the optimal level of water sanitization effort ($u$) over time. The variations (1 to 7) correspond to increasing values of $C_u$, which demonstrate the effectiveness of the control based on cost, with a lower cost (Variation 1) leading to a longer control implementation than the high cost (Variation 7) leading to no control. Note that Variation 3 corresponds to the baseline case shown in Figure 4.

DownLoad: Full-Size Img PowerPoint

Based on our parameter estimation and sensitivity analysis, we wanted to explore the $\beta$ parameter in relation to the optimal control problem. Specifically, from the sensitivity analysis, we know that the number of infected individuals is sensitive to changes in the $\beta$ parameter, which we estimated and is related to our control function. We varied the values of $\beta$ up and down by 0.02% of the estimated value listed in Table 2 in even increments, and the results are shown in Figure 6. We observed a similar behavior as seen in Figure 5, though there is a different spacing to the state variables and control in these simulations. As with different variations of $C_u$, we observed how sensitive the control $u(t)$ is to changes in the value of $\beta$, thus highlighting the need for precise knowledge of the transmission rate. Specifically, between Variations 7 and 1 in Figure 6 for the exposed and infected individuals, the difference is around 100 individuals, which would drastically change an outbreak. Furthermore, we observed that although the evenly spaced variations of $\beta$ led to evenly spaced changes in the control, $u(t)$, they did not result in the same spacing for the state variables $S$, $E$, $I$, and $R$. For instance, the difference between Variations 1 and 2 resulted in close values for the state variables for humans. As with varying $C_u$, when we explored values for $\beta$ below Variation 1 and above Variation 7, we observed the same constant control behavior of either at the cutoff ($u(t)\equiv 0.5$) or staying at zero.

Figure 6.  Numerical simulation results for the model in Eq (2.2) with different variations of $\beta$, including the number of individuals that are susceptible ($S$), exposed ($E$), infected ($I$), and recovered ($R$), and the optimal level of water sanitization effort ($u$) over time. The Variations (1 to 7) correspond to increasing values of $\beta$, which demonstrate the effectiveness of the control based on different values of the transmission rate. When using the lowest $\beta$ in Variation 1, we apply no control due to the reduced transmission, while a high transmission in Variation 7 leads to control implementation throughout. Note that Variation 5 corresponds to the baseline case shown in Figure 4.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

We have formulated a novel model to capture the dynamics of LD for an outbreak. Our model formulation captures the interaction between individuals and a water environment containing the Legionella pathogen. For our model, we calculated the basic reproduction number and formulated an optimal control, which includes the management strategy of sanitization of individual water sources to control the spread of LD. We performed numerical simulations of our base model and optimal control problem, using parameters we estimated from a specific LD outbreak in Murcia, Spain.

From our parameter estimation, we found a baseline fit of our model to an LD outbreak that lasted for about 15 days, thus capturing the basic dynamics of an outbreak in a closed population. We decided to perform a sensitivity analysis using the LHS/PRCC global sensitivity analysis to assess the importance of various model parameters on the disease burden. From the analysis, we found that the number of infected individuals was highly sensitive to the transmission rate $\beta$, incubation period $1/\gamma$, and recovery rate $\lambda$. Meanwhile, the half-saturation rate for the environment, $h_B$, was statistically significant but had a very small PRCC value, which is interesting since both $\beta$ and $h_B$ are related in our model dynamics within Eq (2.1). As our optimal control problem is based on cleaning water from the individual sources, we observe that those parameters related to the pathogen in the environment ($r_B$, $d_B$, and $K_B$) are not all statistically significant. From these analyses, we see the value of better surveillance and data when it comes to LD outbreaks, alongside capturing the overall dynamics of the Legionella pathogen to incorporate the parameters related to the pathogen in the environment into our model. Based on our findings, we suggest exploring methods of decreasing the rate of transmission $\beta$, as it is a key factor that affects the change of infected cases and is one of the parameters that we can reasonably modify through sanitization. Since the primary mode that individuals contract LD is through the inhalation of this waterborne pathogen suspended in aerosolized water, the most effective way of reducing the transmission potential is through reducing the interaction with known contaminated water sources. This led to the idea of implementing sanitization of the individual sources of water in our optimal control problem.

In our optimal control problem, we used our estimated parameters to consider the implementation of sanitization management on an LD outbreak, with the baseline simulation shown in Figure 4. We noted that with sanitization, the number of exposed and infected individuals at the end of the 20-day simulation was reduced significantly. To further study the implications of management strategy accounting for the cost, we varied the cost of water sanitization, $C_u$, to explore its impact on the outcomes, as shown in Figure 5. Additionally, we explored the impact of varying the transmission rate $\beta$ on the management strategy given the sensitivity analysis results, the results of which are summarized in Figure 6. In both scenarios, we noted the drastic changes in the number of exposed and infected individuals with slight changes in $C_u$ or $\beta$. This was as expected, as we found that the transmission rate was highly positively correlated with the number of infected individuals. Alternatively, when we had a low cost of control, we applied control for longer lengths of time, as opposed to when the cost was increased. In general, we observed that as we had a less expensive sanitization method (i.e., low value for the cost $C_u$) or a more difficult outbreak to work with (i.e., high transmission rate $\beta$), we needed to implement our control longer to control the outbreak.

Due to the limitation of surveillance on LD, there is little publicly available high-quality time series data to be used to guide the model formulation. Therefore, the data used for our parameter estimation was limited. Specifically, in the 15-day LD outbreak in Murcia, Spain ^[6], there are 15 data points for us to use to estimate 5 parameters in our model, which was not ideal. Additionally, since the reported data is the number of individuals on the dates where the onset of symptoms occurred, we were restricted to perform the parameter estimation against the cumulative number of infections. Thus, instead of aiming for a perfect fit, we used those estimated parameter values as a baseline poised at starting the mathematical analysis and gaining insight. Because many countries lack appropriate diagnostic, surveillance, and/or reporting systems ^[38], additional computational mathematical modeling and data sciences techniques could be utilized to gain insights of the true incidence of LD from the limited data available. Additionally, since there has been little to no previous effective methods of control for LD, and there is currently no known safe level of Legionella in water systems, we do not have a definite number for $u_{\max}$, which is the maximum feasible amount of sanitization percentage that we think can be cleaned in every source of water. For our model, we elected to use a higher $u_{max}$ in relation to cholera, though this can be further explored with additional knowledge about the feasibility of amount of water sanitized in relation to Legionella ^[24]. Additionaly, there are limitations on the way we model the water environment dynamics; therefore, future work should include capturing the Legionella concentration in the water, $B$. Refining the dynamics for $B$ will aid in better modeling the transmission, and potentially leading to preventative strategies for LD outbreaks ^[8].

For future directions, we are considering a variety of additional factors to calibrate our model. For instance, although Legionella can live in all types of water, they can thrive at 25–42 ℃, with the optimal temperature for its growth being 35 ℃ ^[8]. Datasets with a monthly occurrence indicate that seasonality to be an important epidemiology feature of LD, with more cases reported during the warm season and spike in summer ^[32,39]. With the trend of global warming, it would be important to explore the potential impact of climate change on LD outbreaks. Although LHS/PRCC provides information on the direction (positive or negative) of the association of model factors to the model output (number of infections in our case), it does not distinguish the contribution of single and combined factor interactions. We plan to perform a variance based global sensitivity analysis on these epidemic dynamics with the Sobol method to examine the first-order and total-order effects. While we are able to capture the dynamics of an outbreak in a closed population, our next steps include generalizing and fine-tuning our model for future outbreaks. Some potential sources of calibration include creating an age-structured model that considers individuals over the age of 50 high risk or a sex-based model that incorporates the higher risk for men. Another direction would be investigating the health inequities specific to LD, as it is known that there was a disproportionate number of infections of black individuals compared with white individuals during an LD outbreak in 2018^[40].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

We want to thank Kimi Marimontes (Scripps College 2022) and Kieran Saucedo (Harvey Mudd College 2026) for their contributions to this work in the preliminary stages.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A1.   Details on optimal control formulation

We construct the Hamiltonian as follows:

and obtain our adjoint equations as follows:

Finally, we solve the optimality condition as follows:

where $\xi = S\beta\frac{B}{h_B+B}$.