A comparison and calibration of integer and fractional-order models of COVID-19 with stratified public response

1.   Introduction

On January 25th, 2020, the first case of COVID-19 was confirmed in Toronto, Ontario ^[1]. As of November 2021, SARS-CoV-2, the virus that causes the disease COVID-19, has infected over 1.7 million Canadians ^[2]. Throughout the course of this pandemic, the general public has been burdened with the responsibility of reducing transmission through various public health efforts. Public health officials have successfully lobbied the general public to participate in transmission reduction by employing tools such as face masks, social distance, ventilation upgrades, and vaccinations. It is largely by employing these public health guidelines that further epidemic waves are avoided in local contexts.

Mathematical models of the spread and evolution of COVID-19 have been utilized to help answer questions of public health policy, vaccine deployment, and allocation of treatment resources. Indeed, the successes of compartmental models in epidemiological contexts are hard to overstate. Such mathematical models are incredibly powerful at producing short-term forecasts of the evolution of epidemic waves. More importantly, however, these models create a test-bed by which various public health policy questions can be answered. As such, various situations can be simulated in order to gain insight into the possible effects that changes in the behavior of the populace can create. In either case, whether the model is utilized for prediction or for experimentation, any insights gained by the model are limited by the capacity for the underlying model to describe the spread of the disease within the local population.

Fractional calculus has been employed in various contexts where classical (integer-order) derivatives have met success. Namely, since fractional derivatives are calculated by utilizing global, integral operators, the use of fractional derivatives in differential equations results in a dynamical system capable of capturing the effects of state memory and other non-local information that integer-order methods are incapable of utilizing ^[3,4]. As a result, fractional-order methods have met great success in biological systems in general such as in the modelling of cancers ^[5], modelling the spread of COVID-19 ^[6,7,8,9], and describing the dynamics of HIV ^[10,11]. More specifically, fractional-order methods have found great success in epidemiological compartmental systems in particular ^{[12,13,14,15]}. Fractional derivative operators are not the only derivative operators capable of producing memory effects, for instance, delay differential equations can introduce memory effects as well ^[16]. In principle, a delayed differential operator could be employed as well. However, given the success of fractional operators in the biological literature (and specifically the theoretical epidemiology literature), we focus on fractional derivative operators in this report. From a particular perspective, time series prediction can be accomplished via the use of iterative techniques, machine learning, and deep learning network approaches ^{[17,18,19,20]}. While these are popular methods for time series forecasting, we are focused on increasing the interpolative and predictive strength of phenomenological modeling of epidemics via modifications of classical compartmental models.

From a practical standpoint, there are multiple choices of how to implement a fractional derivative. Of these choices, it is likely the Riemann–Liouville derivative and the Caputo derivative are the most popular. There are certain particulars to consider when choosing between implementing a Riemann-Liouville derivative or a Caputo derivative. For instance, the Riemann–Liouville has certain disadvantages for modeling physical or biological systems. The Riemann–Liouville derivative of a constant is not zero. In addition, if the Riemann-Liouville derivative is used in a dynamical system, then the initial conditions must contain the limit values of Riemann-Liouville derivatives at time $t = 0$. This is concerning when modeling physical or biological systems as the initial conditions become decoupled from physical or biological interpretation. These disadvantages reduce the field of application of the Riemann–Liouville fractional derivative. On the other hand, some of the great advantages of the Caputo fractional derivative are that the Caputo fractional derivative of a constant vanishes, and those dynamical systems under the Caputo fractional operator allow for traditional initial and boundary conditions to be included in the formulation of the problem ^[21].

It is often the case in mathematical biology contexts that the model of the phenomenon being considered contains certain tune-able, unknown parameters that must be calibrated to data in order for the model to produce meaningful descriptions of the physical behavior. The process of determining these unknown parameters is non-trivial, and it is often the case that some of these parameters can not be determined simultaneously, even in the presence of sufficiently large data sets ^{[22,23,24,25]}. These questions of the identifiability of unknown parameters, both structural and practical, are important theoretical considerations for any calibrated model. Here, by performing an identifiability analysis, our previous model of infectious disease in the Canadian province of Ontario can be investigated. Moreover, we compare this model with a modified fractional-order model and demonstrate to accurately describe the dynamics of multiple waves of epidemic infection.

