Nonlinear differential equations with neutral term: Asymptotic behavior of solutions

1.   Introduction

COVID-19 is a contagious disease that can be easily transmitted from one person to another. At the end of 2019, Wuhan, China, was the site of this terrible infection for the first time. Later, this disease gave rise to a worldwide outbreak, which the World Health Organization (WHO) notified in the first quarter of 2020 ^[1]. Almost all of the countries of the world have experienced widespread transmission of this virus. The aforementioned sickness has currently claimed the lives of almost $50$ billion people. COVID-19 has affected more than 500 million people worldwide. We cite ^[2,3] for a detailed explanation of how diseases are contracted and spread. Because the pandemic has affected the entire world, every nation has implemented its own policies. Some nations have implemented lockdowns, while others have issued laws requiring strict safety measures like wearing face masks and avoiding large gatherings, etc. The countries most impacted by the disease are Brazil, Italy, India, Iran, Spain, USA and the UK. Earlier research suggested that the virus originated in animals, but later research showed that it may also transfer from person to person ^[4,5].

Numerous infectious diseases are being investigated for potential treatments, controls, cures, etc. Since epidemiology is an important field of medical science, one of the key subfields from this viewpoint is mathematical biology. In contrast to earlier times, researchers are now becoming increasingly interested in biomathematics and bio-math engineering. In order to give a framework for comprehending the dynamics of diverse infectious diseases, mathematical models have been utilized. The most effective tools for studying different diseases are mathematical models ^[6,7,8]. We can create numerous predictions, regulating strategies, etc ^[9,10,11], for a disease in a community by using mathematical models. Considering these factors, researchers have shown a strong interest in the aforementioned subject, we refer the reader to ^{[12,13,14,15]}. Many mathematical models have recently been developed by researchers for COVID-19 as well; we provide e.g., see ^{[16,17,18,19]}. Major public health difficulties and an economic catastrophe are being faced internationally as a result of the aforementioned disease's recent pandemic status and quick global spread.

In order to theoretically improve existing mathematical models for the aforementioned disease, it is necessary to evaluate the influence of recently developed vaccinations with known high efficacy that have weak transmission-controlling measures. For example, Pearson et al. ^[20] studied COVID-19 immunization in low as well as middle-income countries. they found it to be very cost-effective, according to recent research ^[21], if a vaccination with high efficacy and low cost is available. This will have a big impact on disease eradication. Prior to pharmaceutical interventions like treatment and immunization, non-pharmaceutical precautionary measures like the self-quarantining of confirmed cases, isolation, using of face masks, hand washing, social distancing, lockdowns, avoiding gatherings, and closing schools have been employed. This, however, is not a permanent solution. The best strategy will thus be to develop a suitable immunization and make it broadly accessible in every civilization at a fair price. Investigation of COVID-19, various researchers have developed many models. In order to do this, we have aimed to develop a mathematical model with a class of people who have received vaccinations and qualitatively analyze it while considering the therapeutic choices, immunization of those who are susceptible, and healthcare for infected/hospitalized people. Our model takes into account certain important biological and epidemiological elements of the disease under study, such as the disease's inhibitory effect, mortality rates from infections and natural causes, birth rate, and various vaccination rates.

Here, we build the dynamical model for COVID-19. Additionally, we depict the interactions between several model compartments in flow charts in Figure 1. Consequently, we first explain our model as follows:

In Figure 1, the model is presented as a flowchart. Epidemiological models have been examined using difference equations or classical order derivatives. Some references are cited as ^[22,23]. The aforementioned operators cannot generate the phenomenon's global dynamical behavior because of their local nature ^{[24,25,26,27]}. As a result, academics have now taken the majority of models into account when considering fractional order derivatives ^{[28,29,30,31]}. The growth of the aforementioned field is attributable to its numerous applications to practical issues, particularly in the fields of mathematical biology and dynamical analysis ^{[32,33,34,35,36]}. As far as we are aware, fractional calculus is being used to describe an increasing number of real-world scenarios; we refer the reader to ^{[37,38,39,40,41,42,43]}. Additionally, the idea of fractional calculus has been applied in a number of scientific and technological fields ^[44,45,46]. Because of its global nature and memory-preserving methodology, the theory of fractional calculus has been widely used in the mathematical modeling of numerous diseases for further information, see ^[47,48,49]. The conventional analytical and numerical approaches have been upgraded to handle situations with fractional order derivatives (see references [50,51]). To examine diverse dynamical issues involving fractional order derivatives, the classical decomposition, perturbation, and transformation methods, for instance, have been expanded (some references are cited as ^[52,53]). Traditional numerical techniques, such as RKM tools, Euler and Taylor methods, the Adams-Bashforth method, etc., have been upgraded more and more to cope with fractional order problems (see ^[54,55]). Recently, a novel idea has been identified that describes the crossover behavior of the dynamics, because the states of multiple evolution processes frequently change suddenly. Due to these effects, typical derivatives, whether classical or fractional, cannot be used to appropriately address the relevant multi-phase behaviors. Therefore, scholars ^{[56,57,58,59,60,61]} have suggested this idea, which has the potential to appropriately address the aforementioned behavior to a certain extent. As a result, we also attempt at numerical analysis for our studied mathematical model by looking at multi-phase behaviors under various fractional orders by utilizing some numerical technique based on Newton interpolation polynomials. Discussion and a graphic presentation versus some actual facts are offered.

