A fractional system of delay differential equation with nonsingular kernels in modeling hand-foot-mouth disease

In this article, we examine a computational model to explore the prevalence of a viral infectious disease, namely hand-foot-mouth disease, which is more common in infants and children. The structure of this model consists of six sub-populations along with two delay parameters. Besides, by taking advantage of the Atangana–Baleanu fractional derivative, the ability of the model to justify different situations for the system has been improved. Discussions about the existence of the solution and its uniqueness are also included in the article. Subsequently, an effective numerical scheme has been employed to obtain several meaningful approximate solutions in various scenarios imposed on the problem. The sensitivity analysis of some existing parameters in the model has also been investigated through several numerical simulations. One of the advantages of the fractional derivative used in the model is the use of the concept of memory in maintaining the substantial properties of the understudied phenomena from the origin of time to the desired time. It seems that the tools used in this model are very powerful and can effectively simulate the expected theoretical conditions in the problem, and can also be recommended in modeling other computational models in infectious diseases.


Introduction
Infectious diseases are a branch of science that deals with the diagnosis and treatment of diseases caused by microorganisms. While many infectious diseases such as tuberculosis, plague, leprosy, smallpox, and the flu have existed throughout the history of the world, humans have been fighting microorganisms for centuries. Although some infectious diseases, such as smallpox, have been eradicated through vaccination, new infectious diseases have emerged, such as AIDS and COVID-19. With globalization, global warming, and increasing travel, whether newly emerging viral diseases or resistant bacterial infections, infectious diseases remain at the forefront.
The hand-foot-mouth disease (HFMD) is one of the most common infectious diseases that affect many people around the world. The disease is caused by some viruses like Coxsackievirus A16 and Enterovirus 71. The disease is more common in children under 5 years old. More precisely, the prevalence of the disease is much more common in children under three years of age. Besides, there have been numerous cases of the disease in adults around the world. In most cases, patients have mild symptoms for 7 to 10 days [1]. Some of these symptoms are fever and flu-like symptoms, mouth sores, and skin rash. So far, two main ways of transmitting this contagious disease have been identified: the first way is transmitted through contact with an infected person or contact with tools used by infected people. The disease can also be spread through coughing or sneezing in the air and transmitted through the respiratory transmission to another person.
The first cases of HFMD cases were clinically reported in Canada and New Zealand in 1957. The disease has been reported in many parts of the world, including some parts of Asia, Europe, and the United States. The disease first appeared in Shanghai, China, and spread rapidly to other parts of the country, including Beijing, Shandong, and Jiangxi provinces. According to our information, the epidemic of this disease occurs in cycles of two or three years. The best time to peak the spread of the disease is usually in the summer. Heat and humidity are known as two main factors that aggravate the spread of the disease. The first case of the disease was seen in Thailand in 2003. In a short time, signs of this contagious disease were visible in all cities and provinces of this country [32]. The first national study describing a large number of deaths caused in this country by the disease during an outbreak in 2011 was presented in [46]. HFMD has also been found in many Indian major states [38].
Unfortunately, no specific curative treatment has been found for this viral disease. However, several simple preventive measures such as avoiding direct contact with infected people, quarantining infected children at school, cleaning common utensils regularly, and disinfecting contaminated surfaces have been confirmed as the most effective ways to reduce transmission of the virus to other people in the community.
Due to the high prevalence of this disease in different parts of the world and the importance of identifying and controlling the factors affecting its prevalence, many extensive kinds of research have been conducted by researchers from various scientific aspects. In [26], a new method for the disease prediction using GeoDetector and a Long Short-Term Memory neural network (LSTM) has been proposed. In [40], the authors applied optimal control theory to the HFMD model, including the treatment and vaccination interventions. They presented some control strategies based on minimizing the cost of the intervention and minimizing the number of infected people. Very recently, a fractionalorder model has been utilized to describe the transmission of HFMD in [39]. In this paper, the authors considered two major cases of constant and optimal control. In each case, the existence and uniqueness of positive solutions and the sufficient conditions for the existence and stability of equilibriums were investigated. It is important to mention that their study on the optimality control conditions is based on Pontryagin's maximum principle. Another interesting study on the disease was conducted to investigate the role of air pollution in the prevalence and incidence of the disease in the warm and cold seasons in Wuhan, China [47]. In [42], the authors investigated some strategies to control the disease using a system of ordinary differential with two delay parameters. To estimate and better understand the transmissibility of HFMD, a susceptible-infected-recovered (SIR) model has been utilized in [26]. They claim that the main reason for this choice is that the incubation period of the disease is less than one week. In [13], another new SIR model has been developed to fit the surveillance data containing valuable information on the severity of HFMD in order to accurately estimate the basic reproductive number (R 0 ) of the disease. In [27], a simple SEIR model has been examined to investigate the dynamics of the disease among young children. A discrete SEIADR (susceptible-exposed-infectiousasymptomatic-dead-recovered) epidemic model has been also developed in [34]. To read more articles, please refer to [12,16,17,28,31,33,45,50,51].
Fractional differential calculus has been used in modeling many phenomena in everyday real-world applications [18,[20][21][22][23]25]. Therefore, due to the importance and scope of applications, many efficient numerical methods specific to each of these types of fractional operators have also been introduced [3-8, 10, 11, 15, 19, 24, 30, 35, 37, 43]. However, in [41] the author has conducted a review to point some possible problems and difficulties arising in the construction of fractional-dynamic model analogs of standard models by using fractional calculus. This stems from the fact that some fundamental properties of fractional derivatives (such as the multiplicity principle, the solvability and correspondence principles, and the interpretability principle) violate various known standard rules and properties that are fulfilled for derivatives of integer order.
Information retention (or so-called memory effect) is one of the most basic properties of fractional-order derivatives. This significant property has made these operators one of the most efficient tools in computational models arising in biology. In other words, the basic information of the model from the beginning of time is utilized to characterize the behavior of the phenomenon at any desired time. The use of fractional derivatives has made it possible to employ such a valuable feature. To take advantage of memory-related benefits in the context of fractional-order derivatives, our main objective is to study a nonlinear system of the delayed-fractional model in studying HFMD. For this purpose, the present contribution is structured as follows: Some mathematical background, mainly about recent definitions on fractional calculus model, is formulated in Sect. 2. The proposed model is introduced in Sect. 3. In Sect. 4, we investigate some mathematical frameworks of the model, including the equilibrium points, the basic reproduction number, and the existence and uniqueness of the model's solution. Then, we present some corresponding numerical simulations in Sect. 5. Finally, some conclusions are drawn in the last section of the article.

