Introduction

The outbreak of coronavirus disease 2019 (COVID-19), reported early in Wuhan (China)1 and spread around the world, is creating dramatic and daily changes with profound impacts worldwide. As a consequence the outbreak was declared a pandemic by the World Health Organization (WHO) in March 20202, and by the end of 2020, COVID-19 has infected about 79.2 millions of people in the world, with an approximate cumulative global mortality of 3.2%2. To limit the impact of this deadly virus, a rapid and widespread vaccination of the population is now in place. However, it is established that vaccine are not 100% effective to stop the transmission or infection of COVID-19. In addition, huge numbers of global SARS-CoV-2 infections have led to the emergence of variants, notably Alpha (B.1.1.7 UK), Beta (B.1.351 S. Africa), Gamma (P.1 Brazil), Epsilon (B.1.429 California), Iota (B.1.526 New York), Delta and Kappa (B.1.617.2 and B.1.617.1 India) which make the situation more challenging. In this circumstance to get a complete feature of COVID-19, it is essential to fully understand the key (incubation, recovery and decease) periods.

We already successfully estimated the incubation period of COVID-19 in Canada3. The previous model3 simply allowed for the calculation of the incubation period, while the current model allows for the calculation of all key periods, incubation, recovery and decease. In the present context, we focus on the recovery and decease periods and their correlation with the incubation period. In the current framework, we define the recovery period as the time from the contraction of the coronavirus to recovery, i.e., the incubation period plus the onset time from the symptom to recovery; the latter is the same as the viral shedding of SARS-CoV-2. We describe the decease period in the same way as the recovery period. Understanding the recovery period of disease is very useful information in the struggle against the disease. If the incidence of a disease is remarkably high and the recovery period of the disease is also high then the prevalence of the disease in the country is likely to increase which in turn puts extra health, economic and social burden on this country. Understanding the recovery period of the disease will help governments to plan proper strategies to counter the disease and to organize the requirements such as hospitals, doctors, medical staffs, medical equipment’s, etc. It will also help to implement different social and economic policies which will be essential to fight the disease.

There are several statistical studies4,5,6,7,8,9,10,11,12, based on various samples of patients such as severe, non-severe, ICU, non-ICU, large size, small size, meta-analysis, estimated the recovery time of the current pandemic. In addition to those statistical approaches, there are numerous analytical and computational studies based on mathematical models, involving Ordinary Differential Equations (ODE)13,14,15,16,17,18,19,20,21,22,23 as well as Delay Differential Equations (DDE)24,25,26,27,28,29, to calculate the basic reproduction number \(R_0\) and understand the underlying dynamics of the epidemic. Researchers usually consider single-delay models, occasionally two delays.

To the best of our knowledge, we demonstrate for the first time a substantial compartment based model, with a total 830 partitions, in order to estimate the key (incubation, recovery, decease periods) periods of COVID-19 as well as the bivariate distribution of incubation and recovery periods, and the bivariate distribution of incubation and decease periods. This will be achieved using publicly available database30 of the total number of corona-positive cases, recovery and death toll. There is no scope to verify the database which is the only limitation of the present study. Using the novel model, demonstrated here, we divide the publicly available database into thousands of groups, and these separated classes are the key source for estimating all the key periods. This approach is free from any special type of samples in order to produce the distributions of those periods; it only involves large scale computations for estimating about thousand model parameters. After a single calculation of this method, we can generate the current distributions as well as previous distributions of those periods. In the statistical based approaches, it is usually difficult to consider large incubation, recovery and decease periods if the sample size is small. However, in our approach, we can go well beyond 14 days, the maximum incubation period that we have set in this paper, and beyond the interval 2 weeks to 6 weeks, the range of recovery as well as decease periods that we have considered in the current computations. As of May 23, 2021, the World Health Organization (WHO) had confirmed a total of 1,359,180 cases of COVID-19 in Canada, including 25,231 deaths2. As of May 23, 2021 there are five provinces, out of eleven provinces where we observe significant effect of COVID-19, in Canada with death toll more than 1000 (Fig. 1a), and the recovery and death rates are respectively 96.1% and 0.7% (Fig. 1c). During the first wave of COVID-19 in Canada, January 22, 2020 to July 16, 2020, the recovery and death rates were respectively 66.5% and 8.1% (Fig. 1b). Here, we assume that the recovery and decease periods of COVID-19 remain unchanged i.e., these periods during the first wave and the present time are almost identical, and under this assumption we merely consider the database of first wave, the period before vaccination, for the calculation. In the present context, we do not consider the patient’s gender or age due to lack of the required data.

There are various studies on recovery period, and no result is reported (to the best of our knowledge) on bivariate distributions as mentioned above. The key periods may depend on age31 (median-age/country), hard immunity, public health system, corona testing capacities, daily corona cases, etc. For a better estimation of the key periods for a particular region, we need to study local patients. Data collection is a bottleneck in studying those key periods for COVID-19 or other infectious diseases using clinical survey, and we need a sample of large size for bivariate analysis. However, key periods can easily be estimated using the approach we propose here, the publicly available database along with optimization algorithms.

Figure 1
figure 1

COVID-19 pandemic in Canada: (a) As of May 23, 2021 Canadian out break at-a-glance; the purple, white and cyan digits are represented the number of recovered individuals, death toll and active cases, respectively. (b) Percentage of recovered (purple), deaths (white) and active cases (cyan) in Canada during the first wave, January 22, 2020 to July 16, 2020. (c) Percentage of recovered (purple), deaths (white) and active cases (cyan) in Canada as of May 23, 2021. The image has been generated using Microsoft Paint—Windows 10. Model calculation for Canada during the first wave, January 22, 2020 to July 16, 2020: (d) Estimation of the number of infected individuals. (e) Estimation of the total number of coronavirus cases compared to the available data30. (f) Estimation of the total number of recovered compared to the available data30. (g) Estimation of the total number of deaths compared to the available data30.

Table 1 Comparison of several studies (including the present work) for infectious period along with sample size, mean/median and ranges.
Figure 2
figure 2

Distribution of the recovery and decease periods: Results based on the total recovered cases of the first 177 days during the pandemic in Canada starting from January 22, 2020 i.e., cumulative data as of July 16, 2020. (a) Splitting values of recovered individuals as a function of incubation and recovery periods. (b) Probability density function of the gamma distribution \(\Gamma (\zeta , K, \theta )\) with \(K = 18.62067\) and \(\theta = 1.18892\). The blue bars indicate the densities obtained from the model calculation. (c) Probability density function of the bimodal gamma distribution \(0.9365\, \Gamma (\zeta , K_1, \theta _1) + 0.0635\, \Gamma (\zeta , K_2, \theta _2)\) with \(K_1 = 34.55447\), \(\theta _1 = 0.60847\), \(K_2 = 226.40545\) and \(\theta _2 = 0.17171\). The blue bars indicate the densities obtained from the model calculation. (d) Percentile curves for unimodal and bimodal gamma distributions. (e) Splitting values of the deaths as a function of incubation and decease periods. (f) Probability density function of the gamma distribution \(\Gamma (\eta , K, \theta )\) with \(K = 21.33660\) and \(\theta = 1.03174\). The blue bars indicate the densities obtained from the model calculation. (g) Probability density function of the bimodal gamma distribution \(0.9508\, \Gamma (\eta , K_1, \theta _1) + 0.0492\, \Gamma (\eta , K_2, \theta _2)\) with \(K_1 = 35.00855\), \(\theta _1 = 0.60511\), \(K_2 = 186.11379\) and \(\theta _2 = 0.20636\). The blue bars indicate the densities obtained from the model calculation. (h) Percentile curves for unimodal and bimodal gamma distributions.

