An exploration of modeling approaches for capturing seasonal transmission in stochastic epidemic models

1.   Introduction

Seasonality is a significant factor influencing the transmission dynamics of infectious diseases, with many pathogens exhibiting distinct seasonal patterns ^[1,2,3,4,5]. Climate change, including shifts in temperature, rainfall, and humidity, further intensifies the complex interactions between hosts, parasites, and vectors ^[6,7]. These environmental factors, along with social behaviors and ecological interactions, affect host susceptibility and transmission rates, driving seasonal fluctuations in disease prevalence. This is particularly evident in the transmission of widespread infectious diseases, such as chickenpox, dengue fever, influenza, Lyme disease, malaria, measles, norovirus, pertussis, rotavirus, and West Nile virus ^{[8,9,10,11,12,13,14,15,16,17,18,19]}. Figure 1(a) highlights the key seasonal factors that drive variations in disease transmission. These factors encompass environmental influences such as temperature fluctuations and rainfall patterns, biological interactions between hosts and vectors, and behavioral aspects like social contact patterns. Together, these drivers shape the dynamic nature of transmission rates across seasons.

Soper ^[20] was arguably the first to address seasonality in epidemic models by incorporating the periodic influx and accumulation of susceptibles, influencing epidemic oscillations based on disease introduction timing and varying susceptibility levels. Moreover, Kermack and McKendrick ^[21], and later Anderson and May ^[22], further developed the theory of periodicity in seasonal diseases, emphasizing the role of susceptible dynamics, transmission rates, and demographic factors. Since then, several researchers looking at deterministic and stochastic epidemic model studies have explored the nonlinear effects of periodic parameters, particularly their influence on seasonality in disease transmission ^{[23,24,25,26,27,28,29,30,31,32,33]}.

Several researchers have reviewed the impact of seasonality and temporal heterogeneity in contact rates on infectious disease dynamics ^[2,3,7,34], providing examples from both human and wildlife systems that highlight the complexities of seasonal environmental drivers ^[1,35,36,37]. Other researchers have focused on offering an overview of statistical methods for assessing and quantifying the seasonality of infectious diseases ^[38,39,40], while some have presented mathematical approaches for formulating and simulating stochastic epidemic models ^{[42,43,44,45]}. Although these studies offer valuable insights into seasonality, heterogeneity, and infectious disease modeling, there remains a need for a detailed examination of modeling approaches tailored to seasonal transmission rates, which I aim to address.

Modeling seasonal diseases is challenging, particularly when multiple seasonal drivers are involved. Accurately capturing these effects is crucial for understanding epidemic patterns and developing effective control strategies. Here, I explore core modeling approaches for incorporating seasonality into disease transmission. These approaches include sinusoidal functions, periodic piecewise linear functions, Fourier series expansions, Gaussian functions, and data-driven methods, all demonstrated through relevant real-world examples. Each approach is assessed in terms of complexity, flexibility, relative advantages, and limitations. Moreover, examples from the literature where these approaches have been applied in epidemic models are cited, and their contributions are briefly discussed.

Additionally, we discuss the application of a seasonal transmission rate in a Susceptible-Infected-Recovered (SIR) model. Given that epidemic outcomes can vary widely, even with the same initial conditions, I extend this study to a time-nonhomogeneous continuous-time Markov chain (CTMC) SIR model. This extension aims to better understand the combined effects of variability and periodicity on the risk of a disease outbreak, offering a more accurate representation of the unpredictability in disease transmission, especially in small populations and early outbreak phases. Figure 1(b) offers an overview of the modeling approaches used to capture seasonal variability in transmission rates, including a compartmental diagram of the SIR model.

Furthermore, important outcome measures, such as the basic and instantaneous reproduction numbers, as well as the probabilities of disease extinction and outbreak derived from the branching process approximation (BPA), are briefly reviewed. Numerical simulations are also presented in which the BPA derived from the CTMC SIR model is validated against Monte Carlo simulations to assess the impact of seasonality on disease dynamics. Figure 1(c) summarizes these critical outcomes, highlighting key metrics such as the basic and instantaneous reproduction numbers, the probability of disease outbreak derived from the BPA, as well as the mean time and standard deviation to disease extinction.

The paper proceeds as follows: The seasonal SIR epidemic model is discussed in Section 2, and various modeling approaches are presented in Section 3. Outcome measures, including the reproduction numbers, probability of outbreak, and related numerical examples, are discussed in Section 4. Finally, in Section 5, I provide a summary and discuss future directions for research.

2.   Seasonal SIR epidemic model

In the deterministic SIR model, $S(t)$, $I(t)$, and $R(t)$ represent the number of susceptible, infectious, and recovered individuals, respectively. In the simplest version of the SIR model, only infection and recovery are considered, with a transmission rate of $\beta$ and a recovery rate of $\gamma$, and no births or deaths are included. The seasonal ordinary differential equation (ODE) SIR model, with a total population size of $N$, is presented below, and the corresponding compartmental diagram is shown in Figure 1(b).

Seasonality is introduced into the transmission rate $\beta(t)$, which is a positive, time-dependent, periodic function with period $\omega > 0$, i.e.,

The stochastic SIR model captures the randomness and variability in disease dynamics, particularly when the number of infectious individuals is small or when fluctuations in transmission and other environmental factors significantly affect epidemic outcomes. The continuous-time Markov chain (CTMC) SIR model extends the deterministic framework using the rates from the ODEs in Eq (2.1) to define the infinitesimal transition probabilities, with the state variables represented as discrete random variables:

