Genetic algorithm–optimized SEIR-V: robust parameter estimation for wastewater-based infectious disease modeling across multiple regions

1 Introduction

Since its initial report in Wuhan, China, in December 2019, the COVID-19 pandemic caused by the novel coronavirus (SARS-CoV-2) has evolved into a global health crisis, severely impacting public health systems and the global economy. This unprecedented event underscored the urgent need for advanced modeling approaches to better understand the transmission of respiratory infectious diseases. Developing sophisticated modeling tools is therefore essential, not only to improve preparedness for future outbreaks but also to support data-driven public health interventions during pandemics.

During the COVID-19 pandemic, traditional surveillance systems, primarily dependent on clinical testing and case reporting, revealed significant limitations in timeliness and coverage. In response, wastewater-based surveillance (WBS) has emerged as a transformative public health innovation (Purves et al., 2024), offering a rapid and cost-effective method to monitor and predict epidemic dynamics. By quantifying SARS-CoV-2 viral particles in wastewater, WBS provides an early indicator of community-level infection trends. Multiple studies have confirmed that wastewater viral signals often precede reported clinical cases by several days, making WBS a powerful tool for early warning and outbreak prediction (Wu et al., 2022; Safford et al., 2022).

In recent years, several computational and mechanistic approaches have been applied to improve the predictive performance of wastewater-based epidemiology (WBE). Traditional Susceptible–Infectious–Recovered (SIR) and extended SEIR-Virus (SEIR-V) models (Bivins et al., 2020) have been adapted to incorporate environmental factors such as wastewater temperature, improving their realism and accuracy. Simulating viral loss in sewer systems has further enhanced these models. Machine learning strategies and optimization algorithms, including Particle Swarm Optimization (PSO) (Hart and Halden, 2020), have been employed to fit epidemiological parameters more effectively (Randazzo et al., 2020). Additionally, studies have examined the limitations of exponential series solutions for SIR/SEIR models (Saguti et al., 2021), and a fractional-order vaccine controller based on SEIR dynamics has been proposed to simulate intervention effects. Collectively, these efforts have advanced WBE-driven prediction models toward more accurate epidemic simulations.

WBS offers several well-documented advantages: it captures shedding from asymptomatic and pre-symptomatic individuals (Phan et al., 2023a), requires far fewer resources than mass individual testing, and reflects external influences such as vaccination coverage, mask mandates, and isolation policies (Gondim, 2021; Schmitz et al., 2021; Maged et al., 2023; Amrollahi et al., 2022), often anticipating reported outbreak trends by several days at low operational cost.

Despite its potential, WBS still faces several challenges: (1) the instability of sewage environments, which are sensitive to temperature and other external factors; (2) the uncertain correlation between viral load in wastewater and actual infection counts; and (3) the difficulty of translating wastewater data into actionable epidemiological insights (Fazli et al., 2021; Bertels et al., 2022). To address these limitations, researchers have refined compartmental models such as the Susceptible–Exposed–Infectious–Recovered–Virus (SEIR-V) model (Taghizadeh and Mohammad-Djafari, 2022), which explicitly links viral concentration in wastewater to infection dynamics. However, these models rely heavily on accurate parameter estimation. Traditional optimization techniques, such as least squares, often struggle with the system's non-linearity and high dimensionality, leading to convergence on local optima, poor parameter identifiability, and reduced predictive accuracy. This parameter estimation bottleneck remains one of the greatest barriers to improving WBS-based epidemic modeling (Lai et al., 2023; McMahan et al., 2021).

To overcome these challenges, we propose a Genetic Algorithm–optimized SEIR-V (GA-SEIR-V) framework for robust, multi-parameter optimization within a temperature-aware modeling environment. Unlike conventional methods, Genetic Algorithms (GAs) are well-suited to exploring complex, high-dimensional search spaces and identifying globally optimal parameter sets that reflect regional differences and policy effects. By integrating temperature-dependent viral decay with a global optimization engine, our model provides a more realistic representation of disease transmission dynamics and viral persistence in wastewater (Wang et al., 2024).

We do not claim genetic algorithms or temperature-dependent viral decay as novel in isolation: GA optimization has previously been coupled with compartmental models (Wang et al., 2024), and temperature-driven RNA decay in wastewater is well-established (Weidhaas et al., 2021; Wurtzer et al., 2021; Ahmed et al., 2020). Our contribution is their targeted integration and validation for wastewater-based SEIR-V forecasting, along three concrete axes. First, whereas Wang et al. (2024) couple a GA with a socio-psychological SEIR model calibrated on clinical case counts in a single setting, we calibrate a wastewater-observed SEIR-V directly on viral-load signals and use the GA specifically to resolve the λαβE(0) collinearity that destabilizes least-squares estimation on noisy WBE data. Second, whereas temperature-aware analyses such as Proverbio et al. (2022) retain conventional estimation, we embed a dynamic, sinusoidally-fitted decay rate γ(t) inside a global evolutionary optimizer. Third, we show that this single, unmodified pipeline generalizes across three regions with contrasting epidemic profiles—unimodal, large-scale multi-wave, and Omicron-era multimodal—without per-region restructuring, which we present as practical evidence of robustness rather than a claim of universal superiority.

2 Problem statement

Accurate and reliable epidemic forecasting using wastewater-based epidemiology (WBE) data presents two fundamental and interrelated challenges that this study seeks to address.

The first challenge concerns parameter estimation in compartmental models. Standard SEIR-type models depend on the accurate specification of multiple interdependent parameters including transmission rate (λ), viral shedding rate (β), and fecal load (α), whose values are not fixed but vary across regions, time periods, and public health contexts. Traditional optimization techniques, such as least squares and gradient descent, are ill-suited to this task because they are highly sensitive to initial conditions, prone to converging on local optima, and struggle to navigate the complex, non-linear, and high-dimensional parameter spaces that characterize real-world epidemiological systems. This results in poor parameter identifiability, where multiple parameter combinations may fit the data equally well, and in models that fail to reliably reproduce observed transmission dynamics or generalize across different outbreak settings (Poonia et al., 2022; Phan et al., 2023b).

The second challenge involves the environmental sensitivity of wastewater viral signals. Raw SARS-CoV-2 RNA concentrations measured in wastewater are inherently noisy and subject to significant environmental influence. Among all environmental factors, temperature is the most critical determinant of viral RNA decay rate in wastewater systems. As temperature fluctuates seasonally and regionally, the degradation rate of viral particles changes substantially, directly affecting the detectable viral concentration at the sampling point. When this temperature dependence is ignored—as in most standard SEIR-V implementations—the model is misspecified, leading to systematic discrepancies between simulated viral concentrations and clinically observed infection counts. This limitation is particularly acute when modeling data from regions with distinct seasonal patterns or across multi-year datasets spanning different climatic conditions.

Together, these two challenges—unreliable parameter optimization and environmentally confounded wastewater signals—form the core motivation for developing the GA-SEIR-V framework proposed in this study.

2.1 Motivation