Our manuscript is laid out as follows. In Section 2, we first introduce the integer-order model of SARS-CoV-2 spread in the Canadian province of Ontario that we consider. Then, we extend this model to the context of the fractional Caputo derivative in Section 2.2. Next, we discuss the calibration methodology used to fit the model to the epidemiological data in Section 2.3. We begin the Results and Discussion by presenting and discussing the results of the model calibration process in Section 3.2. We then utilize the results of the fitting procedure in order to perform sensitivity and identifiability analysis post-hoc validating the results of the calibration process. Then, we compare the results of the integer and fractional models in Section 3.3, comparing not only the capacity by which either model can describe the data but also comparing the predicted public response function $\theta(t)$ and its effects on the time-dependent reproduction number ${\mathcal{R}_t}$. Finally, we summarize our conclusions in Section 4.

2.   Materials and methods

2.1.   Integer-Order Mathematical Model

We consider the Distancing-SEIRD model from Eastman et al. ^[26]. This compartmental epidemic model stratifies the population into five distinct subgroups. Individuals in group $S$ represent those susceptible individuals who have not yet been exposed to the virus, individuals in group $E$ represent those who have been exposed to enough viral pathogen as to eventually become infectious themselves after some incubation period, and individuals in group $I$ represent those who are infected with the disease and are infectious to others, while individuals in groups $R$ and $D$ represent those who have either recovered from the disease or have perished mode of COVID-19.

Moreover, it is assumed that the population can further be stratified into two groups: those that actively follow public health guidelines to reduce transmission and those that do not. Rather than re-stratify the existing five compartments into ten compartments total, the authors decide to model the proportion of the population that is reducing transmission as a time-dependent proportion $\theta(t)$. When $\theta(t)$ is near $1$, most of the population is attempting to reduce transmission. In contrast, when $\theta(t)$ is near $0$, most of the population eschews public health guidelines for transmission reduction.

along with initial data

where

is a time-dependent transmission rate determined by $\theta(t)$. Furthermore, in order for the model to be biologically relevant, we consider Equations (2.1) only under non-negative initial data.

Note that while, in principle, following transmission reduction guidelines may not mechanistically cause a change in disease recovery or disease mortality, there still may be a difference in these effects between groups due to the population effects of individuals who follow these behaviors. For instance, an individual who spends less time in public or interacting with friends may recover quicker solely because they spend more time resting as a result or to demographic effects of those more likely to follow transmission reduction guidelines. In any case, as in Eastman et al. ^[26], we make the simplifying assumption that disease-specific parameters are unaffected by the transmission reduction behavior of the host: in effect, we assume that $\alpha_D = \alpha_M$, $\delta_D = \delta_M$, $\gamma_D = \gamma_M$, and $\beta_{DM} = \beta_{MD}$ in Equations (2.1) for simplicity.

2.2.   Fractional-Order Mathematical Model

To model the fractional-order derivative of the epidemic model, we shall introduce a modified fractional differential operator as proposed by Caputo in his work on the theory of viscoelasticity ^[27]. While other fractional derivative operators could be considered ^[28,29,30], there is an abundance of literature in the physical sciences that use the Caputo operator. See for instance ^{[31,32,33,34]}. Indeed the Caputo operator has found great success in the modeling of infectious diseases ^{[35,36,37,38]}.

Definition 1. The Caputo-fractional derivative of order $a\in (k-1, k)$ of $f$ is defined as

for $t > 0$ (where $\Gamma(\cdot)$ refers to Gamma function). In the case $a\in (0, 1)$ we have

2.2.1.   Discretization Options

From a computational point of view, we do not have any closed form for calculating the Caputo derivative Equation (2.4) a-priori. To overcome these issues, we consider a popular numerical discretization of the Caputo derivative operator. The most well-known approximation of Equation (2.4) is given by the so-called L1-formula as in the notation of Lin and Xu ^[39], the L1-discretization of Equation (2.5) over a uniform time mesh $\{0 = t_0, t_1, \ldots, t_{n}\}$ with step size $\Delta t$

where $b_{i} = (i+1)^{1-a} - i^{1-a}$. We define $\zeta_{t}(f)$ as the truncation error, then $\zeta_t(f)$ satisfies the order-estimate

Hence, given $a\in(0, 1)$, the L1 discretization results in at-least linear order accuracy. In order to affect more accurate results, we can consider a higher-order approximation, such as the L1-2 formula presented due to Gao et al. ^[40]

where