The CTMC SIR model is a time-nonhomogeneous Markov process, meaning that the probability of transitioning from one state to another depends on both the current time $t$ and the time elapsed between events. Let $(s, i)$ represent the initial state at time $t$, where $s$ and $i$ are the numbers of susceptible and infected individuals, respectively. Similarly, let $(j, k)$ represent the state at time $t + \Delta t$, with $j$ susceptible and $k$ infected individuals. For a small time interval $\Delta t$, the infinitesimal transition probability from the initial state $(s, i)$ at time $t$ to the state $(j, k)$ at time $t + \Delta t$ is defined as:

where the transitions occur as follows:

Here, $o(\Delta t)$ represents all terms of order $\Delta t$, which become negligible as $\Delta t \to 0$. Since the process is nonhomogeneous, the time between events does not follow a simple exponential distribution, making the standard Gillespie algorithm inapplicable ^[46]. Instead, a Monte Carlo simulation with small time steps $\Delta t$ is employed to numerically simulate sample paths. The events and their corresponding transition rates, which result in changes $\Delta S(t) = S(t+\Delta t) - S(t)$ and $\Delta I(t) = I(t+\Delta t) - I(t)$, are summarized in Table 1.

Both the ODE and CTMC models share the property that $S(t) + I(t) + R(t) = N$, meaning the total population remains constant over time. In both the ODE and CTMC models, when $I(t) = 0$, the system reaches an invariant or absorbing state, indicating that the disease has been eradicated and no further infections can occur, marking the end of the epidemic. Three sample paths for the seasonal CTMC SIR model with $N = 1000$, starting with the initial conditions of a single infected individual, $I(0) = 1$, and $S(0) = 999$, are plotted alongside the solution of the seasonal ODE SIR model in Figure 2.

3.   Modeling approaches

In this section, I explore the different modeling approaches for defining $\beta(t)$ in Eq (2.2), followed by a comparative analysis of the discussed methods, with the transmission rate values demonstrated in the numerical examples aligning with those used in ^[47,48,49].

3.1.   Sinusoidal function

The sinusoidal function is one of the most commonly used methods for incorporating seasonality into mathematical models. It is widely used in epidemiological studies and public health planning for its straightforward representation. Its use dates back to Joseph Fourier's 1822 work, "Théorie analytique de la chaleur ^[50], " which laid the foundation for modeling periodic phenomena. However, its application to infectious disease modeling became more prevalent in the 20th century, particularly in the 1970s for diseases like influenza. Since then, numerous mathematical models have used the sinusoidal function to model seasonal transmission rates ^{[51,52,53,54,55,56,57,58]}. Generally speaking, a sinusoidal transmission rate can be expressed as follows:

where $\beta(t)$ represents the time-varying transmission rate at time $t$, $\bar{\beta}$ is the baseline or average transmission rate around which seasonal fluctuations occur, and $A_{\beta}$ is the amplitude of the sinusoidal function, which determines the extent of seasonal variation. The parameter $\omega$ represents the period of the seasonal cycle (e.g., 12 months for an annual cycle), and $\theta_{\beta}$ is a phase shift parameter that adjusts the timing of the seasonal peak, enabling the model to account for variations in when the peak transmission occurs throughout the year.

One of the major advantages of using a sinusoidal function is its inherent simplicity, which likely contributes to its widespread use in seasonal modeling. Its straightforward nature makes it relatively easy to incorporate into epidemic models, leading to easier interpretation and implementation. Since sinusoidal functions have a very smooth, regular periodic pattern, they are best suited for seasonal diseases that demonstrate clear, regular seasonal patterns. However, the very regularity of sinusoidal functions limits their flexibility in capturing more irregular, sharp, or abrupt seasonal variations.

Sinusoidal functions are frequently used for diseases where seasonal drivers, such as temperature, are closely linked to transmission. This is evident in diseases like influenza and respiratory syncytial virus (RSV), whose prevalence peaks during colder months. In temperate regions, where flu transmission is significantly influenced by temperature, particularly during the winter, a sinusoidal transmission rate, as presented in Eq (3.1), is well-suited for modeling the seasonal dynamics of influenza, as demonstrated in ^[59]. Figure 3 illustrates the total number of Influenza-Like Illness (ILI) cases in New York State for the year 2023, which can be effectively captured using a sinusoidal transmission rate in conjunction with temperature variations throughout the year.

To overcome the limited flexibility of basic sinusoidal functions, the use of Kot-type functions or similar periodic functions, as proposed in ^[7], has gained recent attention in epidemiological modeling. These functions consist of two cosine-modulated terms, making them more complex than simple sinusoidal models, as demonstrated by d'Onofrio et al. in ^[62]. The Kot-type function introduces asymmetric weightings of seasonal regimes, enabling it to capture sharper seasonal variations. The Kot-type function used in ^[62] is given by:

where the new parameter $\epsilon_\beta$ controls the strength of the seasonal forcing, determining how strongly the transmission rate fluctuates due to seasonal effects.

The term inside the fraction in Eq (3.2) modulates the periodic behavior of the transmission rate with two cosine terms, enabling asymmetric seasonal weightings, which makes it more flexible than a simple sinusoidal function. Figure 4 shows a comparison between the transmission rate obtained from the usual sinusoidal function and the transmission rate derived from the periodic Kot-type function.

