Fractional Growth Model with Delay for Recurrent Outbreaks Applied to COVID-19 Data

Abstract

There are a great many epidemiological models that have been implemented to describe COVID-19 data; however, few attempted to reproduce the entire phenomenon due to the complexity of modeling recurrent outbreaks. In this work a fractional growth model with delay is developed that implements the Caputo fractional derivative with $0<\beta \le 1$. Furthermore, in order to preserve the nature of the phenomenon and ensure continuity in the derivatives of the function, a method is proposed to construct an initial condition function to implement in the model with delay. This model is analyzed and generalized to model recurrent outbreaks. The model is applied to fit data of cumulative confirmed cases from Mexico, the United States, and Russia, obtaining excellent fitting corroborated by the coefficient of determination, where $R^2>0.9995$ in all cases. Lastly, as a result of the implementation of the delay effect, the global phenomenon was decomposed into its local parts, allowing for directly comparing each outbreak and its different characteristics.

1. Introduction

Models for growth phenomena characterize dynamics in which the elements they model grow; these elements can be animals, plants, cells, tumors, goods, services, and more recently, individuals infected with COVID-19. These models belong to areas of science such as biology, medicine, economics, physics, and mathematics. These phenomena continue to be a challenge that allows for the creation of increasingly complex models [1,2,3,4,5].

One of the best growth models is the Malthus model [6], where the population growth rate is proportional to the same population, and this model has unbounded growth. Considering the case of growth models with bounded growth, some of the best known models are the Verhulst or logistic model [7] and the Gompertz model [8]; for more models, see [1,2]. These models have sigmoidal behavior and thereby an inflection point, which can be fixed, as in the logistic and Gompertz models, or nonfixed, which is determined from additional parameters of the model.

Previous models were solved by several methods [4,5,9,10,11,12,13,14,15,16,17] and the references therein; these have also been subject to variations such as the incorporation of delay in order to more accurately reflect reality in the phenomena [18,19].

Fractional calculus, although it is as old as the calculus itself, has recently been applied in all areas of knowledge; in particular, for growth models, the order of the fractional derivative is an indicator of how the phenomenon is behaving in relation to classical behavior (with an integer derivative) [4,13,15,16], thus defining the memory effect, where an order of the fractional derivative less than unity is defined as a subdiffusive phenomenon generating an attenuation effect in terms of the characteristics of the phenomenon; if the order of the fractional derivative is greater than unity, a superdiffusive phenomenon is obtained, and characteristics of the phenomenon are accelerated.

Fractional differential equations with delay are a developing area in which results continue to be generalized for linear and nonlinear models, and models with multiple lags [20,21,22,23,24,25,26].

COVID-19 is a recent disease caused by the SARS-CoV-2 virus that was declared as a public health emergency of international importance by the World Health Organization (WHO) on 30 January 2020 [27]. Since then, various models were implemented in order to understand the behavior of this phenomenon, from growth models such as the logistic and Gompertz [28,29,30], models with fractional derivatives [31,32] and models with delay [33,34,35].

Alcántara-López et al. developed a fractional growth model with a nonfixed inflection point in which they fitted data of cumulative confirmed cases by COVID-19 from Mexico, the United States, and Russia, modeling multiple outbreaks and obtaining very good adjustments; however, by modeling multiple outbreaks as a global phenomenon, they implicitly assumed that all outbreaks had already been determined from the first day of the pandemic [32].

The present work proposes a growth model that contains, as in particular cases, the logistic and Gompertz models, in which the Caputo fractional derivative and the delay effect are incorporated. This model is analyzed and implemented to fit COVID-19 data from multiple countries from the start of the pandemic until 8 December 2021, capturing the phenomenon of multiple outbreaks.

This work is organized as follows: Section 2 describes the necessary mathematical tools corresponding to fractional calculus and delay equations; Section 3 develops the proposed growth model, its properties, and behavior for each parameter; Section 4 generalizes the proposed model to model phenomena with multiple outbreaks; Section 5, applies the proposed models to fit the data of confirmed cases by COVID-19 from Mexico, the United States, and Russia; lastly, Section 6, summarizes our conclusions.

2. Mathematical Preliminaries

In this section, the tools used throughout this work pertaining to fractional calculus and delayed fractional differential equations are mentioned.

2.1. Fractional Calculus

There are a large number of definitions of fractional derivatives [36]; however, the Caputo fractional derivative is widely used to model real phenomena because the initial conditions are the same as for classical phenomena, with an integer derivative. For more information on the properties of this and other fractional derivatives, see [37,38].

Definition 1.

Let $-\infty <t_0<b<\infty$. Riemann–Liouville fractional integrals ${}^{RL}{I}_{t_0+}^{\beta }y\left(t\right)$ and ${}^{RL}{I}_{b-}^{\beta }y\left(t\right)$ of order $\beta \in \mathbb{R}$ are defined by

$$\left({}^{RL}{I}_{t_0+}^{\beta }y\right)\left(t\right)=\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_{t_0}^t{(t-s)}^{\beta -1}y\left(s\right)ds,t>t_0,$$

(1)

and

$$\left({}^{RL}{I}_{b-}^{\beta }y\right)\left(t\right)=\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_t^b{(s-t)}^{\beta -1}y\left(s\right)ds,t<b,$$

(2)

where $\mathsf{\Gamma }(·)$ is Euler’s gamma function. These integrals are called left-sided and right-sided fractional integrals.

The Caputo fractional derivative, expressed from the fractional integral, is as follows:

Definition 2.

Let $-\infty <t_0<b<\infty$. Caputo fractional derivatives ${}^C{D}_{t_0+}^{\beta }y\left(t\right)$ and ${}^C{D}_{b-}^{\beta }y\left(t\right)$ of order $\beta \in \mathbb{R}$ are defined by