The SEIR model is a simple, mechanistic, and idealized framework for understanding epidemic spread. Although it provides a fundamental basis for modeling infectious diseases, its predictive performance depends heavily on the accurate specification of parameters. These parameters are not static; rather, they represent complex, dynamic processes influenced by public health interventions, population behavior, and environmental conditions that vary widely across regions. Such variability makes it impossible to define a single, universally applicable parameter set in advance. A critical bottleneck in applying SEIR-based models, particularly to complex data sources such as wastewater surveillance, is parameter optimization. Conventional algorithms often converge to local optima and fail to explore the full parameter space, resulting in solutions that fit portions of the data but fail to capture the underlying transmission dynamics. This limitation greatly reduces the model's reliability and predictive accuracy.

These challenges motivate our adoption of the Genetic Algorithm (GA) approach. We sought a robust optimization strategy capable of effectively exploring high-dimensional parameter landscapes to identify globally optimal solutions where traditional methods fail. By developing the GA-SEIR-V(w) model and validating it across three distinct regions, our goal is to demonstrate a more generalized and accurate modeling approach. This framework not only yields realistic parameter estimates and improved forecasts but also provides deeper insights into the regional drivers of epidemic dynamics.

3 Materials and methods

3.1 SEIR-V model review

The Susceptible–Exposed–Infectious–Recovered–Virus (SEIR-V) model is an extension of the classical SEIR compartmental framework, specifically adapted to incorporate wastewater-based epidemiology (WBE) data. Unlike standard SEIR models that rely exclusively on clinically reported case counts, the SEIR-V model establishes a quantitative mechanistic link between the viral RNA concentration measured in wastewater and the number of actively infected individuals in the community. This linkage makes the SEIR-V framework particularly valuable for early outbreak detection, as wastewater viral signals typically precede clinical case reports by several days (Purves et al., 2024; Wu et al., 2022).

The model partitions the total population into four epidemiological compartments: Susceptible (S), Exposed (E), Infectious (I), and Recovered (R), along with a wastewater viral load variable (V) that serves as an observable proxy for community infection burden. The population dynamics are governed by the following system of ordinary differential equations (ODEs):

S′=-λIS,(1)

E′=λIS-kE,(2)

I′=kE-δI,(3)

V′=αβ(1-γ)I.(4)

The biological interpretation of each compartment and parameter is as follows. Susceptible individuals (S) become exposed upon contact with infectious individuals at a transmission rate λ, which encapsulates the combined effect of contact frequency and infection probability per contact. Exposed individuals (E) progress to the infectious state at rate k, where 1/k represents the mean incubation period. Infectious individuals (I) recover at rate δ, where 1/δ represents the mean infectious duration, and transition into the recovered compartment (R). Under the assumptions of this study, the recovered population will not become a new source of the virus or affect virus transmission; therefore, we do not include the recovered population when constructing the model. The model assumes a closed population operating on a daily time scale, thereby excluding births, deaths, immigration, and emigration.

The wastewater viral load variable (V) accumulates as infectious individuals shed virus-containing fecal matter into the sewer system. The fecal excretion rate α (grams per person per day) and the viral shedding concentration β (viral copies per gram of feces) together determine the daily input of viral RNA into the wastewater. However, viral particles undergo degradation within the wastewater environment due to physical, chemical, and biological factors. The parameter γ represents the fractional daily loss of viral RNA in wastewater, so that only the fraction (1−γ) of shed virus remains detectable at the sampling point (Saththasivam et al., 2021). Solving this system of ODEs provides estimates of the time-varying infectious population I(t) and the corresponding wastewater viral concentration V(t), the latter serving as the primary observable for model calibration.

A critical limitation of the standard SEIR-V formulation, however, is its reliance on a constant viral decay rate γ. In reality, viral RNA degradation in wastewater is strongly temperature-dependent: higher temperatures accelerate RNA hydrolysis and microbial activity, leading to faster decay, while lower temperatures slow degradation and prolong viral detectability (Safford et al., 2022). Treating γ as a fixed constant therefore introduces systematic model misspecification, particularly in datasets that span multiple seasons or geographic regions with distinct climate profiles. This limitation directly motivates the temperature-aware extension described in Section 3.3. A complete summary of all model parameters, their definitions, assigned values, and sources is provided in Table 1.

Table 1. Description and values of parameters used in the SEIR-V and GA-SEIR-V models.

This table summarizes all parameters and variables employed in the SEIR-V model framework, including their definitions, assigned or fitted values, and literature sources.

The symbol “–” indicates that the parameter value was not taken from literature but was estimated from regional data using the genetic algorithm optimization. These parameters include the transmission rate, fecal load, and initial exposed population, which vary by region and were fitted to match the wastewater and case data for each study area.

3.2 GA-SEIR-V framework

Parameter estimation in the SEIR-V model is a non-trivial optimization problem. The parameters λ (transmission rate), α (fecal load), β (viral shedding rate), and E(0) (initial exposed population) are mutually interdependent, meaning that changes in one parameter systematically influence the fitted values of others. This interdependence, combined with the non-linear dynamics of the ODE system, creates a high-dimensional search space with multiple local optima—a landscape that conventional gradient-based methods such as least squares and gradient descent are ill-equipped to navigate reliably. To overcome this limitation, we integrate a Genetic Algorithm (GA) into the SEIR-V framework, creating the GA-SEIR-V model. Once the temperature-dependent decay described in Section 3.3 is incorporated, we refer to the complete model as GA-SEIR-V(w), the framework proposed in this study; this single naming convention is used consistently throughout, with GA-SEIR-V denoting the intermediate (temperature-independent) construction step.

Genetic Algorithms are population-based metaheuristic optimization methods inspired by the principles of biological evolution, specifically natural selection, genetic crossover, and mutation (Holland, 1992). Unlike gradient-based methods that follow a single solution trajectory through the parameter space, GAs maintain a population of candidate solutions that evolve in parallel across generations. This parallel search strategy enables GAs to explore a far broader region of the parameter landscape, reducing the risk of premature convergence to local optima and substantially improving the probability of identifying the globally optimal parameter configuration.

The GA-SEIR-V optimization procedure follows three iterative steps, as illustrated in Figure 1:

Step 1: Population initialization and encoding. Each candidate solution consists of a parameter set {λ, α, β, E(0)} encoded as a genetic string (chromosome). The permissible range for each parameter is defined based on epidemiological priors and literature values (see Table 1). An initial population (100–200 individuals depending on the dataset; see Table 1) is randomly generated within these bounds, providing a diverse starting point that broadly samples the feasible parameter space.

Step 2: Fitness evaluation and selection. Each individual's parameter set is input into the SEIR-V ODE system to simulate the time-varying wastewater viral concentration V(t). The day-to-day incremental change in V(t) represents the model-predicted viral load in wastewater. Model performance is quantified using the Sum of Squared Errors (SSE) between predicted and observed wastewater viral concentrations:

SSEV=∑t=1n(Vtd-Vpred)2(6)

where Vtdis the measured viral concentration in wastewater on day t and V_*pred* is the corresponding model prediction. A lower SSE_*V* indicates a closer correspondence between the model and observed data. Following Darwin's principle of survival of the fittest, parent solutions are chosen by tournament selection (tournament size 2), so that individuals with lower SSE_*V* are preferentially retained for the next generation.

