Mathematical modeling of infectious diseases and the impact of vaccination strategies

1.   Introduction

The field of epidemiological modeling has a rich history, with seminal concepts introduced by Kermack and McKendrick in 1927 ^[1]. Since then, a variety of models have been developed, ranging from deterministic compartmental models to more complex stochastic and network-based approaches. The susceptible-infected-recovered (SIR) model, a specific case of the susceptible-exposed-infected-recovered-susceptible (SEIRS) model, has been widely used in the study of infectious diseases ^[2]. However, the simplicity of these models has led to the development of more complex ones that attempt to account for more details of real-world disease processes. The threshold theorem, a key result in epidemiology, predicts the critical fraction of susceptibles in the population that must be exceeded for an epidemic to occur ^[3].

Recent studies have emphasized the importance of incorporating latent periods in infectious disease models, with a focus on the incidence rate function. Takeuchi ^[4] and Hattaf ^[5] have both explored the use of general incidence rate functions in delayed epidemic models. The study in ^[6] has extended this work by incorporating a dual time delay and Caputo fractional derivative in a viral system model, while Al-Darabsah ^[7] has applied a time-delayed susceptible-vaccinated-exposed-infected-recovered (SVEIR) model with a non-monotone incidence rate to the dynamics of infectious diseases with an imperfect vaccine. Ghosh has also considered the impact of vaccination and delay in SIR models ^[8]. Notably, recent studies have explored general incidence rate functions satisfying mild conditions ^[9,10,11].

Recent work ^[12] has further highlighted the importance of combining vaccination with other therapeutic measures, such as the treatment of infected individuals, to enhance disease control. It demonstrated that the combination of vaccination and treatment can be more effective in reducing the spread of infectious diseases like dengue, compared to vaccination alone, particularly when vaccine efficacy is not optimal. This finding underscores the potential for integrated approaches in epidemiological models to achieve more robust disease control outcomes.

The field of epidemiology has seen significant growth due to the prevalence of infectious diseases, leading to the development of various vaccination methods. Traditional approaches such as susceptible-vaccinated-infected-recovered (SVIR) and SVEIR have been used, with the pulse vaccination strategy (PVS) emerging as an effective method. This strategy has been applied in various models, including SEIR with time delay ^[13], SEIR with latent period ^[14], and SVEIRS with two time delays and saturated incidence ^[15]. These studies have shown that a high vaccination rate, short pulse of vaccination, and long latent period can lead to disease extinction, making PVS a promising tool for disease control. The primary objective of PVS ^[16] is to minimize the number of susceptible individuals, facilitating their direct transition to the recovered population without becoming infected. This strategy involves repeated vaccination actions within a population during periods of rising infection rates, persisting until disease spread ceases.

In this work, we introduce a novel epidemic model that incorporates both impulsive and continuous vaccination strategies aimed at preventing disease occurrence, termed the disease-free state. We will show a unique periodic solution known as the disease-free periodic solution (DFPS). The analysis of the dynamics associated with this solution provides valuable insights into the effectiveness of a fixed and periodic vaccination strategy. Moreover, we explore scenarios where recovered individuals may experience waning immunity, which highlights the importance of ongoing vaccination efforts. To determine whether the disease will persist or be eradicated, we employ the basic reproduction number, $R_0$, utilizing the spectral radius approach. This method, pioneered by Diekmann et al. ^[17], extends our understanding of disease dynamics and informs long-term control strategies. Through the presentation of a hybrid SIRS model incorporating pulse and continuous vaccination strategies, our study contributes to the broader understanding of infectious disease dynamics and offers guidance for effective disease control and prevention measures.

2.   Model formulation and main results

To study and compare the effect of pulse vaccination and continuous vaccination strategies, we introduce the following hybrid model.

When $t \ne k \tau, k \in \mathbb N$, we have the following delay differential equations governing the disease dynamics:

When $t = k \tau, k \in \mathbb N$, we have the following pulse vaccination applied:

