On the fractional SIRD mathematical model and control for the transmission of COVID-19: The first and the second waves of the disease in Iran and Japan.

In this paper, a fractional-order SIRD mathematical model is presented with Caputo derivative for the transmission of COVID-19 between humans. We calculate the steady-states of the system and discuss their stability. We also discuss the existence and uniqueness of a non-negative solution for the system under study. Additionally, we obtain an approximate response by implementing the fractional Euler method. Next, we investigate the first and the second waves of the disease in Iran and Japan; then we give a prediction concerning the second wave of the disease. We display the numerical simulations for different derivative orders in order to evaluate the efficacy of the fractional concept on the system behaviors. We also calculate the optimal control of the system and display its numerical simulations.

Introduction

Coronaviruses, viruses that have a crown, have caused a variety of severe illnesses like Severe Acute Respiratory Syndrome (SARS) and Middle East Respiratory Syndrome (MERS). The disease COVID-19, which was caused by the virus SARS COV-2, was first identified in Wuhan City of China on December 2019. The COVID-19 disease appears in different people with different symptoms. In most people, fever, cough, and the shortness of breath have been the early symptoms; in some others, gastrointestinal symptoms such as anorexia, diarrhoea, and vomiting have been appeared. Additionally, other symptoms such as body aches, the loss of smell and taste senses, and some skin symptoms have been reported. The incidence of COVID-19 in children is lower than adults, and for people with some diseases such as diabetes, high blood pressure, asthma, and autoimmune diseases, it can lead to more serious conditions and sometimes death. Because of the new outbreak of COVID-19, there is no complete information on how it is transmitted, but according to the World Health Organization (WHO), coronavirus is spread through the sneezing and coughing of infected people, and it stays in the air for several hours and enters the body by breathing the polluted air. Also, it sits on the surface and is transmitted to a body by touching the infected environment and then the nose, eyes, or mouth. Due to these transmission kinds of this virus, the suggested ways of preventing the transmission include washing hands with soap, using masks outside the home, not presenting in the polluted environments, and reducing unnecessary presence in the community.

Many mathematicians and biologists have studied the dynamics of realistic systems by using the concept of mathematical modeling (see, for example, [1], [2]). By extending the integer-order derivatives to the fractional-order ones, many researchers have studied the mathematical models of natural phenomena with fractional-order models [3], [4]. Recently, the utilization of fractional operators has been expanded to model real-world systems, an efficient approach which investigates the transmission of diseases as well as their control [5], [6]. With the expansion of COVID-19, researchers have also used mathematical models to simulate the spread of this disease in different environments (see, for example, [7], [8], [9], [10]).

In this research, we study the transmission of COVID-19 by introducing a fractional SIRD model in which the Caputo derivative is used for the system formulation. With the onset of this disease, most countries made a series of measures to control its transmission, which led to a decline in the spread of the disease after the first peak. Iran was the first country that experienced the second peak of COVID-19. We provide a comparative study for the first and the second waves in Iran and the prediction of the disease transmission in the future. Japan was another country that experienced the second wave, while the mortality rate with COVID was low in that country. Thus, we examine the second wave of the disease in Japan and provide a forecast for it.

In the following, Section 2 provides some basic definitions with regard to the fractional calculus. The system of equations for the COVID-19 spread as well as the existence and uniqueness of the solution are presented in Section 3. The equilibrium points and the optimal control of the system are presented in Sections 4, 5 , respectively. Numerical results and the simulations of the first and the second waves as well as the prediction of the second wave in Iran and Japan are discussed in Section 6. Finally, the paper is closed by some conclusions in Section 7.

Fractional calculus

In this section, we first recall some definitions of fractional operators; then we examine the important feature of fractional calculus called memory effect.

Basic definitions

[11]

For , , and , the Caputo fractional derivative is described by in which in the fractional order. Also, for with , the fractional integral of order is defined as

For the fractional-order Caputo operator of order , the Laplace transform is which can be written more simply as

[12]

Let be a function such that and are continuous for . Then for all , there exists some satisfying the condition

Memory effect

