

ORIGINAL ARTICLE 

Year : 2022  Volume
: 12
 Issue : 2  Page : 95107 

Statistical inference of COVID19 outbreak: Delay distribution effect in EQIR modeling of epidemic
Mahnoosh Tajmirriahi^{1}, Zahra Amini^{1}, Rahele Kafieh^{1}, Hossein Rabbani^{1}, Ali Mirzazadeh^{2}, Shaghayegh Haghjooy Javanmard^{3}
^{1} Medical Image and Signal Processing Research Center, School of Advanced Technologies in Medicine, Isfahan University of Medical Sciences, Isfahan, Iran ^{2} Department of Epidemiology and Biostatistics, Institute for Global Health Sciences, University of California, San Francisco, California, United States ^{3} Applied Physiology Research Center, Isfahan Cardiovascular Research Institute, Isfahan University of Medical Sciences, Isfahan, Iran
Date of Submission  18Jun2021 
Date of Decision  02Aug2021 
Date of Acceptance  28Oct2021 
Date of Web Publication  12May2022 
Correspondence Address: Rahele Kafieh Medical Image and Signal Processing Research Center, School of Advanced Technologies in Medicine, Isfahan University of Medical Sciences, Isfahan Iran Zahra Amini Medical Image and Signal Processing Research Center, School of Advanced Technologies in Medicine, Isfahan University of Medical Sciences, Isfahan Iran
Source of Support: None, Conflict of Interest: None
DOI: 10.4103/jmss.jmss_134_21
Background: The world is experiencing another pandemic called COVID19. Several mathematical models have been proposed to examine the impact of health interventions in controlling pandemic growth. Method: In this study, we propose a fractional order distributed delay dynamic system, namely, EQIR model. In order to predict the outbreak, the proposed model incorporates changes in transmission rate, isolation rate, and identification of infected people through time varying deterministic and stochastic parameters. Furthermore, proposed stochastic model considers fluctuations in population behavior and simulates different scenarios of outbreak at the same time. Main novelty of this model is its ability to incorporate changes in transmission rate, latent periods, and rate of quarantine through time varying deterministic and stochastic assumptions. This model can exactly follow the disease trend from its beginning to current situation and predict outbreak future for various situations. Results: Parameters of this model were identified during fitting process to real data of Iran, USA, and South Korea. We calculated the reproduction number using a Laplace transformbased method. Results of numerical simulation verify the effectiveness and accuracy of proposed deterministic and stochastic models in current outbreak. Conclusion: Justifying of parameters of the model emphasizes that, although stricter deterrent interventions can prevent another peak and control the current outbreak, the consecutive screening schemes of COVID19 plays more important role. This means that the more diagnostic tests performed on people, the faster the disease will be controlled.
Keywords: COVID19, EQIR epidemic model, fractional differential equation, stochastic differential equation
How to cite this article: Tajmirriahi M, Amini Z, Kafieh R, Rabbani H, Mirzazadeh A, Javanmard SH. Statistical inference of COVID19 outbreak: Delay distribution effect in EQIR modeling of epidemic. J Med Signals Sens 2022;12:95107 
How to cite this URL: Tajmirriahi M, Amini Z, Kafieh R, Rabbani H, Mirzazadeh A, Javanmard SH. Statistical inference of COVID19 outbreak: Delay distribution effect in EQIR modeling of epidemic. J Med Signals Sens [serial online] 2022 [cited 2022 May 25];12:95107. Available from: https://www.jmssjournal.net/text.asp?2022/12/2/95/345070 
Introduction   
In December 2019, the epidemic was created by a member of coronaviruses, namely, 2019nCoV. Till now, the virus has been spreading worldwide, and millions of people have been infected. Extensive mathematical studies in various fields have been done to find the policies to control the disease in its prevalence. Mathematical models can indicate how infectious diseases spread and what the trend of epidemic is. They can be divided into two general purpose models, forecasting models and mechanistic models. Forecasting models try to fit curves to data using different methods from extrapolations using simple models such as exponential growth^{[1]} to machine learning methods^{[2]} or regressions.^{[3]} These models are simple and are only used in the early days of the outbreak. Mechanistic models can be divided into two groups as follows.
Selfexciting branching point process^{[4]}: These models can easily fit to data and are used to estimate reproduction number and estimation rate. However, they can be used in initial phase of epidemic when the number of infected individuals are small compared to the population size.^{[5]}
Dynamic epidemiological models^{[6]}: In these models, the population is divided into several states or compartments, and different parameters are attributed to each state such as transmission rate and mortality rate. Then, considering specific assumptions for the population, differential equations for each state of the system are obtained and solved. In fact, these models are one of the best ways to explore longterm epidemiologic outcomes and to implement various assumptions about disease and controlling measures.^{[7],[8]} They have been widely used for the COVID19 epidemic because of its simplicity, parametric approach, and good results. Some of these models have contributed improvements in susceptibleinfectedrecovered model with a variable total population^{[9]} and time varying parameters.^{[10],[11],[12]} Some others tried to use susceptibleexposedinfectedrecovered (SEIR) models^{[13],[14]} including control strategies in the SEIR model,^{[15],[16]} and utilizing fractional order of differential equations in SEIR.^{[17],[18]} Meanwhile, some studies have been done to find new dynamic models with more states and relations to achieve more accurate results such as.^{[19],[20],[21],[22]}
However, time frame, randomness, basic assumptions of model, model generality, and amount of data matching are some important factors in evaluating a model, while a simple model may not satisfy all these conditions. Assumptions such as social knowledge, fluctuation in the behavior of populations, and randomness of transmission rate must be inserted in model through fractional stochastic differential equations (SDE).^{[23],[24]} SDE has been widely used in modeling of long range dependent events to indicate the memory of the systems.^{[25],[26]}
In this study inspiring by Chen et al.,^{[27]} we use a dynamic stochastic and deterministic model of fractional order to describe the regional outbreaks of COVID19 and to predict its future trend. In the proposed model, in addition to considering specific states for different parts of the population, latent periods, isolation, and treatment durations are also considered by adding distributed delays to the model. Furthermore, we include memory effects in the model dynamics by means of fractional derivatives. The stochastic nature of population behavior is also considered in this model by assuming the transmission rate as a realization of Brownian motion ^{More Details} stochastic process. All these assumptions lead to solve a fractional ordered distributed delay SDE.^{[28]} We also obtain the basic reproduction number of the epidemic model. After finding parameters of deterministic model from available data through fitting process, results show a good performance of the model in fitting and prediction the trends of epidemic, besides its ability to justify the real behavior of different communities against disease. Stochastic model is used to explore shortterm and longterm trajectories of COVID19 and determines the range of outbreak size.
This paper is organized as follows. After this introduction, we explain deterministic and stochastic model and their parameters in Section 2. In this section, we also discuss about existence and persistence of the model by calculating reproduction number. The numerical results of the model are presented in Section 3. Discussion about the model and conclusions are presented in Sections 4 and 5.
Materials and Methods   
Model description
The important point in spreading trend of COVID19 is transmission from an exposed patient with no symptoms to others which is challenging in SEIR models. We propose an EQIR model,^{[27]} in which the population is divided into four states. The first state is exposed (E) individuals with no symptoms, due to government decisions some of the exposed individuals may be isolated in their latent period. Hence, the second state (Q) contains the quarantined individuals. Third state (I) includes individuals who are confirmed and have no ability of spreading disease. The last state contains removed (R) people (recovered or dead). Since environmental fluctuations can impact transmission rate β, it can be assumed a timevarying function such as exponentially decreasing function.^{[29],[30]} This decreasing trend may be caused because of government interventions such as isolations and public closures at specific time. Here, we assumed following function for transmission rate.
where T is the time of starting interventions and d is decaying coefficient. Another approach to display the time varying nature of transmission rate is considering a random generating process for it by adding a white Gaussian noise to transmission rate β. Hence, we expect extra random curves around the main curve indicating how far the outbreak will be from controllability in short terms. Further assumptions and parameters are adopted for both stochastic and deterministic models. Some important assumptions such as considering extremely large population of susceptible, independency, and homogeneity of infected individuals, and disease equilibrium condition are implicitly considered.
 We assume that it takes averagely 𝜏_{E} days for an exposed individual to have symptoms and be identified. After this latent period, the exposed individual comes to infected state and in this state can no longer transmit the coronavirus to others
 The second assumption is about those exposed individuals who are in the incubation period and are isolated from others. They cannot spread the disease and will be diagnosed during next 𝜏_{Q} days. After 𝜏_{Q} ≤ 𝜏_{E} days, they will be identified as infected ones. We assume q as the isolation rate of currently exposed people
 Infected individuals averagely experience 𝜏_{R} days before complete cure or death. We indicate mortality or recovery rate of this disease by parameter δ. We ignore the natural mortality in our population
 We assume that every exposed, isolated, and treated individual experiences different time period to change his/her states. This can be considered by means of distributed delay functions h_{E}(t), h_{Q}(t), and h_{R}(t) respectively. Where h is a Gaussian probability distribution function. The mean value of this distribution is adopted to parameters 𝜏_{E}, 𝜏_{Q}, 𝜏_{R}. Different variances of each distribution can be determined by fitting the model to real data. Zero variance value, converts h(t) to a Dirac delta function and means that all individuals experiences same latent. These parameters can determine how fast the exposed individuals are identified and isolated
 Ordinary compartmental models implicitly follow a Markov epidemic process, in which the states of individuals are independent from previous steps at each time. This process is socalled memory less. However, when a disease spread within a human population, knowledge of individuals affects their response and it shows an implicit memory. However, the effect of memory of earlier times decreases over time. This decay of longrange memory can be controlled by the order of the fractional derivatives in the corresponding nonlinear fractional differential equations of the system.^{[31]} We assume such a fractional order model and use the RiemannLiouvile definition of fractional derivative of order α as follows:^{[32]}