$$\left({}^C{D}_{t_0+}^{\beta }y\right)\left(t\right)=\left({}^{RL}{I}_{t_0+}^{n-\beta }\frac{d^n}{d{t}^n}y\right)\left(t\right)=\frac{1}{\mathsf{\Gamma }(n-\beta )}{\int }_{t_0}^t{(t-s)}^{n-\beta -1}{y}^{\left(n\right)}\left(s\right)ds,$$

(3)

and

$$\left({}^C{D}_{b-}^{\beta }y\right)\left(t\right)={(-1)}^n\left({}^{RL}{I}_{b-}^{n-\beta }\frac{d^n}{d{t}^n}y\right)\left(t\right)=\frac{{(-1)}^n}{\mathsf{\Gamma }(n-\beta )}{\int }_t^b{(s-t)}^{n-\beta -1}{y}^{\left(n\right)}\left(s\right)ds,$$

(4)

where $n\in \mathbb{N}$ with $n-1<\beta \le n$ and $y^{\left(n\right)}$ is the n-th derivative. These derivatives are called the left-sided and right-sided Caputo fractional derivatives of order β.

Some important properties of the fractional integral and the fractional derivative are listed below.

• Fractional integral applied to a polynomial had the following analytic expression:

  $$\begin{array}{c}\hfill \left({}^{RL}{I}_{t_0+}^{\alpha }{(t-t_0)}^{\beta }\right)\left(t\right)=\frac{\mathsf{\Gamma }(\beta +1)}{\mathsf{\Gamma }(\beta +\alpha +1)}{(t-t_0)}^{\beta +\alpha },\end{array}$$

  (5)

  $$\begin{array}{c}\hfill \left({}^{RL}{I}_{b-}^{\alpha }{(b-t)}^{\beta }\right)\left(t\right)=\frac{\mathsf{\Gamma }(\beta +1)}{\mathsf{\Gamma }(\beta +\alpha +1)}{(b-t)}^{\beta +\alpha }.\end{array}$$

  (6)
• If $\beta =n\in \mathbb{N}$ then $\left({}^C{D}_{t_0+}^{\beta }y\right)\left(t\right)=y^{\left(n\right)}\left(t\right)$ and $\left({}^C{D}_{b-}^{\beta }y\right)\left(t\right)={(-1)}^n{y}^{\left(n\right)}\left(t\right)$.
• The fractional derivative of Caputo and the fractional integral behave as inverse operators, as follows:

  $$\begin{array}{c}\hfill \left({}^C{D}_{t_0+}^{\beta }{}^{RL}{I}_{t_0}^{\beta }y\right)\left(t\right)=y\left(t\right)\mathrm{and}\left({}^C{D}_{b-}^{\beta }{}^{RL}{I}_{b-}^{\beta }y\right)\left(t\right)=y\left(t\right).\end{array}$$
• If the order of the operators is interchanged, one must

  $$\left({}^{RL}{I}_{t_0}^{\beta }{}^C{D}_{t_0+}^{\beta }y\right)\left(t\right)=y\left(t\right)-\sum _{k=0}^{n-1}\frac{y^{\left(k\right)}\left(t_0\right)}{k!}{(t-t_0)}^k$$

  and

  $$\left({}^{RL}{I}_{b-}^{\beta }{}^C{D}_{b-}^{\beta }y\right)\left(t\right)=y\left(t\right)-\sum _{k=0}^{n-1}\frac{{(-1)}^k{y}^{\left(k\right)}\left(b\right)}{k!}{(b-t)}^k.$$

2.2. Delayed Fractional Differential Equations

To model phenomena where the integer derivative is not good enough, the fractional derivative with the properties described above is a useful tool for integrating the memory phenomenon and the delay effect, so these two tools are naturally combined, giving way to more robust models. For more information on delayed differential equations, see [39].

Definition 3.

The fractional delay differential equation (FDDE) is defined as

$$\begin{array}{ccc}\hfill \left({}^C{D}_{t_0+}^{\beta }y\right)\left(t\right)& =& f\left(t,y(t-\tau )\right),t>0,\tau >0,0<\beta \le 1;\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill y\left(t\right)& =& \varphi \left(t\right),-\tau \le t\le 0.\hfill \end{array}$$

(8)

Unlike phenomena without delay, the initial condition is a function, namely, $\varphi$; in many papers, it is handled as a constant. In this work, we propose a way to calculate function $\varphi$ as a function of the phenomenon to be modeled, that is, as a function of f.

In Model (7), if $\tau =0$, then the model returns to a fractional differential equation without delay; the existence and uniqueness of solutions to this model was extensively studied and resolved [40]. For case $\tau >0$, Yang and Cao [21] proved the existence and uniqueness of fractional delay differential equations, where they showed that, if f is continuous, by applying the theory of the fixed point, solutions exist; to ensure the uniqueness of solutions, the Lipschitz condition is added.

3. Fractional Growth Model with Delay

Definition 4.

Logistic and Gompertz growth models are defined through the following differential equations, respectively:

$$\begin{array}{ccc}\hfill \frac{dN}{dt}& =& rN\left(t\right)\left(1-\frac{N\left(t\right)}{N_{\infty }}\right),N\left(0\right)=N_0;\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill \frac{dN}{dt}& =& rN\left(t\right)ln\left(\frac{N_{\infty }}{N\left(t\right)}\right),N\left(0\right)=N_0;\hfill \end{array}$$

(10)

where $N=N\left(t\right)$ represents the size of the population, t is time, $r>0$ is the growth rate, $N_{\infty }>0$ is the maximal number of individuals that the population can sustain or carry the capacity of the environment, and $N_0$ is the initial condition.

These growth models are widely studied and applied to study and fit various phenomena; in particular, having sigmoidal growth, these models have an inflection point given by