The provided model (1.1) can be stated as follows in the sense of the new derivative:

where $^{PWC}{\mathrm{D}}_{+0}$ denotes piecewise classical or Caputo derivative with two sub-intervals in $[0, T]$.

For further clarification, we respectively express the left hand sides of (1.2) as follows:

where ${\mathrm{D}}_{+0}$ and $^{C}{\mathrm{D}}_{+0}^\varphi$ are classical and Caputo derivatives, respectively.

Our manuscript is structured as follows: In the first section, we give introduction. In the Section 2, we give basic results. Stability results are given in the Section 3. Numerical results are given in the Section 4. In the Section 5, the fractional order form of the proposed model is investigated. In the Section 6, we given conclusion.

2.   Basic results

We recollect some fundamental results of fractional calculus.

Definition 2.1. ^[43] Let $\varpi > 0$, then, the non-integer order integral of a function ${\bf{S}}:[0, \infty)\to \mathcal{R}$ is defined as

provided that right-hand side exists. Also, derivative in the Caputo sense is defined as

Lemma 2.2. ^[43] Let ${\mathrm{h}}\in L[0, \infty),$ then, the solution of

is given by

Here, we recollect the piecewise definitions of the fractional order derivative and integral.

Definition 2.3. ^[61] If ${\mathrm{h}}(t)$ is a differentiable function, then the definition of the classical, and fractional piecewise derivative is

such that $^{PWC}_{0}{D}_{t}^\varphi$ is the classical derivative for $0 < t\leq t_1$ and fractional derivative for $t_1 < t\leq t_2$.

Definition 2.4. ^[61] Let ${\mathrm{h}}$ be a differentiable function, then, the classical, and fractional order piecewise integration is given as

such that $^{PF}{{\mathrm{I}}}_{+0}\varphi(t)$ denotes the classical integration for $0 < t\leq t_1$, and Reimann-Liouville integration for $t_1 < t\leq t_2$.

Lemma 2.5. ^[61] Consider the problem with the piecewise fractional order derivative

whose solution is given by

Here, we define the Taylor series from ^[59]; consider

where $\mathscr{F}:[0, \infty)\times {\bf{R}}\rightarrow {\bf{R}}$. The generalized Taylor series is thus given by

The outcomes of fractional calculus are used here. $\mathcal{D}_t^{2\varphi}{\bf{x}} = \mathcal{D}_t^{\varphi}\mathscr{F}(t, {\bf{x}}(t))+\mathscr{F}(t, {\bf{x}}(t))\mathcal{D}_{\bf{x}}^\varphi \mathscr{F}(t, {\bf{x}}(t)).$ By using (2.2), after rearranging the terms, we obtain the general formula shown below.

where

Assume that $\varphi = 1$; the standard method described by (4.2).

3.   Stability analysis

We initially examine our model's feasibility in this part. In order to achieve this, we add each equation in the model (1.1) under the condition that

one has

From (3.1), we have

Then, $\lim_{t\to \infty}\sup \mathcal{N} = \mathcal{N}_0$ if and only if $\lim_{t\to \infty}\sup \mathcal{I} = 0$.

In the system (1.1), the first equation yields

Similarly, the second equation of the model (1.1) yields

From above, we conclude that $\frac{d\mathcal{N}}{dt} < 0$ if $\mathcal{N} > \mathcal{N}_0$. Additionally, we have

Disease free equilibrium points: We analyze whether the system (1.1) has an equilibrium point. $\mathcal{E}_0 = (\mathcal{S}^0, 0, 0, 0) = \left(\frac{B}{d_o}, 0, 0, 0\right)$ is a representation of the model (1.1) disease-free equilibrium.

Endemic equilibrium: From (1.1), lets assume that $\frac{d\mathcal{S}}{dt} = 0, \frac{d\mathcal{E}}{dt} = 0, \frac{d\mathcal{I}}{dt} = 0, \frac{d\mathcal{Q}}{dt} = 0$ and $\frac{d\mathcal{V}}{dt} = 0$; then, by simple calculation one can find the endemic equilibrium to be as follows:

3.1.   Expression for $\mathscr{R}_0$

In epidemiology, the idea of the fundamental reproduction number, or $\mathscr{R}_0$, outlines how illnesses spread and are treated. Both the frequency of the illness in the populace and the best countermeasures to protect the neighborhood's residents from the deadly virus are disclosed by $\mathscr{R}_0$. The most recent method for determining $\mathscr{R}_0$ is as follows. If we set $\chi = (\mathcal{E}, \mathcal{I})$, we have from the system (1.1) the following

where

and

The Jacobian of ${\bf{H}}$ for the disease-free equilibrium is

and for the disease-free equilibrium, the Jacobian of ${\bf{F}}$ is given by

Hence

We have

Hence, the required $\mathscr{R}_0$ is given by

If and only if the reproduction number meets the requirement that $\mathscr{R}_0 < 1$, then we may state that the disease-free equilibrium is stable and the population is disease-free. You can halt the spread of an epidemic by exercising prudence. By sitting the numerical values of the variables above, $\mathscr{R}_0 = 0.98347$. This number shows that the spread of COVID-19 in the locality is well managed.

On the basis of $\mathscr{R}_0$, we established the following theorem.

Theorem 3.1. (ⅰ) If $\mathscr{R}_0 \leq1,$ then positive equilibria for system (1.1) are not possible.

(ⅱ) If $\mathscr{R}_0 > 1$, then there is an endemic equilibrium, often known as a distinct positive (unique) equilibrium $\mathcal{E^*} = (\mathcal{S^*}, \mathcal{I^*}, \mathcal{Q^*}, \mathcal{V^*})$.

Some further results on equilibrium point with global and local stability are presented here for system (1.1).

Theorem 3.2. The disease-free equilibrium $\mathcal{E}^0$ is locally asymptotically stable, if $\mathscr{R}_0 \leq 1$, and unstable if $\mathscr{R}_0 > 1.$

Proof. The Jacobian matrix at $\mathbb{E}^0$ is given below:

From a simple straightforward calculation ${\bf{M}}^0$ has three eigenvalues which are negative; $\Lambda_1 = -(v_\mathcal{S}+d_0+q_\mathcal{S})$, $\Lambda_2 = -(d_0+v_\mathcal{Q})$ and $\Lambda_3 = -d_0$. The next eigenvalues of ${\bf{M}}^0$ are derived from

If $\mathscr{R}^0 > 1$, then,

which means that (3.4) has two roots, i.e., one positive root and one negative root. Therefore, the disease-free equilibrium $\mathbb{E}^0$ has an unstable saddle point. When $\mathscr{R}_0 = 1$, to get the result for the disease-free equilibrium $\mathbb{E}^0$ that is globally asymptotically stable, we consider the function called the Lyapunov function.

Taking the derivative of $L(\mathcal{E}, \mathcal{I})$ with respect to time $t$ gives

Suppose that $\mathcal{I} = 0$ implies that $\frac{dL(\mathcal{E}, \mathcal{I})}{dt} = 0$. Furthermore, when $\mathscr{R}_0 = 1$, which implies that $\mathbb{E}^0$ is globally asymptotically stable in $\Phi$. □

In the next theorem, we study the stability at point $\mathbb{E}^*$.

Theorem 3.3. If $\mathscr{R}_0 > 1$, the endemic equilibrium $\mathbb{E}^*$ is globally asymptotically stable.

Proof. For the system (1.1), the Jacobian matrix is

We see that the eigen values of ${\bf{M}}$' i.e., $\Lambda_1 = -d_0$ and $\Lambda_2 = -(d_0+v_\mathcal{Q})$ are negative. Also, we have

where

Consider the following equation

Using the Routh-Hurwitz theorem, we have that $A_1.A_2-A_3 > 0$ by a simple direct calculation which shows that $\mathbb{E}^*$ is locally asymptotically stable. This completes the proof. □

4.   Numerical simulation and discussion

Using the parameters from Table 2, we conducted the numerical simulations for model (1.1), we describe them in this section. Here, we describe how we simulate our system (1.1) by using the second-order Runge-Kutta (RK2) method as implemented in ^[64]. For this, we took into consideration

Here $\mathcal{x} = (\mathcal{S}, \mathcal{E}, \mathcal{I}, \mathcal{Q}, \mathcal{V})$. The suggested approach allows us to express the relevant iterative formula as

