Generalizing The Fractional SIR Model: Memory Effects And Variant-Dependent Outbreak Dynamics

Classical compartmental epidemic models assume that the instantaneous rate of change of each compartment depends only on the current state of the system—a Markovian, memoryless assumption inherited directly from the ordinary differential equation formalism used since Kermack and McKendrick's (1927) original SIR model. This assumption, while mathematically convenient, is epidemiologically restrictive: it implies that a population's current trajectory toward or away from an outbreak depends only on today's susceptible, infected, and recovered counts, and not at all on the path by which the population arrived at that state. Fractional-order generalizations relax this assumption by replacing the integer-order derivative d/dt with a fractional derivative of order α, which weights the system's entire history, via a power-law kernel, in determining its current rate of change—a mathematical structure originally developed within pure and applied mathematics for viscoelastic and anomalous-diffusion problems, and only more recently imported into epidemic modelling.

This has been applied to COVID-19 Delta-variant data using a fractional SIR formulation, but existing applications typically estimate or assume a single fixed value of α for a given dataset, leaving open the question of how outbreak dynamics change systematically as α varies, and, more specifically, whether variation in α alone could serve as a useful proxy for variant-to-variant differences in transmission behaviour that would otherwise require full re-estimation of a model's transmission and recovery rate parameters for each new variant. This paper addresses that question directly, treating α as a free parameter and studying its effect on peak infection level, outbreak timing, and long-run recovered fraction, using the Delta-variant dataset as a numerical baseline. The remainder of the paper is organized as follows: Section 2 reviews the literature on fractional calculus in epidemic modelling; Section 3 presents the model and numerical method; Section 4 presents results across four values of α; Section 5 discusses the interpretation and limitations of treating α as a variant-differentiating parameter; and Section 6 concludes.

LITERATURE REVIEW

Fractional calculus, though mathematically formalized as early as the late seventeenth century in correspondence between Leibniz and L'Hôpital, saw comparatively little applied use until the twentieth century. Podlubny's (1998) foundational text systematized the theory of fractional differential equations and their numerical solution methods, including the Grünwald-Letnikov discretization used in this paper. Podlubny's work, and the broader applied fractional-calculus literature it helped establish, was initially oriented primarily toward physical and engineering applications—viscoelastic materials, anomalous diffusion processes, and control theory—rather than epidemiology, and it was not until considerably later that fractional-order operators were systematically applied to compartmental disease models.

The application of fractional calculus to epidemic modelling accelerated markedly following the onset of the COVID-19 pandemic, as researchers sought modelling frameworks capable of capturing empirically observed features of COVID-19 case data—such as heavy-tailed waiting-time distributions and apparent sub-exponential early growth in some settings—that classical integer-order models struggled to reproduce without ad hoc modification. Khan and Atangana (2020) modelled the dynamics of the novel coronavirus using a fractional derivative approach, incorporating both fractional-order dynamics and a more detailed compartmental structure distinguishing symptomatic and asymptomatic transmission routes, and demonstrated improved fit to early Wuhan case-count data relative to an integer-order comparison model. Shaikh et al. (2020) developed a fractional-order compartmental model specifically for the Indian COVID-19 outbreak, applying potential control strategies within the fractional framework and analyzing the resulting model's stability properties, representing one of the more directly comparable prior applications to the Delta-variant, India-specific dataset used as a baseline in this paper.

In addition to these disease-specific applications, there has been a large amount of methodological literature related to the numerical solution of fractional differential equations for epidemic applications. Several different numerical schemes are available to approximate the Caputo fractional derivative, including the Grünwald-Letnikov discretization, the Adams-Bashforth-Moulton predictor-corrector scheme, and various spectral collocation schemes, which differ in computational cost, numerical stability, and ease of implementation. The choice of the Grünwald-Letnikov approach for this paper was motivated mainly by conceptual simplicity: the memory-weighting structure is explicit in the summation over the simulated history, and the mechanism by which the fractional order dampens and delays outbreak dynamics is easier to examine and understand than in the predictor-corrector or spectral schemes.

A separate strand of literature has examined the biological or epidemiological interpretation of the fractional order α itself, rather than treating it purely as a curve-fitting device. Proposed interpretations have included heterogeneous population mixing (whereby aggregate population-level dynamics deviate from the homogeneous-mixing assumption implicit in integer-order models), reporting or diagnostic delay (whereby the apparent memory effect partly reflects lags in case ascertainment rather than genuine transmission dynamics), and, more speculatively, genuine biological memory in immune or behavioural response. This paper does not attempt to adjudicate between these competing interpretations, but instead treats α as a parsimonious, empirically fittable summary parameter, in keeping with the majority of the applied fractional-epidemic literature reviewed above, which has generally prioritized goodness-of-fit and forecasting performance over mechanistic interpretability of α specifically.

