**Appendix B. Positiveness of the Delay Model (9)**

We show the positiveness of the solution of the delay model (9). Note that, as per our assumption, *J*(*t*) is a positive function. From (9a), we observe that, if *S*(*t*∗) = 0 at some time point *<sup>t</sup>*∗, then *dS dt* |*t*=*t*<sup>∗</sup> = 0. This implies that *S*(*t*) ≥ 0 for *t* > 0. From (3c) and (3d), we obtain that *R*(*t*) and *D*(*t*) are both increasing functions, and hence, *R*(*t*) and *D*(*t*) also remain positive for all *t*. Next, we prove that *I*(*t*) > 0. Integrating (9c) from 0 to *t*, we obtain

$$\begin{aligned} I(t) &= \, \, \, I(0) + \int\_0^t I(\eta)d\eta - \int\_0^t I(\eta - \tau)d\eta \\ &= \, \, \, I(0) + \int\_0^t I(\eta)d\eta - \int\_{-\tau}^{t-\tau} I(\eta)d\eta. \end{aligned}$$

Since, *J*(*t*) = 0 for *t* < 0, we can write:

$$I(t) := \left. I(0) + \int\_{t-\tau}^{t} I(\eta)d\eta \right| \ge \left. 0 \right.$$

Therefore, *I*(*t*) is positive for all *t*. This shows the positiveness of the solution of System (9).

**Appendix C. Gamma Distributions for Recovery and Death Rates**

**Figure A1.** Probability distribution of recovery (**a**) and death (**b**) as a function of time (in days) after the onset of infection. The red curves show the best-fit gamma distributions using the function "*fitdist(:,'gamma')*" in MATLAB.