The Kot-type functions enable sharper seasonal transitions and can capture more complex, non-standard seasonal patterns, particularly in infectious diseases with irregular seasonality or where human behavior significantly impacts transmission cycles (e.g., school terms, vaccination campaigns). These functions provide more control over the shape and timing of seasonal peaks through asymmetrical adjustments, extending their applicability to a wider range of diseases. However, like any complex model, Kot-type functions require accurate seasonal data to calibrate their parameters effectively. Despite these advantages, standard sinusoidal functions remain widely used due to their ease of use and straightforward interpretation.

3.2.   Periodic piecewise linear function

Instead of using a smooth, continuous function like a sinusoidal curve, the transmission rate can be modeled as a piecewise linear function. This approach accounts for factors such as seasonal weather changes and school terms, dividing the year into distinct periods (e.g., winter, spring, summer, and fall). For example, this could be modeled as:

where $m_1$, $m_2$, $m_3$, and $m_4$ represent the slopes of the linear segments associated with the transmission rate in each season. For instance, $m_1$ may capture the effects of factors such as increased indoor activities, colder temperatures, or changes in host behavior, which can lead to a rise in transmission during winter. This is followed by $m_2$, which reflects the transitional dynamics of spring as weather warms and social behaviors shift. Similarly, $m_3$ may account for seasonal declines or stagnations in transmission during summer due to favorable environmental conditions, such as higher temperatures. In contrast, $m_4$ highlights the transition into fall, marked by cooler weather, school reopenings, or other seasonal influences on transmission trends. Furthermore, the intercepts $c_1$, $c_2$, $c_3$, and $c_4$ define the baseline transmission rate for each season, enabling the model to be effectively calibrated to real-world data. Together, the slopes and intercepts enable the piecewise linear function to flexibly capture and reflect seasonally varying dynamics, accommodating abrupt changes in transmission rates rather than assuming uniform or continuous patterns.

The piece-wise structure of this approach offers better flexibility than sinusoidal functions, enabling it to capture non-smooth seasonal variations and sudden abrupt changes in transmission rates, such as those that might occur during specific periods like school calendars, holidays, or brief intervention periods. However, the sharp spikes, discontinuities, and lack of smoothness in piecewise transmission rates can make the modeling, analysis, and interpretation of seasonal disease dynamics more complex and challenging. Additionally, if calibration with data is required, parameter fitting within each selected time interval can be difficult. To address these limitations, I later discuss the use of spline interpolation methods (e.g., cubic splines, B-splines) and other data-driven approaches to achieve smooth representations of the transmission rate in Section 3.5.

Periodic piecewise linear functions can be an effective approach for modeling seasonal diseases where prevalence shows a sharp, notable spike, with transmission rising or falling dramatically within specific periods. González-Parra et al. ^[63] employed the Multistage Adomian Decomposition Method (MADM) to solve nonlinear differential equations by dividing the time domain into subdomains, producing a piecewise solution that accurately captures periodic transmission rates, especially in highly variable cases. On the other hand, Silva et al. ^[64] used piecewise constant parameters to model the impact of public health policies and human behavior on the transmission dynamics of COVID-19. Other notable applications that incorporate a periodic piecewise function can be found in ^[65,66,67].

Carmona and Gandon ^[47] used a square wave, which alternates abruptly between two fixed values, to study vector-borne disease emergence and enhance Zika virus risk maps. This method effectively models sharp seasonal changes in mosquito activity, providing a realistic framework for predicting and managing Zika outbreaks by capturing the 'high' and 'low' transmission states typical of its seasonal dynamics. Figure 5 illustrates the Zika incidence during the 2016 outbreak in Brazil. The transmission rate fluctuates between a baseline level for most of the year and peaks between February and April, indicating a seasonal increase in transmission risk during these months. This fluctuation is represented by the following equation:

A natural extension of periodic piecewise linear functions is the use of piecewise non-linear functions, which offer greater flexibility than linear segments and are better suited to capturing more complex or irregular seasonal patterns. These functions are particularly useful for diseases where transmission is highly dependent on variable factors such as behavior, interventions, or environmental conditions, providing more tailored fits to real-world data. One of the earliest uses was by Capasso and Serio ^[69], who introduced a non-linear incidence rate with a saturated interaction term for modeling cholera transmission. While less common than piecewise linear or sinusoidal functions, piecewise non-linear functions have been used during the COVID-19 pandemic to model transmission changes due to interventions like lockdowns and vaccination rollouts. For instance, Rohith and Devika ^[70] used piecewise non-linear functions to model the nonlinear incidence rate of the COVID-19 pandemic within a Susceptible-Exposed-Infected-Recovered (SEIR) framework. However, while these models offer increased flexibility, it comes at the cost of more complex parameter fitting and interpretation.

3.3.   Fourier series expansion

In epidemiological modeling, a Fourier series expansion is a useful tool for representing seasonal effects with multiple peaks or varying frequencies. It decomposes the seasonal pattern into a sum of sinusoidal terms, each representing a different harmonic or frequency component. The general form is:

where $\beta_0$ is the baseline transmission rate; $a_n$ and $b_n$ are the Fourier coefficients for the cosine and sine terms, determining the amplitude of each harmonic; $\omega$ is the period; and $N$ is the total number of harmonics included in the Fourier series.