Step 3: Crossover and Mutation. Selected parents undergo blend crossover (BLX-α with α = 0.5), applied with a crossover probability of 0.5, in which each offspring gene is drawn from an interval spanning the two parent values. This continuous recombination combines high-fitness parameter configurations and drives the population toward progressively better-fitting solutions. To maintain genetic diversity and prevent premature convergence, Gaussian mutation (mean 0, standard deviation 1) perturbs each gene with probability 0.1, applied to individuals with a per-individual mutation probability of 0.2. Mutation ensures that the GA continues to explore novel regions of the parameter space throughout the evolutionary process, rather than collapsing prematurely around a suboptimal solution.

Steps 2 and 3 are iterated for a fixed number of generations (25–100 depending on the dataset; Table 1) using the standard generational GA (the eaSimple routine of the DEAP evolutionary-computation library in Python). A hall-of-fame archive retains the best individual encountered across all generations, and the parameter set yielding the lowest SSE_*V* is returned as the optimal solution. This iterative evolutionary process is summarized in Figure 1. The complete GA configuration used for each region is detailed in Table 1.

The SEIR-V ODE system was integrated numerically with SciPy's ode interface (real-valued VODE backend, an adaptive Adams/BDF multistep scheme) advanced at a 1-day output step, matching the daily resolution of the surveillance data; the temperature function T(t) is evaluated at each daily step. All experiments were run with fixed random seeds to ensure reproducibility.

3.3 Temperature-dependent extension [GA-SEIR-V(w)]

Although the GA-SEIR-V framework described in Section 3.2 substantially improves parameter optimization, initial model validation revealed a systematic residual limitation: the model's predictive accuracy degraded during periods of marked seasonal temperature change. This observation pointed to a fundamental biological reality—the degradation rate of SARS-CoV-2 RNA in wastewater is not constant but varies continuously with environmental temperature. Treating γ as a fixed parameter, as in the standard SEIR-V formulation, therefore introduces time-varying model misspecification that GA optimization alone cannot resolve. To address this, we developed a temperature-aware extension of the GA-SEIR-V model, hereafter referred to as GA-SEIR-V(w).

3.3.1 Temperature data and fitting

We first conducted a visual and quantitative analysis of daily temperature records from the three study regions; the fitting time period includes the entire period of the viral outbreak in the study area. Plotting the time–temperature series revealed a clear periodic pattern across all datasets, closely approximating a sinusoidal function (Brozak et al., 2022). We therefore modeled the temporal variation in temperature using the following sinusoidal function:

T(t)=Asin(Bt-C)+D(7)

where T(t) is the estimated temperature (°C) on day t, A is the amplitude of seasonal temperature variation, B is the angular frequency (controlling the period of the seasonal cycle), C is the phase shift (controlling the timing of temperature peaks), and D is the vertical offset representing the mean annual temperature. The four parameters {A, B, C, D} (initialized to {4, 0.2, 0, 15}) were fitted independently for each region using non-linear least squares regression applied to the observed daily temperature records. The resulting fitted values for each region are presented in Table 2.

Table 2. Optimized temperature model parameters for the GA-SEIR-V framework across three study regions.

The table lists the fitted parameter values (A, B, C, and D) obtained from the Genetic Algorithm optimization for California, Greater Boston, and Switzerland. These parameters reflect region-specific differences in seasonal temperature variation used in the GA-SEIR-V simulations.

3.3.2 Integration into the SEIR-V framework

Once the temperature function T(t) was established for each region, we incorporated it into the SEIR-V ODE system by replacing the fixed viral decay parameter γ with a dynamic, temperature-dependent decay function γ(t). Drawing on established relationships between temperature and SARS-CoV-2 RNA decay kinetics reported in the literature (Weidhaas et al., 2021), the decay rate is modeled as an increasing function of temperature, reflecting the empirical observation that higher temperatures accelerate viral RNA degradation in wastewater. By substituting γ(t) into Equation 4, the wastewater viral load equation becomes:

ξ(T) =ln2η0Q10-(T(t)-T0)/10°C,(8)

γ(t)=1-e-ξ(T)tarrive,(9)

V′= αβ(1-γ(t))I.(10)

Following the method used in the literature to calculate the decay function γ(t), we first compute the one-step decay rate ξ(T) of the virus at the corresponding temperature. Here, η₀ denotes the virus half-life at reference temperature T₀, and Q₁₀ represents the rate of change with temperature. In biological systems, Q₁₀ is typically between 2 and 3 (Hart and Halden, 2020; McMahan et al., 2021); in this study, it is set to 2.5 (Reyes et al., 2008). The wastewater treatment (arrival) time t_*arrive* used to compute γ(t) is 18 days.

This formulation allows the model to dynamically adjust the simulated wastewater viral concentration in response to seasonal temperature fluctuations, more accurately representing the real-world environmental conditions under which wastewater samples are collected. The remaining model structure, including the GA optimization procedure described in Section 3.2, is retained without modification. The full GA-SEIR-V(w) workflow—encompassing temperature fitting, dynamic decay rate integration, and GA-driven parameter optimization—is illustrated in Figure 1.

3.3.3 Rationale for temperature as the primary environmental variable

While other environmental factors—including wastewater pH, hydraulic retention time, UV exposure, and flow rate—also influence viral RNA stability, temperature was selected as the primary environmental covariate for two reasons. First, temperature is the most extensively documented determinant of SARS-CoV-2 RNA decay in wastewater systems, with consistent quantitative relationships reported across multiple independent studies (Wurtzer et al., 2021; Ahmed et al., 2020). Second, continuous daily temperature records were available for all three study regions, whereas reliable measurements of other environmental variables were either unavailable or inconsistent across datasets. Future extensions of this framework could incorporate additional environmental covariates as data availability improves, as discussed in Section 5.

3.4 Data sources

Model calibration and validation were performed using publicly available datasets from three geographically and epidemiologically distinct regions: California (United States), Greater Boston (United States), and Switzerland. For each region, three data streams were assembled: (1) daily SARS-CoV-2 RNA concentrations in wastewater, (2) clinically reported daily COVID-19 case counts, and (3) daily ambient temperature records. Using three independent regional datasets with contrasting epidemic profiles and climate conditions was a deliberate design choice intended to evaluate the generalizability and robustness of the GA-SEIR-V(w) framework across diverse real-world settings. Table 3 provides a summary of the key characteristics of each dataset.

Table 3. Summary of wastewater and clinical data used for model evaluation across three study regions.

The table lists the observation period, number of reported infection cases, and corresponding wastewater sample counts for California, Greater Boston, and Switzerland. These datasets were used to calibrate and validate the GA-SEIR-V model under different epidemic and environmental conditions.

3.4.1 California, United States

The California dataset spans September 1, 2020 to November 27, 2022, encompassing 662 days of reported clinical cases obtained from the California Department of Public Health (CDPH) and 78 paired wastewater samples from the WastewaterSCAN network (State Portal of California, 2026). During this observation period, the California COVID-19 epidemic exhibited highly complex, multi-wave dynamics, characterized by immense spatial heterogeneity and substantial environmental noise across a catchment population of approximately 27.7 million. This dataset therefore represents a profoundly large-scale, heterogeneous epidemic profile, providing a stringent stress-test for evaluating model performance under conditions that severely challenge traditional parameter optimization techniques.