Figure 3
figure 3

Bivariate distribution of the incubation and recovery periods: (a) Histogram of the estimated data \(R_{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\) using the model. (b) Fitted bivariate normal distribution. (c) Two-dimensional display of (b); the red region is the highly probable domain for recovery, and x (6.43, 21.91) denotes the center of the region. (d) Fitted nonparametric density estimate with wide 33%; two peaks show that there are two distinguishable high probable regions. (e) Two-dimensional display of (d); two x represent the centers of two high probable regions (3.49, 20.52) and (8.38, 20.35).

Figure 4
figure 4

Bivariate distribution of the incubation and decease periods: (a) Histogram of the estimated data \(D_{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\) using the model. (b) Fitted bivariate normal distribution. (c) Two-dimensional display of (b); the red region is the highly probable domain for decease, and x (6.56, 21.64) denotes the center of the region. (d) Fitted nonparametric density estimate with wide 40%; two peaks show that there are two distinguishable high probable regions. (e) Two-dimensional display of (d); the x represents the center of the high probable red region (8.17, 21.86).

Results

The proposed model assists us to generate new refined recovery and death toll database, \(R_{ik}\) and \(D_{ik}\), by dividing the total recovered individuals and the total number of deaths as of July 16, 2020 into myriad of groups. The new database is the key source for studying all kinds of distributions, reported in the article.

Validation of the proposed model

After estimating the model parameters with sufficiently small values of error functions, we obtain a good agreement (Fig. 1e–g) between the calculated values of the model variables such as total corona-positive cases, number of recovered individuals, etc. and the available data30. The population of the infected group gradually increased until end of April 2020, and thereafter slowed down (Fig. 1d).

Univariate distributions

The groups of recovered individuals \(R_{ik}\), \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\), corresponding to the incubation period (in days) \(\tau _i\), \(1 \le \tau _i \le 14\), and recovery period (in days) \(\zeta _{k}\), \(14 \le \zeta _k \le 42\) can be represented in a matrix form (Fig. 2a). We use the data set \(\sum _{i = 1}^{14} R_{ik}\) for \(\zeta _k = 14, 15, \ldots , 42\) to obtain the frequency distribution for recovery period and the corresponding fitted gamma distributions, unimodal (Fig. 2b) \(\Gamma (\zeta , K_r, \theta _r)\) and bimodal (Fig. 2c) \(0.9365\, \Gamma (\zeta , K_{r1}, \theta _{r1}) + 0.0635\, \Gamma (\zeta , K_{r2}, \theta _{r2})\). Here, the variable \(\zeta\) indicates the recovery period and the parameters \(K_r = 18.62067\), \(\theta _r = 1.18892\), \(K_{r1} = 34.55447\), \(\theta _{r1} = 0.60847\), \(K_{r2} = 226.40545\) and \(\theta _{r2} = 0.17171\) with statistical p value less than 0.01. The mean recovery period we obtain using an unimodal gamma distribution is 22.14 days (95% CI 22.00–22.27); the median of the recovery period is 21.74 days (95% CI 21.61–21.87); the 90th percentile is 28.91 days (95% CI 28.71–29.13); the 95th percentile is 31.20 days (95% CI 30.95–31.45). For a better estimation, we use a bimodal distribution, a linear combination of \(\Gamma (\zeta , K_{r1}, \theta _{r1})\) and \(\Gamma (\zeta , K_{r2}, \theta _{r2})\). The mean of \(\Gamma (\zeta , K_{r1}, \theta _{r1})\) and \(\Gamma (\zeta , K_{r2}, \theta _{r2})\) are 21.02 days (95% CI 20.92–21.12 ) and 38.88 days (95% CI 38.61–39.15), respectively. The percentile curves of unimodal and bimodal gamma distributions show (Fig. 2d) that the median of unimodal and bimodal are the same, although there are slight differences other than the median.

The death toll groups \(D_{ik}\), \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\), and corresponding incubation period (in days)\(\tau _i\), \(1 \le \tau _i \le 14\), and decease period (in days) \(\eta _k\), \(14 \le \eta _k \le 42\) can be represented as a matrix (Fig. 2e). We use the data set \(\sum _{i = 1}^{14} D_{ik}\) for \(\eta _k = 14, 15, \ldots , 42\) to obtain the frequency distribution for decease period and corresponding fitted gamma distributions, unimodal (Fig. 2f) \(\Gamma (\eta , K_d, \theta _d)\) and bimodal (Fig. 2g) \(0.9508\, \Gamma (\eta , K_{d1}, \theta _{d1}) + 0.0492\, \Gamma (\eta , K_{d2}, \theta _{d2})\). Here, the variable \(\eta\) indicates the decease period and the parameters \(K_d = 21.33660\), \(\theta _d = 1.03174\), \(K_{d1} = 35.00855\), \(\theta _{d1} = 0.60511\), \(K_{d2} = 186.11379\) and \(\theta _{d2} = 0.20636\) with statistical p value less than 0.01 for \(\Gamma (\eta , K_d, \theta _d)\), \(\Gamma (\eta , K_{d1}, \theta _{d1})\) and equal to 0.18 for \(\Gamma (\eta , K_{d2}, \theta _{d2})\) . The mean decease period we obtain using an unimodal gamma distribution is 22.01 days (95% CI 21.64–22.39); the median of the decease period is 21.67 days (95% CI 21.31–22.04); the 90th percentile is 28.30 days (95% CI 27.72–28.89); the 95th percentile is 30.39 days (95% CI 29.71–31.10). For better estimation, we use a bimodal distribution, a linear combination of \(\Gamma (\eta , K_{d1}, \theta _{d1})\) and \(\Gamma (\eta , K_{d2}, \theta _{d2})\). The mean of \(\Gamma (\eta , K_{d1}, \theta _{d1})\) and \(\Gamma (\eta , K_{d2}, \theta _{d2})\) are 21.18 days (95% CI 20.90–21.47) and 38.41 days (95% CI 37.41–39.40), respectively. The percentile curves show (Fig. 2h) that the percentiles of unimodal and bimodal distributions are almost the same.

Bivariate distributions