The Fourier series offers versatility in modeling periodic phenomena by adjusting the number of harmonics and coefficients, enabling the capture of multiple seasonal effects, such as yearly, quarterly, or biannual variations. Simpler periodic forms, such as standard sinusoidal functions and the Kot-type functions discussed in Section 3.1—the latter specifically designed to model sharp seasonal transitions—are inherently limited in flexibility, particularly in representing multiple seasonal peaks. In contrast, the sum of sines and cosines with different frequencies, as presented in Eq (3.5), is able to capture complex periodic patterns, including multi-periodic transmission dynamics. For instance, with a small number of harmonics, the series captures simple periodic variations like temperature-driven changes, while a higher $N$ can capture more complex effects, such as multiple peaks due to school vacations or behavioral changes. Figure 6 shows how different values of $N$ change the number of harmonics in the Fourier series. However, having multiple parameters increases complexity, meaning more parameters to estimate and a higher risk of overfitting.

Fourier series expansions are particularly suitable for seasonal diseases with more complex or multi-periodic transmission patterns, or where multiple seasonal drivers, such as rainfall and temperature, influence the transmission rate. The use of Fourier series as an alternative to sinusoidal models was explored by Gavenčiak et al. ^[71], enabling an analysis of how results change when applying Fourier series of varying degrees to capture the complex periodic patterns of SARS-CoV-2 seasonality across Europe. Similarly, Deguen et al. ^[72] applied Fourier series to model the periodic contact rate in a seasonal SEIR model. The coefficients of the Fourier series were estimated by fitting the model to weekly chickenpox incidence data, providing a more accurate representation of seasonal fluctuations in the contact rate throughout the year. Furthermore, Mortensen et al. ^[73] demonstrated how machine learning can enhance predictive accuracy by dynamically adjusting parameters to capture evolving transmission patterns, using a Fourier series to estimate the oscillating infection rate and predict COVID-19 dynamics.

Lyme disease transmission involves complex seasonal dynamics, influenced by the activity of ticks and their hosts, which vary over the year. A Fourier series can capture these multiple seasonal peaks and troughs, providing a nuanced model of transmission patterns. In reality, tick activity gradually increases and decreases with changing temperatures and humidity levels, and the Fourier series is well-suited to this type of gradual variation. Figure 7 shows the weekly fluctuations in Lyme disease cases in the USA for 2022, exhibiting a seasonal pattern that justifies using Fourier series to effectively model and capture these periodic trends.

Extensions of Fourier series expansions in epidemic modeling enable the incorporation of greater complexity to capture intricate seasonal or environmental effects. One such extension involves the use of wavelet transforms, as demonstrated by Kumar et al. in ^[75], where both time and frequency information are captured simultaneously. This method is particularly useful in non-stationary cases where disease transmission patterns change over time. Additionally, the coefficients in the Fourier series expansion can be modeled as functions of time, accommodating evolving seasonality due to interventions, environmental shifts, or behavioral changes, such as changes in vaccine coverage or policy interventions over the years. This makes these approaches more suited to capturing real-world complexities in seasonal disease transmission patterns, especially when patterns are multi-scale or affected by a variety of time-varying factors.

3.4.   Gaussian function

A Gaussian function can be effectively used in epidemic models to capture the smooth and symmetrical rise and fall of transmission rates, which often occur during concentrated periods of seasonal transmission in real-world epidemics. Unlike the linear or piecewise functions discussed earlier, which may introduce abrupt changes, the Gaussian function forms a bell-shaped curve (Figure 8) that peaks at a specific time (e.g., the peak season for a disease) and gradually decreases before and after the peak. The general form of the Gaussian function is given by:

where $A_{\beta}$ is the amplitude representing the peak transmission rate, $\mu$ is the mean indicating the peak month of transmission (center of the bell curve), and $\sigma$ is the standard deviation controlling the spread of the curve and determining how quickly transmission rises and falls.

Much like sinusoidal functions, Gaussian functions are relatively easy to implement and interpret, providing a straightforward approach to capturing smooth, natural, well-defined, and concentrated peaks. Additionally, they have low data requirements, with the major parameters of the Gaussian function, the mean and standard deviation, directly mapping to key epidemiological features that drive seasonality. However, unlike more flexible approaches such as Fourier series, the symmetric, predetermined bell shape of the Gaussian function limits its use for diseases with multiple transmission peaks or irregular, asymmetric seasonal patterns.

Gaussian functions are best suited for seasonal diseases where transmission rates rise and fall within a concentrated period, such as during dengue outbreaks with clear seasonal peaks. For instance, Stolerman et al. ^[76] used Gaussian-like functions within an SIR-network model to represent the time-dependent transmission rates of dengue fever in Rio de Janeiro, incorporating environmental factors such as rainfall that influence mosquito populations. Similarly, Ramírez-Soto et al. ^[77] used a Gaussian function in a SIR-SI model to represent seasonal transmission, enabling the model to capture seasonality observed in dengue outbreaks in Lima, Peru. Several other uses of Gaussian functions in seasonal epidemic models can be found in ^[78,79,80].

To improve the flexibility of the Gaussian function in capturing more complex seasonal transmission patterns, a mixture of multiple Gaussian functions can be used. For example, Setianto et al. ^[80] used a Gaussian function to model the time-dependent transmission rate of COVID-19 as a superposition of multiple Gaussian pulses, effectively capturing dynamic changes over time and accurately representing multiple waves of transmission influenced by external factors such as government interventions, behavioral changes, and vaccinations. Additionally, a skewed (or generalized) Gaussian function can be used to model asymmetrical seasonal patterns, or a seasonal modulation factor can be incorporated into the Gaussian function to adjust its amplitude. For example, an extension of the Gaussian function described in Eq (3.7) was used by Hridoy and Mustaquim in ^[81] to develop an exponentially modulated Gaussian function within a SEIR model, which was employed to capture dengue transmission dynamics in Bangladesh, where the disease is significantly influenced by seasonal variations in rainfall. The functional form of $\beta(t)$ used in ^[81] is represented as follows:

where $\beta(t)$ is expressed as a combination of two components: $\beta_{np}$, which represents the constant baseline transmission rate throughout the year, and $\beta_{p}$, which introduces a seasonal Gaussian peak during the rainy season and logistically decays back to the baseline rate. Here, $A_{\beta}$ denotes the amplitude of seasonal variation, $t_{p}$ indicates the peak time at the midpoint of the rainy season, and $\sigma$ controls the width of the seasonal peak. Figure 9 showcases the exponentially modulated Gaussian function used for modeling the seasonal transmission rate, alongside its correspondence with the average monthly rainfall in Bangladesh for 2023.

The Gaussian function presented in Eq (3.6) defines a single peak centered at $\mu$, with $\sigma$ controlling the spread. In a periodic setting (e.g., $t$ ranging from 0 to a fixed period $\omega$), discontinuities can arise at the boundaries where $t = 0$ and $t = \omega$, depending on the relationship between $\mu$, $\sigma$, and the periodic domain. For instance, when $\mu$ is near the boundary (e.g., close to $t = 0$ or $t = \omega$), the Gaussian curve does not naturally extend across the boundaries, potentially resulting in abrupt transitions and loss of continuity. Similarly, in Eq (3.7), while the logistic term adds flexibility in modeling transitions, the Gaussian component can introduce discontinuities near periodic boundaries unless explicitly handled.

To enforce periodicity and eliminate discontinuities, the Gaussian profiles can be incorporated into a periodic framework by explicitly summing periodic copies of the Gaussian function. The periodic form can be expressed as:

where $n$ indexes the periodic copies and $N$ determines the number of Gaussian functions summed. This ensures smooth transitions and effectively addresses discontinuities at periodic boundaries.

3.5.   Data-driven method

If empirical data is available, a data-driven function can be fitted to reflect the seasonality in the transmission rate. This method directly uses observed data for the transmission rate by integrating real-world factors (e.g., temperature, rainfall, humidity, or vector population changes) instead of relying on a predefined functional form, such as a sinusoidal. The general form of $\beta(t)$ in a data-driven function can be expressed as:

where $f$ is a function that relates the data and parameters to the transmission rate over time, and it can represent a wide range of models, such as machine learning algorithms (e.g., neural networks, random forests), regression models (e.g., generalized additive models), or time-series models (e.g., ARIMA, LSTM); $\text{data}(t)$ represents the set of time-dependent variables that act as seasonal drivers influencing transmission; and $\theta$ is a set of parameters or weights estimated from the data (e.g., regression coefficients, model hyperparameters).

Data-driven methods offer significantly greater flexibility compared to predefined functions, such as those discussed in Eqs (3.1), (3.3), (3.5), and (3.6). These methods can learn patterns directly from empirical data, enabling the capture of complex, non-stationary, irregular, and dynamic seasonal variations, as well as the integration of multiple seasonal drivers. The major advantage of data-driven approaches is their ability to provide superior forecasting performance by learning from historical trends, identifying hidden patterns, and incorporating real-time data. Despite these advantages, data-driven methods typically require large amounts of high-quality data to perform well and may struggle to capture meaningful seasonal patterns if the available data is sparse, noisy, or incomplete. Additionally, machine learning-based data-driven methods can be computationally intensive, prone to overfitting, and difficult to interpret due to the incorporation of multifaceted features and complex interactions ^[84]. These limitations are not necessarily present in predefined functional forms, such as sinusoidal or piecewise functions, which can perform adequately even with minimal data because they rely on simple assumptions about seasonality.

Data-driven methods can practically be applied to any seasonal disease where data is available, but they are particularly well-suited for diseases with seasonality influenced by multiple factors, such as vector-borne diseases like malaria, dengue, and Zika, where transmission is affected by temperature, rainfall, humidity, and mosquito population dynamics. These methods are also ideal when high-quality data is available, supported by advanced public health surveillance systems or detailed environmental and behavioral datasets. Furthermore, data-driven approaches can be the best choice for fine-tuned forecasting and early detection, as they can leverage real-time data to enable quick action in mitigating outbreaks. This capability has been demonstrated by numerous data-driven models during the COVID-19 pandemic ^{[85,86,87,88]}.

Several data-driven functions can be used to model seasonality or other time-dependent behaviors in epidemiological models. With the abundance of data in the modern world, exciting work has been done in recent times using data-driven seasonal models, e.g., ^{[87,89,90,91,92,93,94,95]}. An empirical interpolation method, such as linear interpolation, cubic splines, piecewise polynomials, or kernel smoothing, can be used to define $\beta(t)$ based on observed patterns. For example, Prosper et al. ^[30] applied cubic spline interpolation to map temperature data to model parameters (e.g., transmission rates, mosquito survival rates), effectively incorporating temperature-dependent traits of mosquitoes and parasites into a malaria transmission model.

Beyond interpolation, other data-driven approaches to model seasonality include regression-based models like Generalized Linear Models (GLMs) ^[96], Generalized Additive Models (GAMs) ^[97], and polynomial regression ^[98], which offer flexible alternative for capturing non-linear relationships between transmission rates and seasonal drivers. For instance, Aelenezei et al. ^[99] used linear and logistic regression within an SIR model to estimate time-dependent parameters, such as transmission, infection, and recovery rates, for COVID-19 transmission in Kuwait, driven primarily by policy measures instead of seasonal factors like temperature. Additional spline methods such as B-splines, P-splines, and thin plate splines can provide smooth, continuous approximations to represent seasonal trends without overfitting the data.