3.4.2 Greater Boston, United States

The Greater Boston dataset covers September 30, 2020 to May 24, 2021, comprising 225 days of clinical case data and 245 wastewater samples. This period coincided with the first major autumn–winter COVID-19 wave in the northeastern United States, driven predominantly by early SARS-CoV-2 variants. The epidemic progression in this region follows a clear unimodal pattern—a single rise, peak, and decline—making it the most straightforward dataset for model fitting. This dataset was also used to calibrate the temperature-dependent component of the GA-SEIR-V(w) model, given the availability of high-quality, continuous temperature records for the Greater Boston metropolitan area.

3.4.3 Switzerland

The Switzerland dataset encompasses January 1, 2022 to December 31, 2022, providing 365 days of both clinical and wastewater observations. This period corresponds to the emergence and rapid spread of the Omicron variant (BA.1, BA.2, and subsequent subvariants), which produced multiple overlapping epidemic waves of varying magnitude throughout 2022. Unlike the unimodal Greater Boston dataset, the Swiss data exhibits a multimodal distribution with several distinct peaks, reflecting the complex interplay of variant succession, waning immunity, and evolving public health policies. This dataset presents the most challenging modeling scenario and serves as a critical test of the GA-SEIR-V(w) framework's ability to capture complex, multi-wave epidemic dynamics.

For Greater Boston and Switzerland, all three datasets were used in their entirety for model fitting, with the first 70% of each time series used for model training and the remaining 30% reserved for prospective validation. Due to the specific computational setup for the California dataset, a contiguous 78-day subset was strictly utilized for model training. To ensure data integrity for this large-scale dataset, missing values in the California records were imputed using a Kalman filter, whereas for the other two regions, no imputation was performed for missing data points; gaps in wastewater sampling were handled by interpolating the ODE solution at the observed sampling time points only. This train-test split strategy was applied to ensure fair comparison of predictive performance across all models. To further ensure cross-regional comparability—and to address the differing per-region processing noted above—the supplementary evaluation in Section 4 applies a single unified protocol to every region: calibration on the days up to the split point, followed by the two-window (in-sample/out-of-sample) scoring described there. Regarding missing values, wastewater gaps were linearly interpolated for Switzerland and Luxembourg, whereas the California series was pre-cleaned by Kalman-filter imputation before analysis; no other gap-filling was applied.

3.5 Evaluation metrics

Model performance was assessed from two complementary perspectives: (1) the accuracy of predicted infection counts relative to clinically reported cases, and (2) the precision of predicted epidemic peak timing relative to observed outbreak peaks. Together, these two criteria provide a comprehensive evaluation of both the quantitative fit and the practical epidemiological utility of each model.

3.5.1 Correlation between predicted and reported cases

To quantify the agreement between model-simulated and clinically reported COVID-19 case counts, we employed linear regression analysis. The predicted case counts (ŷ) were regressed against the reported case counts (x) using the following model:

ŷ=b0+b1x,(11)

where b₀ is the intercept and b₁ is the slope, estimated using the ordinary least squares method:

b1=∑i=1n(xi-x̄)(yi-ȳ)∑i=1n(xi-x̄)2, b0=ȳ-b1x̄(12)

where x̄and ȳ denote the sample means of reported and predicted case counts, respectively.

The coefficient of determination (R²) was then computed to quantify the proportion of variance in observed case counts explained by the model:

R2=1-∑i=1n(yi-ŷi)2∑i=1n(yi-ȳ)2(13)

An R² value approaching 1.0 indicates near-perfect explanatory power, while values closer to 0 indicate poor correspondence between predicted and observed trends. Additionally, the Pearson correlation coefficient (r) was calculated to measure the linear association between predicted and reported case counts:

r=∑i=1n(xi-x̄)(yi-ȳ)∑i=1n(xi-x̄)2∑i=1n(yi-ȳ)2(14)

A value of r = 1 indicates perfect positive linear correlation, while r = 0 indicates no linear association. Together, R² and r provide complementary information: R² measures overall model fit, while r captures the directional consistency of predicted trends relative to observed dynamics. In addition, to report error in interpretable case units—as recommended for time-series evaluation—we compute the root-mean-square error (RMSE) and mean absolute error (MAE) between the regression-calibrated model prediction ĉ_*i* = (ŷ_*i*−b₀)/b₁ and the reported case count x_*i*:

RMSE=1n∑i=1n(ĉi-xi)2, MAE=1n∑i=1n|ĉi-xi|.

3.5.2 Epidemic peak timing accuracy

Beyond case count correlation, a critical practical measure of model utility is the accuracy with which the model predicts the timing of epidemic peaks—that is, the day on which the infection rate reaches its maximum. Early and accurate peak prediction is of direct public health relevance, as it enables proactive deployment of interventions before case surges overwhelm healthcare capacity.

Peak timing accuracy was quantified as the absolute difference in days between the model-predicted epidemic peak and the clinically observed peak:

Δtpeak=|tpred-tobs|(15)

where t_*pred* is the day of the model-predicted peak infection count and t_*obs* is the day of the observed peak in reported clinical cases. A smaller Δt_*peak* indicates that the model identifies the outbreak peak closer to the actual event. This metric was computed separately for the GA-SEIR-V(w) model and the traditional least squares–based SEIR-V model, enabling direct comparison of their peak-timing accuracy across the study regions.

3.5.3 Comparative models

To contextualize the performance of the proposed GA-SEIR-V(w) framework, four models were evaluated in parallel across all three datasets. The basic SEIR model uses neither wastewater data nor temperature effects and is fitted by least squares. SEIR-V adds the wastewater compartment under a fixed viral decay rate, again fitted by least squares. SEIR-V(w) further allows the decay rate to vary with temperature. The proposed GA-SEIR-V(w), finally, combines temperature-dependent decay with Genetic Algorithm optimization. This graded sequence isolates the incremental contribution of the wastewater term, the temperature dependence, and the optimizer in turn.

This stepwise comparison is designed to isolate the individual contributions of each model component—wastewater integration, temperature dependence, and GA optimization—to overall predictive performance, allowing a clear attribution of improvement at each modeling stage.

3.6 Analysis of parameters

3.6.1 Temperature parameter

SARS-CoV-2 RNA in wastewater undergoes degradation influenced by several environmental factors, including temperature and propagation time (Saguti et al., 2021). We incorporated the time-varying temperature data (as shown in Table 2 and Figure 1) into the modeling framework to evaluate its effect on model performance. By integrating this temporal variability, the viral degradation rate (γ) becomes a dynamic parameter that changes with both temperature and time. Compared with a constant degradation rate, this temperature-dependent formulation allows the model to better align with observed clinical data, as reflected in a smaller Sum of Squared Errors (SSE). This improvement underscores the critical importance of including temperature variation within the GA-SEIR-V modeling framework to more accurately represent real-world viral persistence in wastewater environments.

3.6.2 SEIR parameters