$$\begin{array}{ccc}\hfill N_{\mathrm{inflection}}=& \frac{N_{\infty }}{2},\mathrm{in}{t}_{\mathrm{inflection}}=& \frac{1}{r}ln\left(\frac{N_{\infty }}{N_0}-1\right);\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill N_{\mathrm{inflection}}=& \frac{N_{\infty }}{e},\mathrm{in}{t}_{\mathrm{inflection}}=& \frac{1}{r}ln\left(ln\left(\frac{N_{\infty }}{N_0}\right)\right);\hfill \end{array}$$

(12)

where these models are shown to belong to the family of growth models with a fixed inflection point.

When considering limit

$$\underset{\mu \to \infty }{lim}\mu x\left(1-x^{1/\mu }\right)=-xln\left(x\right),$$

(13)

it is possible to generalize the right part in Models (9) and (10), namely,

$$f\left(N\left(t\right)\right)=r\mu N\left(t\right)\left(1-{\left(\frac{N\left(t\right)}{N_{\infty }}\right)}^{1/\mu }\right).$$

(14)

Therefore, the fractional growth model with delay is obtained by using the Caputo fractional derivative instead of the integer derivative and applying the delay effect on N in function $f(t,N)$, that is,

$$\begin{array}{ccc}\hfill {}^C{D}_0^{\beta }N\left(t\right)& =& r\mu N(t-\tau )\left(1-{\left(\frac{N(t-\tau )}{N_{\infty }}\right)}^{1/\mu }\right),t>0\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill N\left(t\right)& =& \varphi \left(t\right),t\in [-\tau ,0];\hfill \end{array}$$

(16)

where N, $N_{\infty }$, and t have the same meaning as in Equations (9) and (10); $r=\nu /{\rho }^{\beta -1}$, where $\nu$ is a growth rate, $\rho$ is a reference time introduced in order to maintain dimensional balance in the equation, $\tau$ is constant, and $\mu$ is a dimensionless shape parameter, such that if $\mu =1$ then the model simplifies to the fractional logistic model; if $\mu \to \infty$, then the model is equivalent to the fractional Gompertz model.

Models (15) and (16) with $\tau =0$ were extensively studied in [32], where it was proven that the model belongs to the family of growth models with a nonfixed inflection point, the existence and uniqueness of the solution and its numerical stability were proved, the numerical scheme used was the result of applying the fractional integral in Equation (15), thus converting the differential equation into a Volterra-type equation with the following form:

$$N\left(t\right)=\varphi \left(0\right)+\frac{1}{\mathsf{\Gamma }\left(\beta \right)}{\int }_0^t{(t-s)}^{\beta -1}f\left(N(s-\tau )\right)ds,$$

(17)

where f is as in Equation (14); if $t\in [0,\tau ]$ then $f\left(N\right(s-\tau \left)\right)=f\left(\varphi \right(s\left)\right)$.

3.1. Sensitivity Analysis and Delay Effect

In order to understand the delay effect in the fractional growth model with delay, we show the effect of fractional derivative $\beta$ in combination with delay parameter $\tau$.

Figure 1 shows the effect of delay parameter $\tau$ on the fractional growth model with delay considering classical derivative $\beta =1$, and a constant initial function $\varphi \left(t\right)=\varphi$. Figure 1a shows function $N\left(t\right)$ where, for $\tau =0$, classical behavior was observed; that is, the function had sigmoidal growth, where parameters r and $\mu$ determined the way in which the function reached the inflection point and subsequently limited value $N_{\infty }$.

For $\tau >0$, function $N\left(t\right)$ can be examined in two parts: for $t-\tau \le 0$, function $f(s,N(s-\tau \left)\right)=f(s,\varphi )$ is constant. Applying Equation (5) to Equation (17) shows that $N\left(t\right)$ is a linear function; for $t-\tau >0$; the function behaves similar to the case without delay with initial condition $N_0=\varphi +f\left(\varphi \right)\frac{{\tau }^{\beta }}{\mathsf{\Gamma }(\beta +1)}$, that is, the function is a copy of the case without delay from the initial condition and with a time shift equal to the value of $\tau$.

Figure 1b shows how inflection point $N_{\mathrm{inflection}}$ was the same in all functions and only had one shift in time at which it occurred, $t_{\mathrm{inflection}}$.

Figure 2 shows the effect of the fractional derivative in combination with the delay effect and constant initial function $\varphi \left(t\right)=\varphi$. Figure 2a, with $\beta =1.0$, shows the behavior described in Figure 1; with $t-\tau \le 0$, the function had linear behavior, while with $t>\tau$, the function had classical behavior with initial condition $N_0$. The other cases showed that function $N\left(t\right)$ was strongly influenced in all its characteristics by the value of the fractional derivative; with a lower value of $\beta$, the moment at which the inflection point was reached and the moment at which the function reached the limit value, which was also strongly decreased, were delayed.

Figure 2b shows the consequences of applying the fractional derivative in combination with the delay effect on the derivative of function $N\left(t\right)$. On the one hand, the aforementioned with Figure 2a was confirmed. With $t>\tau$, the inflection point was reduced and delayed the lower value of $\beta$; for $t<\tau$, the function had polynomial behavior with a value of $N\left(t\right)=\varphi +f\left(\varphi \right)\frac{t^{1-\beta }}{\mathsf{\Gamma }(\beta +1)}$. Lastly, at $t=\tau$, function $N^{\prime }\left(t\right)$ had discontinuity caused by the abitrarity of the initial function in combination with the way in which the constants were affected by the fractional integral.

From Morgado et al. [25], the following result could be obtained: the derivative of the function $N\left(t\right)$ can havea discontinuity at $t=\tau$ unless the following equality is satisfied:

$${\left.\frac{dN}{dt}\right|}_{t=\tau }=\frac{d}{dt}{}^{RL}{I}_{t_0}^{\beta }f(\tau ,\varphi \left(0\right)).$$

(18)