Gao et al. ^[40] demonstrated that the truncation error satisfies the following order estimate \(\lVert\zeta_{t}(f)\lVert = O(\Delta t ^{3-a})\). Hence, the L1-2 discretization represents an increase of a single order of magnitude over the L1 discretization alone (and hence for $a\in(0, 1)$, the L1-2 discretization is at least quadratic).

For algebraic simplicity, we first transform the L1-2 discretization from (2.6) to the following form:

in which $c_{0}^{(a)} = a_{0}^{(a)} = 1$ for $i = 1,$ and for $i\geq 2$

For the numerical simulations of the fractional model in this work we consider the approximation of the Caputo fractional operator given by ${}^{C}_0D_{t}^{{a}}f(t_{n})-\zeta_t(f)$ as in Equation (2.7).

2.2.2.   Fractional-Distancing-SEIRD Model

Since the Caputo fractional derivative preserves the integer-order initial data, we can apply the Caputo fractional derivative in place of each derivative operator in Equations (2.1). Moreover, these 5 different applications of the Caputo fractional derivative can all be of unique order $a_S$, $a_E$, $a_I$, $a_R$, and $a_D$.

along with initial data in Equation (2.2) and $\beta(\theta(t))$ is defined as in Equation (2.3). The presence of the $\Delta t^a_i-1/\Gamma(a_i+1)$ factor in front of the derivative operator in Model 2.8 is a result of the fractional-order Taylor expansion, where $\Delta t$ is a characteristic time. In numerical simulations, we took $\Delta t$ to be the time step size of the model. For a full derivation, please see (for instance) ^[17]. We first note that initially in the epidemic model, the majority of the population mass is within the susceptible compartment. Moreover, according to the data from Berry et al. ^[2], the sum of the susceptible and exposed has only decreased in size by around 4% between January 25th, 2020 and, November 1st, 2021. Since reports that the latent period of the disease is on the order of days, we assume that the relative size of the susceptible compartment remains quite large over the course of the time frames being considered ^[41]. Since the relative size of $S$ is the largest amongst the compartments, we suspect that replacing the derivative operator in the $S$ equation with the Caputo fractional derivative operator will have the largest effect on the qualitative dynamics of the system. As a result, we consider $0 < a_S < 1$ with $a_E = a_I = a_R = a_D = 1$. Having only one equation in the system be described as a fractional-order operator and the rest by an integer-order operator is computationally efficient in our approach. In ^[42,43] the authors consider a similar mixed fractional and integer-order scheme. Hence the modified Fractional-Distancing-SEIRD model that we consider is defined as

along with initial data in Equation (2.2) and $\beta(\theta(t))$ is defined as in Equation (2.3).

2.3.   Fitting procedure

To calibrate the models in Equations (2.9) and Equations (2.1), we utilize the data collected by the COVID-19 Canada Open Data Working Group for the Canadian province of Ontario ^[2]. This data-set contains data of the form $\{I_t, R_t, D_t\}$ where $I_t$ represents the cumulative known infections at time $t$, $R_t$ represents the cumulative known recoveries at time $t$, and $D_t$ represents the cumulative known fatalities at time $t$. As in Eastman et al. ^[26] we consider $\theta(t)$ as the simple linear spline. Hence $\theta(t)$ is the linear interpolant of the paired data $\{(0, \theta_0), (\tau_0, \theta_0), (\tau_1, \theta_1), \ldots, (\tau_{20}, \theta_{20})\}$, where $\theta_0$, $\tau_0$, $\ldots$, $\tau_{20}$ are as in Table 2.

To quantify the quality of the fit, we consider a simple weighted sum of squared errors of infection, recovered, and fatality data as in Equation (2.10)

We then minimized this function in both the integer and fractional-order case by solving Equations (2.1) and Equations (2.9) for a given parameter set and calculating the value of Equation (2.10). This cost function value is then minimized via Matlab's implementation of the genetic algorithm ^[44]. We initialized the population used in the genetic algorithm with a stochastic scheme. For each of the parameters $\gamma$, $\alpha$, $\beta_{DD}$, $\beta_{DM}$, and $\beta_{MM}$ in each of the initial population vectors, we initialized the value to be either a randomly selected value or the corresponding value from Eastman et al. ^[26]. Given that those values were selected by the authors for fitting the first wave of SARS-CoV-2 in Ontario, we would expect similar values to apply in this case. All other parameters were initially randomly selected. In both the fractional and integer-order case, we ran the genetic algorithm on a population of 1500 vectors for a maximum of 30,000 generations stopping early if the average relative fitness value has not changed by more than $10^{-6}$ for at least 100 generations. Finally, the fits are polished via Nelder-Mead's simplex algorithm to ensure convergence.