where function is Γ() Gamma function and n is an integer. Considering above assumptions, the continuous time of deterministic EQIR model is constructed as:
where α_{i}, i = 1 …4 can be obtained by parameter identification. In order to simulate this model, we use discrete time implementation of the model in two cases deterministic and stochastic models, both of them are presented in Appendices A and B.
Optimization problem
Each model has several parameters which must be carefully determined. Some of them must be extracted from epidemiological texts and others must be identified by fitting models to the available data. In second case, we encounter an optimization problem. For this purpose, Let θ denotes vector of parameters of specific state which must be estimated, then our parameters identification problem has the following form:
where ∁_{obj} denotes the available data of that component and ║║_{2} is related to the 2norm. In order to solve this problem, we use LevenbergMarquadt (LM) method^{[37]} for deterministic model. This method uses an iterative algorithm to solve nonlinear least square problems. In stochastic disease modeling, x_{t} is real number of infectious individuals which may be controlled by parameters θ such as transmission rate, y_{t} is reported data, and the goal of modeling is determining the posteriori distribution of parameters that is p(θ│y_{t}).^{[38]} An ordinary approach to estimate distribution of parameters is pseudomarginal Markov Chain Monte Carlo (PMCMC) method.^{[39]} In this method, samples of desired probability density function p(θ│y_{t}) are generated by using a proposal function which usually considered a normal distribution. A sampler generates some samples (particles) from proposed distribution. Then, the acceptance ratio of each sample is calculated and decision is made to keep or reject that sample. After lots of iterations, remaining particles represent p(θ│y_{t}).
Reproduction number
By means of the same method used in Kiskinov, Zahariev and ÖZalp, Demirci's study,^{[40],[41]} one can show that there exists a unique, nonnegative solution for fractional differential equations depicted in Equation 2. For analysis the stability of the model, we adopt a Laplace transform approach to show the prevalence of individual in each state explained in.^{[42]} For this purpose, we describe the basic reproduction number, typically denoted with R_{0}, of the model. For a deterministic epidemic model, if R_{0}<1 the disease dies out and if R_{0}<1 the disease persists and outbreak will happen, meaning that there exists an epidemic equilibria, in addition to the diseasefree equilibria. For a stochastic model, this number indicates existence of unique stationary diseasefree distribution or Epidemic distribution in states. There are several methods to estimate R_{0}.^{[43]} However, in an ordinary SEIR model, R_{0} = βμ_{1} where μ_{1} is average infectious latent period. In our proposed model, we generalize the SEIR model by first assuming that latent periods follow a Gaussian distribution and second by using fractional order derivatives. It is shown^{[41],[44]} that the fractional order of model does not affect R_{0}, however, distribution of latent periods affects it, and controlled reproduction number is defined as R_{c} = βμ¯_{1}.^{[42]} Where β¯ indicates reduction in β due to system outflows occurred before exposed state. is defined as follows:
where ꞷ denotes the outflow of state, F¯_{1} is survival function and L [.] denotes Laplace transform. For f_{I} ~ N(𝜏_{I}, σ_{I}) we have:
[Figure 1] shows the block diagram of the proposed model. It is worthy to mention that COVID19 can spread during latent period and transmission rate is considered for exposed state and β¯ = β What about infectious latent period? We assumed that individuals in state E are isolated with rate q and the out flow of this state is q. Isolated individuals are no longer infectious; they may recover or die after a latent period of distribution h_{Q}(t) or they may become symptomatic and escape to state I with outflow rate of q. Furthermore, Exposed individuals may become symptomatic and identified after a latent period of distribution h_{E}(t) and will recover or die after a removal latent period of distribution h_{R}(t). Thus, R_{c} can be estimated as follows:
where μ_{R} = 𝜏_{R}, θ = pr(isolatedexposed) and μ¯_{E}, μ¯_{Q} defend as follows:
Parameter R_{c} can be decreased by either increasing q or decreasing 𝜏_{E} and 𝜏_{Q} and this can be done by identifying exposed individuals via screening schemes. Parameter θ can be considered equal 0.5 and means for every exposed individual the probability of being isolated or not is the same.
In stochastic case, we assumed that infectious rate β follows a normal distribution β ~ N(β_{0}, σ), that is R_{cs} also has a normal distribution as follows:
It is clear that any change in transmission rate can directly be interpreted as change in the reproduction number.
Results   
In this section, we provide numerical solutions to illustrate main theoretical results on the proposed deterministic and stochastic models. For this purpose, the model is implemented in MATLAB 2017a. Starting points are considered to be E (0) = 2, I (0) = 0, Q(0) = 0, R(0) = 0. Parameters of these models and their corresponding descriptions are shown in [Table 1].
Model fit was done using cumulative confirmed, dead, and recovered regional data available on World Health Organization documents.^{[45]} Based on the model, we used the LM and PMCMC methods to solve the optimization problem and identify parameters. To show that our model can adapt to changes in parameters over time, we fitted it with accumulated data of Iran, USA, and South Korea with two scenarios of short and long term of predicting data. Identified parameters of models for three selected countries are illustrated in [Table 2].
In the left panel of [Figure 2], we illustrated the sample distribution of transmission rate acquired from PMCMC method for three countries up to August [Table 2] for more clarity. According to this figure, posteriori distributions of transmission rate follow Gaussian distribution due to random walk assumption of transmission rate. In the right panel, we provided a plot of R_{c} with respect to parameters of (10), (11). As revealed in this figure, in the stochastic assumption, R_{c} has a normal distribution with mean value equals to the R_{c} of deterministic model and a variance depending on the variance of β.  Figure 2: Left: Normal distributions of transmission rate in three countries up to August. Right: Variation of Rc with respect to parameters depicted in Equations (10), (11) for σ_{q} = 2,𝜏_{q} = 5,σ_{E} = 2,𝜏_{E} = 8,μ_{R} = 11,β_{0} = 0.27, q = 0.45,σ = 0.04
Click here to view 
Results of deterministic model
Short term prediction
We started from the accumulated data of South Korea (from February 20 to May 5). This country was very successful in controlling the spread of the disease in early stages [Figure 3]. Besides low transmission rate, delay parameters (𝜏_{E}, 𝜏_{Q}, σ_{E}) were also very low which indicated that the exposed people have been identified and isolated very soon. The variance of this identification time, confirmed that approximately all exposed people were identified after at most 11 days. We then fitted the model to data of USA (from March 1 to May 5) and Iran (from February 20 to May 5) which their confirmed cases have been growing rapidly [Figure 4]. The spread rate of disease and delay parameters have higher values in the USA and Iran than South Korea's which indicate weaker performance to identify exposed cases. Indeed, approximately all exposed people were identified after at most 16 days [Table 2]. We assumed that, for both countries, government interventions have been applied 10 days after beginning of the epidemic (T = 10). According to prediction results, if these interventions continued, the COVID19 pandemic could have been controlled in Iran by June and in the USA by August.  Figure 3: Numerical results of data fitting in deterministic model for South Korea. (a) Cumulative confirmed, (b) daily confirmed, (c) cumulative dead, (d) daily dead cases. Continuous line is prediction result and grey dots correspond to real data
Click here to view 
 Figure 4: Numerical results of data fitting in deterministic model; (a), (b), (c), (i), (j), (k) for IRAN. (d), (e), (f), (l), (m), (n) for USA. Continuous line is prediction result and gray dots correspond to real data