Equation (18) shows that, to solve a model that includes both a fractional derivative and the delay effect, it is not enough to observe the continuity of the function; the smoothness of this solution also needs to be observed. Therefore, an arbitrary initial function is not recommended, much less a constant initial function that generates as many problems as those already discussed.

3.2. Initial Function Construction

There is a need for an initial function and to fit data to fractional models with delays that are continuous and whose derivative is also continuous. In this section we propose a method in which, starting from a specific initial condition $N\left(0\right)=N_0$, an initial function $\varphi \left(t\right)$ can be constructed for $t\in [-\tau ,0]$ that complies with these requirements.

The idea behind this proposal consists of dividing the differential equation with delay into two parts: in the first, it is solved:

$$\begin{array}{ccc}\hfill {}^C{D}_{0-}^{\beta }N\left(t\right)& =& r\mu N\left(t\right)\left(1-{\left(\frac{N\left(t\right)}{N_{\infty }}\right)}^{1/\mu }\right),t\in [-\tau ,0]\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill N\left(0\right)& =& N_0;\hfill \end{array}$$

(20)

where ${}^C{D}_{0-}^{\beta }N\left(t\right)$ is the left fractional derivative of Caputo. The obtained function satisfies that it and its derivative are both continuous. We then define $\varphi \left(t\right)\equiv N\left(t\right)$ for $t\in [-\tau ,0]$ and Models (15) and (16) are solved with the obtained function.

Figure 3 shows the effect of delay parameter $\tau$ in the fractional growth model with delay, where the initial function was calculated from a punctual condition, and dashed vertical lines were added at $t=\tau$ for $\tau =0,50,100,200$ and 300. Figure 3a with $\tau =0$ shows the previously mentioned classical behavior; in cases with $\tau >0$ a displacement of function $N\left(t\right)$ with $\tau =0$ was observed for $t>\tau$; $t\le \tau$ showed how function $N\left(t\right)$ was “completed” while keeping the nature of function f (Equation (14)). Figure 3b confirms these claims; on the one hand, the function, in addition to being continuous, had a continuous derivative at all times; for cases $\tau >0$, the function was shown to be a displacement of function $N\left(t\right)$ with $\tau =0$, where the value of the displacement was equal to $\tau$. Likewise, by having a continuous derivative at all times, it was possible to discuss the inflection point that, like the function with $\tau =0$, was shown to have the same value $N_{\mathrm{inflection}}$, but it shifted a distance $\tau$.

Figure 4 shows the solution of the fractional growth model with nonzero delay, $\tau >0$, varying the order of the fractional derivative and its integer derivative $N^{\prime }\left(t\right)$. For $t>\tau$ the shown behavior is consistent with that shown in Figure 2, that is, for $\beta =1.0$, the function had exponential growth until it reached the inflection point, which became slower until limit value $N_{\infty }$ was reached; for $\beta <1.0$, the growth of the function suffers the memory effect expected by the fractional derivative, causing a delay in the appearance of the inflection point, its value, and the necessary time to reach the limit value, which was also affected and decreased by the value of fractional derivative $\beta$.

For $t\le \tau$, Figure 4a shows how the behavior of the function holds for all orders of the fractional derivative. On the other hand, Figure 4b shows how function $N\left(t\right)$ was not only continuous, but also differentiable with a continuous derivative for all values of the order of fractional derivative $\beta$.

4. Fractional Growth Model with Delay for Recurrent Outbreaks

Growth models, including the fractional growth model with delay, have the limitation of describing a single sigmoidal growth, thereby showing a single inflection point. However, phenomena such as the cumulative confirmed cases of COVID-19 showed behavior of multiple instances of sigmoidal growth that are reflected in multiple inflection points, i.e., cases of maximal confirmed daily cases.

Therefore, in this section, the fractional growth model with delay, Equation (15), is generalized to include the behavior of multiple instances of sigmoidal growth, where from a punctual initial condition, the method described in Section 3.2 was applied to construct initial function $\varphi \left(t\right)$.

We assumed a growth phenomenon with k number of sigmoidal growth, that is, k outbreaks. Since the number of cumulative confirmed individuals in a time t is the result of the sum of cumulative confirmed individuals by already developed outbreaks and the developing outbreak, the principle of superposition is applicable. Therefore, the global phenomenon was modeled as the sum of each local model, that is,

$$N\left(t\right)=\sum _{j=1}^k{N}_j\left(t\right),$$

(21)

where $N_j\left(t\right)$ for $j=0,\cdots ,k$ satisfies Equation (15) with its respective parameters.

Considering Equation (21), the growth model that incorporates the fractional derivative and delay is

$${}^C{D}_0^{\beta }N\left(t\right)=\sum _{j=1}^k{}^C{D}_0^{{\beta }_j}{N}_j\left(t\right)=\sum _{j=1}^k{r}_j{\mu }_j{N}_j(t-{\tau }_j)\left(1-{\left(\frac{N_j(t-{\tau }_j)}{N_{j,\infty }}\right)}^{1/{\mu }_j}\right),t>0.$$

(22)

In general, if the method proposed in Section 3.2 was not considered, the initial condition for the model shown in Equation (22) would be required to be of form

$$N_j\left(t\right)={\varphi }_j\left(t\right)\mathrm{for}t\in [-{\tau }_j,0]\mathrm{and}j=1,\cdots ,k;$$

(23)

that is, for each outbreak of the phenomenon, an initial function is required that is defined and is continuous at all times until the delay effect is covered, that is, for $t-{\tau }_j\le 0$.

However, applying the method to construct initial function $\varphi \left(t\right)$, a punctual initial condition is sufficient for each outbreak, i.e.,

$$N_j\left({\tau }_j\right)=N_{j,0}\mathrm{for}j=1,\cdots ,k;$$

(24)

where $N_{j,0}$ is a constant.

The solution to Model (22) with initial punctual conditions, Equation (24), exists and is unique. For this, the right-hand side of the equality in Equation (22) is continuous and conducting some algebraic operations to show that it satisfies the Lipschitz condition, which does not require much effort but was omitted due to space issues.