Fractional calculus is one of the most powerful tools in the modeling of real-world phenomena with memory effect. The effect of memory in fractional derivatives has been studied by many researchers (see, for example, [13], [14], [15]). In the mathematical models with classical differential operators, the responses do not depend on the previous instance, so these models do not have memory. However, the fractional mathematical models include memory as their inherent property. Consider an integrable function in and time instances such that . If , we have If , then the first integral is canceled, and we obtain The above equations show that in the classical case , the value of depends only on , but in the fractional case , the value of depends on and the previous range . Thus, the fractional integral in depends on what happened before ; this shows the effect of memory on the fractional process, while the integral in the classical case has nothing with the range before .

The nature of coronavirus transmission shows that the patient has no specific symptoms in the first few days of infection, but he/she is the carrier of the disease. Therefore, the population of infected people now affect the number infected individuals in the next few days. Hence, the population of infected people at time is affected by the number of infected persons in the previous days, and the chain continues so that the number of infected people at time is influenced by the number of infected people from the time to , indicating that the coronavirus transmission is a system with memory; therefore, it is more appropriate to model this phenomenon by fractional-order derivatives.

Model formulation

In this section, we provide a Caputo-type fractional SIRD model for the transmission of COVID-19. Generally, the SIRD models consider the population to be constant, but here we take into account the natural birth and death rates. We also consider that the total includes four sub-categories: susceptible persons , infected people , recovered individuals , and died people with this disease so that . The diagram of the model is shown in Fig. 1. Also, the following compartmental model describes the spread of COVID-19 where is a nonnegative initial state vector. The parameters of the model are also the per capita birth rate , the recovery rate , the transmission rate , the virus-induced average fatality rate , and the per capita natural death rate .

Fig. 1

Diagram of the COVID-19 transmission.

Fractional-order systems are related to nonlocal effects, history, and memory, the features which exist in different biological systems. The integer-order system (1) suffers from the lack of internal memory effects regarding the COVID-19 transmission, so we are going to modify the aforesaid ordinary model to a fractional-order one. To this end, we replace the ordinary derivatives with the Caputo-type fractional operators. Also, an auxiliary parameter is used to solve the dimensional mismatching problem [16]. Thus, the fractional-order model for the spread of COVID-19 is given as follows where and .

Feasibility region

Suppose that ; then we state that is positively invariant by the next lemma.

With respect to the system (2) , the set is positively invariant.

To prove this lemma, first we show that the solutions are non-negative by using the method presented in [17]. We have which are all nonnegative. Consider the auxiliary fractional differential system where , and is the initial state vector. We will show that for all , the solutions of the auxiliary system are nonnegative. If not, there exists a time instance at which the solutions are negative. Consider Then , and one of the components is zero. Let . Since and is continuous, it implies that for some . Using Theorem 2.1, we obtain , which shows that is nonnegative. Similarly, it can be shown that other components are also nonnegative, which is a contradiction. Therefore, we obtain that for all as ; thus, the solutions are nonnegative. In the following, we first add the first three equations of the system (2). So Applying the Laplace transform results in where are the initial size of susceptible, infected and recovered population, respectively. With some calculations, we get Thus, if , then for . If , since for all , then . Consequently, with respect to the system (2), the set is positively invariant. □

Existence and uniqueness

To prove that the existence of a unique solution for the system (2), we first consider (2) in the form The integral operator is applied to the both sides, which results In the next theorem, we prove the Lipschitz and the contraction conditions for the kernels –.

The Lipschitz condition and the contraction are established for if the following relation holds

According to the values of at and , one can write By considering , where is a bounded function, we obtain Therefore, the Lipschitz condition holds for , and if , then is a contraction. □

Similarly, we can show that satisfy the Lipschitz condition, i.e. where , , , and are bounded. The kernels are also contraction if we have for . It should be noted that if the system (2) has an equilibrium point at the origin, then one can use the method proposed in [18] to prove the existence and uniqueness analysis.

Now suppose the following recursive formulas with respect to the system (2) with the initial state vector . We take the norm of as By the Lipschitz condition (4), we get Similarly, we obtain Then we can obtain In the next theorem, we prove the existence of the solution for the system (2).

The fractional-order transmission model of COVID-19 (2) has a solution if there exists a time instant such that

By considering a recursive technique and using Eqs. (5), (7), we have This shows that the system has a continuous solution. We claim that a solution of the system (2) is constructed by the above functions. For this purpose, consider So We repeat the method; then At , we get Then . Similarly, we have . □