Click here to view 
[Figure 5] shows epidemic trajectories in the presence and absence of interventions in the USA and Iran. We plotted the effective reproduction number over time reflecting the impact of transmission rate reduction. It can be seen that although interventions effected the trend of epidemic and limited its spread, but does not decrease very much. We tried to find an appropriate decaying coefficient to predict if epidemic stops sooner. For d = 0.008 reproductive number decreases more rapidly and equals to one late in June and the epidemic completely goes to extinction in the USA and Iran. Comparing with current decreasing trend of transmission rate, either stricter rules must be enforced or patient screening schemes must be implemented continuously.  Figure 5: (a and c) Epidemic trajectories derived from the decaying and nondecaying transmission rates during the first 100 days of the epidemic for USA and Iran respectively. (b and d) Effective reproduction number over time reflecting the impact of transmission rate reduction for two decaying factors
Click here to view 
Longterm prediction
In early stage, South Korea was so successful in handling the COVID19 pandemic, but this country experienced the second wave of infection from August 2020. We updated model's parameters and depicted the results in [Figure 6]. According to this figure, this pandemic will be controlled in South Korea by the end of 2020.  Figure 6: Numerical results of data fitting in deterministic model for South Korea.(a) Cumulative confirmed, (b) daily confirmed, (c) cumulative dead, (d) Daily dead cases. Continuous line is prediction result and grey dots correspond to real data
Click here to view 
Furthermore, in Iran and the USA, social limitations decreased in June when public places, jobs, administrative activities, etc., were reopened. Contact rate of individuals increased, and consequently, transmission rate was intensified, while isolation rate decreased. This caused spreading the disease into new populations or resurging in places that let down their guard too soon. We applied these points in our model by changing parameters of model from June to September (second peak) and after September (third peak) fitting to new data in Iran and the USA. Since readjustments of parameters have been done at a certain time (like a step function) and the derivative function of cumulative data has been used to indicate daily infected and death persons, discontinuities have been created in the daily curves (like Dirac delta functions) according to [Figure 7]c, [Figure 7]f, [Figure 7]k, and [Figure 7]n. In order to precisely follow the outbreak trend, we have detected the amount of restrictions and cumulative data in mentioned countries to adjust the decaying coefficient, d, in (1). For interventions leading the reduction of β, we decreased decaying coefficient at the time point of starting interventions, and we increased decaying coefficient when the restrictive interventions decreased. We have done the precise adjustment by a new fitting of updated cumulative data to the model and have predicted outbreak trend.  Figure 7: Numerical results of data fitting in deterministic model; (a), (b), (c), (i), (j), (k) for Iran. (d), (e), (f), (l), (m), (n) for USA. Continuous line is prediction result and gray dots correspond to real data
Click here to view 
We found that, as it was predicted, transmission rate and delay parameters have increased while isolation rate has decreased [Table 2] and [Figure 7]. According to prediction results, if these conditions continued, COVID19 pandemic will be continued in Iran and the USA by the end of 2020. Hence, unfortunately, there will be a deplorable situation. This indicates that stricter and more effective preventive laws and more diagnostic tests should be considered to avoid more disasters. We have considered a decreasing function to determine the decreasing nature of the transmission rate. Indeed, when countermeasures such as voluntary social distancing by the population, short term interventions such as global local lockdowns and general quarantine occur, its effect slowly decays the transmission rate for a period of time. Further abrupt changes in these interventions can be incorporated by adjusting the value of decaying coefficient at certain times. However, if restrictions were not sufficient from the beginning time (T) or were not observed properly, their effect in reducing the transfer rate may disappear sooner, and amounts of underestimation may occur in the prediction when the model is used for longterm predictions. This may necessitate more increase in decaying coefficients.
Results of stochastic model
In order to simulate the stochastic model, we used all fitted parameters of deterministic model except β. Parameter β identified in the deterministic model, was used as an initial guess of transmission rate for the PMCMC method to estimate variance of transmission rate [[Figure 2] (left panel)]. Documentation and implementation of MATLAB example of this method are available on.^{[46]} Since transmission rate follows the random walk, general trend of epidemic is preserved during ensemble realization of the model. These fluctuations in transmission rate can model daily changes in the number of infected cases better than deterministic model
Results of stochastic model are shown in [Figure 8] and [Figure 9]. Corresponding to [Figure 8] (shortterm prediction), if health restrictions had been continued, probability of outbreak with greater size would have been very low and the outbreak would have been going to be slowly disappeared by the end of June in Iran and by the end of August in the USA.  Figure 8: Short term stochastic epidemic realizations of the proposed model in South Korea (a and c), Iran (e and g), USA (I and k).The black dots correspond to the actual outbreak trajectory and the cyan blue lines correspond to 500 stochastic realizations of cumulative and daily number of infected cases respectively. Red line indicates median. The corresponding distribution of outbreak sizes are shown for South Korea in (b and d), Iran in (f and h), and USA in (j and l). The vertical dashed line indicates the predicted COVID19outbreak size
Click here to view 
 Figure 9: Long term stochastic epidemic realizations of the proposed model in South Korea (a and c), Iran (e and g), USA (I and k). The black dots correspond to the actual outbreak trajectory and the cyan blue lines correspond to 500 stochastic realizations of cumulative and daily number of infected cases respectively. Red line indicates median. The corresponding distribution of outbreak sizes are shown for South Korea in (b and d), Iran in (f and h), and USA in (j and l). The vertical dashed line indicates the predicted COVID19 outbreak size