5. Applications to COVID-19 Data

Recently, the pandemic caused by COVID-19 has impacted the whole world, infecting a large number of people, and differently impacting each country by testing their sanitary and other measures to avoid contagion. In particular, the duration of the pandemic has allowed for recurrent outbreaks, where each outbreak must be independently analyzed since the creation of various currently existing vaccines.

In this section, the fractional growth model with delay for recurrent outbreaks was applied to fit cumulative confirmed COVID-19 cases data from multiple countries, namely, Mexico, the United States, and Russia. The effectiveness of the model describes growth phenomena with multiple instances of sigmoidal growth, resulting in good coincidence between the inflection point of each outbreak (i.e., the peak of the outbreak) and moments of maximal contagion. The outbreaks were also compared with each other, allowing for inferences among the behaviors of each outbreak in each country.

The used numerical scheme was the result of applying the fractional integral to the fractional derivative of Caputo, Property 4 in Section 2.1, that is, Models (22) and (23) were converted into integral Equation (17) and a quadrature was used.

To obtain the parameters that fit function $N\left(t\right)$ to the cumulative confirmed cases, optimization function curve_fit from Python was used.

This function uses nonlinear least squares to fit some function to data. In order to ensure good fitting and the lowest possible error, function $N^{\prime }\left(t\right)$ was adjusted to data that resulted from applying a moving average of 7 days to the daily confirmed cases. Lastly, to confirm the effectiveness of the obtained parameters, coefficient of determination $R^2$ was calculated between function $N\left(t\right)$ and cumulative confirmed cases.

5.1. Mexico Data

In Mexico, confirmed COVID-19 cases began on 28 February 2020; since then, information regarding infections, suspects, negatives, and deaths has been constantly updated [27,41].

Figure 5 shows confirmed COVID-19 cases from 28 February 2020 to 8 December 2021, where accumulated confirmed cases are shown in red, and daily cases are shown in blue.

Due to the long duration of the pandemic, cumulative confirmed cases data showed behavior of multiple instances of sigmoidal growth; in particular, when observing daily cases, 3 outbreaks were clearly shown.

Figure 6 shows the fitting to data of confirmed cases in Mexico by the multiple sigmoidal growth model with $k=3$ outbreaks. Figure 6a shows function $N\left(t\right)$, Equation (22), in comparison with the cumulative confirmed cases on the logarithmic scale; indeed, a very good fit was shown by the $N\left(t\right)$ function to the point that they were visually indistinguishable. Figure 6b shows the derivative of the function, $N^{\prime }\left(t\right)$, in comparison with the daily cases; although the data showed some dispersion, function $N^{\prime }\left(t\right)$ followed the behavior of the data, as could be observed during the periods of increased and decreased daily cases.

Table 1. Parameters obtained by the fractional growth model with delay for the adjustment of data of cumulative confirmed cases of Mexico including confidence intervals at 95% for each parameter.

Parameters	${\mathit{\tau }}_{\mathit{j}}$ [days]	${\mathit{N}}_{\mathit{j},0}$ [cccc] 1	${\mathit{\beta }}_{\mathit{j}}$	${\mathit{r}}_{\mathit{j}}$ [days ${}^{-\mathit{\beta }}$]	${\mathit{\mu }}_{\mathit{j}}$	${\mathit{N}}_{\mathit{j},\mathit{\infty }}$ [cccc] 1
Optimal value	0	$1.16\times {10}^0$	$9.20\times {10}^{-1}$	$2.24\times {10}^{-2}$	$1.58\times {10}^8$	$1.08\times {10}^6$	First Sprout $j=1$
Lower bound	0	$1.02\times {10}^0$	$9.06\times {10}^{-1}$	$2.14\times {10}^{-2}$	$1.54\times {10}^8$	$1.03\times {10}^6$
Upper bound	0	$1.31\times {10}^0$	$9.33\times {10}^{-1}$	$2.34\times {10}^{-2}$	$1.62\times {10}^8$	$1.12\times {10}^6$
Optimal value	215	$2.03\times {10}^4$	$9.37\times {10}^{-1}$	$5.04\times {10}^{-2}$	$1.21\times {10}^0$	$1.31\times {10}^6$	Second Sprout $j=2$
Lower bound	215	$9.08\times {10}^3$	$7.10\times {10}^{-1}$	$1.00\times {10}^{-5}$	$8.53\times {10}^{-1}$	$1.25\times {10}^6$
Upper bound	215	$3.16\times {10}^4$	$1.00\times {10}^0$	$1.21\times {10}^{-1}$	$1.56\times {10}^0$	$1.36\times {10}^6$
Optimal value	440	$1.06\times {10}^4$	$9.89\times {10}^{-1}$	$4.56\times {10}^{-2}$	$1.36\times {10}^0$	$1.43\times {10}^6$	Third Sprout $j=3$
Lower bound	440	$5.17\times {10}^3$	$6.33\times {10}^{-1}$	$1.00\times {10}^{-5}$	$1.05\times {10}^0$	$1.39\times {10}^6$
Upper bound	440	$1.60\times {10}^4$	$1.00\times {10}^0$	$1.55\times {10}^{-1}$	$1.68\times {10}^0$	$1.48\times {10}^6$

¹ cccc = Cumulative Confirmed COVID-19 cases.

shows the obtained parameters by fitting the fractional growth model with delay, considering $k=3$ outbreaks, to data of cumulative confirmed cases for Mexico; results are shown in Figure 6.