Based on prior studies, we defined the incubation period, the transition from exposed individuals to infected individuals, as 3 days, represented by k = 1/3. The recovery period for infected individuals was set to 8 days, represented by δ = 1/8 (Safford et al., 2022; Hart and Halden, 2020; Killingley et al., 2022). The initial infectious count was reconstructed from the first wastewater measurement as I(0) = V(0)/(αβ(1−γ(0))) rather than set to zero, so that the observed initial signal is honored; the cumulative viral concentration (V) was initialized from the first available data point and the exposed seed E(0) was fitted. The viral decay rate (γ) varies dynamically with temperature fluctuations in wastewater. Consequently, the key parameters requiring optimization via Genetic Algorithm (GA) are α (fecal load), β (viral shedding rate), and E(0) (initial exposed population).

3.6.3 Genetic algorithm parameters

Within the GA optimization framework, the parameters λ, α, β, and E(0) are jointly encoded as a single individual (chromosome). The population size (100–200) and number of generations (25–100) were set per region (Table 1), ensuring a sufficient search space. Tournament selection (size 2), blend crossover (BLX-α, α = 0.5, applied with probability 0.5), and Gaussian mutation (per-gene probability 0.1, per-individual probability 0.2) were used to preserve population diversity and prevent premature convergence. This configuration maintains the genetic diversity essential for effective evolutionary search and convergence toward a globally optimal solution. A detailed summary of all model parameters is presented in Table 1.

3.7 Model assumptions and their scope

For transparency and to delimit the domain of validity of the framework, we state the modeling assumptions explicitly and note their principal implications. The model treats the population as closed, neglecting births, deaths, and migration; this is reasonable over a single epidemic wave but is a simplification for the 662-day California record, where demographic turnover and mobility occur. It further treats the population as homogeneous and well mixed at the aggregated sewershed level, without resolving spatial heterogeneity or age structure, which is a stronger idealization for large, diverse catchments such as California and Switzerland than for compact ones. The recovered compartment R is omitted on the assumption that recovered individuals neither transmit nor shed detectable virus; this ignores waning immunity and reinfection, which are most relevant to the Omicron-dominated Switzerland dataset and are a likely contributor to the peak-timing offset reported in the Discussion. The shedding rate β is fixed at a literature-derived value because the wastewater signal depends only on the product αβ (Equation 4); fixing β resolves the α–β collinearity and renders α identifiable, at the cost of not capturing inter-individual and variant-specific shedding variation. The initial infectious count is not set to zero but reconstructed from the first wastewater measurement, I(0) = V(0)/[αβ(1–γ(0))], with the exposed seed E(0) fitted, so that the non-zero initial wastewater signal is honored. Temperature, finally, is the sole environmental covariate; pH, UV exposure, flow rate, and hydraulic retention time are omitted owing to inconsistent data availability across regions. The most influential of these assumptions—the fixed β and the omission of waning immunity—are revisited in the Discussion, and the effect of the fixed shedding rate is quantified by the sensitivity analysis reported in the Results.

The most influential of these assumptions (fixed β and omission of waning immunity) are revisited in the Discussion, and the effect of the fixed shedding rate β is quantified by the sensitivity analysis reported in the Results.

3.8 From wastewater concentration to reported cases

The framework operates in two stages. First, the GA calibrates the model by minimizing SSE_*V* (Equation 6) between the simulated and observed wastewater viral series, so that the wastewater signal drives parameter estimation. Second, the calibrated model yields the latent transmission dynamics, from which the predicted number of new clinical cases on day t is taken as the daily influx into the infectious compartment, kE(t) (equivalently, the daily increment of cumulative infections). These predicted daily cases are compared with reported clinical case counts by the linear-regression analysis of Section 3.5. This calibrate-on-V, evaluate-on-cases design provides an out-of-signal check: the clinical case series is never used during optimisation.

The decay function γ(t) (Equations 8, 9) is evaluated using an effective arrival time t_*arrive*, representing the cumulative time over which shed virus is exposed to temperature-dependent decay before measurement. We adopt t_*arrive* = 18 days; we note that this exceeds typical hydraulic residence times in a single sewer network and should therefore be read as a lumped effective parameter (aggregating in-sewer transport, storage, and pre-measurement holding) rather than a literal pipe residence time. Its influence is a natural subject for the sensitivity analysis noted above.

3.9 Well-posedness, stability, and numerical considerations

3.9.1 Positivity and invariant region

The non-negative orthant is forward-invariant: on each boundary face the corresponding derivative is non-negative (e.g., S′|S=0=0, E′|E=0=λIS≥0, I′|I=0=kE≥0), so trajectories starting from non-negative data remain non-negative. Since S′+E′+I′ = −δI ≤ 0, the total N = S+E+I is non-increasing and bounded by N(0); hence all state variables remain in the compact set Ω = {(S, E, I)≥0:S+E+I ≤ N(0)}.

3.9.2 Existence and uniqueness

The right-hand side of Equations 1–4 is polynomial in the state variables and therefore continuously differentiable, hence locally Lipschitz, on Ω. The temperature-dependent decay enters only through T(t) = Asin(Bt−C)+D, which is continuous and bounded, giving a continuous, bounded γ(t)∈[0, 1); the extended right-hand side remains Lipschitz in the state variables and continuous in t. By the Picard–Lindelöf theorem a unique solution exists locally, and because trajectories stay in the compact invariant set Ω, this solution extends uniquely to all t ≥ 0. The temperature term thus does not compromise well-posedness.

3.9.3 Basic reproduction number and equilibrium stability

The disease-free equilibrium is (S, E, I) = (S₀, 0, 0). Applying the next-generation-matrix method to the infected compartments (E, I) gives new-infection and transition matrices F=[0λS000], V=[k0-kδ], so that the spectral radius of FV-1 yields:

R0=λS0δ.(16)

The disease-free equilibrium is locally asymptotically stable when R₀ < 1 and unstable when R₀ > 1. Because the model omits vital dynamics and S(t) is monotonically non-increasing, there is no interior endemic equilibrium; the epidemic exhibits a single burnout per wave, consistent with the compartmental structure adopted here.

3.9.4 Numerical integration

The ODE system is integrated with SciPy's ode interface using the real-valued VODE backend, an adaptive-order, adaptive-step Adams/BDF multistep scheme with built-in local error control (default relative/absolute tolerances); the BDF option provides stability for stiff regimes. The solver is sampled at a 1-day step to match the daily surveillance resolution, and T(t) is evaluated at each step.

3.9.5 GA convergence and complexity

As a stochastic global search over a non-convex landscape, the GA carries no deterministic guarantee of reaching the global optimum; however, elitism via the hall-of-fame archive makes the best-so-far fitness monotonically non-increasing across generations, and empirically the fitness plateaus within the generation budget used here. Rather than the arbitrary SSE < 0.2 criterion, we therefore run a fixed generation budget (Table 1) and select the best archived individual. The cost per run is O(GPC_*ODE*), where G is the number of generations, P the population size, and C_*ODE* the cost of one ODE integration.

3.9.6 Parameter identifiability

The wastewater observable scales with the product αβ (Equation 4), so α and β are not jointly identifiable from V alone; we fix β to break this collinearity, leaving α practically identifiable. The transmission rate λ is constrained by the growth phase and E(0) by the initial signal level, so the reduced set {λ, α, E(0)} is identifiable in practice given a single wavefront, a point we return to when interpreting multi-wave regimes.