The system can be described using the flow diagram presented in Figure 1, where the arrows show the direction of the single individual to become susceptible, then infected, and in the end recovered. However, in this research work, reinfection is possible because immunity wanes over time, allowing recovered individuals to become susceptible to the disease again. Another thing to note is the vaccination parameters $\theta$ and $\lambda$ that point from susceptible to the recovered compartment by making the individuals immune to the disease. Other parameters are standard and their definitions are explained in the table.

We note several differences between the model above and classical models. While many traditional approaches use autonomous systems, our model employs delay differential equations. Specifically, we chose a delay formulation over the SIRS model's latent class to better represent the incubation period with a uniform distribution, in contrast to the exponential distribution used in the SIRS model. Moreover, in our model, the nonlinear term gets multiplied by $e^{-\mu \omega}$. This term represents the probability that an infected individual remains infectious after a delay $\omega$, which might correspond to the incubation period or the time between initial infection and becoming infectious to others. Similar approaches were considered earlier, see e.g., ^[18,19,20]. This is similar to the exponential term in Nicholson's blowfly equation, which describes the survival probability of blowflies as their population increases ^[21].

We note that similar models, without continuous vaccination, were studied earlier. In ^[19], authors discuss a SEIRS epidemic model with pulse vaccination, time delay, and varying total population size. The study highlighted the effectiveness of pulse vaccination in disease control and eradication, emphasizing that a short interpulse time or a high pulse vaccination intensity can lead to disease elimination. Although the authors did not directly calculate the basic reproduction number, they estimated both upper ($R^*$) and lower ($R_*$) bounds. These estimates helped determine the conditions for disease persistence: values of $R^* > 1$ indicate that the disease is likely to persist in the population, while values of $R_* < 1$ suggest that the disease will likely fade out. A later study ^[20] significantly improves ^[19] by establishing an exact threshold parameter for determining disease eradication and uniform persistence. By deriving the basic reproduction number $R_0$ and analyzing the global attractivity of the disease-free periodic solution, the paper enhances the understanding of disease dynamics and control strategies. Building upon the work of ^[20], authors in ^[22] provided a more comprehensive analysis of global stability in the context of pulse vaccination and nonlinear incidence rates, contributing to the generalization and advancement of mathematical models in epidemiology.

3.   Threshold dynamics

Note that the impulsive delay systems (2.1) and (2.2) are a generalized version of the work in ^[20] with the inclusion of continuous vaccination and a continuation of the study in ^[22] with vaccination policy involved. Thus, one can easily show that the total size of the whole population satisfies the following inequalities:

In other words, $N(t)$ is bounded by a constant $\dfrac{{\Lambda}}{\mu}$. Thus, it is sufficient to analyze the solutions in a bounded and biologically feasible set:

Therefore, the set $\Omega$ is positively invariant with respect to the initial system and the model is well-posed.

Arguing as in ^[22], one can show that there exists a unique globally asymptotically stable disease-free periodic solution $E_0 (t) = (S^*(t), 0, \frac{{\Lambda}}{\mu} - S^*(t))$ of the systems (2.1) and (2.2) where $S^*(t)$ is given by:

where $X = \dfrac{{\Lambda}(\mu + \alpha)}{\mu (\mu + \lambda + \alpha)}.$ Thus, it can be concluded that the recovered compartment oscillates in the same way as the susceptible part with the same period $\tau$.

3.1.   Implicit formula for the basic reproduction number

The basic reproduction number cannot be calculated using the next-generation matrix method because the system is not continuous. Therefore, we use the linearization technique at $E_0(t)$, which approximates $I$ in (2.1) by replacing the non-linear term with $i(t) : = \beta e^{-\mu \omega}S^*(t-\omega)I(t-\omega)$. This leads to the differential equation:

Following similar steps as in ^[22], one can show that $i(t)$ satisfies the following equation:

where $i_0(t) = \beta e^{-\mu \omega}S^*(t-\omega) I(0) e^{-(\mu+d+\gamma)(t-\omega)}$ and

Let us define the linear operator $L: \mathbb{R} \rightarrow \mathbb{R}$:

Note that $K(t, x)$ can be expressed as $K(t, x) = c(t)G(x)$, where $c(t) = \beta e^{-\mu \omega}S^*(t-\omega)$,

It follows from ^[23] that the basic reproduction number can be represented as the spectral radius $r(L)$ of the linear operator $L$. Thus, to find the basic reproduction number $R_0 = r(L)$, the following eigenvalue problem can be solved:

In other words, one can solve

Using the similar steps as in ^[22] one derives the following expression for $R_0$ from the last equation:

It is easy to see that when $\omega = 0$, the above expression is simplified as follows:

We are ready to state the results on the threshold dynamics of the model based on the basic reproduction number $R_0$. The proof will be omitted as it is the same as in Theorems 4.1 and 4.2 of ^[22].

Theorem 1. If $R_0 < 1$, then the disease-free periodic solution $E_0(t)$ of the systems (2.1) and (2.2) are globally attractive. If $R_0 > 1$, then there exist $\xi > 0$ such that $I(t)$ satisfies $\liminf_{t\to \infty}I(t) \ge \xi$ for every positive solution of the systems (2.1) and (2.2). That is, the disease will persist.

Figure 2 illustrates the scenario where global stability is achieved, with parameters set to the values specified in Table 4.

4.   Comparison of various vaccination strategies

In this section, we outline three vaccination strategies: periodic impulsive vaccination (pulse), continuous vaccination (continuous), and a hybrid approach (hybrid). These strategies are evaluated based on their effectiveness in reducing mortality during infectious disease outbreaks.

The pulse vaccination strategy administers vaccines at specific intervals. The key parameters are the pulse period $\tau$, determining the frequency of vaccination events, and the pulse rate $\theta$, representing the proportion of the population vaccinated during each pulse. This approach is effective when resources are limited and vaccines are deployed in bursts.

The continuous vaccination strategy involves vaccinating at a constant rate $\lambda$ over time, gradually increasing population immunity. This method is suited for large-scale programs with a steady vaccine supply.

The hybrid strategy combines pulse and continuous approaches. A portion of the population is vaccinated continuously at rate $\lambda$, with additional pulses at intervals $\tau$ and pulse rate $\theta$. This strategy balances continuous coverage with periodic boosts, ideal for dynamic environments requiring both steady and rapid responses.

To simulate disease dynamics and compare results, we use the parameter values specified in Table 4 for calibration. These values are selected to closely approximate the COVID-19 parameters for Kazakhstan ^[22]. To compare vaccination strategies, we use the percentage of deaths averted as the dependent variable and control for time by setting it to 365 days in our simulations. This is calculated by comparing the total number of deaths in a scenario without vaccination to the total number of deaths with the implemented vaccination strategy. Specifically, the percentage of deaths averted is determined using the formula:

This calculation quantifies the impact of the vaccination strategy on reducing mortality by comparing it to a baseline scenario where no vaccination is applied. The contour plot, depicted in Figure 3, illustrates the percentage of deaths averted as a function of the pulse vaccination period ($\tau$) and pulse vaccination intensity ($\theta$). The color gradient, ranging from cool to warm hues, indicates the proportion of deaths averted, with specific contour lines highlighting the 25, 50 and 75% thresholds.

Figure 4 visually represents how the continuous vaccination rate ($\lambda$) impacts the percentage of deaths averted in a population over a specified period.

To facilitate meaningful comparisons of vaccination strategies, we maintain a consistent framework by fixing the vaccination capacity to vary from 10% to 100% of the total population. While the correlation between vaccination and reduced deaths may seem intuitive, it is not immediately clear which strategy is more effective in averting deaths when vaccination capacity is taken as a control variable. Since the amount of vaccination and its cost are highly correlated, fixing the vaccination capacity effectively means treating cost as a control variable.