For each outbreak, delay parameters $\tau =$ 0, 215, and 440 proposed by inspecting the data were considered; the proposal to the delay parameter did not disturb the fitting since, when modifying delay value $\tau$, the punctual initial condition value associated with the outbreak was modified. Values of $r_j$ indicated the growth rate for each shoot $j=1,2,3$; by Equation (15), if ${\mu }_j\approx 1$, the j-th outbreak had logistic-type behavior; if ${\mu }_j\gg 1$, then the j-th outbreak had Gompertz-type behavior; $N_{j,\infty }$ indicated the maximal number of individuals to be infected in the j-th outbreak; lastly, ${\beta }_j$ reflects the memory effect contained in the j-th outbreak, where ${\beta }_j=1$ indicates classical behavior or no-memory effect. If ${\beta }_j<1$, the memory effect allows for a delay in the characteristics of the phenomenon as described in Section 3.1; the effect is greater when the value of the parameter is lower.

In order to show the uncertainty of each parameter, Table 1 shows, in addition to the found fit values, the lower and upper bounds of each parameter considering an uncertainty level of $95\%$; these confidence intervals were calculated from Student’s t distribution.

Table 4 reflects the effectiveness of the fitting by the model, and it shows dates corresponding to the days of maximal contagion in each outbreak and dates corresponding to the inflection points of function $N\left(t\right)$. Coefficient of determination $R^2>0.9999$ was obtained, which was an excellent fit by the model; furthermore, the difference between the date of maximal recorded contagion and the date corresponding to the inflection point of each outbreak was 13, 11, and 6 days, respectively, which reaffirmed the accuracy of the model.

Figure 7 shows the derivative of function $N_j$ with $j=1,2,3$ shifted to the left by a value of $\tau$; this facilitates highlighting the advantages of applying the delay effect to a multiple sigmoidal growth model. Therefore, outbreaks can be directly compared. The first outbreak, despite having a lower growth rate r, had a higher value of $\mu$, which indicates that it had a higher contagion rate $r\mu$ of order $1e6$, which is contrary to the second and third outbreaks, where the contagion rate was of the order of $1e-1$; however, the value of fractional derivative $\beta$ was the lowest compared to those of subsequent outbreaks, which can be attributed to implemented sanitary measures.

5.2. US Data

In the United States, confirmed COVID-19 cases began on 22 January 2020. Figure 8 shows data for confirmed COVID-19 cases from 22 January 2020 to 8 December 2021, where cumulative cases are shown in red, while daily cases are shown in blue.

For data fitting, the data from 29 February 2020 were considered as, before this date, cumulative data did not have exponential behavior; therefore, a growth model could not be applied. Considering these data, behavior of multiple instances of sigmoidal growth was clearly observed; in particular, $k=5$ outbreaks were considered.

Figure 9 shows the results of fitting to the confirmed cases from the United States by the multiple sigmoidal growth model considering $k=5$ outbreaks. Figure 9a shows function $N\left(t\right)$, Equation (22), superimposed on cumulative data cases on a logarithmic scale, showing an excellent fit, where practically, from the beginning of time, the function and data are indistinguishable. Figure 9b shows its derivative, that is, $N^{\prime }\left(t\right)$ compared to daily cases on the same dates where, despite the dispersion of the data, function $N^{\prime }\left(t\right)$ followed the nature of the data reproducing increases and decreases in daily cases, and dates of maximal contagion.

Table 2. Parameters obtained by the fractional growth model with delay for the adjustment of data of cumulative confirmed cases of the United States, including confidence intervals at $95\%$ for each parameter.

Parameters	${\mathit{\tau }}_{\mathit{j}}$ [days]	${\mathit{N}}_{\mathit{j},0}$ [cccc] 1	${\mathit{\beta }}_{\mathit{j}}$	${\mathit{r}}_{\mathit{j}}$ [days ${}^{-\mathit{\beta }}$]	${\mathit{\mu }}_{\mathit{j}}$	${\mathit{N}}_{\mathit{j},\mathit{\infty }}$ [cccc] 1
Optimal value	0	$9.74\times {10}^{-1}$	$7.17\times {10}^{-1}$	$8.73\times {10}^{-2}$	$1.08\times {10}^8$	$2.76\times {10}^6$	First Outbreak $j=1$
Lower bound	0	$8.41\times {10}^{-1}$	$6.81\times {10}^{-1}$	$7.48\times {10}^{-2}$	$1.03\times {10}^8$	$2.34\times {10}^6$
Upper bound	0	$1.11\times {10}^0$	$7.52\times {10}^{-1}$	$9.97\times {10}^{-2}$	$1.13\times {10}^8$	$3.18\times {10}^6$
Optimal value	80	$9.70\times {10}^4$	$1.0\times {10}^0$	$3.80\times {10}^{-2}$	$1.70\times {10}^0$	$5.27\times {10}^6$	Second Outbreak $j=2$
Lower bound	80	$1.0\times {10}^0$	$1.0\times {10}^{-1}$	$1.00\times {10}^{-5}$	$1.00\times {10}^{-5}$	$3.27\times {10}^6$
Upper bound	80	$3.28\times {10}^5$	$1.00\times {10}^0$	$2.63\times {10}^{-1}$	$5.00\times {10}^0$	$7.27\times {10}^6$
Optimal value	200	$3.77\times {10}^5$	$9.64\times {10}^{-1}$	$4.45\times {10}^{-2}$	$1.31\times {10}^0$	$2.42\times {10}^7$	Third Outbreak $j=3$
Lower bound	200	$1.00\times {10}^0$	$7.30\times {10}^{-1}$	$1.00\times {10}^{-5}$	$5.70\times {10}^{-1}$	$2.24\times {10}^7$
Upper bound	200	$8.02\times {10}^5$	$1.00\times {10}^0$	$1.08\times {10}^{-1}$	$2.04\times {10}^0$	$2.60\times {10}^7$
Optimal value	380	$1.76\times {10}^5$	$9.84\times {10}^{-1}$	$4.01\times {10}^{-2}$	$4.39\times {10}^5$	$4.39\times {10}^6$	Fourth Outbreak $j=4$
Lower bound	380	$9.26\times {10}^3$	$7.80\times {10}^{-1}$	$1.00\times {10}^{-5}$	$1.00\times {10}^{-5}$	$3.57\times {10}^6$
Upper bound	380	$3.43\times {10}^5$	$1.00\times {10}^0$	$9.53\times {10}^{-2}$	$8.67\times {10}^6$	$5.21\times {10}^6$
Optimal value	480	$1.00\times {10}^4$	$9.16\times {10}^{-1}$	$6.40\times {10}^{-2}$	$7.73\times {10}^0$	$1.29\times {10}^7$	Fifth Outbreak $j=5$
Lower bound	480	$1.0\times {10}^0$	$1.00\times {10}^{-1}$	$1.00\times {10}^{-5}$	$1.00\times {10}^{-5}$	$1.15\times {10}^7$
Upper bound	480	$3.32\times {10}^4$	$1.0\times {10}^0$	$9.60\times {10}^{-1}$	$1.94\times {10}^1$	$1.44\times {10}^7$