4 Results

4.1 Overview of model performance

The graphical results (Figures 2–4) illustrate the modeling and forecasting performance of the proposed framework. In each figure (left), the portion before the green dashed line represents the wastewater data sequence used for model training, while the green dashed line marks the predicted outbreak peak. The red dots correspond to measured viral loads in wastewater. In each figure (middle), the black solid line shows simulated COVID-19 case estimates generated by the model, while the red curve depicts the clinically reported cases. The blue and red dashed lines indicate the respective peak times of both curves. Each figure (right) presents the correlation between simulated and observed case numbers. The best-fitting model parameters are summarized in Table 4.

Table 4. Best-fitting parameter estimates and single-run performance for the GA-SEIR-V(w) model across the three primary regions.

The table lists the virus transmission rate (λ), fecal load (α), viral shedding rate (β), initial exposed population E(0), the Sum of Squared Errors (SSE), and the calibration-window coefficient of determination (R²). These are single best-run values; the run-to-run distribution and a fairly re-tuned least-squares baseline are reported in Table 5.

Across the three regions the GA recovers well-identified parameter sets and reproduces the calibration-window case trajectory closely. As the following subsections make explicit, however, these single-run calibration fits should be read together with the run-to-run distribution, the out-of-sample analysis, and the fair baseline comparison, which qualify how far the strong in-window fit translates into genuine forecasting skill.

4.2 Region-specific results

The region-by-region descriptions below report the single best-run fit within each calibration window (the values in Table 4). The corresponding run-to-run distributions, a fairly re-optimized least-squares baseline, and out-of-sample behavior are given by the unified analysis in the subsections that follow; the two should be read together, as the single-run figures here are illustrative rather than a claim of consistent superiority.

4.2.1 California

In the California model, the calibration window coincided with complex, multi-wave transmission in a large-scale, heterogeneous population. Within this window the GA model fitted the observed trend well, whereas a same-setting least squares fit was markedly poorer. The GA model achieved a calibration-window R² of 0.7979, against 0.4373 for that least squares fit. This illustrates the GA's better in-window fit on this difficult record; its out-of-sample behavior, and a fairly re-optimized least-squares baseline, are examined below.

4.2.2 Greater Boston

The Greater Boston record covers the early pandemic phase and follows a clean, largely single-wave profile. The estimated transmission rate was λ = 8.85 × 10⁻⁸, and the observed peak around January 2021 is recovered by the GA to within a few days. The GA model reached a calibration-window R² of 0.8922; on this simple record a well-tuned least squares fit is essentially as good (0.8748 for the same-setting fit here, and even higher for the fairly re-optimized baseline in Table 5). Greater Boston is therefore best read as a case where global optimization is not needed rather than one where it dominates.

Table 5. Two-window evaluation (mean over 30 independent GA runs).

For each region we report the calibration-window (in-sample) and extrapolation-window (out-of-sample) R² for the GA and for a fairly re-optimized least-squares (LS) baseline (20 restarts, identical protocol). Strong calibration fit does not guarantee out-of-sample skill on the long multi-wave records.

4.2.3 Switzerland

The Switzerland record spans 2022 and is dominated by overlapping Omicron subvariant waves. The GA placed the epidemic peak in January 2022, slightly earlier than observed, with a calibration-window R² of 0.6887, against 0.3083 for the same-setting least squares fit (see Table 5 for a fairly re-optimized baseline). As with the other regions, this single-run in-window comparison should be read alongside the distribution and the fair baseline; the 1-month peak offset reflects the difficulty a single-compartment SEIR structure has with the rapid immune-evasion and reinfection dynamics of Omicron.

4.3 In-sample fit, out-of-sample prediction, and comparison with least squares

To separate calibration quality from genuine forecasting skill, we evaluate every region under a unified two-window protocol over 30 independent GA runs: an in-sample calibration window (used for fitting) and an out-of-sample extrapolation window (never seen during fitting). For a fair comparison, the least-squares baseline is re-optimized from 20 random restarts under the identical protocol, rather than taken from a single initialization. Table 5 reports the mean R² in each window for the GA and this fair least-squares baseline.

Within the calibration window the GA reproduces the originally reported fit—for California, R² = 0.78, matching the 0.7979 reported earlier—and, on the two large or complex records where gradient-based least squares is prone to poor local optima, it is far more reliable (California 0.78 vs. 0.20; Luxembourg 0.35 vs. 0.25). On the single-wave Greater Boston record a well-tuned least squares is already excellent and the GA merely matches it (0.96 vs. 0.97), and on the short Switzerland calibration window neither method yields a meaningful in-sample R². Out of sample, however, long multi-wave records remain hard for both methods, and predictive skill beyond the calibration window is limited: the GA's advantage lies in robust parameter estimation for complex, heterogeneous data, not in universal predictive superiority.

Finally, cross-correlation between the daily wastewater viral-load signal and reported cases shows near-synchronous alignment at the sampling resolution of these datasets (peak correlation at a lag of 0 days for Greater Boston and Switzerland, −3 days for Luxembourg, and −10 days for California). In these records the wastewater signal tracks, rather than substantially leads, aggregated case reports; establishing a longer lead time would require finer-grained or less-aggregated case data.

4.4 Comparison with other optimizers and data-driven baselines

To place the Genetic Algorithm in context, we compared it against three further parameter-estimation strategies within the identical SEIR-V(w) model—Differential Evolution (DE) (Storn and Price, 1997), Particle Swarm Optimization (PSO) (Kennedy and Eberhart, 1995), and Bayesian estimation by Markov-chain Monte Carlo (MCMC) (Hastings, 1970)—and against three purely data-driven baselines that predict cases directly from the wastewater signal without any mechanism [linear regression, XGBoost (Chen and Guestrin, 2016), and an LSTM (Hochreiter and Schmidhuber, 1997)]. All methods use the same two-window protocol; Table 6 reports out-of-sample (extrapolation) R².

Table 6. Out-of-sample (extrapolation-window) R² for every method under the unified protocol.

Top block: three global optimizers and least squares, all estimating the same SEIR-V(w) parameters. Bottom block: data-driven models predicting cases directly from the wastewater signal (no mechanism); these reach in-sample R²≈1.0 but do not generalize. Negative values are worse than predicting the mean. “–” denotes a run that failed to converge.

Bold values indicate the highest R-squared value (best result) for each region in the out-of-sample extrapolation window, showing which method achieved superior prediction performance for that specific region and evaluation metric.