The accuracy of $\beta(t)$ from Eq (3.8) is heavily dependent on the availability and reliability of the corresponding data. The choice of the data-driven function significantly influences how well the model captures the actual seasonal variation and smooths the data. Consequently, parameter estimation techniques such as Least Squares Estimation (LSE) ^[100], Maximum Likelihood Estimation (MLE) ^[101], or Bayesian Estimation ^[102] can be complementary methods used to fine-tune data-driven functions for modeling seasonality in epidemiological models, enabling the function's parameters to be adjusted to best match the empirical data. Figure 10 shows an example of using linear interpolation in conjunction with least-squares optimization to fit $\beta(t)$ to the observed dengue incidence data from Figure 9 (b).

For further details on an evaluation framework that combines various features, error measures, and ranking schemes for forecast evaluation, see ^[103]. For a side-by-side comparison of results from a stochastic agent-based model and a structured metapopulation stochastic model, refer to ^[104]. Additionally, novel concepts in machine learning, such as identifying treatment effects and explaining model outputs, are discussed in ^[105]. For a discussion of machine learning models, hybrid methods that combine mechanistic models with statistical flexibility, and the challenges of deploying these models in real-world forecasting, see ^[106].

3.6.   Comparison of modeling approaches

A comparative analysis of various modeling approaches, focusing on key factors such as complexity, data requirements, flexibility, relative advantages and limitations, as well as application examples, is presented in Table 2. The evaluation uses terms like "Low, " "Moderate, " "High, " and "Very High" to indicate the relative demands of each method.

4.   Outcome measures

In this section, I discuss essential outcome measures such as the basic, seasonal, and instantaneous reproduction numbers, as well as the probability of disease outbreak using the branching process approximation from the CTMC SIR epidemic model.

4.1.   Reproduction numbers

Reproduction numbers are considered one of the most critical metrics for understanding the spread of infectious diseases in epidemic models. Among the different reproduction numbers, I focus on the basic reproduction number ($\mathcal{R}_0$) and the instantaneous reproduction number ($\mathcal{R}(t_0)$).

For the SIR model, the disease-free equilibrium is given by $S = N$, and $I = R = 0$. The basic reproduction number is simply the ratio of the transmission rate, $\beta$, to the recovery rate, $\gamma$ ^[107], i.e.,

$\mathcal{R}_0$ represents the average number of secondary infections produced by a single infectious individual in a completely susceptible population and can be interpreted as: if $\mathcal{R}_0 > 1$, the disease is expected to spread within the population; if $\mathcal{R}_0 = 1$, the disease will remain stable, neither growing nor declining; and if $\mathcal{R}_0 < 1$, the disease is expected to die out.

For seasonal epidemic models with periodic parameters, the reproduction number differs between a seasonal environment and a constant environment, where all parameters are constant ^[108]. Wang and Zhao's approach, which uses Floquet theory ^[109], can be applied to calculate the seasonal reproduction number. For models with more complex structures, such as the SEIR model, which includes an additional exposed stage, or other models with multiple infectious stages, the next generation matrix (NGM) is the standard method for computing the basic reproduction number ^[110]. However, for simpler models with a single infection stage, such as the SIR model, the basic reproduction number remains unchanged between seasonal and constant environments. The basic reproduction number for the seasonal SIR model with periodic $\beta(t)$, as described in Eq (2.2), is given by:

To better capture the overall ability of a disease to spread in a periodically changing environment, a time-dependent reproduction number, $\mathcal{R}(t_0)$, referred to as the "instantaneous reproduction number, " is defined as

For a comprehensive review on mathematical methods for determining $\mathcal{R}_0$ in ODE disease models and its use in guiding control strategies, see ^[111,112].

4.2.   Probability of disease outbreak

Although the reproduction numbers computed from deterministic ODE SIR models are important metrics for guiding public health responses to infectious disease outbreaks, the use of a branching process approximation (BPA) to estimate the probability of disease outbreak is an essential tool for assessing real-world uncertainties that deterministic models cannot capture.

I briefly review the BPA applied to the CTMC SIR epidemic model ^[28,42,113], with a focus on the infected stage at the disease-free equilibrium (DFE), as depicted within the dotted curve in the compartment diagram in Figure 1. When the number of infected individuals is small, the BPA assumes negligible interactions between them, where each infected individual independently generates secondary cases.

To compute the seasonal probability of disease extinction (i.e., no major outbreak), the BPA is applied at the DFE, where $S = N$. Let

represent the probability of transitioning from $i$ infected individuals at time $t_0$ to $j$ infected individuals at time $t$.

The asymptotic periodic probability of extinction is derived by solving the Backward Kolmogorov Differential Equation (BKDE), applied to the BPA of the CTMC SIR model ^[42]. These backward equations, simpler for computing the probability of extinction than the forward equations, help assess the likelihood of reaching an absorbing state, which marks the end of an epidemic. However, our primary focus is on the stochastic dynamics at the onset of an epidemic, when most of the population remains susceptible. The BKDE for the infected state $I$ is formulated using the transition probabilities $p_{i, j}(t_0, t)$:

where $\tau \in [t_0, t_0 + \Delta t_0]$. By subtracting $p_{i, j}(t_0 + \Delta t_0, t)$, dividing by $\Delta t_0$, and taking the limit as $\Delta t_0 \to 0^+$, we arrive at the BKDE:

Assuming $i = 1$ and $j = 0$, with $p_{1, 0} \equiv p_{1, 0}(t_0, t)$, and applying the simplifying assumptions $p_{2, 0} = (p_{1, 0})^2$ and $p_{0, 0} = 1$, we obtain the BKDE for the probability of extinction:

Simplifying further, we get:

The function inside the square brackets in Eq (4.4) represents an "infinitesimal" probability generating function (pgf), also called the "offspring" pgf ^[42,113]. The asymptotic periodic probability of extinction can be numerically calculated by solving Eq (4.4) for the BPA of the CTMC SIR model. The probability of extinction, $P_{\text{extinction}}(1, t_0)$, is given by:

The probability of extinction depends on both the number of infected individuals introduced, $I(t_0) = i$, and the time of introduction, $t_0$. Due to the independence of infected individuals, the asymptotic probability of disease extinction given $I(t_0) = i$ is:

Finally, the probability of a disease outbreak is:

Bacaër and Ait Dads introduced the use of a multitype branching process approximation for predicting seasonal outbreaks in epidemic models and later confirmed that the basic reproduction number ($\mathcal{R}_0$) from nonautonomous ODE models with periodic coefficients acts as a threshold for disease extinction ^[114]. Similarly, Ball examined the relationship between $\mathcal{R}_0$ and the probability of disease extinction in seasonal epidemic models ^[115]. When $\mathcal{R}_0 \leq 1$, the probability of disease extinction is one, i.e., $\mathcal{P}_{\text{extinction}}(i, t_0) = 1$. However, when $\mathcal{R}_0 > 1$, the probability of extinction becomes periodic and less than one. This finding implies that seasonality plays a critical role in determining the long-term persistence or extinction of an epidemic, highlighting that an $\mathcal{R}_0 > 1$ is not always sufficient for the disease to establish itself when seasonality is pronounced.

Figure 11 displays the seasonal transmission rate $\beta(t_0)$, the instantaneous reproduction number $\mathcal{R}(t_0)$, the probability of disease outbreak $\mathcal{P}_{\text{outbreak}}(t_0)$, as well as the mean and standard deviation (SD) of the time to extinction, where $\beta(t_0)$ is modeled using Eq (3.1).

As shown in Figure 11, the graphs of $\beta(t_0)$, $\mathcal{R}(t_0)$, and $\mathcal{P}_{\text{outbreak}}(t_0)$ exhibit similar shapes. However, while the extrema of $\beta(t_0)$ and $\mathcal{R}(t_0)$ coincide exactly, the extrema of $\mathcal{P}_{\text{outbreak}}(t_0)$ occur earlier in time. This phenomenon, sometimes referred to as the "winter is coming" effect, has been observed across epidemic models, such as SEIR and vector-host models ^{[28,45,47,58]}. The shifts in the extrema of the probabilities of outbreak, $\mathcal{P}_{\text{outbreak}}(t_0)$, are prominently highlighted by Carmona and Gandon ^[47] in scenarios with high transmission rates. They observed that the timing of pathogen introduction relative to fluctuations in transmission rates plays a pivotal role. Notably, the probability of an outbreak can become vanishingly small even during high transmission periods if the pathogen is introduced shortly before a low transmission season. This timing creates a demographic trap, subjecting the pathogen to unfavorable conditions that drastically reduce its likelihood of survival, leading to disease extinction. Additionally, a Monte Carlo simulation is employed to compute the sample paths of the CTMC models. The results from the simulations show good agreement with the branching process approximation, as illustrated by the simulation circles in Figure 11.

5.   Discussion

Seasonality plays a pivotal role in the epidemiology of infectious diseases, significantly influencing disease dynamics and informing the development of effective mitigation strategies. In this paper, I provide an overview of various modeling approaches for incorporating seasonality into transmission rates in epidemic models, ranging from simple sinusoidal functions to more advanced data-driven techniques. Each approach offers distinct advantages depending on the complexity of seasonal patterns, the availability of data, and the specific characteristics of the disease under investigation. While this manuscript serves as a review, it is not intended to be a meta-analysis or an exhaustive survey of all literature on seasonality and infectious disease modeling. Instead, it provides a comprehensive exploration of key approaches for capturing seasonal transmission rates, supported by relevant examples from the literature to illustrate their applications and limitations.

Predefined functions, such as sinusoidal and piecewise linear, are popular due to their simplicity and ease of use. These functions perform well when modeling diseases with clear seasonal patterns. However, their simplicity can limit their ability to capture more complex or irregular seasonal trends. Moreover, Gaussian functions can be better suited for modeling single, concentrated peaks, while Fourier series expansions can be more effective at capturing multiple seasonal peaks, making them particularly valuable for diseases with more complex seasonal dynamics.

Extensions that address the limitations of predefined functions have also been discussed in this paper. For instance, the use of Kot-type functions can overcome the restricted control over the shape and timing of seasonal peaks commonly associated with traditional sinusoidal functions, as demonstrated in Figure 4. Similarly, extending a traditional Gaussian function through the use of an exponentially modulated Gaussian function, as shown in Figure 9(a), can be a valuable enhancement. This extension helps to overcome the limitations of the standard Gaussian function in capturing sharp, narrow peaks, enabling greater flexibility in representing rapid changes in seasonal patterns.

In contrast to predefined functions, data-driven methods offer greater flexibility as they can adapt to the complexities of real-world disease patterns, particularly when multiple factors influence transmission. Specifically, methods such as empirical interpolation methods, generalized additive models (GAMs), machine learning algorithms, and time-series models excel in capturing intricate seasonal variations driven by multiple variables. These approaches are especially relevant for vector-borne diseases, where transmission is often influenced by interrelated factors, such as environmental changes, host behavior, and vector dynamics.