MODEL AND NUMERICAL METHOD

The base (integer-order) fractional SIR model is defined by dS/dt = -βSI/N + ηR, dI/dt = βSI/N - γI, dR/dt = γI - ηR, where β is the transmission rate, γ the recovery rate, and η a waning-immunity rate returning recovered individuals to the susceptible class. The fractional generalization replaces each d/dt with the Caputo fractional derivative of order α, discretized here using the Grünwald-Letnikov (GL) approximation, D^α f(t_k) ≈ h^{-α} Σ_{j=0}^{k} c_j f(t_{k-j}), where the binomial coefficients c_j are computed recursively as c_0 = 1, c_j = c_{j-1}(1 - (α+1)/j). This yields an explicit time-stepping scheme in which the state at each step depends on the full simulated history, weighted by the coefficients c_j—the discrete manifestation of the model's memory. As α approaches 1, the coefficients c_j for j ≥ 1 approach zero, recovering the classical, memoryless integer-order scheme as a limiting special case; as α decreases toward 0, the coefficients decay more slowly with j, meaning that increasingly distant past states exert non-negligible influence on the current rate of change.

Parameters were taken from the Delta-variant dataset (susceptible, infected, and recovered counts normalized to population fractions; transmission rate β = 1, recovery rate γ = 0.07298, giving a classical reproduction number R0 = β/γ ≈ 13.7), with a small waning-immunity rate η = 0.01 introduced to allow long-run recirculation. The model was solved for α ∈ {1.0, 0.9, 0.8, 0.7} over a 60-day horizon with step size h = 0.5 days, using a direct Python implementation of the Grünwald-Letnikov recursion described above rather than a general-purpose fractional-differential-equation solver library, in order to maintain full transparency and reproducibility of the specific numerical scheme applied.

RESULTS

Note. Left: infected-class trajectories for four fractional orders α. Right: peak infected fraction as a function of α.

Figure 1: Infected-Class Trajectories and Peak Infected Fraction Across Fractional Orders α

Fractional order α	Peak infected fraction	Time of peak (days)	Final recovered fraction
1.0 (classical)	0.766	19	0.822
0.9	0.392	28	0.421
0.8	0.144	54	0.116
0.7	0.0006	60 (still rising)	0.0003

Table 1: Outbreak Severity and Timing Across Fractional Orders α

The classical integer-order case (α = 1) reproduces standard SIR outbreak behaviour: a sharp early peak around day 19, with roughly 82% of the population eventually infected. As α is reduced, the outbreak both flattens and delays substantially—at α = 0.9, peak infection roughly halves and shifts nine days later; at α = 0.8, peak infection falls by an order of magnitude and shifts to day 54; and at α = 0.7, the outbreak fails to take off at all within the simulated horizon, with infected fraction remaining below 0.1% throughout. This pattern is monotonic across all four values of α examined, with no non-monotonicity or reversal observed at any point in the tested range, suggesting a robust and consistent dampening relationship between fractional order and outbreak severity over this parameter range.

The magnitude of the effect is worth emphasizing: a reduction in α from 1.0 to 0.7—a 30% reduction in the fractional order—corresponds to more than a thousand-fold reduction in peak infected fraction (from 0.766 to 0.0006) within the simulated horizon. This non-linear, highly sensitive relationship between α and outbreak severity indicates that even modest uncertainty in the estimated value of α, when fitting the model to real variant-specific data, could translate into very large uncertainty in projected outbreak severity—a practical concern for any application of this framework to real-time variant forecasting, discussed further below.

DISCUSSION

The systematic dampening effect of decreasing α is consistent with the interpretation of fractional order as encoding a form of population-level 'memory' or sub-diffusive behaviour in transmission—effectively, a slower, more history-dependent accumulation of new infections that suppresses the rapid exponential takeoff characteristic of the classical Markovian SIR model. This suggests a potential use of α as a compact, single-parameter proxy for variant-to-variant differences in effective transmissibility: rather than re-estimating β and γ separately for each new variant, a shift in α alone can reproduce qualitatively similar dampening or acceleration effects. This is not a claim that α has a direct biological interpretation equivalent to a specific transmission mechanism; rather, as discussed in Section 2, it is offered as a parsimonious modelling device that a specific variant's case data could, in principle, be used to calibrate, in keeping with the predominant fit-oriented rather than mechanism-oriented use of α in the applied fractional-epidemic literature reviewed above.