Click here to view 
However, due to removal of restrictions, we changed the parameters of model as we described in section 3.1.1 and observed that if current conditions remain, the pandemic will last at least by the end of 2020 in mentioned countries, corresponds to [Figure 8] (long term prediction), and more and more spreading peaks will occur. In [Figure 8] and [Figure 9] the corresponding distributions of outbreak size of cumulative [b, f, j] and daily [d, h, l] numbers are also represented for three countries. According to [Figure 9], the distribution of outbreak size tends to have a large cumulative and daily frequencies of small sizes meanwhile a very heavier tail than its corresponding distribution in [Figure 8]. It indicates that final size of outbreak will dramatically increase compared to the scenario that could happen in shortterm prediction. It confirms that in the absence of comprehensive controlling efforts, the transmission rate has increased rapidly, resulting in a sizable epidemic.
Evaluation of the deterministic model
In order to evaluate the model in tracking the trend of the outbreak, we gathered new data which was not used in parameter estimation of the proposed EQIR model. We calculated the normalized mean square error (MSE) between predicted, and true number of daily confirmed infected and dead cases as follows:
where x(t), x̂(t) denote the true and x, predicted time series respectively and , indicates the average values of corresponding time series. Results of MSE calculated for T = 30 (from October 20,2020, to November 20,2020) are reported in [Table 3] for more details. According to this table, normalized MSE of the proposed model is very low, and the proposed model has good performance in prediction of new cases. It must be mentioned that errors by day are probably greater than the cumulative errors. This can be visually observed through cumulative and daily plots as well [Figure 6] and [Figure 7].
Forecasting the longterm trend of epidemic in Iran
We have adjusted parameters of the proposed dynamic model such as transmission rate, and isolation rate at different time points due to the changes occurred in the behavior of Iranian population resulting in the growth of infection rate in this country and caused several consecutive peaks in the outbreak trend. We provided the numerical results of data fitting of confirmed infected cases for the proposed EQIR deterministic and stochastic models from February 20,2020 to August 01,2021. [Figure 10] shows epidemic trajectories in Iran since the beginning of the outbreak. We also provided prediction of number of patients in Iran for the next 2 weeks in [Table 4]. According to the prediction results, if these conditions continued, the 5^{th} peak of the COVID19 pandemic will be continued in Iran till September 2021.  Figure 10: Numerical results of data fitting in deterministic and stochastic model in Iran; (a) cumulative infected cases, (b) daily infected cases
Click here to view 
 Table 4: Number of predicted patients and dead cases in Iran from 01 August 2021 to 14 August 2021