In our methodology, we define the percentage of deaths averted as a function of the vaccination parameters and use Python's minimize function to find the optimal parameter values. The minimize function employs the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm as the default optimization method. BFGS employs gradient information to find local minima efficiently without computing the Hessian matrix ^[24]. This method helps in identifying the optimal configurations for the vaccination strategies by minimizing the defined function.

The tables offer valuable insights into both the optimal parameters for vaccination deployment (Table 2) and the resulting impact on mortality rates (Table 3). Across diverse vaccination rates and strategies, the data unveils nuanced relationships between parameter configurations and mortality outcomes. Particularly striking is the positive correlation between higher vaccination rates and increased percentages of deaths averted, highlighting the pivotal role of vaccination coverage in mortality mitigation. Moreover, variations in vaccination strategies, including pulse, continuous, and hybrid approaches, yield comparable effectiveness levels, with subtle performance discrepancies observed at varying vaccination rates.

Comparing the percentage of deaths averted under different vaccination strategies—pulse, continuous, and hybrid—across various vaccination rates unveils notable effectiveness differences. For example, at a 50% vaccination capacity, the pulse strategy averts approximately 45.61% of deaths, while the continuous and hybrid strategies avert around 45.18 and 45.69%, respectively, indicating relatively similar effectiveness levels among these strategies. Furthermore, delving into the optimal timing and frequency of periodic pulse vaccinations, exemplified by parameters such as ($\tau$, $\theta$), highlights insights into maximizing the death avertion rate. For instance, for a 50% vaccination rate, an optimal configuration might be ($\tau = 23$, $\theta = 0.0460$), shedding light on how parameter variations impact vaccination campaign scheduling and mortality outcomes. Additionally, examining how these strategies influence vaccine coverage rates over time, particularly at critical vaccination rates like 50%, provides valuable insights into the dynamic interaction between vaccination strategies and their efficacy in mortality mitigation.

The regression analysis in Figure 5 shows strong linear relationships between the parameters and the percentage of deaths averted. For pulse vaccination with a fixed period $\tau = 23$, an $R^2$ of 0.997 indicates a nearly perfect linear relationship between vaccination capacity and deaths averted. A 10% increase in coverage results in a 9.2% rise in deaths averted (p-value = $1.37\times10^{-11}$). The continuous vaccination rate $\lambda$ also has a robust relationship, with a slope of 22374.39 and an $R^2$ of 0.9876, indicating that small increases in $\lambda$ can significantly improve outcomes (p-value = $6.47\times10^{-9}$). Similarly, the intensity of pulse vaccination $\theta$ has a near-perfect linear relationship with a slope of 998.5 and an $R^2$ of 0.998, highlighting its importance in reducing mortality (p-value = $1.71\times10^{-10}$). These results emphasize the effectiveness of both continuous and pulse vaccination strategies in saving lives, with optimal parameter adjustments offering substantial public health benefits.

5.   Sensitivity analysis

The objective of this analysis is to determine the sensitivity of different parameters in the SIRS model to the final outcomes of infected, vaccinated, and deceased. Moreover, we also study the sensitivity on the basic reproduction number. We use the partial rank correlation coefficient (PRCC) to quantify the relationships between input parameters and the model outcomes, providing insights into the significance and impact of each parameter on the disease dynamics.

5.1.   Sensitivity analysis of parameters across different compartments

A total of 4000 parameter sets were generated using Latin hypercube sampling from the parameter ranges provided in Table 4, and the model was simulated for 365 days for each set. PRCC values were computed between each parameter and the final outcomes of the compartments to measure the sensitivity and the strength of association. PRCC values close to +1 or $-1$ indicate a strong positive or negative correlation, respectively, while values near 0 indicate weak or no correlation. The PRCC values for each parameter with respect to the final outcomes in the compartments $I$, $V$, and $D$ were computed and plotted. Statistical significance was assessed using $p$-values, the significance of the PRCC values was indicated by asterisks, with *** for $p < 0.001$, ** for $p < 0.01$, and * for $p < 0.05$.