Some basic preliminaries on fractional operators
Employing new definitions in mathematics always makes extensive progress in modeling various phenomena around us in the world. One of these new and powerful topics is the use of fractional calculus concepts. In this section, we will have an overview of some of known basic definitions corresponding to fractional calculus that are used in what follows.
For this purpose, we present the definition of Riemann-Liouville fractional integral and derivative operators RL I α and RL D α of order α > 0, respectively. Then, we define the Caputo fractional derivative C D α of order α > 0 as a modification of the Riemann-Liouville fractional derivative. Further, we present some properties of the Mittag-Leffler function E α,β (z), α, β > 0 (see [48]).

Definition 1 For a given integrable function S(t), the fractional integral operator in the
Riemann-Liouville sense C D -α of order α > 0 is given by [36] where (·) denotes the well-known gamma function.
Definition 2 For a given function S(t) in C[0, T], the fractional derivative operator in the Riemann-Liouville sense C I -α of order α > 0 is given by The following two important propositions are consequences of the above definitions: Definitions (2) and (3) are different from each other, and the relation between the two types of fractional derivatives is as follows: The Caputo derivative has the main advantage that the initial condition of the corresponding problem has the same value as the ordinary differential equation. Moreover, for a constant-valued function, the Caputo derivative is zero.

Definition 4 The Mittag-Leffler function of two parameters is given by
where α, β > 0, C denotes the complex plane.

Definition 6
The AB-Caputo integral operator of order α is defined as (6) and (5), the following union is established [2]:

The delayed-fractional version of the model
In recent years, modeling real-world problems using the fractional delay differential equations (FDDEs) has attracted much attention among mathematicians, physicists, and engineers. Due to the wide application of these equations, the theoretical and practical aspects of this category of equations have been extensively studied by researchers in many research articles such as [42,44,49].
In the light of these facts and after employing the AB-Caputo fractional derivative AB D α and defining two-time delays τ 1 and τ 2 in the model presented in [42], we arrive at the following FDDE for the spread of the disease: where λ = β(I c (tτ 2 ) + μV ca (tτ 2 ))/N (t), and subject to given initial conditions The effective populations in this model are categorized into six state variables as follows: the susceptible subpopulation S(t), the treated and vaccinated subpopulation T ν (t), the clinically infectious subpopulation I c (t), the recovered subpopulation R(t), the vaccinated subpopulation V ν (t), and finally, the vaccinated carrier subpopulation V ca (t). We symbolize the sum of all subpopulations by N (t). In this model, the transmission rate is denoted by β. Moreover, b is the birth rate, η 2 represents the transforming rate from I c (t) to R(t). Also, η 2 is utilized to describe the vaccinated carrier rate, ω is the recovery rate, φ is the rate of protection loss, τ 1 and τ 2 are two delay parameters. Moreover, μ denotes the natural death rate, treating and vaccinating rate of susceptible people is shown by ρ 1 , culling of clinical infective, and vaccinated carrier rate is denoted by δ. Finally, the rates corresponding to the treating vaccinated, vaccinating susceptible, and recovery of treated people are , ρ 2 , and η 4 , respectively [29].
In consecutive subsections, we analyze some theoretical aspects of the fractional model outlined in (8).
-The endemic equilibrium of the model is , , , A detailed survey of the stability of these equilibrium points can be found in the reference [42].