Although models that incorporate data-driven methods are powerful, they are not without challenges. These approaches rely heavily on large, high-quality datasets. In the absence of reliable data, they may struggle to identify meaningful patterns or risk overfitting to noise, which can reduce their predictive accuracy. Therefore, incorporating parameter estimation techniques such as LSE, MLE, and Bayesian estimation can add significant value to data-driven models, helping to improve their robustness and accuracy. Furthermore, machine learning-based models can be computationally intensive and often lack interpretability, making it difficult to understand the underlying mechanisms driving their predictions. Despite these limitations, the flexibility and adaptability of data-driven methods make them valuable tools in epidemiology, particularly for forecasting and early detection of disease outbreaks.

The comparative analysis of the modeling approaches in Table 2 underscores that each method has its own unique strengths and limitations. Ultimately, the selection of an appropriate model depends on several factors, including the characteristics of the disease being studied, the quality and availability of data, and the research objectives. These considerations are crucial in determining whether a simpler, predefined function or a more flexible, data-driven approach will be most effective in capturing the seasonal dynamics of disease transmission.

Additionally, I discuss the seasonal transmission rate in the context of a seasonal SIR model and extend the deterministic model to a time-nonhomogeneous Markov process, specifically a continuous-time Markov chain (CTMC) SIR model. Key outcome measures, such as the basic reproduction number ($\mathcal{R}_0$), the instantaneous reproduction number ($\mathcal{R}(t_0)$), and the probability of disease outbreak $\mathcal{P}_{\text{outbreak}}(t_0)$, are also addressed. We derive the probability of an outbreak using the branching process approximation, based on the CTMC SIR model. These tools can provide valuable insights for both modelers and policymakers, aiding in the planning and implementation of effective interventions. As demonstrated by the "winter is coming" effect in Figure 11, the extrema of the $\mathcal{P}_{\text{outbreak}}(t_0)$ occur earlier than those predicted by the seasonal transmission rate, $\beta(t)$. This suggests that disease mitigation strategies or control interventions may need to be implemented in advance of the peaks predicted by the seasonal driver.

While I focus on fundamental methods, readers should be aware that other techniques are available, which may also be relevant depending on the modeling context. For instance, predefined functional forms, such as Fourier or Gaussian functions, can be used to define seasonal transmission rates in conjunction with parameters estimated from the data, as demonstrated in ^[73] and ^[81], respectively. Besides, adding a noise term, as in ^[116], can account for environmental variability and random changes (e.g., weather anomalies), or using a hybrid approach that combines a sinusoidal function with step functions ^[117] or covariate-based adjustments ^[118], can improve accuracy.

Although I demonstrate the use of a seasonal transmission rate in a stochastic SIR epidemic model, these approaches are broadly applicable to other epidemic settings, such as SEIR or vector-host models with multiple patches and age structures. Furthermore, beyond modeling seasonality in transmission rates, these approaches can be extended to other seasonal parameters, such as seasonal recovery rates or vector control interventions. Therefore, researchers could explore the incorporation of different compartments, structures, and interventions in a stochastic seasonal setting to gain deeper insights into disease dynamics.

On the other hand, incorporating seasonal forcing can lead to more complex disease dynamics, such as multiple stable cycles and chaos ^{[119,120,121,122]}. While these complex behaviors are not widely discussed in the literature due to the unpredictability associated with chaos, which limits the practical applicability of such models, the ability to model both stable cycles and chaotic dynamics has significant implications for public health strategies. A promising future direction would be to explore how chaotic dynamics can be harnessed for practical disease forecasting and control.

Complex interactions between environmental noise, demographic factors, and disease transmission require improved statistical approaches and long-term datasets to fully capture their effects on seasonal disease dynamics. Additionally, human contact patterns remain difficult to model and need further exploration within the context of seasonal epidemics. Researchers should prioritize the development of more robust models that integrate both environmental and human behavioral data. A deeper understanding of these factors will enhance the design of control strategies that align with seasonal disease dynamics, leading to more effective interventions.

Given that many infectious diseases are significantly influenced by seasonal patterns, incorporating these variations into epidemiological models is crucial for improving the accuracy of outbreak predictions and peak intensities, as well as for designing timely interventions. Time-heterogeneity, driven by social and environmental factors, remains a highly productive area of research in mathematical and computational epidemiology. By summarizing key approaches for modeling seasonality and evaluating their strengths and limitations, this study serves as a launching point for researchers and educators, offering handy guidance for the informed parameterization of future mathematical models. Moving forward, it is imperative to prioritize seasonality in infectious disease modeling. As Grassly and Fraser ^[3] aptly noted, "Seasonality has moved from the center ground of early infectious disease epidemiology to the periphery; a distraction often ignored or assumed of minor importance." Re-centering this perspective can lead to improved prediction and prevention of seasonal infectious disease outbreaks.

Use of AI tools declaration

The author declares Artificial Intelligence (AI) tools isn't used in the creation of this article.

Acknowledgments

The author expresses heartfelt gratitude to his Ph.D. advisor, Dr. Angela Peace, for her insightful feedback and guidance throughout this research. The author is also deeply appreciative of the Associate Editor and anonymous reviewers for their thoughtful suggestions. Additionally, the author confirms that no external funding was received to support this study.

Conflict of interest

The author declares that there are no conflicts of interest.