key: cord-325862-rohhvq4h authors: Zhang, Yong; Yu, Xiangnan; Sun, HongGuang; Tick, Geoffrey R.; Wei, Wei; Jin, Bin title: Applicability of time fractional derivative models for simulating the dynamics and mitigation scenarios of COVID-19 date: 2020-06-04 journal: Chaos Solitons Fractals DOI: 10.1016/j.chaos.2020.109959 sha: doc_id: 325862 cord_uid: rohhvq4h Fractional calculus provides a promising tool for modeling fractional dynamics in computational biology, and this study tests the applicability of fractional-derivative equations (FDEs) for modeling the dynamics and mitigation scenarios of the novel coronavirus for the first time. The coronavirus disease 2019 (COVID-19) pandemic radically impacts our lives, while the evolution dynamics of COVID-19 remain obscure. A time-dependent Susceptible, Exposed, Infectious, and Recovered (SEIR) model was proposed and applied to fit and then predict the time series of COVID-19 evolution observed over the last three months (up to 3/22/2020) in China. The model results revealed that 1) the transmission, infection and recovery dynamics follow the integral-order SEIR model with significant spatiotemporal variations in the recovery rate, likely due to the continuous improvement of screening techniques and public hospital systems, as well as full city lockdowns in China, and 2) the evolution of number of deaths follows the time FDE, likely due to the time memory in the death toll. The validated SEIR model was then applied to predict COVID-19 evolution in the United States, Italy, Japan, and South Korea. In addition, a time FDE model based on the random walk particle tracking scheme, analogous to a mixing-limited bimolecular reaction model, was developed to evaluate non-pharmaceutical strategies to mitigate COVID-19 spread. Preliminary tests using the FDE model showed that self-quarantine may not be as efficient as strict social distancing in slowing COVID-19 spread. Therefore, caution is needed when applying FDEs to model the coronavirus outbreak, since specific COVID-19 kinetics may not exhibit nonlocal behavior. Particularly, the spread of COVID-19 may be affected by the rapid improvement of health care systems which may remove the memory impact in COVID-19 dynamics (resulting in a short-tailed recovery curve), while the death toll and mitigation of COVID-19 can be captured by the time FDEs due to the nonlocal, memory impact in fatality and human activities. Fractional calculus can provide a promising tool in modeling biological phenomena, as reviewed recently by Lonescu et al. [1] . For example, fractional-derivative equation (FDE) models have been applied to capture complex dynamics in biological tissues [2] , tumor growth [3] , DNA sequencing [4] , drug uptake [5] , and salmonella bacterial infection in animal herds [6] . Most recently, the FDE models have been applied to model the Pine Wilt Disease [7] , the human respiratory syncytial virus (HRSV) disease [8] , the harmonic oscillator with a position-dependent mass [9] , human liver using Caputo-Fabrizio fractional derivative [10, 11] , and tumor-immune surveillance [12] . Motivated by these successful applications, this study tests whether the FDEs can be applied to model the dynamics and mitigation scenarios of coronavirus, an emerging and critical research area that has not been focused by the fractional calculus community before. The novel coronavirus disease 2019 (COVID-19) outbreak, a respiratory illness that started (first detected) in late December 2019, is a pandemic infecting >336,000 people in more than 140 countries with the average fatality rate of 4.4% globally (data up to 3/22/2020) [13] . The COVID-19 pandemic is infiltrating almost every aspect of life, damaging global economy, and altering both man-made and natural environments. Urgent actions have been taken but further effective and efficient strategies are promptly needed to confront this global challenge. To address this challenge and promptly guide the next efforts, it is critical to model the COVID-19 outbreak. Mathematical models are among the necessary tools to quantify the COVID-19 and, therefore, testing the applicability of such FDE models under this new global pandemic is the primary objective motivating this study. This study aims to model the COVID-19 evolution dynamics (i.e., transmission, infection, recovery, and death evolution) for representative countries with apparent coronavirus cases, including China, the United States (U.S.), Italy, Japan, and South Korea using mathematical models, most specifically the FDE models described previously. It should be noted that the COVID-19 spread in these countries experienced different starting (initiation and detection) times, an important fact to consider when applying these models. For example, China had passed the peak of the coronavirus outbreak, finally reaching a milestone with no new local infections on 3/19/2020 (79 days from the onset, 12/31/2019, in Wuhan, China), while the U.S. coronavirus cases soared past 10,000 on the same day. This study will apply the core characteristics of the COVID-19 outbreak obtained in China to estimate the COVD-19 spread in the U.S., as well as other countries where the number of affected people has not yet reached its number of peak cases. In addition, this new pandemic may last for a relatively longer time than expected [14] . No vaccine against SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) is currently available [15] . Indeed, a vaccine for prevention and infection control may not be ready before March 2021 for the COVID-19, considering a minimum of 4~5 weeks for trials and at least 1 year for safety evaluation and final deployment. Efficient strategies are therefore needed to mitigate the COVID-19 outbreak. Possible non-pharmaceutical scenarios, such as isolation of cases and contact tracing, can be evaluated using mathematical models [16] in order to identify the most efficient mitigation strategies going forward. This is another major motivation and the secondary task of this study using the FDE models. To address the questions mentioned above, this study is organized as follows. Section 2 proposes an updated SEIR model with an FDE-based component for COVID-19, where "S", "E", "I", and "R" stand for Susceptible, Exposed, Infectious, and Recovered people, respectively [17] . This model is then applied to fit and predict the COVID-19 spread in various provinces and major cities in China, resulting in abundant datasets to derive the core characteristics of the COVID-19 6 dynamics in evolution. Section 3 predicts the spread of COVID-19 in the other countries using the knowledge that is gained from the China case study. Section 4 proposes a fully Lagrangian approach with the time FDE to model the spatiotemporal evolution of COVID-19, and then applies it to evaluate non-pharmaceutical scenarios to mitigate the virus spread. Section 5 reports the main conclusions, specifically the feasibility of FDEs for capturing COVID-19 evolution. In the appendix, the impacts of the non-singular kernel and fractional derivative type on model results are further discussed for readers particularly interested in these fractional calculus techniques. The main contributions of this work, therefore, include 1) the first application of FDEs in modeling the evolution of the COVID-19 death toll, 2) an updated SEIR model with a transient recovery rate to better capture the dynamics of COVID-19 pandemic within China and for other countries, and 3) a particle-tracking approach based on stochastic bimolecular reaction theory to evaluate the mitigation of the spread of the COVID-19 outbreak. Several mathematical models have been applied for epidemic analysis of COVID-19. The most widely used one, to date, is the well-known SEIR model. For example, Peng et al. [18] (using data up to 2/16/2020) proposed a generalized SEIR model to successfully estimate the key epidemic parameters of COVID-19 in China, and predicted the inflection point and ending time of confirmed COVID-19 cases. The SEIR model was also applied by Li et al. [19] (using the observed data up to 2/6/2020) to compare the effect of city lockdowns on the transmission dynamics for different cities in China. The SIR model with time-dependent transmission and recovering rates was used by Chen et al. [20] (using data up to 2/20/2020) to analyze and predict the number of confirmed cases of COVID-19 in China. The SIR model was extended by Wang et al. [21] (using data up to 2/12/2020) to incorporate various time-varying quarantine protocols for assessing interventions on the COVID-19 epidemic in China. The SEIR model and its modifications were also successfully applied by others [22] [23] [24] [25] , mostly for assessing the early spreading of COVID-19 in China. Previous applications of the popular SEIR model, however, may contain high uncertainty since they had limited data access for only a short period of the COVID-19 outbreak. As will be shown below, the COVID-19 dynamics have changed dramatically in the last three months, likely due to the improvement/adjustment of screening/testing techniques, public hospital system capabilities, and the government's control policies for contagious diseases. An updated, much better version of the SEIR model is therefore needed and can be reliably built now for China as more detailed and complete datasets are available that includes the coronavirus outbreak data (i.e., # of people infected) well beyond the peak in China. Other models have also been used for specific purposes related to COVID-19, including the global metapopulation disease transmission model to predict the impact of travel limitations on the epidemic spread [26] , the transmission model for risk assessment [27] , and the synthetic contact matrix model for quantifying reproductive ratios for COVID-19 [28] . To the best of our knowledge, the classical stochastic epidemic models, such as the discrete or continuous time Markov chain model and the stochastic differential equation model, have not yet been applied for COVID-19 spread scenarios. A preliminary stochastic model built upon the FDE for evaluating COVID-19 spread in a city will be developed and applied in Section 4. This section focuses on the deterministic model. The classical SEIR model considers constant parameters, which may not adequately capture time-dependent dynamics of epidemic (the COVID-19) spread. Hence, the SEIR model is updated in Section 2.1, and then tested in Section 2.2. The classical SEIR model containing four populations (S, E, I, and R) takes the form [29] : where is the stock of susceptible population; is the number of persons exposed to or in the latent period of the disease; is the stock of infected population; is the stock of recovered population; r denotes the number of susceptible people whom the infected people contact daily; N is the sum of all the four groups of people ( ( ) ( ) ( ) ( ) (which is a constant representing the constancy of population N); p is the constant rate of infection (i.e., representing the probability for the infected people to transform the susceptible people into infected ones); is the constant rate for the exposed person to be transformed into an infected one; and is the constant recovery rate (defining the speed for the infected person to be cured or expired). To allow for possible time-sensitive rates and time nonlocal-dependency for COVID-19 evolution, we revise model (1) as: where represents the number of deaths (which is one component of I); is the rate for the healthy, susceptible person to be transferred to an infected one from exposed persons (note that COVID-19 patients in the incubation period might be contagious too); and is the number of healthy, susceptible people that are contacted by exposed people daily. Now, the infection rates can change with time, and the infected persons are removed from the risk of infection via a time-dependent rate term of ( ). If the recovered individuals can return to a susceptible status due to, for example, a loss of immunity, then the partial differential equation (PDE) (2d) for the time rate of change of requires one more (sink) term: , where is the rate of the recovered individuals returning to a susceptible status. We replace the integer-order derivatives in the classical SEIR (1) by the fractional-order derivatives in (2) , to capture the possible nonlocal impact, such as the memory impact or any apparent delay, in the COVID-19 outbreak. Herein we use the death evolution (2e) as an example. Hence, the FDE (2e) containing the death probability of (while the other patients are cured) is modified as: which is the Caputo fractional derivative [30, 31] with order ( ) (note that all the other indexes , , , and in model (2) are in the same range of (0, 1]). The Caputo fractional derivate was selected here because the initial condition for Caputo derivative takes the same form as that for the integer-order PDE. Particularly, the Caputo derivative allows the utilization of initial values (of the integer order) with known physical interpretations. When the order =1, model (2e) reduces to the classical integer-order PDE for the death evolution. The fractional PDE (2e) is used here for two reasons. First, the evolution of deaths and cures may be characterized by a random process, due to the fact that the exact time for the (recovered) person to be initially infected is unknown (i.e., the patient that died or was cured today may have been diagnosed yesterday or last week). Second, some patients may not be treated in time after being infected, making the death toll evolve with a time memory. Therefore, we extend the classical mass-balance equation of death cases using the FDE to characterize the random property and memory impact embedded in the temporal evolution of mortality. The fractional PDE (2e) and its classical version will be compared below using real (observed) data. We apply the SEIR model (2) does not exhibit any late-time tailing behavior, which supports the application of the integer-order, time-local model (2c). The fast decline of the late time "current infected population" is most likely due to the improvement of health care facilities, which tends to accelerate the recovery of infected people and remove the possible memory impact on the disease recovery evolution. Other complex SEIR models were also applied to model COVID-19 spread, such as the one proposed by Tang et al. [32] which has 12 rates/probabilities and 8 groups of people in SEIR. Numerical results show that, compared with the SEIR model (2), the complex model proposed by Tang et al. [32] (with solutions shown by the dotted lines in Fig. 1a ) accurately fits the observed data at the early stage, but then overpredicts the spread of COVID-19 observed after 2/12/2020. The dynamics of transmission for COVID-19, especially the recovery rate, therefore, changed over time in China, likely due to the time-dependent conditions (i.e. improvements/adjustment)s in medical care as mentioned above. A dynamic time-local SEIR model, therefore, may be preferred for modeling COVID-19 spread in China. In addition, compared to the FDE (2e), the best-fit solution using the classical PDE for death evolution (see the black, dotted line in Fig. 1a ) slightly overestimates the late-time growth of mortality. The actual death toll in Hubei province grew slower than that estimated by a constant rate model, indicating that the memory impact may affect the late-time dynamics of death and can be better captured by the FDE (2e). The best-fit solutions using model (2) fit the evolution of the infected and recovered populations well for the data recorded from Hubei province and three large cities closely related (such as in transportation or economic cooperation) to Wuhan, China ( Fig. 1) . The model was also shown to predict well the observed time series of COVID-19 spread from 3/8/2020 to 3/22/2020 for most places, except for Shanghai City (Fig. 1d) . This is likely due to the number of overseas COVID-19 cases imported into Shanghai, whereby the number of cases were observed to increase rapidly after 3/3/2020, causing inconsistency in the affected population and failure of the model. Shanghai Pudong International Airport, one of the two airports located in Shanghai City, is the eighth-busiest airport in the world and the busiest international gateway of mainland China. When excluding the coronavirus cases imported from overseas, model (2) predicts the COVID-19 spread evolution data in Shanghai (Fig. 2) . Therefore, model (2) works well for various places in China that may not be experience a significant influx of imported cases, as external sources can easily break the internal evolution, especially the asymptotic status, of COVID-19 in China. The resultant time-dependent recovery rate ( ) is depicted in Fig. 3 , where the rate fitted by the latest observation data point within the fitting period (i.e., 3/7/2020) remains stable in the following prediction period. The best-fit recovery rate is the highest for Shanghai (except for the impulse of ( ) for Wenzhou as discussed below), which is expected since Shanghai has the best public health system of all of these cities. Contrarily, Hubei shows the lowest recovery rate, likely due to its delayed response and the relatively limited public health capability at the beginning of the outbreak compared with Shanghai. Wenzhou exhibits an impulse in the infection and recovery dynamics of coronavirus (Fig. 1b) , different from that observed for the other places studied herein. On February 27, 2020, the number of Wenzhou's infected people suddenly declined, combined with the sudden increase of the number of people cured. This abrupt change can be effectively captured by the SEIR model (2) with an impulse in the recovery rate ( ). This recovery impulse is most likely due to a new hospital, the No. 2 Affiliated Hospital of Wenzhou Medical University, recently built in this city in early February, which significantly improved the public health system. The first discharged cases of coronavirus from this hospital appeared in late February 2020, resulting in the sudden increase in the total recovery rate. In addition, a relatively large number of people working in Wuhan returned to Wenzhou in late January, and it appears that the improved efficient screening process successfully identified the number of infected cases. The new cases were ~29,000 from 1/24/2020 to 1/31/2020 in Wenzhou (with an average of 3,600 new patients per day), who were then immediately centralized for treatment. It appears that this fast response helped to alleviate the spread of coronavirus in Wenzhou. The best-fit parameters of model (2) are listed in Table 1 , and the initial values for each group of people are listed in Table 2 . We reveal three behaviors in model parameters. First, the best-fit "S"-shaped recovery rate ( ) (Fig. 3) can be described by the sigmoid function ( ) ( ) (where a, b, and c are factors), showing that the recovery rate increases exponentially before reaching a stable condition. This increase is likely due to the fact that healthcare facility capabilities improved with time (with an accelerating rate) before reaching their asymptote or maximum capacity. Second, the rates and probabilities (r, , p, and ) affecting the COVID-19 transmission/infection evolution slightly change in space and remain stable for a given site ( Table 1 ). The small spatial fluctuation of these rates may be due to the similar strategies implemented by We also introduce an index C to quantify the infection severity of COVID-19 at different places: where is the maximum number of cumulative infected people at the given site. A smaller C represents a greater infection severity of coronavirus. There is a power-law relationship between the regional population N and the maximum cumulated number of infected people (Fig. 4) . This empirical formula may be used to approximate the largest cumulative number of infections, which will be applied in the subsequent section for predicting the COVID-19 evolution outside of China, where the coronavirus infection has not yet reached its peak number of cases. Different countries are applying different modes to slow the the spread of the COVID-19. In the next sub-sections, we discuss several representative countries and then fit/predict the virus spread there. To To decrease the acceleration of COVID-19 spread, Italy's mode is now similar to China: lockdown of the full population. The predicted COVID-19 spread in Italy is plotted in Fig. 6a . Although Italy has followed China's mode of national isolation, the number of infected people increased rapidly from 3/14/2020 to 3/19/2020 (~495 new cases per day). To account for the delayed national quarantine compared with China, we decrease the C index (while also increasing the upper limit of the cumulative infection). The COVID-19 evolution prediction results show that there may be a turning point in the next two weeks when the current infected cases begin to decline. We also separate the death toll from the number of recovered cases. South Korea's mode of combating the spread of COVID-19 is through fast detection and tracking of the disease. South Korea is using efficient mobile diagnoses tests and accurate tracing of infected cases to maintain a low death rate even with a large infected population. The mobile method can test 20,000 people per day (the maximum capability on 3/12/2020), and apps for cell phones and/or credit cards can accurately track the routes of infected people with the help of local government (without an invasion of privacy), so that warnings can be immediately delivered to the general population to obviate the places with high risk. The current infected population may have passed its peak number of cases around 3/20/2020, and the prediction shows that the COVID-19 outbreak may be well controlled in ~35 days from 3/22/2020 (Fig. 6b) . Japan's mode of combating the spread of COVID-19 is compatible with that of the U.S., in addition to other changes such as enhanced education/outreach and rapid treatment of the infected cases. Specific policies include social distancing (which might be a key barrier to the spread of the novel coronavirus), personal hygiene, and quarantine of the infected cases. The current data and modeling results (Fig. 6c) show that Japan has an efficient way so far to limit the maximum population infected and slow the spread of COVID-19, while this outbreak may last for a while. The model-predicted COVID-19 spread in the U.S., using the fitted infection and recovery rates from China (Fig. 5) , reveals the impact of one possible mitigation scenario for COVID-19: coronavirus lockdowns, which have now been implemented by some states in the U.S. such as California and New York. Other non-pharmaceutical options can and should also be evaluated using mathematical models, considering the recent surge of infected cases in the U.S. When the number of infected persons is initially small compared to that of susceptible people, the infected and susceptible people are not well mixed and hence the system is not homogeneous. Under such conditions, a stochastic model is needed, as the deterministic, continuum models (such as the SEIR model) assume well-mixing of components for a homogeneous system [33, 34] . Hence, this section develops and applies a stochastic model to evaluate non-pharmaceutical scenarios for mitigating the COVID-19 with a small number of initial infections. The random walk based stochastic model for COVID-19 spread is analogous to a mixing-limited bimolecular reaction-based mechanism/condition [35] . When a reactant A particle (representing a susceptible individual A) meets a reactant B particle (representing an infectious person B), a chemical reaction may occur if the collision energy is large enough to break the chemical bond (meaning that the susceptible person A may be infected if satisfying additional criteria such as A and B are close enough, and B touches his/her face after receiving coronavirus from A). Therefore, the condition of A being infected is not deterministic but rather a random, probability-controlled process. This probability is related to various factors, such as the duration that A and B are in contact, the infectivity rate, and the distance between the two people, which may be characterized parsimoniously by the interaction radius R that controls the number of reactant pairs (susceptible + infectious) in a potential reaction (infection) [35] . Hence, the core of the random walk stochastic model to simulate COVID-19 spread is to define the interaction radius The analogous development and similarities between bimolecular reactions and the SIR model can also be seen from their governing equations. The time-dependent SIR model takes the form [36] : where ( ) and ( ) denote the transmission rate and recovery rate at time t, respectively. The rate equations for irreversible bimolecular reaction take the form [35] : where , , and denote the concentrations of A, B, and C, respectively; and ( ) is the forward kinetic coefficient of reaction. Equations (5a) and (5b) are functionally similar to equations (4a) and (4b), respectively, if the recovery rate ( ) . Therefore, following the argument in Zhang et al. [35] and Lu et al. [37] , we derive analytically the interaction radius R for the SIR model (4): where V denotes the volume of the domain, is the time step in random walk particle tracking, is the mass (or weight) carried by each A particle, is the initial concentration of A (which can be assumed to be the normalized value 1 here), and denotes the initial number of susceptible people. The movement of A, B, and C can be described by the following time FDE [35] : where denotes A, B, and C, respectively; denotes the fractional capacity coefficient (which controls the ratio between the immobile and mobile population); denotes the mean moving speed; and denotes the macrodispersion coefficient. After defining the interaction radius R, the particle tracking scheme proposed by Zhang et al. [35] and Lu et al. [37] with particle trajectories defined by the time FDE (7) can be applied to model the transmission of coronavirus between the susceptible and infectious people. In addition to pharmaceutical strategies including vaccine and therapeutic drug development, and herd immunity that may either take a while or have a high risk, non-pharmaceutical scenarios can be tested. Several particle-tracking based stochastic models were proposed recently [38] to evaluate non-pharmaceutical scenarios to mitigate coronavirus spreading in a city. Here we evaluate three related scenarios (described below) using the stochastic model proposed above. In the stochastic model, we assume that 10 days after being infected, the person will be removed because of being cured or expired (dead). This is because the median disease incubation period has been estimated to be 5.1 days [39] . For simplicity purposes, the interaction radius R (6) remains constant, since the constant interaction radius was found to be able to efficiently capture the temporal variation of effective reaction rates in mixing controlled reaction experiments or simulations [35, 37] . The initial number of A and B particles is 10,000 and 4, respectively. The Lagrangian solutions of the COVID-19 outbreak for the three scenarios are depicted in Fig. 7 . Modeling results for scenarios 1, 2, and 3 show a peak in the curve of newly infected people at time t=28, 65, and 32, respectively, demonstrating that the virus spreads the fastest for the scenario without mitigation constraints (i.e., scenario 1, where the number of the total infected people or cases increases by one order of magnitude every 10 days in the rising limb), as expected. However, the value for this peak (scenario 1, =198 people) is lower than that for scenario 3 (=267 people), although the total number of the infected people for scenario 1 (9,336) is slightly larger than scenario 3 (9,322). This may be due to a greater separation of infection cases for the higher number of initial coronavirus carriers in scenario 1, which causes a lower and relatively flatter COVID-19 evolution peak compared to that of scenario 3. Scenario 2 has the lowest peak value (=121 people) and the most-delayed peak in the curve for new cases, and the total infection time is almost doubled compared to the other two scenarios, indicating that people living with strict social distancing may also suffer from a much longer period of COVID-19 threat. It is also noteworthy that the overall trend of the solution (simulation) of scenario 1 (initial surge without special constraints, Fig. 7a ) is similar to that for Italy which had a delayed response to the COVID-19 outbreak initially (Fig. 6a) , and scenario 3 solution (a lower peak value and a longer duration due to social distancing, Fig. 7b ) is similar to that observed in Japan where social distancing actions have been implemented (Fig. 6c) . The simulated particle plumes plotted in Fig. 8 reveal the subtle discrepancy between the three mitigation scenarios. Scenario 1 assumes that four initial cases were initially located on the right side of the city, while the whole population (10,000 susceptible persons) was distributed randomly in the 11 domain (Fig. 8a) . The trajectory of each person is assumed to follow (two-dimensional) Brownian motion with retention (described by Eq. (7)), to capture the random vector for each displacement and the random waiting time between two consecutive motions (described by the time fractional derivative term in Eq. (7)). The virus moved quickly from east to west (Figs. 8b and 8c) , spreading over the entire city before all the infected people were cured or expired (dead) at time t=69 (Fig. 8d) . A total of 664 susceptible people (6.6% of the total population) distributed randomly around the city were never infected. Scenario 2 assumes that social distancing can reduce the infection probability, which can be characterized by a smaller reaction rate or a smaller interaction radius in our Lagrangian approach. The virus was spreading at a much lower rate from west to east than that in scenario 1 (Figs. 8e~8g), reaching a stable condition (i.e., # of cases neither increasing or decreasing) at a later time (t=125) and leaving more susceptible people unaffected (2,734 total, which is 27.3% of the population). Therefore, social distancing is effective in limiting the spread of coronavirus among people. Note, however, this scenario assumes that every person in this city strictly maintains social distancing; otherwise a surge of infections may occur the same way as that shown in scenario 1. Scenario 3 assumes self-quarantine. Notably, not all of the infected people can be effectively quarantined due to the following: 1) people can be infected without coronavirus symptoms; 2) people in the incubation period can transmit the infection; and 3) limited health care facilities and capabilities for the large influx of patients. For example, according to Imai et al. [40] , many infected people could not be appropriately screened initially in Wuhan City, China. Under this condition, we assume that 50% of people infected and diagnosed are immediately quarantined, while the remaining infected people (Fig. 8i) can still cause the spread of coronavirus (Figs. 8j~8l). Self-quarantine, therefore, may not be as effective as maintained social distancing. Fractional calculus provides a useful tool to modeling complex dynamics in biology, and this study extended the FDE for modeling the coronavirus outbreak. The COVID-19 is a pandemic Third, a stochastic model based on the Lagrangian scheme for the time FDEs, analogous to a mixing-limited reaction mechanism model, showed that self-quarantine may not be as effective as strict social distancing, since not all the infected people can be diagnosed and immediately quarantined. While strict social distancing can apparently slow COVID-19 spread, the pandemic may last longer. This is another case that fractional calculus may be used to explore COVID-19 outbreak. Therefore, one of the main contributions of this study is to extend the application of FDEs to model dynamics and mitigation scenarios of the coronavirus spread. This appendix quantifies the potential impacts of the non-singular kernel and the type of the fractional derivative on the COVID-19 death toll simulation. First, the nonsingular time-fractional definitions also provde promising modeling tools for real-world fractional dynamics. For example, the Atagana-Baleanu fractional derivative is defined as [41] : where ( ) ∑ ( ) represents the single-parameter Mittag-Leffler function. We employ this definition to extend the SEIR model and then use the finite difference method to simulate the evolution of the COVID-19 death toll, which leads to the following solver: (Fig. 9) , which is consistent with the conclusion in the previous study (see [42] ). The resultant death roll peak is lower than that simulated by the conventional kernel function, but it does not mean that the kernel used in the Caputo fractional derivative is no longer valid. Indeed, the solution of the Caputo fractional derivative can capture well the obsvered death peak (Fig. 1a) . Applications of the nonsingular kernel functions deserve further study in the future. Second, the Riemann-Liouville fractional derivative is defined as [43] : The Caputo fractional derivative listed in Section 2.1 and the Riemann-Liouville fractional derivative (9) are related by [44] : where , and the operators and represent the Caputo and Riemann-Liouville fractional derivatives, respectively. ☒ 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. Credit author statement: The role of fractional calculus in modeling biological phenomena: A review Fractional calculus models of complex dynamics in biological tissues Modeling with fractional difference equations Fractional dynamics in DNA Simulation of drug uptake in a two compartmental fractional model for a biological system On fractional SIRC model with salmonella bacterial infection Semi-analytical study of Pine Wilt Disease model with convex rate under Caputo-Febrizio fractional order derivative A new fractional HRSV model and its optimal control: A non-singular operator approach The fractional features of a harmonic oscillator with position-dependent mass A new study on the mathematical modelling of human liver with Caputo-Fabrizio fractional derivative New aspects of time fractional optimal control problems within operators with nonsingular kernel A new fractional model and optimal control of a tumor-immune surveillance with non-singular derivative operator Johns Hopkins Coronavirus Resource Center A mathematical model for the novel coronavirus epidemic in Wuhan, China COVID-19 -New insights on a rapidly changing epidemic Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts A contribution to the mathematical theory of epidemics Epidemic analysis of COVID-19 in China by dynamical modeling The lockdown of Hubei province causing different transmission dynamics of the novel coronavirus (2019-ncov) in Wuhan and Beijing. medRxiv A time-dependent SIR model for COVID-19 An epidemiological forecast model and software assessing interventions on COVID-19 epidemic in China Effectiveness of control strategies for Coronavirus Disease 2019: a SEIR dynamic modeling study SEIR Transmission dynamics model of 2019 nCoV coronavirus with considering the weak infectious ability and changes in latency duration Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study Preliminary prediction of the basic reproduction number of the Wuhan novel coronavirus 2019-nCoV The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak Risk assessment of novel coronavirus COVID-19 outbreaks outside China Estimation of country-level basic reproductive ratios for novel Coronavirus (COVID-19) using synthetic contact matrices Linear model of dissipation whose q is almost frequency independent-II Impact of absorbing and reflective boundaries on fractional derivative modes: Quantification, evaluation and application Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions An introduction to stochastic processes with applications to biology Discrete and continuous SIS epidemic models: A unifying approach Evaluation and linking of effective parameters in particle-based models and continuum models for mixing-limited bimolecular reactions Networks: An Introduction Lagrangian simulation of multi-step and rate-limited chemical reactions in multi-dimensional porous media Why outbreaks like coronavirus spread exponentially and how to "flatten the curve The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application Report 2: Estimating the potential total number of novel Coronavirus cases in Wuhan City New Fractional Derivatives with Nonlocal and Non-Singular Kernel: Theory and Application to Heat Transfer Model Time fractional derivative model with Mittag-Leffler function kernel for describing anomalous diffusion: Analytical solution in bounded-domain and model comparison Fractional integrals and derivatives Stochastic models for fractional calculus Resources. Bin Jin: Supervision, Funding acquisition