In order to prove that the solution of the system (2) is unique, assume that another solutions , , , and exist for the system (2). We can write Taking the norm, we get By the Lipschitz condition (4), we obtain Then

If the following inequality holds, then the solution of the spread of COVID-19 model (2) is unique

From the above condition, we obtain With this inequality and (7), we conclude that , so we obtain . Following the same manner, we can derive , , and . □

Reproduction number and equilibrium points

The first equilibrium point showing a disease-free state is the point in which , so it is obtained by simplifying the equations as where is the total population. To find the other steady-states, since does not appear in the first three equations of the system under study, we remove the fourth equation and form the following algebraic system We solve the above equations and obtain the equilibrium point such that Here we use the next generation method to find the basic reproduction number  [19]. First, we consider the first three equations of the system (2) in a compact form as follows where and The values of the Jacobin matrices of and at the equilibrium point are Then the basic reproduction number is computed from where is the next generation matrix. By performing simple algebraic calculations, the value of is obtained as follows The reproduction number is epidemiologically important and indicates the ability to transmit the infection and the continuation of the disease.

Stability analysis of equilibrium points

The Jacobian matrix of the fractional-order model (2) is computed by The Jacobian matrix at the equilibrium point is also obtained as In the following, we first recall the definition of marginal stability.

A system is marginally stable if all the poles on the imaginary axis are distinct, and all the remaining poles have negative real parts.

We state the necessary conditions for marginal stability by the following theorem.

The condition is sufficient for the marginal stability of the equilibrium point .

Following the relation , we obtain the characteristic equation at the point whose roots (eigenvalues) are , , and . If , then So . Since , , and , then the system (2) is marginally stable at . □

Sensitivity analysis of

To investigate the sensitivity of , we use the method used by Zafar et al. [3]. Since all parameters of the model (2) are positive, we have As can be seen, is sensitive to changes in the model parameters such that increases with and decreases with .

Optimal control

Currently, there is no vaccine to control the transmission of COVID-19, so the spread of the disease should be controlled by following health instructions and creating social distance. Meanwhile, creating social distance plays an important role in controlling the spread of COVID-19. We denote the social distance by where represents the social gathering. If we have social distance, then it is , and if there is not any social distance, then it will be . The aim here is to decrease the number of infected people by creating social distance while reducing the cost of implementing this control strategy. For this end, a fractional control system is considered as below If we remove the social distance, then it becomes ; thus, the SIRD model (2) is recovered from the control system (8). In addition, represents the control space for .

For fractional-order systems, the optimal control theory is based on the Pontryagin’s maximum principle [20]. Thus, we consider the performance index for the control system (8) as follows where are weighting coefficients corresponding to the number of infected individuals and the cost of control, respectively. We minimize the cost functional by the following theorem.

Let be a measurable control function for all ; then the optimal control minimizing the functional subject to the system (8) is obtained as follows

We consider the Hamiltonian function as follows where , are the co-state variables with , and they satisfy which is a fractional system of right Caputo derivative equations. By the Pontryagin principle, we obtain the following optimality condition Then according to the sign of , we consider the optimal control as follows:

If , then .

If , then .

If , then .

Therefore, by solving the boundary value problem (8), (10) using the presented method in [21], the optimal control is derived. □

Simulation and numerical findings

Hereinafter, we employ the fractional Euler method to obtain the approximate solution of the transmission model (2) [22].

Numerical method

Let the model (2) be considered in the compact form where , is a vector function which is continuous and satisfies the Lipschitz condition, and is the initial condition. We apply the fractional integral operator to Eq. (11); hence, we derive Let and where and is a natural number. Denote the approximation of at by ; then the fractional Euler method [22] provides where the coefficients are computed by the following equation Note that Theorem 3.1 in [22] proved the stability of the presented method. In addition, the solutions of the model (2) are calculated from where , , , and .

Numerical simulation

Case I: Iran