For the integer-order case, Equations (2.1) were solved via Runge-Kutte of order 4/5. For the fractional-order case, Equations (2.9) were solved via the L1-2 operator in Equation (2.7) for the $S$ compartment and the Runge-Kutte 4/5 operator in the remaining compartments. Both were solved over the time mesh $(t_0, t_1, \ldots, t_n) = (0, 1, \ldots, 376)$.

2.3.1.   Sensitivity and Identifiability Analysis Procedures

The model calibration method, as outlined in Section 2.3, contains a large number of unknown parameters. In this section, we outline the process by which we verify that the fitting procedure can leverage sufficient information in order for the method to uniquely identify optimizing parameters. As such, we borrow the definition of practical identifiability from Miao et al. ^[45]

Definition 2 (^[45]). A dynamical system parameterized by $\xi$ is practically identifiable if $\xi$ can be uniquely determined from the measurable system output $\{y(t_i)\}$.

There are various methods by which one can determine the set of practically identifiable parameters for a given model. Moreover, these different methods are not guaranteed to produce the same set of identifiable parameters. For instance, in ^[22], López et al. consider such practical identifiability analysis methods such as the variance method, SVD method, QR method, and Monte-Carlo method. Moreover, the authors note that while the various methods rarely produce the same number of identifiable parameters, combinations of methods are often able to produce an identifiable subset of appropriate dimensions.

For a given output signal $y$, the sensitivity matrix is defined as

from which the Fisher information matrix is defined $F = \chi^T\, \chi$.

We follow the SVD and QR methods from ^[22] (though these methods are not unique to these authors, see for instance, ^{[23,24,45,46,47]}). Both the QR method and the SVD method depend upon the sensitivity matrix from Equation (2.11). In particular, both methods fail if the Fisher information matrix is singular. Given the output signal

There are certain choices of weights $W_I$, $W_R$, and $W_D$ that result in duplicate columns in $\chi$ (and hence, a singular Fisher information matrix).

Lemma 3. Let $(S, E, I, R, D)$ be a solution to Equations (2.1), then

Proof. First, let $S_\alpha = \frac{\partial S}{\partial \alpha}$, $E_\alpha = \frac{\partial E}{\partial \alpha}$, and $I_\alpha = \frac{\partial I}{\partial \alpha}$ then

Similarly,

If we define $S_\delta$, $E_\delta$, and $I_\delta$ analogously, then by a similar argument we find

Now the system

has a unique solution and so $S_\alpha = S_\delta$, $I_\alpha = I_\delta$, and $E_\alpha = E_\delta$

Lemma 4. Let $(S, E, I, R, D)$ be a solution to Equations (2.1). For $y = \frac{w_I\, x_v\, I+w_R\, R+w_D\, D}{w_I+w_R+w_D}$ then $\frac{\partial y}{\partial\alpha} = \frac{\partial y}{\partial\delta}$, if and only if, $w_R = w_D$.

Proof. Let $\hat{w}_I = \frac{w_I}{w_I+w_R+w_D}$ and let $\hat{w}_R$ and $\hat{w}_D$ be defined analogously. Then,

similarly,

Hence, by Lemma 3,

from which the result follows.

Hence, we consider the SVD and QR method on the model output in Equation (2.12) under the additional assumption $W_R\neq W_D$. Moreover, we notice without loss of generality; we can take any weight in Equation (2.12) to be unitary. As a result, we set $W_R = 1$ and consider $W_D\neq 1$ with $W_I$ arbitrary.

The sensitivity matrix in Equation (2.11) can only be known exactly for systems that are simple enough to yield an analytic solution to the differential equations. For the purposes of our model, we approximate $\chi$ by performing a central difference derivative with step size $h$. When solving the differential systems under the perturbed parameters, we use a relative error tolerance of $h^2$ in order to ensure the numerical stability of the scheme ^[23].

3.   Results and discussion

3.1.   Practical identifiability analysis