Bar plots illustrate the PRCC values for each parameter across all compartments. Positive and negative correlations were color-coded, and significant correlations were marked. A color bar was included to represent the magnitude of PRCC values. Scatterplots were generated to visualize the relationship between each parameter and the final values of the compartments. These plots aid in understanding the nature of correlations and the spread of data points.

PRCC values for each parameter affecting compartment $I$ are depicted in Figure 6. Additionally, a correlation scatterplot is provided in Figure 8 to visualize the sensitivity of the parameters on $I(365)$. From the flow diagram (Figure 1), we know that the parameters $\beta$, $\mu$, $d$, and $\gamma$ have a direct effect on $I$.

The analysis shows that while the continuous vaccination rate ($\lambda$) negatively impacts $I$, the rate of losing immunity ($\alpha$) has a positive effect. Moreover, increasing the pulse vaccination period also leads to an increase in $I$.

The dynamics of PRCC values for $I(t)$ are shown in Figure 8. An interesting observation is that, during the initial stages until approximately day 225, the incubation period $\omega$ has a negative effect on $I(t)$. However, after day 225, this relationship becomes positive.

While we do not consider death as a separate compartment in the model, it can still be tracked using the equation:

The PRCC values for each parameter affecting the death compartment $D$ are depicted in Figure 6. Additionally, a correlation scatterplot is provided in Figure 11 to visualize the sensitivity of the parameters on $D(365)$.

From Figure 6, it is evident that the transmission rate ($\beta$) has a very strong positive correlation with $D$. In contrast, the vaccination rate ($\lambda$), recovery rate ($\gamma$), and incubation period ($\omega$) all have moderate negative effects on the death toll.

The dynamics of PRCC values for $D(t)$ are shown in Figure 11. We note that from day 300 onwards, the PRCC values start converging.

To study the sensitivity of the parameters to the vaccinated population, we introduce the vaccination compartment $V$. The rate of change of $V$ is given by:

depending on whether time $t$ is divisible by $\tau$ or not, respectively. From Figure 12, we see that all the parameters have a significant impact on $V$. The scatterplot in Figure 13 provides a visual representation of these correlations, illustrating the sensitivity of the parameters on the vaccination compartment.

As before, we provide the PRCC dynamics of $V(t)$ on the parameters in Figure 14.

5.2.   Numerical approximation of $R_0$ and sensitivity analysis

The basic reproduction number depends on two pulse vaccination coefficients such as $\tau$ and $\theta$, therefore, the numerical approach is considered to check the dependence on these two parameters, which is the discretization method.

The discretization approach is a procedure of transforming the continuous functions or expressions into discrete intervals to calculte them using computer software. So the finite set of values are used to get an answer close to the original system's solution.

From Eq (3.2):

Making a change of variable $t-x = s$:

Dividing the integral into two parts and discretizing the second part:

The integral and summation sign can be swapped to notice that:

The result is as follows:

If $i < k$:

If $i \geq k$:

where

The plots in Figure 15 depict the convergence behavior of the basic reproduction number ($R_0$) approximation method across different numbers of subintervals ($N$). Notably, it is observed that convergence stabilizes within $N = 100$ steps, indicating that this value provides a reliable approximation. Therefore, for subsequent estimations of $R_0$, we adopt $N = 100$ as it ensures accurate results while maintaining computational efficiency.

We consider parameter ranges as outlined in Table 4 and employ a random selection of 2000 parameters using Latin hypercube sampling to compute PRCC for $R_0$ against each parameter.