Iran was the first country that experienced the second wave among the others getting involved in the COVID-19 epidemic. To simulate the first and the second waves in Iran, we used the reported data for the infected cases from the world meter website [23]. The model (2) describes an endemic model (with vital dynamics) such that the births are added at the rate into the class and deaths reduce the classes at the rate . Life expectancy in Iran is almost 73 years, so years and (the balance of births and natural deaths). Also, the average infectious period is days, and the mortality rate with this disease is . Additionally, the initial conditions for the first peak (started from February 22th) are , , , and . The auxiliary parameter is also considered as for all cases in the sequel. According to the reported cases by WHO for the first wave of infected people in Iran, we used the fmincon package of minimization and obtained for the model (2) with . Since we used the fractional-order model for the fitting purpose, the obtained cannot be the correct value for the integer-order case. Hence, we fitted the integer-order model separately and obtained . In Fig. 2, we plotted the results for the fractional transmission with and the integer transmission () in addition to the available real data. Also, Table 1 reported the errors related to each model. The results indicate that the fractional-order model has better performance and less errors. To show the influence of the fractional order, we plotted the state variables for the optimal order and the other orders in Fig. 3, Fig. 4. The figures indicate that the derivative order does not influence the behavior of the functions, although the obtained plots are different for various orders.

Fig. 2

The plots of fractional-order model () and integer-order model () for the infected people in Iran (the first wave).

Table 1

Absolute error (AE) and relative error (RE) for in Iran (the first wave).

Model	υ	AE	RE
Integer	1	10.9786	0.043051
Fractional	0.994	10.2456	0.042314

Fig. 3

The dynamics of susceptible and infected groups in Iran for the optimal order and the other orders .

Fig. 4

The dynamics of recovered and died groups in Iran for the optimal order and the other orders .

With the outbreak of COVID-19 in Iran, the government enforced quarantine rules until April 18th, so the number of infected people declined rapidly as shown in Fig. 2. This figure also depicts that with the continuation of this practice, the disease could be controlled in days. However, the quarantine rules were lifted on April 18th; thus, by reopening of jobs, the second-wave of COVID-19 began on May 5th. In the second wave, according to the published information, the mortality rate was about , and with the fitting technique, the amounts of and for the fractional- and the integer-order models were obtained, respectively. Fig. 5 shows the fractional-order results with and the classic integer responses with along with the reported cases of infected people from May 5th; Table 2 also indicates the errors for both cases. The results indicate that the errors of the fractional-order responses are lower, a fact which confirms its superiority. Fig. 6, Fig. 7 portray a forecasting, according to the second wave, for the spread of COVID-19 in Iran. According to these figures, it takes about eight months to pass the second wave, and about 70,000 people will die from this disease. To check the fractional order effect on the dynamics the model (2), the state variables were also plotted for the optimal order and the other orders in Fig. 6, Fig. 7. The plots display that the derivative order does not influence the general functions’ behavior, but the achieved values are different for various fractional orders, and the differences increase over the time.

Fig. 5

The plots of fractional-order model () and integer-order model () for the infected people in Iran (the second wave).

Table 2

Absolute error (AE) and relative error (RE) for in Iran (the second wave).

Model	υ	AE	RE
Integer	1	7.3249	0.03234
Fractional	0.9952	7.03468	0.03106

Fig. 6

The prediction of susceptible and infected groups in Iran for the optimal order and the other orders .

Fig. 7

The prediction of recovered and died people in Iran for the optimal order and the other orders .

Case II: Japan

We consider the second wave of COVID-19 transmission in Japan from June 30th. According to the available data, the mortality rate with respect to the COVID-19 transmission for the second wave in Japan was very low and about . Life expectancy in Japan is almost 85 years, so years. Also, the other parameter is the average infectious period days, and the initial conditions are , , , and . According to the reported cases by WHO for the infected people in the second wave in Japan, we used the fmincon package of minimization and obtained for the fractional-order model with and for the integer-order model with . In Fig. 8, we plotted the results for the fractional and the integer derivatives along with the reported cases from June 30th to July 20th. As can be seen, both the integer- and fractional-order models with and , respectively, follow the real data very well, but Table 3 shows that the fractional-order model has less errors. Fig. 8 also forecasts the second wave of COVID-19 transmission in Japan, which will take about months to subside. Fig. 9, Fig. 10 show the results of the model in Japan, which indicates that the treatment process is well progressed, and eventually, by the end of the second wave, the number of deaths will have reached to nearly 1150 persons. Also, Fig. 9, Fig. 10 display the results for different orders of fractional derivative; as can be seen, a small change in the derivative order makes a big difference in the resulting values, while the general functions’ behaviors are the same.