The identifiability analysis procedure outlined in Section 2.3.1 is dependent upon the choice of weights $(W_I, W_R, W_D)$ used in the weighting of the model output function in Equation (2.12). As observed at the end of Section 2.3.1: to avoid duplicate columns in the sensitivity matrix in Equation (2.11) and to avoid duplication of the model output function in Equation (2.12), we consider $W_R = 1$, $W_D\neq 1$, with $W_I$ arbitrary. We performed the SVD and QR method across a lattice of weight choices where $W_I\in\{1, 2, \ldots, 100\}$ and $W_D\in\{2, 3, \ldots, 100\}$. After this search, we selected a (non-unique) weight choice that gave the largest set of identifiable parameters in $W_I = 100$, $W_R = 1$, $W_D = 10$. We also note that this choice of weighting overemphasizes the importance of accurately fitting the active infections and ranks fitting the fatalities as more important than fitting the recoveries. While this is not the only choice of weights one could consider, we believe it represents a reasonable decision to fixate on the importance of tracking active cases and accurately reporting fatalities.

We present the results of the identifiability analysis in Table 1. Notably, for the integer-order case, there are two commonly unidentifiable parameters (namely $x_v$ and $\theta_{20}$). Eastman et al. ^[26] had formulated $x_v$ to represent the "visible proportion of the infected population." Notably, this phenomenological factor was introduced to account for those individuals with active cases that have not been tested but nonetheless are seeding secondary cases. Given the difficulty of measuring this in another capacity, we simply set $x_v = 1$ for the purposes of performing the fitting procedure. Now, $\theta_{20}$ represents the final interpolating point in the public response function being considered. Given that $\tau_{20}$ corresponds with early February 2021, a time period wherein Ontario was in province-wide lockdown; we set $\theta_{20} = 1$ to force a maximal reduction in transmission by the general populace. Similarly, for the fractional-order case, we note that $x_v$, $\theta_{20}$, and $a_S$ are the common unidentifiable parameters. As such we consider $x_v = \theta_{20} = 1$ as in the integer-order case. For $a_S$, we consider three possibilities: $a_S = 0.7$, $a_S = 0.8$, and $a_S = 0.9$ and cross-compare the results in Section 3.3. We chose these values of $a_S$ between $0.7$ and $1$ as has been observed elsewhere in the literature (for instance, a value of $a_S = 0.725$ was found optimal in Rajagopal et al. ^[48] and various values of fractional-order between $0.85$ and $1$ were considered in Khan et al. ^[49] for similar compartmental models of COVID-19).

3.2.   Model calibration results

We present the result of the integer-order fitting procedure in Figure 1. In particular, this fit was obtained by setting $x_v = 1 = \theta_{20}$ with $\theta_0 = 0$, as in Table 2. The values of the parameters discovered by the fitting process used to construct this figure can be found in Table 3.

Next, we present the fits of the fractional-order fitting procedure. These fits were likewise obtained by setting $x_v = 1 = \theta_{20}$ with $\theta_0 = 0$. However, we also considered varying $a_S$ amongst $a_S = 0.7$, $a_S = 0.8$, and $a_S = 0.9$. The results of the fit with $a_S = 0.7$ are presented in Figure 2. As the plots are visually quite similar between the fractional-order cases with any deviations being difficult to see with the naked eye, we relegate plotting the results with $a_S = 0.8$ and the results with $a_S = 0.9$ to the supplement. The values of the parameters used to acquire these fits are reported in the appropriate columns of Table 3 for all three fractional-order fits.

3.3.   Comparisons between Integer and fractional-order methods

In qualitatively assessing the fits, we note the integer-order model struggles to hug both waves considered as evidenced by how consistently the infections line falls below the data between April 2020 and Sept 2020 in Figure 1. Moreover, our fitting method is biased towards fitting infections more closely when the infections are largest, as was the case in the wave of the pandemic in the months of November 2020–February 2021. Despite this, the fractional-order fits visually fit the infection curves to the data in both waves much more closely while still maintaining the appropriate curvature in the recovered curves. All models have difficulty fully matching all three curves at once, as evidenced by the gap between the fatalities curves and the data. By tuning the weights, different qualities of fit can be achieved for different purposes. For the purposes of this manuscript, we considered fitting infections as the most important.