Click here to view 
Discussion   
In this study, we proposed a dynamic EQIR model of COVID19 spread. Statedependent structure of this model makes it possible to consider different descriptive parameters. We considered 14 parameters for this model among them up to10 parameters were fitted to real data using optimization methods. Fitting process can be simply done because of the possibility of parameter description in real world. Furthermore, the model can be easily adapted with changing the interventions of human societies due to reality of parameters. For instance, effects of isolation and screening schemes of exposed people are expressed by isolation rate (q) and different latent periods 𝜏_{E}, 𝜏_{Q} respectively. The effect of contact rate decreasing activities such as home quarantine and closure of public places, which leads to lower transmission rates, is considered by an exponentially decreasing function of time for transmission rate. In addition, fractional order of differential equations determines the effect of society's memory in the face of disease. This degree of freedom allows the model to fit better to data as we observed for South Korea whose fractional order of infected state was 0.8 where for Iran and the USA this value was near one. We obtained the relation between parameters of model and reproduction number and clearly justified it through parameters of model. Results of model demonstrated that besides increasing isolation rate, decreasing 𝜏_{E}, 𝜏_{Q} by means of identifying exposed people, help societies effectively control the outbreak.
We also implemented a stochastic model by adding a Gaussian noise to transmission rate of deterministic model. By considering such SDEs, we could insert model some of realworld fluctuations behavioral trend of disease. Results of this model can help us to estimate maximum intensive care units needed per day in the worst condition. In addition, the stochastic model shows that if the rate of spread has limited changes, outbreak size will be small.
We simulated short and longterm scenarios to predict outbreak size for both models and justified predictions due to changing in parameters in accordance to behavior of communities. Using different scenarios, this model can continuously monitor the trajectory of the disease and help control it.
Conclusion   
In this paper, we proposed a mathematical fractional order distributed delay dynamic model for spreading the trend of COVID19 in two deterministic and stochastic versions. It is a flexible model that can observe and predict outbreak of diseases such as COVID19, which can spread during their incubation period. In addition, multiple model parameters with physical descriptions allow us to apply these changes to the model at any time and forecast the future outbreak of COVID19. However, some higher degrees of contacts such as transmission within healthcare institutions and households are not considered in this model. Moreover, Due to unknown safety period obtained after infecting by COVID19 and the possibility of renewed infectious after immunity, this model must be developed by adding new states and routes for accurate prediction of disease outbreak.
Financial support and sponsorship
This work was supported in part by the ViceChancellery for Research and Technology, Isfahan University of Medical Sciences, under Grant 199067.
Conflicts of interest
There are no conflicts of interest.
Appendix A: Discrete deterministic model
For discretization of the model we use the backward Euler scheme.^{[33],[34]} Using this discretization method, N time intervals t_{0}< t_{0}<...< t_{N= T} with a uniform time step size are defined. For approximating the fractional derivative at time t = t_{n}, the shifted Grunewald formula^{[35]} is used as follows
Now, let λ = u^{α} β_{n} where β_{n} is discrete version of time varying function introduced in Equation 1. Let γ = u^{α} q where q is isolation rate. Let v = u^{α} δ where δ is considered as mortality or recovery rate. Thus, discrete representation of model is illustrated as follows.
where variables E_{n}, I_{n}, Q_{n} are numerical solutions of E(t_{n}), I(t_{n}), Q(t_{n})_{} and h (.) as explained in assumptions and defined as follows
Appendix B: Discrete stochastic model
For stochastic model we consider that β is time varying function following the random walk. Modeling time varying nature of the transmission rate as a random walk, implicitly declare that the amount of transmission rate in one day is not much different from previous day. Thus, transmission rate is considered as stochastic differential equation as follows:
where W_{t} denotes the wiener process^{[36]} and is σ fluctuation of transmission rate. If we apply Ito's lemma to the process, we get:
Now, let where and is assumed to be equal to u. ξ_{n} is Gaussian random variable N(0,1) Let where denotes the k^{th} Gaussian random variable N(0,1) Then, Model's differential equations can be written as follows (backward Euler approximation)
where variables E_{n}, I_{n}, Q_{n} are numerical solutions of E(t_{n}), I(t_{n}), Q(t_{n}) and h(.) as explained in assumptions and defined as follows
References   
1.  Li Q, Guan X, Wu P, Wang X, Zhou L, Tong Y, et al. Early transmission dynamics in Wuhan, China, of novel coronavirusinfected pneumonia. N Engl J Med 2020;382:1199207. 
2.  Kafieh R, Arian R, Saeedizadeh N, Amini Z, Serej ND, Minaee S, et al. Haghjooy Javanmard S. COVID19 in Iran: forecasting pandemic using deep learning. Computational and mathematical methods in medicine. 2021 Feb 25;2021. 
3.  Pirouz B, Shaffiee Haghshenas S, Shaffiee Haghshenas S, Piro P. Investigating a serious challenge in the sustainable development process: Analysis of confirmed cases of COVID19 (new type of coronavirus) through a binary classification using artificial intelligence and regression analysis. Sustainability 2020;12:2427. 
4.  Schoenberg FP, Hoffmann M, Harrigan RJ. A recursive point process model for infectious diseases. Ann Inst Stat Math 2019;71:127187. 
5.  You C, Deng Y, Hu W, Sun J, Lin Q, Zhou F, et al. Estimation of the timevarying reproduction number of COVID19 outbreak in China. Int J Hyg Environ Health 2020;228:113555. 
6.  Hethcote HW. Three basic epidemiological models. In: Applied Mathematical Ecology. Springer, Berlin, Heidelberg; 1989. p. 11944. 
7.  Giordano G, Blanchini F, Bruno R, Colaneri P, Di Filippo A, Di Matteo A, et al. Modelling the COVID19 epidemic and implementation of populationwide interventions in Italy. Nat Med 2020;26:85560. 
8.  López L, Rodo X. A modified SEIR model to predict the COVID19 outbreak in Spain and Italy: simulating control scenarios and multiscale epidemics. Results in Physics. 2021 Feb 1;21:103746. 
9.  Cooper I, Mondal A, Antonopoulos CG. A SIR model assumption for the spread of COVID19 in different communities. Chaos Solitons Fractals 2020;139:110057. 
10.  Law KB, Peariasamy KM, Gill BS, Singh S, Sundram BM, Rajendran K, et al. Tracking the early depleting transmission dynamics of COVID19 with a timevarying SIR model. Sci Rep 2020;10:21721. 
11.  Chen YC, Lu PE, Chang CS, Liu TH. A timedependent SIR model for COVID19 with undetectable infected persons. IEEE Trans Netw Sci Eng 2020;7:327994. 
12.  Marinov TT, Marinova RS. Dynamics of COVID19 using inverse problem for coefficient identification in SIR epidemic models. Chaos Solitons Fractals X 2020;5:100041. 
13.  Rǎdulescu A, Williams C, Cavanagh K. Management strategies in a SEIRtype model of COVID 19 community spread. Sci Rep 2020;10:21256. 
14.  Grimm V, Mengel F, Schmidt M. Extensions of the SEIR model for the analysis of tailored social distancing and tracing approaches to cope with COVID19. Sci Rep 2021;11:4214. 
15.  He S, Peng Y, Sun K. SEIR modeling of the COVID19 and its dynamics. Nonlinear Dyn 2020;101:166780. 
16.  Yang Z, Zeng Z, Wang K, Wong SS, Liang W, Zanin M, et al. Modified SEIR and AI prediction of the epidemics trend of COVID19 in China under public health interventions. J Thorac Dis 2020;12:16574. 
17.  Alrabaiah H, Arfan M, Shah K, Mahariq I, Ullah A. A comparative study of spreading of novel corona virus disease by ussing fractional order modified SEIR model. Alexandria Eng J 2021;60:57385. 
18.  Iwata K, Miyakoshi C. A simulation on potential secondary spread of novel coronavirus in an exported country using a stochastic epidemic SEIR model. J Clin Med 2020;9:E944. 
19.  Liu Z, Magal P, Seydi O, Webb G. A COVID19 epidemic model with latency period. Infect Dis Model 2020;5:32337. 
20.  Mahajan A, Sivadas NA, Solanki R. An epidemic model SIPHERD and its application for prediction of the spread of COVID19 infection in India. Chaos Solitons Fractals 2020;140:110156. 
21.  Korolev I. Identification and estimation of the SEIRD epidemic model for COVID19. J Econom 2021;220:6385. 
22.  Rajendran S, Jayagopal P. Accessing COVID19 epidemic outbreak in Tamilnadu and the impact of lockdown through epidemiological models and dynamic systems. Measurement (Lond) 2021;169:108432. 
23.  Lekone PE, Finkenstädt BF. Statistical inference in a stochastic epidemic SEIR model with control intervention: Ebola as a case study. Biometrics 2006;62:11707. 
24.  Britton T, Pardoux E, Ball F, Laredo C, Sirl D, Tran V.C. Stochastic epidemic models with inference. Berlin: Springer; 2019. p. 187. 
25.  Tajmirriahi M, Amini Z. Modeling of seizure and seizurefree EEG signals based on stochastic differential equations. Chaos Solitons Fractals 2021;150:111104. 
26.  Tajmirriahi M, Amini Z, Hamidi A, Zam A, Rabbani H. Modeling of retinal optical coherence tomography based on stochastic differential equations: Application to denoising. IEEE Trans Med Imaging 2021;40:212941. 
27.  Chen Y, Cheng J, Jiang Y, Liu K. A time delay dynamical model for outbreak of 2019nCoV and the parameter identification. J Inverse Ill Posed Probl 2020;28:24350. 
28.  Alofi A, Cao J, Elaiw A, AlMazrooei A. Delaydependent stability criterion of Caputo fractional neural networks with distributed delay. Discrete Dynamics in Nature and Society. 2014 Jan 1;2014. 
29.  Boatto S, Bonnet C, Cazelles B, Mazenc F. SIR model with time dependent infectivity parameter: Approximating the epidemic attractor and the importance of the initial phase. HalInria 2018. 
30.  Liu Z, Magal P, Seydi O, Webb G. Predicting the cumulative number of cases for the COVID19 epidemic in China from early data. Mathematical Biosciences and Engineering: MBE. 2020;17:304051. 
31.  Saeedian M, Khalighi M, AzimiTafreshi N, Jafari GR, Ausloos M. Memory effects on epidemic evolution: The susceptibleinfectedrecovered epidemic model. Phys Rev E 2017;95:022409. 
32.  Agarwal R, Belmekki M, Benchohra M. A survey on semilinear differential equations and inclusions involving RiemannLiouville fractional derivative. Adv Differ Equ 2009;2009:147. 
33.  Stetter HJ. Analysis of discretization methods for ordinary differential equations. BerlinHeidelbergNew York: Springer; 1973 May. 
34.  Liu J, Peng B, Zhang T. Effect of discretization on dynamical behavior of SEIR and SIR models with nonlinear incidence. Appl Math Lett 2015;39:606. 
35.  Dimitrov Y. Numerical approximations for fractional differential equations. arXiv preprint arXiv:1311.3935. 2013 Nov 15. 
36.  Durrett R, Durrett R. Brownian motion and martingales in analysis. In: Wadsworth Advanced Books & Software California. 1984. 
37.  Moré JJ. The LevenbergMarquardt algorithm: implementation and theory. InNumerical analysis: Springer, Berlin, Heidelberg: 1978. pp. 10516. 
38.  Endo A, van Leeuwen E, Baguelin M. Introduction to particle Markovchain Monte Carlo for disease dynamics modellers. Epidemics 2019;29:100363. 
39.  Andrieu C, Roberts GO. The pseudomarginal approach for efficient Monte Carlo computations. Ann Stat 2009;37:697725. 
40.  Kiskinov H, Zahariev A. On fractional systems with RiemannLiouville derivatives and distributed delayschoice of initial conditions, existence and uniqueness of the solutions. Eur Phys J Spec Top 2017;226:347387. 
41.  ÖZalp N, Demirci E. A fractional order SEIR model with vertical transmission. Math Comput Model 2011;54:16. 
42.  Yan P, Chowell G. Quantitative methods for investigating infectious disease outbreaks: Cham, Switzerland: Springer, 2019. 
43.  van den Driessche P. Reproduction numbers of infectious disease models. Infect Dis Model 2017;2:288303. 
44.  Farman M, Saleem MU, Ahmad A, Ahmad MO. Analysis and numerical solution of SEIR epidemic model of measles with noninteger time fractional derivatives by using Laplace Adomian Decomposition Method. Ain Shams Eng J 2018;9:33917. 
45.  
46.  
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6], [Figure 7], [Figure 8], [Figure 9], [Figure 10]
[Table 1], [Table 2], [Table 3], [Table 4]