Fig. 8

The prediction of infected people in Japan by fractional-order model () and integer-order model () (the second wave).

Table 3

Absolute error (AE) and relative error (RE) for in Japan (the second wave).

Model	υ	AE	RE
Integer	1	4.9104	0.02201
Fractional	0.9961	4.5753	0.02135

Fig. 9

The plots of susceptible and infected people in Japan for the optimal order and the other orders .

Fig. 10

The plots of recovered and died people in Japan for the optimal order and the other orders .

Optimal control

Here we use some numerical simulations to investigate the effect of the suggested control strategy, social distance, on the outbreak of COVID-19. To investigate the effect of social distance on the second wave in Iran and Japan, we used the parameters presented in the previous section. Also, we consider , the weighting constants , and the final intervention times and days for Iran and Japan, respectively. To find the optimal control, we combined our suggested numerical method (Section 6.1) and an iterative algorithm proposed in the study [24]. Indeed, we implemented the proposed method in Section 6.1 forward and backward in time for the state and co-state equations, respectively. Note that the convergence and the stability of this technique were developed in [24]. In addition, the utilized scheme was employed in the previous literature for some practical cases such as the solution of general fractional optimal control problems [25], and the optimal control of a human respiratory syncytial virus infection [26]. Fig. 11 shows the population of susceptible and infected people in Iran in the presence and the absence of the aforesaid control strategy. It is apparent that the social distance reduces the number of infected people. Fig. 12 plots the variation in the number of susceptible and infected people in the presence and the absence of the control strategy in Japan. Our goal was to reduce the number of infected people; the results confirm that the number of infected people (peak of the disease) decreased, and since the initial number of infected people was small, this wave was ended faster, and the spread of the disease was controlled by the social distance strategy.

Fig. 11

The plots of susceptible and infected people in Iran in the presence and the absence of control.

Fig. 12

The plots susceptible and infected people in Japan in the presence and the absence of control.

Conclusion

This research presented a fractional SIRD model in the Caputo sense for the spread of COVID-19 between humans. The equilibrium points, the basic reproduction number, and the feasibility region of the system (2) were determined. Also, we proved that a nonnegative solution exists for the model and checked the stability of the disease-free steady-state. An optimal control for the system under study was calculated considering the social distancing strategy. In the numerical part, the fractional Euler method provided the approximate solutions of the system (2). Some numerical simulations were presented to investigate the first and the second waves of the disease in Iran and Japan. With the application of quarantine laws and the closure of jobs in Iran until April 18th, the first wave came to an end by May 2th. However, by reopening the jobs on April 18th, the value of was increased from to , and the second wave was started after two weeks. The second wave was continued for months, and about 70,000 people were died at the end of the second wave. Also, comparing the results between the integer-order model, the fractional-order case, and the reported data showed that the fractional-order one follows the real data better. By plotting the results for different fractional orders, the effect of the fractional order on the results was investigated, which differed only in the resulting values and did not have any special effect on the general functions’ behavior. Next, the second wave of the disease was simulated in Japan, but the amount of produced was less than that in Iran, so we predicted that the second wave in Japan would take about months. The number of infected cases was increased to 12,000, and the number of deaths from this disease was increased to 1200, which was much less than that in the second wave in Iran. This can be due to two reasons: first, because the value of was low, the number of infected cases was slowly increased giving more time to the treatment department; second, the treatment system in Japan might be more successful. Additionally, some numerical simulations for the optimal control strategy in Iran and Japan were presented. The results showed that by creating social distance and observing basic health issues, the number of patients in Iran were decreased; also, if the desired measures continued, the number of patients would reach zero. In Japan, where the initial number of infected people was lower, social distance reduced both the peak and the duration of the disease, so the prevalence of the illness was controlled.

The fractional differential system preserves the historical memory of the system, so in the second wave we should use the data from the beginning of the outbreak, but is different in the two waves. Therefore, with the proposed model, we could not use the primary data for the second wave. Therefore, young researchers can examine how and by which method they can overcome this weakness.

Availability of data and materials

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.