As we can see from Table 3, the fractional-order fitting procedures all converged to the same form of the public response function $\theta(t)$. Moreover, Figure 3 shows a comparison of the exact values of the function $\theta(t)$ between the integer-order model and the fractional-order model. As we discussed in Section 3.3, the similarities between the fractional-order columns of Table 1 is an unsurprising result due to the practical unidentifiability of $a_S$. In this case, some parameters remain constant when the unidentifiable parameter is varied. Indeed, the results of Table 3 suggest that the unidentifiability of $a_S$ is due to a correlation between $a_S$ and $\beta_{DD}$ (or, equivalently, $\beta_{DD}$) in the fractional-order model.

Moreover, the fractional-order fitting procedures all discovered values of $\beta_{DD}$ that were very close to those of $\beta_{DM}$ (with changes only appearing in the fifth significant figure at the $10^{-14}$ level), suggesting that, at least in the fractional case, the model may be simplified by assuming $\beta_{DD} = \beta_{DM}$. We also note that the values of $\gamma$, $\alpha$, and $\delta$ varied only trivially for different values of $a_S$. For comparison with the integer-order method, we note that $\alpha$ and $\delta$ are very similar values between all four methods corresponding to a value of $1/(\alpha+\delta)\approx 12$ days between all four methods suggesting that an individual is infectious for, on average, 12 days. Note that due to the way recovered tests are processed in Ontario (namely, cases without any medical follow-up are marked as "recovered" if 14 days have passed since symptom onset), this value is artificially biased towards a value of 14. One change between the integer-order fit and the fractional-order fits can be found in the value of $\gamma$. For the fractional-order fits, this value suggests a mean incubation period of approximately one week, while for the integer-order fits, this value suggests a mean incubation period of approximately $2$-$3$ days.

While the exact values of the public response function $\theta(t)$ differ between the integer-order model and the fractional-order model, many of the qualitative features of the function are maintained, as is evidenced in Figure 3. In particular, considering the fractional-order $\theta(t)$ between April 2020 and August 2020, the value seems to oscillate around some baseline values. This suggests that the public response function could be replaced either with a constant value on this domain or a constant value with appropriately scaled sinusoidal oscillations. A similar phenomenon can be observed for the integer-order model in the same time frame, albeit with a different baseline through a similar oscillatory frequency. Both methods predict a drop in the public response function heading out of Summer 2020 into Autumn 2020, potentially due to the return of school. In the fractional-order case, the assumption that $\theta_{20} = 1$ appears to be more realistic than in the integer-order case, as evidenced by $\theta_{19}$ being near $1$ in the former but not the latter.

3.4.   Comparisons for prediction

For the purposes of ascertaining which of the models has the strongest predictive strength, we included additional historical data corresponding to an additional 42 days of time series data. This additional data was not used for fitting but only for evaluating the predictive power of the fit. To that end, the model was fit on the first 376 data points and evaluated on the remaining 42 data points. Now, the model requires a particular form of $\theta(t)$ for these additional data values. For the purposes of this comparison, we considered two such future $\theta(t)$ schedules: one where $\theta(t > 376) = \theta(376)$ (i.e., $\theta(t)$ was fixed at $\theta = 1$ for $t\ge 376$) and another where $\theta(t)$ was fixed at $\theta = 1$ for half of the new time points $(376\le t\le 397)$ and $\theta = 0$ for half of the new time points $(397\le t)$. The motivation for these $\theta$ schedules is as follows: for the first schedule, we simply projected the last fit value of $\theta$ naively forward, and for the second schedule, we recognize that a stay-at-home order in Ontario was set to be lifted by mid-February 2021 as a result we would predict a reduced $\theta$ value corresponding to greater degrees of social mixing. For both of these schedules, we judged the quality of the prediction by considering the value under the cost function (2.10) for only the new time points. In both experiments, the fractional model corresponding to $\alpha_s = 0.8$ performed the best, and so it is the only one we report. For the first $\theta$, schedule $\chi_{\rm frac}\approx 0.26\, \chi_{\rm int}$, and for the second $\chi_{\rm frac}\approx 1.02\, \chi_{\rm int}$. This suggests that both the integer and fractional-order models are both similarly sensitive to the form of $\theta(t)$. As a result, the model presented is only appropriate for predicting over short time intervals (as sensitivities to the particular form of $\theta(t)$ prescribed for prediction begin to dominate shortly thereafter). However, in the first schedule considered where $\theta(t)$ is changing, the fractional-order model appears to do a better job at qualitatively describing the dynamics of the system. The results of these predictions are summarized in Figure 4. Note the similarities in the bottom plot suggesting that both methods perform similarly poorly when the prescribed form of $\theta(t)$ is a poor match to the data. Moreover, note that the presence of the so-called memory effect of fractional operators could be responsible for the increased predictive power for short-term predictions when the form of $\theta(t)$ is not perturbed. However, we also note that it is evident that the trajectories will soon begin to diverge from the data if the projection is carried forward in the top left figure in Figure 4.