A limitation is that the waning-immunity parameter η and the initial conditions were held fixed across all four values of α; a fuller analysis would examine the joint sensitivity of outbreak severity to α together with these other parameters, since it is plausible that the strength of the α-driven dampening effect documented here itself depends on the specific values of β, γ, and η used—a joint sensitivity analysis of the kind conducted for a carrying-capacity SEIR model examining parameters K and β jointly would be a natural extension of the present single-parameter sweep. A second limitation, flagged already in Section 2, is that this paper does not fit α directly to variant-specific case-count data for several distinct variants, but instead explores its effect across a plausible range of values using Delta-variant parameters as a fixed baseline; a direct comparative fit across, for instance, Alpha, Delta, and Omicron case-count trajectories would provide a more direct empirical test of whether α varies systematically across variants in the manner this paper's exploratory analysis suggests it plausibly could.

A third consideration concerns numerical stability and interpretability at the lower end of the tested α range: the near-complete suppression of outbreak growth observed at α = 0.7 within the 60-day simulation horizon raises the question of whether this reflects a genuine epidemiological prediction (a variant so strongly dampened by memory effects that it fails to establish sustained transmission) or an artefact of the specific numerical horizon and parameter combination examined. Extending the simulation horizon beyond 60 days, or systematically varying the initial infected seed fraction, would help clarify whether the α = 0.7 case represents a genuine sub-critical regime or merely a substantially delayed, rather than entirely suppressed, outbreak—a distinction with materially different public health implications that the present single-horizon analysis cannot fully resolve.

A further methodological consideration concerns the computational cost of the Grünwald-Letnikov scheme relative to classical integer-order integration. Because the GL recursion at each time step requires summing over the entire simulated history, its computational cost grows quadratically with the number of time steps, in contrast to the constant per-step cost of a classical Runge-Kutta integrator. For the 60-day, 120-step simulations conducted in this paper this cost is negligible, but for longer-horizon simulations—for instance, modelling a full multi-year, multi-variant epidemic trajectory—this quadratic scaling could become a meaningful practical constraint, motivating the use of short-memory or fast approximate GL schemes documented elsewhere in the fractional-calculus numerical methods literature, which truncate the summation to a fixed recent window rather than the full history, trading a small amount of accuracy for substantially improved computational scaling. Such schemes were not required for the present study but would become relevant for any direct extension to multi-year comparative variant fitting.

The comparative magnitude of the dampening effect documented here also invites a cautious methodological note regarding model identifiability: because a reduction in α produces qualitatively similar outbreak-flattening behaviour to a reduction in the classical transmission rate β, or an increase in the classical recovery rate γ, it is plausible that α could be at least partially confounded with these classical parameters when all three are estimated jointly from a single variant's case-count trajectory. Distinguishing a genuine memory effect from a simple re-scaling of the classical rate parameters would likely require either independent biological evidence for memory-driven transmission behaviour, or a joint fitting exercise across multiple variants in which β and γ are constrained to remain relatively stable (reflecting shared underlying disease biology) while α is allowed to vary—an identification strategy left for future empirical work building on the exploratory, single-baseline analysis presented here.

CONCLUSION

This paper generalized a single-variant fractional SIR model into a family of models indexed by fractional order α, and showed numerically that decreasing α produces a substantial, monotonic dampening and delay of outbreak severity, with peak infection falling by more than three orders of magnitude across the range α ∈ [0.7, 1.0]. This offers a candidate mechanism—memory-order variation—for representing variant-to-variant differences in transmission dynamics within a fractional-calculus framework, complementing existing single-variant fractional SIR applications reviewed in Section 2, and suggests that direct empirical fitting of α across multiple named variants, combined with careful attention to its potential confounding with classical transmission and recovery parameters, is a promising, currently underexplored direction for future work in this area.

REFERENCES

1. Kermack, W. O., & McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society A, 115(772), 700–721.
2. Khan, M. A., & Atangana, A. (2020). Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative. Alexandria Engineering Journal, 59(4), 2379–2389.
3. Podlubny, I. (1998). Fractional differential equations. Academic Press.
4. Shaikh, A. S., Shaikh, I. N., & Nisar, K. S. (2020). A mathematical model of COVID-19 using fractional derivative: Outbreak in India with dynamics of transmission and control. Advances in Difference Equations, 2020, Article 373.