Here, ${\bf{K}}_1 = h\mathscr{F}(t_i, {\bf{x}}_i(t_i))$, also, $h = t_{i+1}-t_i$. We simulated our results by using the numerical approach described in (4.2) and use the numbers from Table 2 in order to evaluate with our built-in model (1.1). We use the results of the numerical simulations to demonstrate our theoretical results. The parameter values and their descriptions for the numerical simulations are listed in Table 2. However, some parameter values that are unknown have been assumed for the sake of illustration.

We applied 300 days of the actual data from Pakistan. The overall population of the nation is approximately equal to $222.60$ million, according to the data in ^[65,67]. Nearly $47%$ of the nation's population was fully vaccinated during the past nine months, while $57%$ was partially vaccinated. In order to analyze the dynamics of disease propagation and the impact of the vaccination on its cure, we simulated the outcomes for $300$ days as shown in Figures 2–6, respectively by using the numerical approach described by (4.2) of RK2 type. Here, we describe the model (1.1).

Figures 2 through 6 depict the spread of the disease over a 300-day period at the specified speeds. The exposed class and the susceptible class as shown by the actual data that were employed. The infection spread relatively quickly in less than 50 days. As a result, there were also more persons under quarantine. More people had full or partial vaccinations during this time. We can notice that there were more persons who had received vaccinations.

5.   Some numerical treatment via fractional calculus

Here, we develop a numerical scheme based on RK2 method for the model (1.2) as expressed in terms of (1.3). We have developed a numerical method for the two sub-intervals of $[0, T]$. We express the system (1.2) in terms of 1.3 by using ${\bf{(x)}} = (\mathcal{S}, \mathcal{E}, \mathcal{I}, \mathcal{Q}, \mathcal{V})$ as follows:

We may write the aforesaid system (5.1) by using (4.2) for classical order derivative when $0 < t < t_1$ and (2.3) with the (2.4) for the fractional order Caputo derivative with $t_1 < t < T$. Hence, we can write the first equation of the system (5.1) as follows:

where $h = t_{i+1}-t_i$, and

Similarly, the equations for the other compartments can be written as follows:

and

Here, we plot the approximated results for two sets of fractional orders in two cases.

Case Ⅰ: When $0 < \varphi \leq 0.6$, the corresponding solutions exhibiting crossover behavior were obtained as in Figures 7–11 for different fractional orders.

Case Ⅱ: Now we present the numerical results corresponding to different fractional orders such that $0.6 < \varphi < 1$ in Figures 12–16, respectively which show the crossover behaviors of different compartments.

Here, we discuss the results of modeling of the piecewise derivative shown in Figures 7–16. We observed a the sudden change after 150 days. So, for various fractional orders in $[0, 50)$ and $[50,300]$, we have plotted the numerical findings for systems (1.3) and (5.1). It is possible to see the declines and rises in several compartments with obvious crossover behavior. This novel idea makes a significant contribution to the understanding of how a system that experiences sudden changes in its dynamics behaves. Additionally, the fractional order derivative is a global operator that furthers our comprehension of the phenomenon. In Figure 17, we have compared some real data of infected cases with the simulated results at given fractional orders. Data applied here were taken from Pakistan from January 1, 2021 to January 30, 2021 (see ^[65]). We see that the simulated results were very similar to the real data which confirm the authentication of the fractional order model that we have considered.

Also, we have compared the simulated data for the vaccinated class for the period of 11 days i.e., March 1, 2022 to March 11, 2022 (see ^[67]) in Figure 18.

6.   Conclusions

Here, a dynamical system for the COVID-19 disease has been examined. The traditional method for determining the conditions for both local and global stability, as well as for computing $\mathscr{R}_0$ at equilibrium locations, is to employ classical analysis. Additionally, the classical model scenario was numerically simulated. We then developed a strong algorithm to carry out the numerical analysis of the suggested model under the conditions of the new concept, extending the RK2 approach to piece-wise derivative applications. Using actual data from Pakistan for the first 300 days by splitting into two sub intervals, we have carried out a number of numerical simulations for various fractional orders. We have found that fractional calculus that involves applying in terms of piecewise derivative is more effective at clearly indicating the rapid changes in various model compartments. Crossover behavior is what we call such behavior. In the aforementioned process, a system experiences a rapid change in its state of dynamics and runs its crossover behavior is terminated. The aforementioned idea can be further refined to encompass more complex dynamical systems with various types of fractional differential operators. Also, in the future, we will combine stochastic type and fractals-fractional differential type operators in piecewise form for various infectious disease models.

Use of AI tools declaration

The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.

Acknowledgements

This work was supported and funded by the Deanship of Scientific Research at Imam Mohammad Ibn Saud Islamic University (IMSIU) (Grant Number IMSIU- RG23124).

Conflict of interest

The authors declare that they have no competing interest.