**4. Comparison of Models (3) and (9) and SIR (4)**

We compare the results obtained from the distributed model (3), the delay model (9), and the conventional SIR model (4). In [31], we showed that Model (3) can describe the epidemic spread more precisely, and it can exactly capture the recovery and death dynamics by using suitable distributed recovery and death rates. However, the main constraint in using the distributed model (3) is the availability of distributed recovery and death rates. Instead, the average duration of the disease is easier to determine. We show that, in the delay model (9), if we take *τ* as the mean of the recovery and death distribution functions involved in the distributed model (3), then the delay model (9) gives a close solution to the distributed model (3).

To estimate the recovery and death distributions in the model (3), we use the function *fitdist(:,gamma)* in MATLAB. This function is used to fit a vector of data *X* = (*x*1, *x*2, ··· , *xn*) by a gamma distribution of the form <sup>1</sup> *<sup>b</sup>a*Γ(*a*) *<sup>x</sup>a*−1*e*−*x*/*b*, where *<sup>a</sup>* and *<sup>b</sup>* are the shape and scale parameters. This function gives the maximum likelihood estimators of *a* and *b* for the gamma distribution, which are the solutions of the simultaneous equations

$$\log \mathfrak{a} - \Psi(\mathfrak{a}) = \log \left( \tilde{\mathcal{X}} / \left( \prod\_{i=1}^{n} \mathfrak{x}\_i \right)^{1/n} \right),$$

$$\hat{b} = \bar{\mathcal{X}} / \hat{\mathfrak{a}},$$

where *X*¯ is the sample mean of the data *X* and Ψ is the digamma function given by

$$\Psi(\mathbf{x}) = \Gamma'(\mathbf{x})/\Gamma(\mathbf{x}).$$

The function *fitdist(:,gamma)* estimates the shape and scale parameters with the 95% confidence interval. Using this statistical technique, we estimated the recovery and death distributions for the model (3) for the data used in [31] and given by the formulas

$$r(t) = \frac{p\_0}{b\_1^{a\_1} \Gamma(a\_1)} t^{a\_1 - 1} e^{-\frac{t}{b\_1}}, \qquad d(t) = \frac{(1 - p\_0)}{b\_2^{a\_2} \Gamma(a\_2)} t^{a\_2 - 1} e^{-\frac{t}{b\_2}}.$$

(Figure A1 in Appendix C), where the estimated parameter values are given in Table 1. Here, *p*<sup>0</sup> is the survival probability, and its estimated value, from the data, is *p*<sup>0</sup> = 0.97 [34].


**Table 1.** Estimated parameter values (gamma distributions).

The estimated average time to recovery is 17.85 days and to death is 13.19 days. The value of survival probability *p*<sup>0</sup> is 0.97 [34], that is, out of 100 infected individuals, 97 infected will recover. This estimate matches with most of the COVID-19 epidemic data from various countries [34,35]. Thus, we can take average disease duration *τ* as 17.7 days. The corresponding value in the SIR model is *<sup>r</sup>*<sup>0</sup> + *<sup>d</sup>*<sup>0</sup> ≈ 1/17.7 days<sup>−</sup>1.

It is essential to mention here that, instead of a gamma distribution, we can use the Erlang distribution *Ek* <sup>=</sup> *<sup>λ</sup>kxk*−1*e*−*λ<sup>x</sup>* (*k*−1)! , *<sup>x</sup>*, *<sup>λ</sup>* <sup>≥</sup> 0 and *<sup>k</sup>* <sup>∈</sup> <sup>N</sup>. Since, for lower values of *<sup>k</sup>*, e.g., *E*1, *E*<sup>2</sup> distributions, with a reasonable choice of *λ*, give a significant probability of recovery or death at the beginning of infection, we should look for an Erlang distribution with higher values of *k*.

We found that the Erlang-8 distribution

$$E\_8 = \frac{\lambda^8 x^{\mathcal{T}} e^{-\lambda x}}{(\mathcal{T})!},$$

with *λ* = 0.4545, and the Erlang-6 distribution

$$E\_6 = \frac{\lambda^6 x^5 e^{-\lambda x}}{(5)!} \text{ }$$

with *λ* = 0.4548, closely match with the gamma distributions for recovery and death rates respectively, as estimated above. The estimated values of the shape parameters *a*<sup>1</sup> and *a*<sup>2</sup> associated with the gamma distributions can be approximated as *a*<sup>1</sup> ≈ 8 and *a*<sup>2</sup> ≈ 6 for the recovery and death rate distributions, respectively. Thus, one can use the *E*<sup>8</sup> and *E*<sup>6</sup> distributions instead of the gamma distributions for recovery and death rates. However, the result will be similar in both cases, and we continue our numerical simulation with the gamma distributions.

Though the parameters of the three models correspond to each other, the distributed model (3) and the delay model (9) both give far different dynamics as compared to the SIR model (4), whereas the distributed model (3) and the delay model (9) give reasonably close dynamics to each other (Figure 4).

**Figure 4.** Comparison of the solutions of the distributed model (3) (blue curve), the delay model (9) (green curve), and the conventional SIR model (4) (red curve). (**a**) Susceptible *S*(*t*), (**b**) infected *I*(*t*), (**c**) recovered *R*(*t*), and (**d**) dead *D*(*t*). The values of parameters: *N* = 105, *I*<sup>0</sup> = 1, *β* = 0.2, for the SIR model *r*<sup>0</sup> + *d*<sup>0</sup> = 1/17.7; for the delay model *τ* = 17.7; for the distributed model (3): *a*<sup>1</sup> = 8.06275, *b*<sup>1</sup> = 2.21407, *a*<sup>2</sup> = 6.00014, *b*<sup>2</sup> = 2.19887, *p*<sup>0</sup> = 0.97.

The comparison of the final size of epidemic *Sf* , maximal number of infected *Im*, and the time to the maximal number of infected *tm* of the distributed model (3) with the gamma distribution, the delay model (9), and the SIR model (4) is shown in Figure 5 for different values of the transmission rate *β*. As before, the maximal numbers of infected individuals *Im* in the models (3) and (9) are sufficiently close, but they are much higher than for the SIR model (4). The times to maximum infected *tm* in the models (3) and (9) are reasonably close, but less than for the SIR model (4). The final size of the epidemic *Sf* is more or less the same for all three models. Similar properties are observed if the gamma distribution is replaced by the Erlang distribution with the corresponding parameters.

This result indicates the relevance of the delay model. If we do not have sufficient individual-level data to estimate the recovery and death distributions, but we have an approximate value of disease duration, then we can describe the epidemic progression in a sufficiently precise way.

**Figure 5.** Comparison of the maximal number of infected individuals *Im* (**a**), the time to reach the maximal number *tm* (**b**), and the final size of the epidemic *Sf* (**c**) for the distributed model (3) (blue curves), the delay model (9) (green curves), and the SIR model (4) (red curves) for different values of *β*. The values of parameters: *N* = 105, *I*<sup>0</sup> = 1, for the SIR model *r*<sup>0</sup> + *d*<sup>0</sup> = 1/17.7; for the delay model *τ* = 17.7; for the distributed model (3): *a*<sup>1</sup> = 8.06275, *b*<sup>1</sup> = 2.21407, *a*<sup>2</sup> = 6.00014, *b*<sup>2</sup> = 2.19887, *p*<sup>0</sup> = 0.97.