¹ cccc = Cumulative Confirmed COVID-19 cases.

shows obtained parameters by fitting the fractional growth model with delay, considering $k=5$ outbreaks, to the data of cumulative confirmed cases from the United States; results are shown in Figure 9.

For each outbreak, the used delay parameters were ${\tau }_j=0$, 80, 200, 380, and 480, proposed by inspection; considering these delay values, initial punctual conditions $N_{j,0}$ were obtained for the j-th outbreak with $j=1,\cdots ,5$; if the value of ${\tau }_j$ was modified, the value of $N_{j,0}$ was also directly modified.

In order to show the uncertainty of each parameter, Table 2 shows, in addition to the found fit values, the lower and upper bounds of each parameter considering an uncertainty level of $95\%$; these confidence intervals were calculated from Student’s t distribution.

Table 4 shows coefficient of determination $R^2>0.9996$, which indicates the effectiveness of the fitting by the model. Additionally, dates of maximal recorded contagion and dates corresponding to the inflection points of each outbreak are shown, showing that the difference between these dates was 6, 13, 19, 1, and 2 days, respectively, thus reaffirming the precision of the model.

Figure 10, similar to Figure 7, shows functions $N_j^{\prime }\left(t\right)$ for $j=1,\cdots ,5$ in order to compare the different outbreaks obtained in the fitting. In this comparison, the first and fourth outbreaks had Gompertz-type behavior since ${\mu }_j\gg 1$, while other outbreaks had logistic-type behavior ${\mu }_j\approx 1$. Furthermore, in combination with the growth rate of each outbreak, $r_j$, outbreaks with a Gompertz-type behavior had a higher contagion rate, $r_j{\mu }_j$, which can be seen by observing how quickly it arrived to the inflection point. Lastly, the order of the fractional derivative, ${\beta }_j$ showed that the second outbreak had classical behavior, ${\beta }_2=1$; the first outbreak obtained the lowest value in this parameter, implying that, in this outbreak, the memory effect was greater, allowing for a delay in all the characteristics of the phenomenon, decreasing this effect for the fifth, third, and fourth outbreaks, in this order.

5.3. Russia Data

In Russia, the data of confirmed cases by COVID-19 started from 31 January 2020. Figure 11 shows these data until 8 December 2021, where cumulative confirmed cases are shown in red, while daily confirmed cases are shown in blue.

For the fitting of the model, the data from 4 March 2020 to 8 December 2021 were considered, since before the proposed date, the data did not show exponential growth; therefore, a growth model could not be applied. Considering these data, behavior of multiple instances of sigmoidal growth was observed; in particular, $k=4$ outbreaks were considered.

Figure 12 shows the fitting to data of confirmed cases in Russia by the fractional growth model with delay for recurrent outbreaks considering $k=4$ outbreaks. Figure 12a shows the result of the model in comparison with cumulative cases on a logarithmic scale, showing excellent fit by the model. Figure 12b shows the comparison with the derivative of function $N\left(t\right)$, that is, the $N^{\prime }\left(t\right)$ function was compared with the reported daily cases.Data from Russia showed the least dispersion compared to data from Mexico and US, and function $N^{\prime }\left(t\right)$ followed the behavior of the data, reproducing increases and decreases in daily cases, and the maximal contagion values of each outbreak.

Table 3. Parameters obtained by the fractional growth model with delay for the adjustment of data of cumulative confirmed cases of Russia including confidence intervals at $95\%$ for each parameter.