3.5.   Time-dependent reproduction number

Modeling the spread of an epidemic allows one to calculate metrics that are useful for targeting control of the disease spread. Other than standard epidemiological data like active infections, recoveries, and deaths, another important metric to track is that of the reproduction number. Famously, the basic reproduction number (denoted ${\mathcal{R}_0}$) represents the number of expected secondary cases that would be seeded by a single infectious individual in an otherwise wholly susceptible population. This number is a function not only of the disease being considered but also represents behavioral quirks of the citizenry under investigation. As the pandemic has progressed, behavior patterns of the general public have changed to embrace transmission reduction protocols through the implementation of remote work, maintenance of social distance, the wearing of face masks, etc. As a result, even if the virus itself does not change, the number of secondary cases can change due to changing public behavior. Hence, while ${\mathcal{R}_0}$ is an important number to calculate, the time-dependent effective reproduction number ${\mathcal{R}_t}$ is often of greater interest. This is especially true later in a pandemic as the decrease in transmission due to the boost in immunity from previous infection and vaccination is captured by ${\mathcal{R}_t}$ but not by ${\mathcal{R}_0}$. In fact, ${\mathcal{R}_t}$ is often argued to be the most important number to track in order to manage an epidemic outbreak. Broadly, ${\mathcal{R}_t}$ represents the expected number of secondary cases seeded by a single infectious individual at time $t$ during the pandemic. If ${\mathcal{R}_t} < 1$ is maintained for a prolonged period of time, then one can expect that transmission of the virus is slowing. If the actions of the public keep ${\mathcal{R}_t}$ below $1$ for long enough, then the virus will eventually die out in the population. In contrast, if ${\mathcal{R}_t} > 1$, then the transmission is increasing within the population.

Here we used a Bayesian method that was developed by Bettencourt and Riberio ^[50] to estimate the effective reproduction number ${\mathcal{R}_t}$ from each of the trajectories of the four fit procedures. We then plot the mean value of these four procedures along with the union of the 95% confidence intervals of each procedure. This visualization is presented in Figure 5.

Following Bettencourt and Riberio, we assumed that new infection cases arise according to a Poisson distribution with a mean value equal to $0.25$. This corresponds to a serial interval of $4$ days for COVID-19 ^[51]. Further, it is assumed that the distribution of ${\mathcal{R}_t}$ is a Gaussian centered around ${\mathcal{R}_{t-1}}$ with a standard deviation $\sigma = 0.1$. As the visualization in Figure 5 demonstrates, ${\mathcal{R}_t}$ remained above the bifurcation value of $1$ for most of the Spring of 2020 and the Autumn of 2020 before falling below 1 briefly in the early months of 2021.

4.   Conclusion

We considered a modified model of SARS-CoV-2 spread in the Canadian province of Ontario. Under this expanded model, we demonstrated that a large number of the parameters were practically identifiable from the data. We then fit these parameter values using a combination of a genetic algorithm and a simplex method for the identifiable parameters and assigned reasonable choices to the unidentifiable parameters. Moreover, we expand our model to accurately describe the dynamics of multiple waves of epidemic infection with a modified fractional-order model. Thus, the fitting process was repeated for an integer-order model as well as three fractional-order models using the discretized Caputo operator. While the Caputo operator has been observed to be beneficial in mathematical biological systems, the results of the fitting procedures between these methods are qualitatively quite similar, though a modest decrease in cost function value of approximately 27% was observed by considering the fractional-order operator over the integer-order one. Future research should consider the effects of vaccination in our model. Moreover, our proposed model could be extended to account for hospitalized, quarantined, or home isolation individuals. With respect to the fractional-order model, comparing the effects of different definitions of the fractional derivative operator on the physical behavior of the model trajectories would be of interest.

Acknowledgments

Financial support by the Natural Sciences and Engineering Research Council of Canada (NSERC)(MK) is gratefully acknowledged.

Conflict of interest

The authors declare there is no conflict of interest.