To analyze the bivariate distribution, we use the software Statgraphics32, based on the statistical package R. Using the elements \(R_{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\), we obtain a bivariate histogram (Fig. 3a) for the incubation and recovery periods. There are two peaks at the points (3, 19), i.e., for \(\tau _i =3\) and \(\zeta _k = 19\), and (8, 20), i.e., for \(\tau _i =8\) and \(\zeta _k = 20\), corresponding to the high densities of recovered individuals. We estimate the histogram using a bivariate normal distribution \({\mathscr {N}}(m_{\tau }^{(r)}, m_{\zeta }, \sigma _{\tau }^{(r)}, \sigma _{\zeta }, \rho _{\tau \zeta }\)) (Fig. 3b) where the variables \(\tau\) and \(\zeta\) represent the incubation and recovery periods, respectively. The mean \(m_{\tau }^{(r)}\) and standard deviation \(\sigma _{\tau }^{(r)}\) of the incubation period are 6.43 (95% CI 6.27–6.59) and 3.06 (95% CI 2.96–3.18), respectively; the mean \(m_{\zeta }\) and standard deviation \(\sigma _{\zeta }\) of the recovery period are 21.91 (95% CI 21.63–22.18) and 5.33 (95% CI 5.14–5.53), respectively; the correlation between incubation and recovery periods \(\rho _{\tau \zeta }\) is − 0.11. The two dimensional representation of the bivariate normal distribution (Fig. 3c) shows that the highly probable recovery region (red in the figure) is a nested domain of \(\tau = 6.43\) and \(\zeta = 21.91\). To precisely analyze the highly probable region, we estimate the histogram (Fig. 3a) using a nonparametric density function with a width of 33%, low and high percentage give a more local and global estimation, respectively, and we obtain a distribution with two peaks (Fig. 3d). Two distinguishable peaks indicate that there are two separate highly probable regions surrounding the points \(\tau = 3.49\), \(\zeta = 20.52\) and \(\tau = 8.38\), \(\zeta = 20.35\) (Fig. 3e). The bivariate mixture distribution analysis shows that we can estimate the histogram of \(R_{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\) using a combination of two bivariate normal distributions, \(0.94 {\mathscr {N}}(m_{\tau }^{(r1)}, m_{\zeta }^{(1)}, \sigma _{\tau }^{(r1)}, \sigma _{\zeta }^{(1)}, \rho _{\tau \zeta }^{(1)}) + 0.06 {\mathscr {N}}(m_{\tau }^{(r2)}, m_{\zeta }^{(2)}, \sigma _{\tau }^{(r2)}, \sigma _{\zeta }^{(2)}, \rho _{\tau \zeta }^{(2)})\) where the superscript 1 (resp. 2) represents the parameters for the first (resp. second) component. The parameters of the first component are \(m_{\tau }^{(r1)} = 6.43, m_{\zeta }^{(1)} = 20.83, \sigma _{\tau }^{(r1)} = 3.05, \sigma _{\zeta }^{(1)} = 3.36, \rho _{\tau \zeta }^{(1)} = -0.16\), and those of the second component are \(m_{\tau }^{(r2)} = 6.42, m_{\zeta }^{(2)} = 37.83, \sigma _{\tau }^{(r2)} = 3.25, \sigma _{\zeta }^{(2)} = 3.45, \rho _{\tau \zeta }^{(2)} = -0.47\).

Using the elements \(D_{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\), we obtain a bivariate histogram (Fig. 4a) for the incubation and decease periods. There are two peaks at the points (3, 22), i.e., for \(\tau _i =3\) and \(\eta _k = 22\), and (9, 23), i.e., for \(\tau _i =9\) and \(\eta _k = 23\), corresponding to the high densities of deaths. We estimate the histogram using a bivariate normal distribution \({\mathscr {N}}(m_{\tau }^{(d)}, m_{\eta }, \sigma _{\tau }^{(d)}, \sigma _{\eta }, \rho _{\tau \eta }\)) (Fig. 4b) where the variables \(\tau\) and \(\eta\) represent the incubation and decease periods, respectively. The mean \(m_{\tau }^{(d)}\) and standard deviation \(\sigma _{\tau }^{(d)}\) of the incubation period are 6.56 (95% CI 6.36–6.76) and 3.00 (95% CI 2.86–3.15), respectively; the mean \(m_{\eta }\) and standard deviation \(\sigma _{\eta }\) of the decease period are 21.64 (95% CI 21.33–21.94) and 4.43 (95% CI 4.23–4.65), respectively; the correlation between incubation and decease periods \(\rho _{\tau \eta }\) is − 0.008. The two dimensional representation of the bivariate normal distribution (Fig. 4c) shows that the highly probable decease region (red in the figure) is a nested domain of \(\tau = 6.56\) and \(\eta = 21.64\). To precisely analyze the highly probable regions, we estimate the histogram (Fig. 4a) using a nonparametric density function with a width of 40%, low and high percentage give a more local and global estimation, respectively, and obtain a distribution with two peaks (Fig. 4d), one in the high probability region (red in figure) and another one in the second high probability region (yellow in figure). The highly probable region is surrounding the point \(\tau = 8.17\), \(\eta = 21.86\) (Fig. 4e). The bivariate mixture distribution analysis shows that we can estimate the histogram of \(D_{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\) using a combination of two bivariate normal distributions, \(0.97 {\mathscr {N}}(m_{\tau }^{(d1)}, m_{\eta }^{(1)}, \sigma _{\tau }^{(d1)}, \sigma _{\eta }^{(1)}, \rho _{\tau \eta }^{(1)}) + 0.03 {\mathscr {N}}(m_{\tau }^{(d2)}, m_{\eta }^{(2)}, \sigma _{\tau }^{(d2)}, \sigma _{\eta }^{(2)}, \rho _{\tau \eta }^{(2)})\) where the superscript 1 (resp. 2) represents the parameters for first (resp. second) component. The parameters of the first component are \(m_{\tau }^{(d1)} = 6.58, m_{\eta }^{(1)} = 21.19, \sigma _{\tau }^{(d1)} = 2.98, \sigma _{\eta }^{(1)} = 3.49, \rho _{\tau \eta }^{(1)} = 0.03\), and those of the second component are \(m_{\tau }^{(d2)} = 5.81, m_{\eta }^{(2)} = 38.61, \sigma _{\tau }^{(d2)} = 3.33, \sigma _{\eta }^{(2)} = 2.15, \rho _{\tau \eta }^{(2)} = -0.29\).

Onset time from symptom to recovery

Using the fact that \(\tau + \theta = \zeta\) and the property of expectation \(E(T + \Theta ) = E(T) + E(\Theta )\), we calculate the mean Onset Time from Symptom to Recovery (OTSR) \(E(\Theta )\) (Table 1), where \(\theta\) is the variable corresponding to \(\theta _{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\); T and \(\Theta\) are the random variables corresponding the incubation period and OTSR, respectively. There is a good agreement between the calculated values, mean of OTSR, short OTSR and long OTSR, with the reported works (Table 1) of earlier studies. However, these calculated values do not show excellent concordance with some other studies, because we consider all recovery cases, mild to moderate, severe, hospitalized (ICU, non-ICU), non hospitalized, in Canada. For example, Voinsky et al.4 reported a study with a sample of 5769 patients, not including severe COVID-19 cases. In fact, they mentioned that severe cases were reported to be discharged from the hospital on average 8 days longer than mild to moderate patients requiring hospitalization.

Discussion

In the present context, we estimate the recovery as well as decease periods using a novel compartment based model and publicly available database. Here, we consider a maximum length of the incubation period of 14 days, and the ranges of the recovery and decease periods are from 2 to 6 weeks. However, in our method, we can go well beyond all those ranges; the longer ranges simply require a long computational time. Notice that our method could apply the proposed model to estimate key periods for any infectious disease, as along as similar data are available. The proposed model is not a prediction model, and hence does not depend on any database. The subsection ‘Validation of the proposed model’ in the ‘Results’ section is presented only to justify the validation of our model with COVID-19 database of Canada. Calculating the incubation and recovery periods for other counties is naturally possible using the proposed model.

The multi-group database \(R_{ik}\), \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\), generated from the model, is the key source to compute all types of distribution of the recovery period, univariate, bimodal and bivariate. The bimodal gamma distribution of the recovery period, \(0.9365\, \Gamma (\zeta , K_{r1}, \theta _{r1}) + 0.0635\, \Gamma (\zeta , K_{r2}, \theta _{r2})\), demonstrates that the recovery period of 93.65% recovered individuals obeys the distribution \(\Gamma (\zeta , K_{r1}, \theta _{r1})\), and that of 6.35% recovered individuals obeys the distribution \(\Gamma (\zeta , K_{r2}, \theta _{r2})\). Thus, there are two groups of recovered individuals with short recovery period, 21.02 days (on average), and long recovery period, 38.88 days (on average). The characteristics of those two groups may depend on age, underlying health condition, immunity, etc. The database of numerous groups \(D_{ik}\), \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\), generated from the model, is the key source to compute all types of distribution of the decease period, univariate, bimodal and bivariate. The bimodal gamma distribution of the decease period, \(0.9508\, \Gamma (\zeta , K_{d1}, \theta _{d1}) + 0.0492\, \Gamma (\zeta , K_{d2}, \theta _{d2})\), demonstrates that the decease period of 95.08% deaths obeys the distribution \(\Gamma (\zeta , K_{d1}, \theta _{d1})\), and that of 4.92% deaths obeys the distribution \(\Gamma (\zeta , K_{d2}, \theta _{d2})\). Thus, there are two groups of deaths with short decease period, 21.18 days (on average), and long decease period, 38.41 days (on average). The characteristics of those two groups may depend on age, underlying health condition, immunity, etc. The calculated results employing the proposed model show that the recovery and decease periods are the same. It seems that the survival period of the coronavirus is the same as that of human, in the form of immunity.

The bivariate normal distribution of incubation and recovery periods indicates a recovery window of \(4.82 \le \tau \le 8.49\) and \(19.27 \le \zeta \le 25.72\) as the highly probable domain for recovery. The bivariate normal distribution of incubation and decease periods indicates a decease window of \(4.55 \le \tau \le 8.45\) and \(19.35 \le \eta \le 24.85\) as the highly probable domain of deaths. The study shows that the recovery and decease windows almost coincide within these key periods. To determine precisely the recovery as well as the decease windows, we use nonparametric distributions. Under the nonparametric analysis we identify two recovery windows, \(2.27 \le \tau \le 4.38\), \(18.41 \le \zeta \le 22.79\) and \(6.42 \le \tau \le 9.63\), \(17.81 \le \zeta \le 23.93\) , and one decease window, \(6.34 \le \tau \le 9.41\), \(20.17 \le \eta \le 23.69\). Nonparametric analysis provides some discrepancy between the recovery and decease windows.

The bivariate mixed distribution, \(0.94 {\mathscr {N}}(m_{\tau }^{(r1)}, m_{\zeta }^{(1)}, \sigma _{\tau }^{(r1)}, \sigma _{\zeta }^{(1)}, \rho _{\tau \zeta }^{(1)}) + 0.06 {\mathscr {N}}(m_{\tau }^{(r2)}, m_{\zeta }^{(2)}, \sigma _{\tau }^{(r2)}, \sigma _{\zeta }^{(2)}, \rho _{\tau \zeta }^{(2)})\), of the incubation and recovery periods demonstrates that 94% recovered individuals obey \({\mathscr {N}}(m_{\tau }^{(r1)}, m_{\zeta }^{(1)}, \sigma _{\tau }^{(r1)}, \sigma _{\zeta }^{(1)}, \rho _{\tau \zeta }^{(1)})\), the bivariate normal distribution, with recovery window \(m_{\tau }^{(r1)} - \sigma _{\tau }^{(r1)} \le \tau \le m_{\tau }^{(r1)} + \sigma _{\tau }^{(r1)}\), \(m_{\zeta }^{(1)} - \sigma _{\zeta }^{(1)} \le \zeta \le m_{\zeta }^{(1)} + \sigma _{\zeta }^{(1)}\) and 6% recovered individuals obey the bivariate normal distribution \({\mathscr {N}}(m_{\tau }^{(r2)}, m_{\zeta }^{(2)}, \sigma _{\tau }^{(r2)}, \sigma _{\zeta }^{(2)}, \rho _{\tau \zeta }^{(2)})\) with recovery window \(m_{\tau }^{(r2)} - \sigma _{\tau }^{(r2)} \le \tau \le m_{\tau }^{(r2)} + \sigma _{\tau }^{(r2)}\), \(m_{\zeta }^{(2)} - \sigma _{\zeta }^{(2)} \le \zeta \le m_{\zeta }^{(2)} + \sigma _{\zeta }^{(2)}\). The bivariate mixed distribution, \(0.97 {\mathscr {N}}(m_{\tau }^{(d1)}, m_{\eta }^{(1)}, \sigma _{\tau }^{(d1)}, \sigma _{\eta }^{(1)}, \rho _{\tau \eta }^{(1)}) + 0.03 {\mathscr {N}}(m_{\tau }^{(d2)}, m_{\eta }^{(2)}, \sigma _{\tau }^{(d2)}, \sigma _{\eta }^{(2)}, \rho _{\tau \eta }^{(2)})\), of the incubation and decease periods demonstrates that 97% deaths obey the bivariate normal distribution \({\mathscr {N}}(m_{\tau }^{(d1)}, m_{\eta }^{(1)}, \sigma _{\tau }^{(d1)}, \sigma _{\eta }^{(1)}, \rho _{\tau \eta }^{(1)})\) with decease window \(m_{\tau }^{(d1)} - \sigma _{\tau }^{(d1)} \le \tau \le m_{\tau }^{(d1)} + \sigma _{\tau }^{(d1)}\), \(m_{\eta }^{(1)} - \sigma _{\eta }^{(1)} \le \eta \le m_{\eta }^{(1)} + \sigma _{\eta }^{(1)}\) and 3% deaths obey the bivariate normal distribution \({\mathscr {N}}(m_{\tau }^{(d2)}, m_{\eta }^{(2)}, \sigma _{\tau }^{(d2)}, \sigma _{\eta }^{(2)}, \rho _{\tau \eta }^{(2)})\) with decease window \(m_{\tau }^{(d2)} - \sigma _{\tau }^{(d2)} \le \tau \le m_{\tau }^{(d2)} + \sigma _{\tau }^{(d2)}\), \(m_{\eta }^{(2)} - \sigma _{\eta }^{(2)} \le \eta \le m_{\eta }^{(2)} + \sigma _{\eta }^{(2)}\).

Conclusions

We have developed a novel compartment based model to divide the publicly available database of total confirmed cases, recovered cases, and number of deaths into numerous subgroups to obtain the distributions of the recovered and decease periods. The outcomes of this study can be divided into three categories; these are univariate, univariate (bimodal) and bivariate distributions. We obtain mean recovery and decease periods from the univariate distribution. We observe two groups of recovered individuals as well as deaths: a short recovery (decease) period and a long recovery (decease) period from the univariate (bimodal) distribution. From the bivariate analysis, we investigate the correlation between the incubation and recovery periods as well as the correlation between incubation and decease periods. The model itself and the procedure to solve it, are the core of this work, and it can be applied to any infectious disease in any region. We obtain the distributions of the key periods from the population, considering all types of cases (non-hospitalized, non-ICU, ICU) of recovered individuals and deaths, which is naturally better than any sample-dependent result. In this approach, we do not need any clinical survey; the publicly available data on confirmed cases, recovery and death toll, are sufficient to analyze the univariate and bivariate distributions. The current model can be extended to study age-based key periods, but for this purpose we need an age dependent database. The monotonic iteration scheme, introduced for better estimation, can be applied to numerical analysis problems.

Methods

In this section, we introduce a compartment based infectious disease model including a large number of partitions, Lockdown, Susceptible, Removed, Infected, fourteen compartments of Confirmed cases, hundreds of compartments of Recovered and Deaths. The model is constructed as a set of coupled delay differential equations involving few thousands of variables and parameters, and will be used, not as a prediction tool, but (1) for constructing the myriad groups of recovered individuals and death tools and (2) estimating accurately the recovery and decease periods. This model will however have to be parameterized and validated using existing data, in order to justify its accuracy and its application in the proposed methodology.

The model

Modeling the spread of pandemics is an essential tool for projecting its outcome. By estimating important epidemiological parameters using the available database and optimization techniques, we can make predictions of different intervention scenarios. Compartment based model, where the population of a region is distributed into several population groups, such as susceptible, infected, total cases, etc., is a simple but useful tool to demonstrate the panorama of an epidemic.

The proposed model is an extension of our previous work3, including a very large number of compartments of recovered and deaths individuals; the schematic diagram of the model is presented in Fig. 5a. The following are the underlying principles of the present model.

  • The total population is constant (neglecting the migrations, births and unrelated deaths) and initially every individual is assumed susceptible to contract the disease.

  • The disease is spread through the direct (face-to-face meeting) or indirect (through air current, common used or delivery items like door handles, grocery products) contact of susceptible individuals with the infected individuals.

  • The quarantined area or the compartment for corona cases contains only members of the infected population who are tested corona-positive.

  • The virus kills a part of the people it infects; the survivors represent the recovered group.

  • There is a non-pharmaceutical policy (stay at home), commonly known as lockdown, to stop the spread of the disease.

  • The group of asymptomatic patients is a part of infected individuals, and the never-tested recovered asymptomatic patients can be removed from the infected group. If an asymptomatic patient dies, it is counted after investigation.

Based on the above principles, we consider the following compartments:

  • Lockdown (insusceptible) (L): the group of persons who are keeping themselves safe.

  • Susceptible (S): the group of individuals who can be infected.

  • Infected (I): the group of people who are spreading the contiguous disease.

  • Removed (V): the group of recovered asymptomatic patients without testing.

  • Confirmed cases (C): the group of individuals who tested corona-positive.

  • Recovered (R): the group of recovered individuals who tested corona-positive.

  • Deaths (D): the group of deaths individuals who tested corona-positive.

In the present context, we assume that there is no overlap between these two compartments, infected (I) and confirmed cases (C). In other words, tested corona-positive individuals are assumed to be unable to substantially spread the disease due to isolation and are immune to re-infection after recovery33. The aim of the present work is to estimate the distribution of the recovery and decease periods of COVID-19. In this goal, we split the compartment C into J subcomponents \(C_1, \ldots , C_J\), the compartment R into \(J \times M\) subcomponents \(R_{ik}\) for \(i = 1, \ldots , J\) and \(k = 1, \ldots , M\) and the compartment D into \(J \times M\) subcomponents \(D_{ik}\) for \(i = 1, \ldots , J\) and \(k = 1, \ldots , M\) where

$$\begin{aligned} C(t)= & {} \sum _{i=1}^J C_i(t) \,\,\, \text {or}\,\,\, C^{(m)} = \sum _{i=1}^J C_i^{(m)}, \end{aligned}$$
(1)
$$\begin{aligned} R(t)= & {} \sum _{i=1}^J \sum _{k=1}^M R_{ik}(t) \,\,\, \text {or}\,\,\, R^{(m)} = \sum _{i=1}^J \sum _{k=1}^M R_{ik}^{(m)}, \end{aligned}$$
(2)
$$\begin{aligned} D(t)= & {} \sum _{i=1}^J \sum _{k=1}^M D_{ik}(t) \,\,\, \text {or}\,\,\, D^{(m)} = \sum _{i=1}^J \sum _{k=1}^M D_{ik}^{(m)}. \end{aligned}$$
(3)

In (1), (2) and (3) m represents the time index, and \(C_i^{(m)}\), \(R_{ik}^{(m)}\) and \(D_{ik}^{(m)}\) represents the total corona-positive cases corresponding the incubation period \(\tau _i\), recovered individuals corresponding the incubation period \(\tau _i\) and onset time \(\theta _{ik}\) i.e., recovery period \(\zeta _k = \tau _i + \theta _{ik}\) and death toll corresponding the incubation period \(\tau _i\) and onset time \(\mu _{ik}\) i.e., decease period \(\eta _k = \tau _i + \mu _{ik}\), respectively, presented in Fig. 5a.

The time-dependent model is the following set of coupled delay differential equations, for \(i=1,\ldots ,J\):

$$\begin{aligned} \left\{ \begin{array}{lcl} \frac{dS}{dt} &{}=&{} -\beta \frac{SI}{N} - \alpha S + \nu L \,\ \\ \frac{dI}{dt} &{}=&{} \beta \frac{SI}{N} -\gamma I - \beta \sum _{i = 1}^J \delta _i \frac{ S(t -\tau _i)I(t -\tau _i)}{N} \,\ \\ \frac{dV}{dt} &{}=&{} \gamma I \,\ \\ \frac{dC_i}{dt} &{}=&{} \delta _i\beta \frac{ S(t -\tau _i)I(t -\tau _i)}{N} - \beta \sum _{k = 1}^M\lambda _{ik}\delta _i\frac{S(t-\tau _i-\theta _{ik})I(t-\tau _i-\theta _{ik})}{N} \,\ \\ &{} &{}- \beta \sum _{k = 1}^M\kappa _{ik}\delta _i\frac{S(t-\tau _i-\mu _{ik})I(t-\tau _i-\mu _{ik})}{N} \,\ \\ \frac{dR_{ik}}{dt} &{}=&{} \lambda _{ik}\delta _i\beta \frac{S(t-\tau _i-\theta _{ik})I(t-\tau _i-\theta _{ik})}{N} \,\ \\ \frac{dD_{ik}}{dt} &{}=&{} \kappa _{ik}\delta _i\beta \frac{S(t-\tau _i-\mu _{ik})I(t-\tau _i-\mu _{ik})}{N} \,\ \\ \frac{dL}{dt} &{}=&{} \alpha S - \nu L \,\ \end{array} \right. \end{aligned}$$
(4)

where the real positive modeling parameters \(\alpha\), \(\beta\), \(\gamma\) \(\delta _i\), \(\lambda _{ik}\), \(\kappa _{ik}\) and \(\nu\) are the rate of lockdown, the rate of infection, the rate of recovery from the asymptomatic group, the rate of tested corona-positive corresponding the incubation period \(\tau _i\), the rate of recovery corresponding the recovery period \(\zeta _k\), the rate of decease corresponding the decease period \(\eta _k\) and the rate of transit from lockdown compartment to susceptible compartment, respectively. The variables \(S(t - \tau _i)\) and \(I(t - \tau _i)\) denote the cumulative data of \((t -\tau _i)\) days, i.e., total number of suspected and infected individuals of \((t - \tau _i)\) days. The factors \(\delta _i\beta S(t - \tau _i)I(t - \tau _i)/N\), \(\lambda _{ik}\delta _i\beta S(t - \tau _i - \theta _{ik})I(t - \tau _i - \theta _{ik})/N\), \(\kappa _{ik}\delta _i\beta S(t - \tau _i - \mu _{ik})I(t - \tau _i - \mu _{ik})/N\) convey the rate of individuals who were infected \(\tau _i\) days ago, the rate of individuals who were infected \(\tau _i + \theta _{ik}\) days ago and recovered, the rate of individuals who were infected \(\tau _i + \mu _{ik}\) days ago and died, respectively. It follows from (4), that for any t

$$\begin{aligned} L(t) + S(t) + I(t) + V(t) + C(t) + R(t) + D(t) = N \, , \end{aligned}$$
(5)

where N (constant) is the total population size. We can define a group of new variable \(T_i\) for \(i = 1, \ldots , J\) such that

$$\begin{aligned} T_i = C_i + \sum _{k = 1}^{M} R_{ik} + \sum _{k = 1}^{M} D_{ik}, \end{aligned}$$
(6)

and

$$\begin{aligned} T(t) = \sum _{i=1}^J T_i(t) \,\,\, \text {or}\,\,\, T^{(m)} = \sum _{i=1}^J T_i^{(m)}, \end{aligned}$$
(7)

where T, total confirmed cases, is the group of individuals who tested corona positive (active cases + recovered + deaths). From Eq.(4) we can generate three different sets of coupled delay differential equations for \(i=1,\ldots ,J\) and \(k=1,\ldots ,M\)

$$\begin{aligned}&\left\{ \begin{array}{lcl} \frac{dS}{dt} &{}=&{} -\beta \frac{SI}{N} - \alpha S + \nu L \,\ \\ \frac{dI}{dt} &{}=&{} \beta \frac{SI}{N} -\gamma I - \beta \sum _{i = 1}^J \delta _i\frac{ S(t -\tau _i)I(t -\tau _i)}{N} \,\ \\ \frac{dT_i}{dt} &{}=&{} \delta _i\beta \frac{ S(t -\tau _i)I(t -\tau _i)}{N} \,\ \\ \frac{dL}{dt} &{}=&{} \alpha S - \nu L \,\ \end{array} \right. \end{aligned}$$
(8)
$$\begin{aligned}&\left\{ \begin{array}{lcl} \frac{dS}{dt} &{}=&{} -\beta \frac{SI}{N} - \alpha S + \nu L \,\ \\ \frac{dI}{dt} &{}=&{} \beta \frac{SI}{N} -\gamma I - \beta \sum _{i = 1}^J \delta _i \frac{ S(t -\tau _i)I(t -\tau _i)}{N} \,\ \\ \frac{dR_{ik}}{dt} &{}=&{} \lambda _{ik}\delta _i\beta \frac{S(t-\tau _i-\theta _{ik})I(t-\tau _i-\theta _{ik})}{N} \,\ \\ \frac{dL}{dt} &{}=&{} \alpha S - \nu L \,\ \end{array} \right. \end{aligned}$$
(9)

and

$$\begin{aligned} \left\{ \begin{array}{lcl} \frac{dS}{dt} &{}=&{} -\beta \frac{SI}{N} - \alpha S + \nu L \,\ \\ \frac{dI}{dt} &{}=&{} \beta \frac{SI}{N} -\gamma I - \beta \sum _{i = 1}^J \delta _i \frac{ S(t -\tau _i)I(t -\tau _i)}{N} \,\ \\ \frac{dD_{ik}}{dt} &{}=&{} \kappa _{ik}\delta _i\beta \frac{S(t-\tau _i-\mu _{ik})I(t-\tau _i-\mu _{ik})}{N} \,\ \\ \frac{dL}{dt} &{}=&{} \alpha S - \nu L \,\ \end{array} \right. \end{aligned}$$
(10)

where Eqs. (8), (9) and (10) can be used to calculate incubation period3, recovery period and decease period, respectively. In the present context, we focus on recovery as well as decease periods. We solve Eqs. (8), (9) and (10) using matlab inner-embedded function dde23 with particular sets of model parameters. To solve the initial value problem, in the interval \([t_0, t_1]\), we consider \(L(t_0)\), \(S(t_0)\), \(I(t_0)\), \(T(t_0)\), \(R(t_0)\) and \(D(t_0)\) as follows:

$$\begin{aligned} \left\{ \begin{array}{lcl} L(t_0) &{}=&{} 0 \, ,\\ S(t_0) &{}=&{} N - L(t_0) - I(t_0) - V(t_0) - T(t_0)\, , \\ I(t_0) &{}=&{} q \, , \\ V(t_0) &{}=&{} 0 \, , \\ T(t_0) &{}=&{} {\widetilde{T}}(t_0) \, , \\ R(t_0) &{}=&{} {\widetilde{R}}(t_0) \, , \\ D(t_0) &{}=&{} {\widetilde{D}}(t_0) \, , \end{array} \right. \end{aligned}$$
(11)

where \({\widetilde{T}}(t_0)\), \({\widetilde{R}}(t_0)\) and \({\widetilde{D}}(t_0)\) are the available data at time \(t_0\), and q is the initial value adjusting parameters. Initially, there are no lockdown individual and no removed individuals from the asymptomatic group so that we can consider \(L(t_0)\) = 0 and \(V(t_0)\) = 0. It follows from (7) and (11)

$$\begin{aligned} \sum _{i = 1}^{J} T_i(t_0) = T(t_0) = {\tilde{T}}(t_0) . \end{aligned}$$
(12)

In the present context \({\tilde{T}}(t_0) = 0\), since there were no corona-positive cases reported on January 22, 2020. As a consequence, we also take \(T_i(t_0) = 0\) for \(i = 1, 2, \ldots , J\), and the similar assumptions are valid for \(R_{ik}(t_0)\) and \(D_{ik}(t_0)\) i.e., \(R_{ik}(t_0) = 0\) and \(D_{ik}(t_0) = 0\) for \(i = 1, 2, \ldots , J\) and \(k = 1, 2, \ldots , M\).

Figure 5
figure 5

Model, methodology and estimated values of the parameters: (a) Schematic diagram of the present compartmental based model, total 830 compartments. Here \(\Lambda _{ik} = \lambda _{ik}\delta _i \beta S(t - \tau _i -\theta _{ik})I(t - \tau _i -\theta _{ik})/N\) and \(K_{ik} = \kappa _{ik}\delta _i \beta S(t - \tau _i -\mu _{ik})I(t - \tau _i -\mu _{ik})/N\) for \(i = 1, 2, \ldots , J\) and \(k = 1, 2, \ldots , M\). We consider \(J = 14\) and \(M = 29\). (b) Bubble diagram of the foundation of the present work, splitting publicly available database, total cases (T), recovered individuals (R) and death toll (D), into myriad groups. (c) Sketch of the Monotonic Iteration Scheme (MIS); for ‘recovery’ calculation \(\Delta _i = \{ \lambda _{ik} | k = 1, 2, \ldots , 29\}\) and for ‘decease’ calculation \(\Delta _i = \{ \kappa _{ik} | k = 1, 2, \ldots , 29\}\) and \(i = 1, 2, \ldots , 14\). (d) Sketch of the optimization scheme for the primary, \({\mathbf {P}}_0\), and secondary, \({\mathbf {P}}_1\) and \({\mathbf {P}}_2\), parameters. \({\mathbf {P}}_{1/2}\) indicates either \({\mathbf {P}}_1\) or \({\mathbf {P}}_2\). (e) Estimated values of the primary parameters. (f) Estimated values of the secondary parameters, upper panel: \(\lambda _{ik}\) and lower panel: \(\kappa _{ik}\). (g) Iteration verses error function in MIS, upper panel: estimating \(\lambda _{ik}\) and lower panel: estimating \(\kappa _{ik}\).

Parameter estimation

We focus on the exponential growth phase of the COVID-19 epidemic in Canada; one can use this approach to estimate the incubation, recovery period and decease periods for any region affected by this infectious disease. The time resolved (daily updated) database30 provides the number of total corona-positive cases, the number of recovered individuals and the death toll. We define two groups of model parameters: primary parameters, the parameters involved in Eq. (8) i.e., q, \(\alpha\), \(\beta\), \(\gamma\), \(\delta _i\) for \(i = 1, 2, \ldots , J\) and \(\nu\), and secondary parameters, the parameters involved in Eqs. (9) and (10) other than the primary parameters i.e., \(\lambda _{ik}\) and \(\kappa _{ik}\) for \(i = 1, 2, \ldots , J\) and \(k = 1, 2, \ldots , M\). We use the estimated values of the primary parameters to optimize the secondary parameters. The optimal values of the primary parameters \({\mathbf {P}}_0 = (q, \alpha , \beta , \delta _1(t), \ldots , \delta _J(t), \nu )^T\), q is the initial value of I(t), is obtained by minimizing the error function \(Er({\mathbf {P}}_0)\), defined as

$$\begin{aligned} Er({\mathbf {P}}_0) = \frac{1}{K}\sqrt{ \sum _{m = 1}^{K} (T^{(m)}({\mathbf {P}}_0) - {\widetilde{T}}^{(m)})^2 }\,\,\, , \end{aligned}$$
(13)

where \({\widetilde{T}}^{(m)}\) is the available data of total corona-positive cases on the particular mth day, and \(T^{(m)}\) is the calculated results obtained from the system (8). The integer K, used in (13), is the size of the data set. Due to the complexity of the error function, the minimization using the matlab function fminsearch requires a very large number of iterations. We use the similar error functions \(Er(\mathbf {P_1})\) and \(Er(\mathbf {P_2})\) to optimize the secondary parameters \(\mathbf {P_1} = (\lambda _{11}, \ldots , \lambda _{JM})^T\) and \(\mathbf {P_2} = (\kappa _{11}, \ldots , \kappa _{JM})^T\), defined as

$$\begin{aligned} Er(\mathbf {P_1}) = \frac{1}{K}\sqrt{ \sum _{m = 1}^{K} (R^{(m)}({\mathbf {P}}_0^{\text {op}},\mathbf {P_1}) - {\widetilde{R}}^{(m)})^2 }\,\,\, , \end{aligned}$$
(14)

and

$$\begin{aligned} Er(\mathbf {P_2}) = \frac{1}{K}\sqrt{ \sum _{m = 1}^{K} (D^{(m)}({\mathbf {P}}_0^{\text {op}}, \mathbf {P_2}) - {\widetilde{D}}^{(m)})^2 }\,\,\, , \end{aligned}$$
(15)

where \({\mathbf {P}}_0^{\text {op}}\) is the estimated values of \({\mathbf {P}}_0\); \({\widetilde{R}}^{(m)}\) and \({\widetilde{D}}^{(m)}\) are the available data of total number of recovered individuals and total number of death toll; \(R^{(m)}\) and \(D^{(m)}\) are the calculated results obtained from the Eqs. (9) and (10), respectively.

Numerical experiment

In this section, we present a detailed description of the computational procedure for the proposed model. On 23 January 2020, a 56-year old man admitted to Toronto hospital emergency department in Toronto with a new onset of fever and nonproductive cough, and returning from Wuhan, China, the day prior34,35. It is believed this is the first confirmed case of 2019-nCoV in Canada, and according to the government report, the novel coronavirus arrived on the Canadian coast on January 25, 2020, first reported case. The above information suggests that the start date of the current pandemic in Canada is possibly to be January 22, 2020. Additionally, some research studies reported that the estimation of the incubation period is from 2 to 14 days, and recovery as well as decease period of COVID-19 is from 2 to 6 weeks2,36. As a consequence, in the present study we consider \(J = 14\) i.e, 14 delays for the incubation period, and \(M = {42-14+1} = 29\) i.e, 29 delays for the recovery as well as decease periods. Here we consider a calculation of 177 days, from January 22, 2020 to July 16, 2020, duration of the first wave in Canada. The purpose of the model is to separate the publicly available database T, R and D into myriad groups \(T_i\), \(R_{ik}\) and \(D_{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\) (Fig. 5b).

Then the local minimum computed by the optimization algorithm depends on the initial values of the parameters: for q, \(\alpha\), \(\beta\), \(\nu\) we consider any positive random number less than unity, where as a choice of \({\varvec{\delta }}=(\delta _1, \ldots , \delta _{14})^T\) is tricky. For this purpose, we consider a vector of 14 positive random numbers \({\varvec{\delta }}\) such that \(\delta _1< \cdots< \delta _4 < \delta _5> \delta _6> \cdots > \delta _{14}\) and \(\sum _{i = 1}^{14} \delta _i = 0.9\). We observe, from numerous numerical experiments, the renormalization factor 0.9 works perfectly for the computation. The estimated values of the primary parameters \({\mathbf {P}}_0^{\text {op}}\) are presented in Fig. 5e, and the value of the error function \(Er({\mathbf {P}}_0)\) = 41.64. The estimated values of the primary parameters are related to Eq. (8), the set of coupled delay differential equations, and Eq. (13), the error function. Using the estimated values of the primary parameters, we optimize the secondary parameters \(\lambda _{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\) related to Eqs. (9) and (14). The choice of the initial values of \(\lambda _{ik}\) is such that for any fixed i, \(i = 1,2, \ldots , 14\), the first fourteen \(\lambda _{ik}\)s i.e., \(\{ \lambda _{i1}, \lambda _{i2}, \ldots , \lambda _{i14} \}\) are in ascending order, and the rest i.e., \(\{ \lambda _{i15}, \lambda _{i16}, \ldots , \lambda _{i29} \}\) are in descending order; and \(\sum _{k = 1}^{29} \lambda _{ik} = 0.72\). After optimization, we obtain the value of the error function \(Er(\mathbf {P_1}^{\text {op}})\) = 236.47. Using the estimated values of the primary parameters, we optimize the secondary parameters \(\kappa _{ik}\) for \(i = 1, 2, \ldots , 14\) and \(k = 1, 2, \ldots , 29\) related to Eqs. (10) and (15). The choice of the initial values of \(\kappa _{ik}\) is such that for any fixed i, \(i = 1,2, \ldots , 14\), the first fourteen \(\kappa _{ik}\)s i.e., \(\{ \kappa _{i1}, \kappa _{i2}, \ldots , \kappa _{i14} \}\) are in ascending order, and the rest i.e., \(\{ \kappa _{i15}, \kappa _{i16}, \ldots , \kappa _{i29} \}\) are in descending order; and \(\sum _{k = 1}^{29} \kappa _{ik} = 0.1\). After optimization, we obtain the value of the error function \(Er({\mathbf {P}}_2^{\text {op}})\) = 52.82. The values of the error functions \(Er({\mathbf {P}}_1^{\text {op}})\) and \(Er({\mathbf {P}}_2^{\text {op}})\) are not sufficiently small. To overcome that difficulties, here, we introduce a Monotonic Iteration Scheme (MIS).

Monotonic iteration scheme

To optimize the parameters \(\xi _{ik}\) for \(i = 1, 2, \ldots , J\) and \(k = 1, 2, \ldots , M\), in the present context \(J = 14\) and \(M = 29\), we use a MIS. However, the MIS can be applied for other numerical/optimization problems with any finite integer values of J and M. The schematic diagram of MIS is presented in Fig. 5c, and consists of the following steps.

  • Step 1 We decompose the parametric domain \(\Delta = \{ \xi _{ik} | i = 1, 2, \ldots , 14; k = 1, 2, \ldots , 29 \}\) into 14 subdomains \(\Delta _i = \{ \xi _{ik} | k = 1, 2, \ldots , 29 \}\) so \(\Delta = \{ \Delta _1, \Delta _2, \Delta _3, \ldots , \Delta _{14} \}\).

  • Step 2 We optimize the subdomain \(\Delta _1\) and consider the other parameters \(\Delta _2, \Delta _3, \ldots , \Delta _{14}\), as constants. After first iteration, we get estimated parameters \(\Delta _1^{\text {op}}\); the entire parametric domain is \(\Delta ^{(1)} = \{ \Delta _1^{\text {op}}, \Delta _2, \Delta _3, \ldots , \Delta _{14} \}\), and the error function \(Er(\Delta ^{(1)})\).

  • Step 3 In the second iteration, we optimize the subdomain \(\Delta _2\) and keeping the other subdomains of \(\Delta ^{(1)}\) unchanged. After second iteration, we get estimated parameters \(\Delta _2^{\text {op}}\); the entire parametric domain is \(\Delta ^{(2)} = \{ \Delta _1^{\text {op}}, \Delta _2^{\text {op}}, \Delta _3, \ldots , \Delta _{14} \}\), and the error function \(Er(\Delta ^{(2)})\).

  • Step 4 Repeated the same procedure discussed in \(\mathbf{Step} 3\).

The optimization of the subdomain \(\Delta _2\), demonstrated in Step3, is related to minimizing the error function such that \(Er(\Delta ^{(1)}) \ge Er(\Delta ^{(2)})\); the equality sign holds for \(\Delta _2^{\text {op}} = \Delta _2\). The error function of the \((n+1)\)th iteration, \(Er(\Delta ^{(n+1)})\) cannot be greater than that of nth iteration, \(Er(\Delta ^{(n)})\), because of this characteristic of the error function we define the approach as MIS. The flow chat of the optimization scheme is presented in Fig. 5d. The upper and lower panels of Fig. 5f show the estimated values of the secondary parameters \(\lambda _{ik}\) and \(\kappa _{ik}\), respectively, obtained from the MIS. The upper and lower panels of Fig. 5g show the values of the error functions \(E_R\), using MIS to optimize \(\lambda _{ik}\), and \(E_D\), using MIS to optimize \(\kappa _{ik}\), for fourteen iteration steps and \(Er({\mathbf {P}}_1^{\text {opMIS}}) = 104.07\), \(Er({\mathbf {P}}_2^{\text {opMIS}}) = 26.66\) where \({\mathbf {P}}_1^{\text {opMIS}}\) and \({\mathbf {P}}_2^{\text {opMIS}}\) are the estimated values of \({\mathbf {P}}_1^{\text {op}}\) and \({\mathbf {P}}_2^{\text {op}}\), respectively, using MIS. Figure 5g shows that the MIS works efficiently to get better estimations.