The basic reproduction number
Employing the next generation matrix method [14], the basic reproduction number of model (8) is given by [42]

Existence of the solution
To investigate and prove the existence of a solution for the model, we can first apply the AB-Caputo integral operator (6) to the sides of model (8) to give the following equalities: Defining N(F(t)) = [N 1 (F(t)), N 2 (F(t)), . . . , N 6 (F(t))], and F 0 = [S(0), T ν (0), I c (0), R(0), V ν (0), V ca (0)], Equation (13) can be considered as follows: The following iterative process can be defined: Subtracting two consecutive terms gives For simplicity, we set n (t) = F n (t) -F n-1 (t). Hence, one has Thus, we get n (t) = F n (t) -F n-1 (t) , Whenever N satisfies the Lipschitz condition with respect to F, thus Hence the following inequality will be obtained: Now, replacing n-1 (t) by its value, we get Also, we have And finally, By choosing we are now in a position to define the following sequence for F(t): where n (t) → 0 when n → ∞. Thus, Now, we can write Taking the standard norm on both sides of the above equation, we conclude that If n → ∞, the right-hand side of the equation tends to zero. So, one gets And that was what we were trying to prove in this subsection.

The uniqueness of the solution
In this part, we are looking for the proof of the uniqueness of the solution related to the model. To do this, let us consider that model admits two solutions F(t) and N(t). Then we can write . . .

A numerical method
In recent years, various approximate methods have been used to solve the system of fractional-order differential equations. In each of these methods, a specific idea is used to discretize the problem and approximate the solution of the system.
In order to express the method, we first need to consider a fractional delay differential equation given by To determine the numerical method, we first consider the domain discretization as follows: Keeping these notations in mind, we obviously have R(t j ) = ψ(t j ), j = -k, -k + 1, . . . , -1, 0. In other words, it can be written The main assumption of the numerical method is that we have approximated the value as the solution function in points R(t j ), j = -k, -k + 1, . . . -1, 0, 1, . . . , n, and now we are looking for the value of R(t n+1 ). The next idea used is to apply the definition of the integral operator (6), which results in the following relation: Setting t = t n+1 in (29) gives Now, we utilize the product trapezoidal quadrature formula to approximate (30). Then, the corrector approximation for R(t n+1 ) is obtained as follows: where k 1 = ceil(τ 1 / ), k 2 = ceil(τ 2 / ), and The unknown term R P (t n+1 ) that appears in (31) is also calculated using the product rectangle rule as follows: where