Parameters	${\mathit{\tau }}_{\mathit{j}}$ [days]	${\mathit{N}}_{\mathit{j},0}$ [cccc] 1	${\mathit{\beta }}_{\mathit{j}}$	${\mathit{r}}_{\mathit{j}}$ [days ${}^{-\mathit{\beta }}$]	${\mathit{\mu }}_{\mathit{j}}$	${\mathit{N}}_{\mathit{j},\mathit{\infty }}$ [cccc] 1
Optimal value	0	$1.01\times {10}^0$	$7.81\times {10}^{-1}$	$5.17\times {10}^{-2}$	$1.63\times {10}^8$	$1.32\times {10}^6$	First Outbreak $j=1$
Lower bound	0	$9.88\times {10}^{-1}$	$7.75\times {10}^{-1}$	$5.03\times {10}^{-2}$	$1.62\times {10}^8$	$1.27\times {10}^6$
Upper bound	0	$1.94\times {10}^0$	$7.87\times {10}^{-1}$	$5.30\times {10}^{-2}$	$1.65\times {10}^8$	$1.36\times {10}^6$
Optimal value	180	$3.36\times {10}^4$	$9.86\times {10}^{-1}$	$2.57\times {10}^{-2}$	$3.02\times {10}^0$	$3.61\times {10}^6$	Second Outbreak $j=2$
Lower bound	180	$2.09\times {10}^4$	$8.20\times {10}^{-1}$	$5.71\times {10}^{-4}$	$2.04\times {10}^0$	$3.84\times {10}^6$
Upper bound	180	$4.63\times {10}^4$	$1.00\times {10}^0$	$5.08\times {10}^{-2}$	$4.01\times {10}^0$	$3.74\times {10}^6$
Optimal value	420	$1.02\times {10}^5$	$9.89\times {10}^{-1}$	$3.56\times {10}^{-2}$	$1.28\times {10}^0$	$2.60\times {10}^6$	Third Outbreak $j=3$
Lower bound	420	$5.02\times {10}^4$	$9.22\times {10}^{-1}$	$1.71\times {10}^{-2}$	$2.30\times {10}^{-1}$	$2.11\times {10}^6$
Upper bound	420	$1.53\times {10}^5$	$1.00\times {10}^0$	$5.40\times {10}^{-2}$	$2.33\times {10}^0$	$3.10\times {10}^6$
Optimal value	550	$8.23\times {10}^4$	$9.90\times {10}^{-1}$	$2.36\times {10}^{-2}$	$2.29\times {10}^1$	$4.31\times {10}^6$	Fourth Outbreak $j=4$
Lower bound	550	$1.0\times {10}^0$	$9.43\times {10}^{-1}$	$1.61\times {10}^{-2}$	$1.00\times {10}^{-5}$	$3.87\times {10}^6$
Upper bound	550	$2.14\times {10}^5$	$1.00\times {10}^0$	$3.11\times {10}^{-2}$	$1.84\times {10}^2$	$4.74\times {10}^6$

¹ cccc = Cumulative Confirmed COVID-19 cases.

shows obtained parameters by fitting the fractional growth model with delay, considering $k=4$ outbreaks, to data of cumulative confirmed cases for Mexico; results are shown in Figure 12.

Table 4. Statistical results of the fit of the fractional growth model with delay for recurrent outbreaks.

Country	${\mathit{R}}^2$	Forecast Peak	Real Peak
Mexico	$0.9999$	19 July 2020, 11 January 2021, and 13 August 2021	1 August 2020, 22 January 2021, and 19 August 2021
US	$0.9996$	13 April 2020, 1 August 2020, 20 December 2020, 10 April 2021, and 31 August 2021	7 April 2020, 19 July 2020, 8 January 2021, 11 April 2021, 29 August 2021
Russia	$0.9999$	17 May 2020, 10 December 2020, 18 July 2021, and 6 November 2021	11 May 2020, 24 December 2020, 9 July 2021, and 6 November 2021

shows the result of the determination coefficient, where $R^2>0.9999$, reflecting the excellent effectiveness of the model; this, added to the good fit by the model to data from Mexico and the US, shows the efficacy of the model to describe phenomena with multiple instances of sigmoidal growth. Additionally, dates of maximal contagion and dates corresponding to the inflection points of each outbreak are shown, which show that the difference between these dates was 6, 14, 9, and 0 days, respectively, reaffirming the precision of the model.

Figure 13 shows the $N_j^{\prime }\left(t\right)$ functions for $j=1,\cdots ,4$ shifted to the left by a value of ${\tau }_j$, respectively, showing the advantages of implementing the delay effect in a model with multiple instances of sigmoidal growth in order to directly compare the different outbreaks that compose the global phenomenon and different characteristics.

Obtained parameters by the fitting, shown in Table 3, indicate that the first outbreak had Gompertz-type behavior, ${\mu }_1\gg 1$; the other outbreaks had logistic-type behavior, ${\mu }_j\approx 1$; this, in combination with $r_j$, implied hat the first outbreak had a higher contagion rate $r_j{\mu }_j$, followed by the fourth, second, and third outbreaks. Lastly, the order of the fractional derivative of each outbreak, ${\beta }_j$ shows that the first outbreak was the one that had the greatest memory effect, allowing for there to be a delay in the various characteristics of the phenomenon, ${\beta }_1=0.781$; for the following outbreaks, a similar value of ${\beta }_j\approx 0.99$ was obtained, which indicated that, on the one hand, the memory effect was almost zero (since the value of ${\beta }_j$ was very close to 1), and on the other hand, it could reflect the consistency of sanitary measures implemented in Russia.

In order to show the uncertainty of each parameter; Table 3 shows, in addition to the found fit values, the lower and upper bounds of each parameter considering an uncertainty level of $95\%$; these confidence intervals were calculated from Student’s t distribution.

Lastly, Figure 12b shows that, between the third and fourth outbreaks, there was no large decrease in confirmed daily cases, unlike decreases in previous outbreaks or decreases between outbreaks in Mexico and US. Considering the date of this event, it can be associated with the impact of the implementation of vaccines applied to this country, evidencing a possibility of improvement in this matter.

6. Conclusions

A fractional growth model with delay was proposed where the Caputo fractional derivative was implemented with $0<\beta \le 1$. The proposed model generalizes Gompertz and logistic growth models, and it belongs to the class of models with a nonfixed inflection point. In order to ensure the smoothness of the function, that is, the continuity of the derivatives of the function, a method was proposed to construct an initial condition function for the model with delay from a punctual initial condition, allowing for the nature of the growth model to be preserved. The model developed with the method to construct the initial function was generalized to model phenomena with recurrent outbreaks. This model was implemented to describe data of cumulative confirmed cases of COVID-19 from Mexico, the United States, and Russia from the beginning of the pandemic in each country until 8 December 2021; the fitting showed great precision, corroborated by respective calculated correlation coefficients, where $R^2>0.9995$ in all cases. Likewise, the derivative of the function, that is, $N^{\prime }\left(t\right)$, was compared with the daily confirmed cases of each country, and in addition to an excellent fitting, the model could reproduce moments of maximal contagion through the inflection points of each outbreak. Lastly, as a result of the implementation of the delay effect, it was possible to decompose the global phenomenon into its local parts, allowing for directly comparing each outbreak and its respective characteristics.