The sensitivity analysis of the basic reproduction number ($R_0$) provides valuable insights into the dynamics of the modeled system. This analysis not only highlights the system's responsiveness to changes in parameters but also underscores the inherent uncertainties in the model. To address sensitivity, we adopted the approach from ^[25], which offers robust methods for managing and analyzing uncertainty in complex systems. We refer to ^[26] for other approaches to address uncertainty. Notably, the transmission rate ($\beta$) exhibited a strong positive correlation (PRCC Value: 0.730) with $R_0$. Conversely, the continuous vaccination rate ($\lambda$), displayed a moderate negative correlation (PRCC Value: $-0.297$), suggesting its influence in reducing the system's basic reproduction number. Parameter ($\alpha$), representing the rate of losing immunity, demonstrated a moderate positive correlation (PRCC Value: 0.280), implying its contribution to the system's dynamics. Parameters $d$ and $\gamma$, representing the disease-induced death rate and recovery rate, respectively, exhibited moderate negative correlations (PRCC Values: $-0.222$ and $-0.241$), highlighting their impact on the $R_0$. Additionally, parameter $\omega$, representing the latent period of the disease, displayed a weaker negative correlation (PRCC Value: $-0.103$), indicating its limited influence on the system. The period between two pulse vaccinations ($\tau$) showed weaker correlations (PRCC Values: $-0.103$). The pulse vaccination intensity $\theta$ showed insignificant impact (P-value: 0.266). Comparatively, although both continuous vaccination and pulse vaccination strategies contribute to disease control efforts, the sensitivity analysis suggests that increasing the continuous vaccination rate has a more direct and significant impact on reducing the basic reproduction number. The absence of a similar effect of pulse vaccination stems from its dependence on two interrelated parameters, $\tau$ and $\theta$, which should not be analyzed in isolation.

6.   Conclusions

This study evaluated various vaccination strategies, including pulse, continuous, and hybrid methods, to determine their effectiveness in controlling infectious diseases and reducing mortality. Key findings revealed that continuous vaccination has a strong negative correlation with the basic reproduction number ($R_0$), making it a robust method for reducing transmission rates and mitigating mortality. The transmission rate ($\beta$) showed a strong positive correlation with $R_0$, while the continuous vaccination rate ($\lambda$) displayed a moderate negative correlation. At a 50% vaccination capacity, the pulse strategy averted approximately 45.61% of deaths, while continuous and hybrid strategies averted around 45.18 and 45.69%, respectively. Sensitivity analysis indicated that continuous vaccination has a more direct impact on reducing $R_0$ compared to pulse vaccination. Although pulse vaccination is effective, it requires precise scheduling for optimal results. Ultimately, the study highlights the importance of high vaccination coverage and optimal strategy selection in public health policies to effectively combat infectious diseases.

In this work, we fixed the vaccination capacity as a control variable, which indirectly accounts for the associated cost. Other studies have explored how different vaccination strategies affect economic costs (see e.g., ^[28]), and similar analyses could be integrated into our model to further understand the economic implications of vaccination strategies.

One limitation of our study is the assumption that all susceptible individuals will eventually be vaccinated. In reality, as seen with COVID-19, a significant portion of the population may choose to remain unvaccinated, as discussed by Ledder ^[27]. Future work could incorporate more realistic vaccination strategies into the model to address this issue. Another limitation of this study is that impulse vaccination strategies are rarely used in real-world scenarios, and continuous vaccination rates ($\lambda$) typically vary over time, which complicates the validation of our results with real-life data. Another possible avenue for future work is to investigate the dependence of vaccination strategies on time. In this study, we fixed the time period to 365 days and treated it as a control variable rather than a dependent one. Future research could explore whether the optimal vaccination strategy varies with different time periods.

Use of AI tools declaration

The authors declare that, in minor parts AI tools were used to proofread the text.

Acknowledgments

The second author acknowledges research funding by the Science Committee of the Ministry of Science and Higher Education of the Republic of Kazakhstan (Grant No. AP19676669).

Conflict of interest

The authors declare that there is no conflict of interest.