Two clear messages emerge. First, within the mechanistic model the global and Bayesian estimators are essentially interchangeable (GA ≈ DE ≈ PSO ≈ MCMC), and all of them decisively beat gradient-based least squares only on the large, complex California record (calibration-window R²≈0.740.78 vs. 0.20); elsewhere the estimators are close. The advantage therefore stems from using a global or Bayesian estimator rather than from the Genetic Algorithm specifically. Second, the data-driven baselines fit the calibration window almost perfectly (XGBoost and LSTM reach in-sample R²≈0.991.00) but fail to extrapolate: their out-of-sample R² is negative—worse than predicting the mean—on every multi-wave record (e.g., XGBoost −7.3 on Switzerland), and only the short single-wave Greater Boston series is forecastable (≈0.70). The mechanistic SEIR-V(w) model, by contrast, never collapses to strongly negative R²: its physical constraints trade some in-sample fit for markedly more stable extrapolation (Switzerland 0.53 vs. ML's −2 to −7). This is the central justification for a mechanism-based approach with robust parameter estimation over black-box fitting.

Table 6 collects the comparison, and Figure 5 visualizes it: panel (a) shows the four global optimizers landing at nearly the same calibration R², and panel (b) shows the data-driven baselines fitting the training window well but turning sharply negative out of sample.

4.5 Robustness checks

4.5.1 Sensitivity to the fixed shedding rate

Because the wastewater equation enters only through the product αβ, the fecal-load and shedding-rate parameters are not separately identifiable. Sweeping β over a 16-fold range and re-fitting confirms this: the calibration-window R² stays essentially flat while the fitted α scales inversely as 1/β, which justifies the fixed-β modeling choice while treating α as the free shedding-load parameter (Figure 6). The same collinearity makes the results robust to implementation choices: re-running all four regions with the sinusoidal temperature fit of Equation 7, the decay form of Equation 8, and a uniform incubation rate (1/k = 3 days) reproduced every two-window R² to two decimal places, because the genetic algorithm absorbs such changes into the fitted α without affecting fit or prediction.

4.5.2 Peak-timing accuracy

When the epidemic peak falls within or near the calibration window—as in Greater Boston—the model recovers the peak day to within about 5 days; when the dominant peak lies far beyond the short calibration window, as in California, Switzerland, and Luxembourg, peak timing cannot be recovered reliably. Accurate peak-timing prediction is thus contingent on the peak falling within or near the calibration window, consistent with the limited out-of-sample skill reported above. Moreover, most of these records are multi-wave: a prominence-based detector identifies two waves in California and Greater Boston, three in Switzerland, and four in Luxembourg. A single-wave SEIR-V structure can represent only one such peak, so we report the dominant-peak timing here with the explicit caveat that secondary waves lie structurally outside the present model; reproducing a full multi-wave sequence would require the immunity-waning or variant-structured extensions discussed in Section 5.

4.5.3 Expanding-window cross-validation

Finally, we probe whether the model is useful when parameters are updated as new data arrive. In an expanding-window cross-validation, each fold re-fits on all data up to a given day and validates on the next 60 days. Rolling short-horizon skill is substantially higher than the one-shot long-range extrapolation reported above (Figure 7), indicating that the framework is operationally useful for short-range tracking with periodic re-calibration, even though a single long-range forecast from a short calibration window is unreliable.

4.6 Fractional-order extension

Fractional-order compartmental models replace ordinary derivatives with Caputo fractional operators and can, in principle, encode memory effects and non-exponential growth and decay. We implemented a Caputo fractional-order version of the SEIR-V(w) system (Grünwald–Letnikov scheme, with q = 1 recovering the integer-order Euler update) and added the fractional order q to the GA search. Table 7 reports the selected order and the resulting calibration and extrapolation R² against the integer-order model. A fractional-order SEIR formulation has been applied directly to wastewater-based COVID-19 data (Meetei et al., 2023), and fractional dynamic formulations of COVID-19 spread have been proposed more broadly (Sioofy Khoojine et al., 2022).

Table 7. Fractional-order (Caputo) vs. integer-order SEIR-V(w).

For each region we report the fractional order q selected by the GA and the calibration- and extrapolation-window R² for the integer-order (q = 1) and fractional-order models. A memory effect (q < 1) is detected everywhere, but it does not systematically improve fit or prediction.

Across all four regions the GA selected a fractional order below one (q between 0.73 and 0.88), confirming a detectable memory effect in the data. However, the fractional order did not systematically improve either calibration fit or out-of-sample prediction over the integer-order model: the calibration fit was in most cases slightly worse (e.g., California 0.30 vs. 0.73) and out-of-sample prediction was essentially unchanged. We therefore retain the integer-order formulation as an adequate and more parsimonious choice for these datasets, while flagging the detected memory effect as a concrete target for future refinement.

5 Discussion

This study developed and tested a Genetic Algorithm–optimized SEIR-V model with temperature-dependent viral decay [GA-SEIR-V(w)] for wastewater-based epidemic modeling, using three primary regions—California, Greater Boston, and Switzerland—together with a fourth, long multimodal record (Luxembourg) for the robustness analyses. The results indicate that the GA-SEIR-V(w) framework is most beneficial for complex, multi-wave outbreaks—where it improved on traditional least squares–based estimation—while remaining competitive with, though not uniformly superior to, a well-tuned least squares fit on simpler single-wave data, consistent with the run-to-run distribution and fair-baseline comparison reported in Section 4. This delineation supports the value of integrating evolutionary computation with environmentally-aware epidemiological modeling in the regimes where conventional estimation is most fragile.

The most fundamental contribution of this work is the demonstration that parameter estimation—rather than model structure—is the primary bottleneck limiting the predictive accuracy of SEIR-based wastewater epidemiology models. This finding is consistent with observations reported by Phan et al. (2023a), who identified parameter identifiability as a central challenge in SEIR-V modeling using wastewater data, and by McMahan et al. (2021), who noted that conventional optimization methods frequently produce unstable parameter estimates when applied to noisy wastewater surveillance data. By replacing least squares optimization with a Genetic Algorithm, the GA-SEIR-V(w) model effectively navigates the high-dimensional, non-convex parameter landscape that characterizes real-world epidemic systems, yielding parameter configurations that more faithfully reproduce observed transmission dynamics. This is particularly evident in the California dataset, which spans 662 days and encompasses a highly complex, large-scale surveillance environment. Here least squares reached an R² of only 0.4373 (SIR-V baseline), missing most of the structure buried under population heterogeneity and environmental noise, whereas the GA reached 0.7979 within the calibration window—an 82.4% relative gain. Our unified in-sample/out-of-sample analysis (Section 4) corroborates this calibration-window advantage on the large, complex California record—R² = 0.78 for the GA vs. 0.20 for a fairly re-tuned least squares baseline over 30 seeds—while making explicit that out-of-sample forecasting beyond the calibration window is limited for all methods on this multi-year record. The GA's contribution is therefore best understood as more robust parameter estimation on difficult data, rather than a guarantee of long-range predictive accuracy.

The integration of temperature-dependent viral decay represents the second major contribution of this framework. The degradation of SARS-CoV-2 RNA in wastewater is well established as a temperature-sensitive process, with multiple studies reporting substantially faster RNA decay at higher ambient temperatures (Weidhaas et al., 2021; Wurtzer et al., 2021). By modeling the viral decay rate γ as a dynamic, sinusoidally varying function of time rather than a fixed constant, the GA-SEIR-V(w) model more accurately captures the seasonal fluctuation in wastewater viral signals. This is especially important for datasets spanning multiple seasons, such as California (over 26 months) and Switzerland (12 months), where a fixed decay rate assumption would introduce systematic bias into both the fitted parameters and the model predictions. Proverbio et al. (2022) similarly found that incorporating environmental covariates into SEIR-V frameworks improved model fit for longer observation windows, further supporting the rationale for our temperature-aware extension.

The four regions also show where the framework generalizes and where it breaks down. In Greater Boston, where the epidemic followed a clean unimodal profile, the GA-SEIR-V(w) model achieved the highest R² of 0.8922, indicating excellent agreement between simulated and observed case trajectories. This strong performance on unimodal data is consistent with findings by Bivins et al. (2020) and Wu et al. (2022), who demonstrated that mechanistic compartmental models perform most reliably under stable, single-wave epidemic conditions where the underlying transmission dynamics are well-approximated by standard SEIR assumptions. In contrast, the Switzerland dataset, characterized by multiple overlapping Omicron subvariant waves throughout 2022, presented a substantially more complex modeling challenge. The GA-SEIR-V(w) model achieved a calibration-window R² = 0.6887 for Switzerland (against 0.3083 for a same-setting least squares fit; see Section 4 for a fairly re-optimized baseline), yet it placed the epidemic peak approximately 1 month earlier than observed. This temporal offset likely reflects the inherent limitation of single-compartment SEIR frameworks in capturing the rapid immune evasion and reinfection dynamics that characterized Omicron transmission, effects that are not explicitly represented in the standard SEIR structure. Future model extensions incorporating immunity-waning compartments or variant-specific transmission parameters could address this limitation.

The peak-timing analysis (Section 4) shows that this benefit is conditional. When the epidemic peak falls within or near the calibration window—as in Greater Boston—the model recovers the peak day to within about 5 days; when the dominant peak lies far beyond the short calibration window, as in California, Switzerland, and Luxembourg, peak timing cannot be recovered reliably. Where it does apply, early and accurate peak prediction is operationally valuable: it enables health authorities to pre-position resources, activate surge-capacity protocols, and implement targeted non-pharmaceutical interventions before case surges overwhelm clinical infrastructure. Safford et al. (2022) emphasized that the true value of wastewater-based surveillance lies not merely in detecting viral presence but in providing actionable early-warning signals—a capability the framework supports within its calibration horizon.

The genetic algorithm is only one of several population-based global searches that suit this problem, so it is fair to ask whether its particular mechanics matter. Particle Swarm Optimization (PSO), also applied to wastewater epidemiology (Hart and Halden, 2020), shares the population-based global-search property, and GA's crossover operator—which recombines high-fitness parameter sets and has no direct analog in PSO—is often argued to aid exploration. When we tested this directly inside the identical SEIR-V(w) model (Section 4), however, GA, DE, PSO, and MCMC were essentially interchangeable: what mattered was using a global search at all rather than gradient-based least squares, not which global search. We therefore make no claim of an advantage specific to GA; its role here is as a reliable, well-understood member of that family.

Several limitations of this study warrant acknowledgment. First, the analysis is based on three regional datasets, which, while carefully selected to represent contrasting epidemic profiles, may not fully capture the diversity of wastewater surveillance contexts globally, particularly in low- and middle-income settings where wastewater infrastructure, sampling frequency, and data quality differ substantially from those of the study regions. Second, the SEIR-V(w) model assumes a homogeneous, well-mixed population and does not account for spatial heterogeneity, age-structured contact patterns, or the effects of vaccination on susceptibility and shedding dynamics—all of which are known to influence epidemic trajectories, especially during later pandemic phases. Third, while temperature was identified as the primary environmental determinant of viral decay in wastewater and was therefore selected as the sole environmental covariate, other factors including pH, UV radiation, hydraulic retention time, and dilution effects from rainfall also contribute to wastewater viral signal variability (Bertels et al., 2022). Incorporating these additional covariates, where data are available, would likely further improve model accuracy. Finally, the fixed β (viral shedding rate) assumption, while grounded in literature values, may not fully account for inter-individual and variant-specific variation in fecal shedding, which has been reported to span several orders of magnitude across infected individuals (Schmitz et al., 2021).

A further methodological consideration concerns the choice of integer-order dynamics. Fractional-order compartmental models, which replace ordinary derivatives with fractional operators, can encode memory effects and non-exponential growth and decay that have been argued to characterize noisy environmental surveillance data (Rezapour et al., 2020; Trisilowati et al., 2023). Rather than leaving this to speculation, we implemented a Caputo fractional-order version of the SEIR-V(w) system and added the fractional order q to the GA search (Section 4). The GA selected q < 1 in all four regions (0.730.88), confirming a detectable memory effect; however, the fractional order did not systematically improve either calibration fit or out-of-sample prediction over the integer-order model. We therefore retain the integer-order formulation—for which existence, uniqueness, positivity, and disease-free stability are established in Section 3—as an adequate and more parsimonious choice for these datasets, while flagging the detected memory effect as a concrete target for future refinement (for example, region- or compartment-specific fractional orders).

Despite these limitations, the GA-SEIR-V(w) framework represents a meaningful advance in the parameter-estimation stage of wastewater-based epidemic modeling. Its ability to obtain robust, well-identified fits across four regions—spanning stable single-wave, large-scale multi-wave, and long multimodal profiles—with a single unmodified pipeline is rarely demonstrated in existing WBE studies. Reliable long-horizon forecasting, by contrast, remains an open problem that none of the compared methods solved on the longer multi-wave records; closing the gap between strong in-sample calibration and dependable out-of-sample prediction—for instance through immunity-waning or variant-structured extensions, or data assimilation that updates parameters as new observations arrive—is the natural next step. As wastewater surveillance expands globally, coupling evolutionary parameter estimation with such extensions offers a scalable pathway toward more dependable early-warning systems.

6 Conclusion

This study presents an enhanced Genetic Algorithm–optimized SEIR-V model (GA-SEIR-V) that integrates temperature-dependent viral decay to improve the accuracy of wastewater-based epidemic forecasting. By coupling the SEIR framework with an evolutionary optimization approach, the model effectively addresses key limitations of traditional parameter estimation methods, particularly their susceptibility to local optima and inability to capture non-linear system dynamics. Through validation across multiple regional datasets, the GA-SEIR-V model showed its strongest advantage on complex, multi-wave outbreaks, where conventional least squares is prone to poor local optima; on simpler single-wave data a well-tuned least squares fit was competitive. The results indicate that GA-driven parameter optimization, combined with temperature-aware decay, primarily improves robustness and parameter identifiability—and reduces reliance on hand-tuned initial values—rather than uniformly maximizing accuracy across all settings. The proposed framework offers a robust, adaptive, and generalizable approach for epidemic modeling based on wastewater-based epidemiology (WBE). Beyond its predictive advantages, it provides valuable insights into how regional environmental conditions and public health policies shape viral transmission dynamics. In future work, this framework can be extended to include additional environmental factors such as pH, humidity, and wastewater flow dynamics, as well as machine learning-based hybrid architectures that further enhance predictive accuracy. The integration of real-time wastewater surveillance with GA-optimized epidemiological modeling holds significant promise for developing early-warning systems that strengthen public health preparedness and response during future pandemics.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Sahil Jain, Jaypee Institute of Information Technology, India

Reviewed by: Swarna Kanchan, Marshall University, United States

Javad Damirchi, Semnan University, Iran