Simulation results and discussion
In this section, the numerical simulations corresponding to the proposed model along with the attached approximate method are examined. In all graphs, the following values are used as the default of the constants, unless we intentionally change some of these values in the graphs ourselves. These changes are expressed in each graph.
The effect of fractional-order parameter (α) on the disease evolution, in the absence of delay parameters τ 1 = τ 2 = 0, is shown in Fig. 1. It is seen that when α increases, then the susceptible and clinically infected classes decrease, and treated and vaccinated, recovered, vaccinated, and vaccinated carrier classes increase. In this case, all model solutions converge to the internal equilibrium point of the system. The effect of the treating of vaccinated people rate ( ) on the disease prevalence, in the absence of delay parameters τ 1 = τ 2 = 0, is presented in Fig. 2. In some diagrams of Figure 1 The impact of α on the results for τ 1 = τ 2 = 0 Figure 2 The impact of on the results for τ 1 = τ 2 = 0 and α = 0.95 this figure, it can be seen that a change in the value of the parameter can cause very significant changes in the overall condition of the system. In a way, by taking some values, the behavior of some solutions will be completely different from others in general. For example, by taking = 0.8, the graphs related to the two solutions clinically infected classes and vaccinated carrier classes take on a damping state and tend to zero.
Further, the effect of the rate of vaccination (ρ) on the disease evolution, in the absence of delay parameters τ 1 = τ 2 = 0, is displayed in Fig. 3. In this figure it is clear that by taking ρ = 0.6, 0.7, and 0.8 the graphs related to the two solutions clinically infected classes and vaccinated carrier classes take on a damping state and tend to zero. In other words, in these cases, the system moves to the disease-free equilibrium point. These conditions are fully consistent with the theoretical arguments presented in this paper and reference [42].
The effect of the rate of vaccination (b) on the disease dynamics, in the absence of delay parameters τ 1 = τ 2 = 0, has been plotted in Fig. 4. The diagrams in this figure clearly show that the curves corresponding to b = 0.1 always create the lowest values in all possible To investigate the significant effect of the rate of culling of clinical infective and vaccinated carrier people (δ) on the disease dynamics, in the absence of delay parameters τ 1 = τ 2 = 0, we have presented Fig. 5. A distinctive feature of the diagrams in this figure is that the final values for system solutions always have an alignment correlation to the parameter δ. In other words, by increasing the value of parameter δ, the equilibrium points of the system are placed in higher levels. Figure 6 exhibits the effect of the rate of delay parameter (τ 2 ) on the disease dynamics of the model when we take τ 1 = 3. As we have seen in this figure, by increasing the delay parameter τ 2 , the solutions S(t), T ν , V ν (t) have similar trends. The same is true also for I c (t), R(t), V ca (t). Also, by increasing the value of the parameter τ 2 , the system shows more oscillating behavior. The impact of ρ 2 on the results for τ 1 = 10, τ 2 = 100, and α = 0.95 Figure 7 demonstrates the effect of the rate of delay parameter (τ 1 ) on the disease variation of the model when we take τ 2 = 15. The occurrence of oscillating properties in this figure is very obvious. In addition, by increasing the value of the delay parameter τ 1 , the system reveals more severe instability behavior. Figure 8 demonstrates the effect of the rate of vaccinating susceptible people (ρ 2 ) on the disease dynamics for delay parameters τ 1 = 10, τ 2 = 100. In this figure, the high sensitivity of the model's response behavior to the parameter ρ 2 can be seen. For example, by considering ρ 2 = 0.1, the behavior of the system is quite volatile and unstable. While for the bigger values of this parameter, the system's response show more stable behavior.
The effect of the rate of culling of clinical infective and vaccinated carrier people (δ) on the disease dynamics, for delay parameters τ 1 = 5, τ 2 = 10, has been displayed in Fig. 9.
The influence of the transmission rate (β) on the disease dynamics, for delay parameters τ 1 = 3, τ 2 = 10, is investigated in Fig. 10. It can be seen that for the β = 0.9, the system tends toward the disease-free equilibrium B 1 , and for the other considered values to the endemic Figure 9 The impact of δ on the results for τ 1 = 5, τ 2 = 10, and α = 0.95

Figure 10
The impact of β on the results for τ 1 = 3, τ 2 = 10, and α = 0.95 equilibrium point of B 2 . Therefore, β can play a very important and constructive role in the behavior and dynamics of the disease.

Conclusion
Humans have been exposed to infectious diseases throughout history and have always tried to be able to control the spread of the disease and then find an effective treatment process to treat the disease. One way to control infectious diseases is to study the dynamic systems that describe the extent and prevalence of the disease. On the other hand, applying the modern concepts presented in the study of these dynamic systems can lead to dramatic advances in the study and control of this disease. In this article, we examined the prevalence of hand-foot-mouth disease in a certain population that has been described by a fractional system of ordinary differential equations. The model is constructed using the Atangana-Baleanu fractional derivative and two constant parameters to apply the time delay in the solutions. In order to take advantage of the concept of memory in the evolu-tion of the model, we have employed a well-known derivative with fractional order as well as two delay parameters in the model. Several numerical simulations revealed the effects of the fractional parameter, and the existence of two delay parameters in the disease dynamic. Through some performed experiments, we also investigated the sensitivity analysis of the model to some numerical parameters of the model. In some cases, it was observed that with a slight change in some of these parameters, very fundamental changes in the behavior of system responses occur. Another noteworthy aspect of this paper is the high degree of influence of the model on the delay parameters. For some specific values for the parameters in the system, unstable chaotic behavior can be observed in the system responses. The results of the model presented in this paper provide the ability to investigate the effects of time delay parameters, prophylactic vaccination, reactive vaccination, prophylactic treatment, and reactionreflecting parameters on disease outbreak in a population. The process presented in this article can be applied to other infectious disease models.