A computational approach to identifiability analysis for a model of the propagation and control of COVID-19 in Chile

Abstract

A computational approach is adapted to analyze the parameter identifiability of a compartmental model. The model is intended to describe the progression of the COVID-19 pandemic in Chile during the initial phase in early 2020 when government declared quarantine measures. The computational approach to analyze the structural and practical identifiability is applied in two parts, one for synthetic data and another for some Chilean regional data. The first part defines the identifiable parameter sets when these recover the true parameters used to create the synthetic data. The second part compares the results derived from synthetic data, estimating the identifiable parameter sets from regional Chilean epidemic data. Experiments provide evidence of the loss of identifiability if some initial conditions are estimated, the period of time used to fit is before the peak, and if a significant proportion of the population is involved in quarantine periods.

Keywords:

• COVID-19 model
• dynamical quarantines
• parameter estimation
• identifiability of parameters
• basic reproduction number
• simulated annealing

1. Introduction

1.1. Scope

The COVID-19 pandemic has received overwhelming news coverage worldwide. We recall that the coronavirus disease 2019 (COVID-19), caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), was declared a global pandemic by the World Health Organization (WHO) on March 11, 2020 [11, 64]. This highly contagious unprecedented virus has impacted government and public institutions, strained the health care systems, restricted people in their homes, and caused country-wide lockdowns resulting in a global economic crisis [20, 32, 61]. The impact of COVID-19 at the time of revision of this work (July 5, 2023) amounts to roughly 767 million cases moreover, 6.95 million deaths worldwide, of which 2.96 million have been reported in the Americas.

The morbidity, mortality, and economic indicators associated with the COVID-19 pandemic point to a devastating picture for Latin American countries. High poverty rates, poor leadership, high prevalence of underlying health conditions such as obesity and diabetes, and frail healthcare systems have exacerbated the impact of the novel coronavirus in this region [9, 23, 30, 45]. In Chile, the Ministry of Health reported the first COVID-19 case on March 3, 2020, becoming the fifth country in Latin America after Brazil, Mexico, Ecuador and Argentina to report COVID-19 cases. Soon after, the Chilean government put in place a number of interventions including the closure of all daycares, schools, and universities (March 16), border controls, telework recommendations (March 18), closure of non-essential businesses (March 19), national night curfews (March 22), as well as targeted lockdowns at the level of municipalities since March 26, 2020. As of July 6, 2023, Chile has accumulated 5288437 cases and 61586 deaths [40]. It is also worth noting that Chile has tested at a higher rate than any other Latin American country [58]. Fortunately, Chile's mass COVID-19 vaccination campaign to immunize about 80% of the total population started in February 2021, and more than 90% of that target population has been fully vaccinated as of July 6, 2023.

One way to understand and explain these epidemic phenomena is by applying mathematical models and surveillance data, and the COVID-19 outbreak has been studied with different focuses, differing in their uses, their mathematical form, and objectives [1, 6, 31, 43, 44, 47]. Since the first case reported of the COVID-19 disease, and the multiple models that have surged to respond to the emergency and its dynamics [53, 54], we can understand their contributions, as well as their limitations [2, 4, 62]. The mathematical models utilized to confront the coronavirus 2019 disease improve situational awareness, assess epidemiological characteristics, and inform the base of evidence for prevention strategies. Specifically, for models derived from the famous SEIR model, one must guarantee that information recovered by them is highly accurate and that one confidently understands the biological mechanisms to lay the ground for developing optimal strategies for public health. The successful application of such epidemic models depends upon our ability to estimate key transmission parameters. These estimates are subject to two principal sources of uncertainty: the noise of the data, and the assumptions built into the model. Both difficulties give rise to the problem of non-identifiability of parameters. That is, the information recovered by applying the model to available datasets is not representative of the dynamics that we desire to interpret. This problem has been studied intensively and there are various approaches to handle non-identifiability. For instance, in [16, 22, 29, 48, 57, 69] guidelines and criteria are proposed to avoid such identifiability difficulties. Ignoring these factors can result in misleading inferences with potentially incorrect public health policy decisions [12, 14, 51].

We are inspired by the dynamics of COVID-19 that occurred in Chile in the initial phase of 2020 to explore these factors in estimating parameters. Then, we propose a compartmental model with nonlinear differential equations for the Chilean dynamics, and we use it to distinguish the noise from the data and the structure. With this in mind, we developed an identifiability study employing a computational approach such as that described in [51]. The identifiability analysis, as explained in [50], determines how well the parameters of a model can be estimated from experimental data, so it is crucial for interpreting and determining confidence in model parameter values. Structural identifiability focuses on a mathematical model and the practical identifiability on data. Hence a model with structurally identifiable parameters may still, be nonidentifiable in practice due to a lack of information gleaned from real data. Thus, it is essential to understand both types of identifiability when applying an epidemic model to a specific situation [27, 39, 51]. Therefore, we applied a computational approach to investigate the identifiability of the parameters from the compartmental model mentioned. This idea motivates an alternative parameter identifiability analysis to a complex compartment model, where other techniques can be rigorous. Specifically, we employ synthetic data from our model using various parameter sets. This parameter estimation is carried out with a defined optimization process with the simulated annealing method. Fitting all data points generated in the simulated epidemic curve, we aim to recover the parameters assumed to create synthetic data. If this approach is successful, we may confirm structural identifiability, which we aim to validate by adding a measurement error structure through a bootstrapping process [14, 51]. For the case of practical identifiability, we only use a limited number of data points of the epidemic during the initial phase (e.g. pre-peak). We employed our sets of identifiable parameters to fit regional Chilean data to compare the conclusions from the previous analysis based on synthetic data.

1.2. Related work

We first recall that introductions to compartmental epidemiological models are provided in [3, 5, 18, 35, 60, 66]. This class of models goes back to the well-known work by Kermack and McKendrick [28]. Within a compartmental model the population is subdivided into at least two epidemiological compartments (e.g. susceptible and infected). The rates of progression between the compartments, as well as the incidence and possibly birth and death of individuals, need to be specified, which leads to a dynamical system, that is a system of coupled, nonlinear ordinary differential equations (ODEs) that describe the progression of the epidemic.

In many cases, the compartmental approach formulated as a dynamical system provides indications of the basic reproduction number $R_0$, where this parameter plays the role of a threshold value for the system's dynamics. This role of $R_0$ as a ‘threshold value’ differs from the original epidemiological significance, namely that $R_0$ gives the number of secondary cases one infectious individual will produce in a population consisting only of susceptible individuals [35, p. 21].

The conclusions obtained by a compartmental model are derived from parameter estimations and allow describing the dynamics present in a data set. Then, the identifiability analysis is crucial to guaranteeing that the parameter estimates are adequate. Specifically, the identifiability analysis addresses whether the parameters of a compartmental model can be uniquely identified. There are different strategies for the identifiability, for example, for models expressed by linear ODEs [38, 65], nonlinear ODEs [39], or partial differential equations [50]. The type of data to fit also plays an essential role. For example, one needs to select between incidence or prevalence curves [12, 19, 56]. In the present theory two types of identifiability are considered: structural and practical. The first involves the model and its structure, and the study can be centred on theoretical or computational techniques (e.g. [13, 16, 19, 25, 33, 39, 50, 57]), and the second involves the data (e.g. [22, 27, 29, 56, 57]). The analysis is separate because a model can be structurally identifiable, but real-world data employed to fit the model results may lead to unidentifiable parameters for reasons related to the characteristics of the data available.

We extend the computational approach developed by Rossa and Chowell [51] to assess parameter identifiability in a compartmental model defined by nonlinear ODEs. Then, the computational approach that combines the optimization process determined by the simulated annealing method and a bootstrap process allows us to assess the identifiability of the corresponding parameters. In particular, structural identifiability is assessed by employing synthetic data with known parameters and the practical one with the same synthetic data, but the model is fitted only to the initial phases of the data curve (e.g. pre-peak periods). In both computations we expect to recover the set of assumed parameters with small uncertainty and error. These properties ensure that the model is indeed associated with identifiable parameters. Finally, we apply the identifiability analysis to various sets of parameters to estimate. We vary their elements in increasing numbers and measure the ability of the model to estimate the parameter sets in each case by fitting the data curves generated with assumed parameters. This procedure can be understood as a sensitivity analysis.

In the application to Chilean data, we consider the parameter sets identifiable from our computational approach (with synthetic data) and fit some regional Chilean data. We hope to confirm the conclusions from the identifiability study with our application to the Chilean case.

In our study, a parameter set is identifiable if it satisfies some criteria defined in Section 2.5, which are an adaptation from the definition presented in [51], considering analyzing the confidence intervals and the mean square error (MSE).

2. Methods

2.1. Compartmental model

We apply our methodology to a complex compartmental transmission model inspired by the strategies implemented by the Chilean government to mitigate the COVID-19 emergency in the initial phase of the outbreak. For this purpose, we consider a single population and adopt a simplified version of the model described in [15], combined with the description of individuals in quarantine proposed in [26]. Individuals within the model are classified as susceptible (S), in home quarantine (Q), latent ($E_1$), partially infectious but not yet symptomatic ($E_2$), asymptomatic (A), infectious and will not be tested ($I_{\mathrm{n}}$), infectious and will be tested and isolated ($I_{\mathrm{s}}$), hospitalized/isolated infectious (J), recovered (R), and deceased (D). Constant population size is assumed, i.e. (1) $N:=S+Q+E_1+E_2+A+I_{\mathrm{n}}+I_{\mathrm{s}}+J+R+D=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$(1) Note that we stipulate one single class of asymptomatic individuals. At the same time, in [15] a distinction is made between asymptomatic individuals who will not be tested and those who are asymptomatic and will be tested and isolated.

It is assumed that home-quarantined individuals do not have contact with infected individuals, due to severe travel restrictions and rigorous supervision by their local communities, do not have contact with infected individuals. The parameters p and $1/\lambda$ represent the percentage of individuals in quarantine and the quarantine duration, respectively. Therefore, class Q removes susceptible individuals from the infection dynamics when $p\ne 0$ or $\lambda \ne 0$, and there is no quarantine when $p=\lambda =0$.

Five classes can contribute to new infections, namely $E_2$, A, $I_{\mathrm{n}}$, $I_{\mathrm{s}}$, and J. Each class has a relative transmission linked with the parameters $q_{\mathrm{e}}$, $q_{\mathrm{a}}$, q and h. These parameters are the relative transmissibilities from partially infectious population $E_2$, asymptomatic population A, asymptomatic/infected population that wear PPE, and the infected population isolated or in hospitals J. Consequently, the product hJ models the contribution of hospitalized individuals to transmission.

If we denote by $r_{X\to Y}$ the rate at which individuals move from class X to class Y, then susceptible individuals move to $E_1$ at the rate $r_{S\to E_1}=\frac{\text{β}}{N}(q_{\mathrm{e}}{E}_2+q_{\mathrm{a}}qA+q{I}_{\mathrm{n}}+q{I}_{\mathrm{s}}+hJ),$ where β denotes the overall transmission rate. Individuals in $E_1$ progress to $E_2$ at rate ${\kappa }_1$. Individuals from $E_2$ are partially infectious, with relative transmissibility $q_{\mathrm{e}}$, and progress at rate ${\kappa }_2$, where a proportion ${\rho }_{\mathrm{a}}$ become asymptomatic and partially infectious (relative transmissibility $q_{\mathrm{a}}$), and $1-{\rho }_{\mathrm{a}}$ become fully infectious. Among the proportion $1-{\rho }_{\mathrm{a}}$ that become fully infectious, ${\rho }_{\mathrm{s}}$ will be tested, while $1-{\rho }_{\mathrm{s}}$ will be undetected. Asymptomatic individuals who are not tested and symptomatic individuals wear personal protective equipment (PPE) (such as wearing masks in public, handwashing, etc.) and thus have relative transmissibility q, which quantifies the effectiveness of those protective behaviours. Individuals within classes A and $I_{\mathrm{n}}$ (who are not tested) recover at rate ${\gamma }_1$. Those tested (class $I_{\mathrm{s}}$) will progress to the hospitalized and isolated class at diagnosis rate α. Relative transmission within hospitals and isolated places may occur, but our analysis assumes perfect isolation. Hospitalized and isolated indiividuals progress to the recovered class at rate ${\gamma }_2$ or to the deceased class at rate δ. Therefore, the model is defined by the following system of non-linear ODEs for a single population, where all variables in capital letters are functions of time t, i.e. it is understood that $S=S\left(t\right)$, $E_1=E_1\left(t\right)$ and so on, and a prime denotes the derivative, that is $S^{\prime }\equiv \mathrm{d}S/\mathrm{d}t$, etc.: (2) $\begin{array}{rl}{S}^{\prime }& =-\frac{\text{β}S}{N}(q_{\mathrm{e}}{E}_2+q_{\mathrm{a}}qA+q{I}_{\mathrm{n}}+q{I}_{\mathrm{s}}+hJ)-pS+\lambda Q,\\ Q^{\prime }& =pS-\lambda Q,\\ E_1^{\prime }& =\frac{\text{β}S}{N}(q_{\mathrm{e}}{E}_2+q_{\mathrm{a}}qA+q{I}_{\mathrm{n}}+q{I}_{\mathrm{s}}+hJ)-{\kappa }_1{E}_1,\\ E_2^{\prime }& ={\kappa }_1{E}_1-{\kappa }_2{E}_2,\\ A^{\prime }& ={\kappa }_2{\rho }_{\mathrm{a}}{E}_2-{\gamma }_{1}A,\\ I_{\mathrm{n}}^{\prime }& ={\kappa }_2(1-{\rho }_{\mathrm{a}})(1-{\rho }_{\mathrm{s}})E_2-{\gamma }_1{I}_{\mathrm{n}},\\ I_{\mathrm{s}}^{\prime }& ={\kappa }_2(1-{\rho }_{\mathrm{a}}){\rho }_{\mathrm{s}}{E}_2-\alpha I_{\mathrm{s}},\\ J^{\prime }& =\alpha I_{\mathrm{s}}-({\gamma }_2+\delta )J,\\ R^{\prime }& ={\gamma }_1(A+I_{\mathrm{n}})+{\gamma }_{2}J,\\ D^{\prime }& =\delta J,\\ C^{\prime }& =\alpha I_{\mathrm{s}}.\end{array}$(2) The auxiliary variable C tracks the cumulative number of diagnosed/reported cases from the start of the outbreak, and $C^{\prime }$ represents the new diagnosed cases. The quantity C is not a state of the system of equations but simply a class to track the cumulative incidences cases; thus, individuals from population are not moving to class C. The initial values of these state variables are denoted by $S\left(0\right),Q\left(0\right),E_1\left(0\right),\dots ,C\left(0\right)$, where $S\left(0\right)=N-Q\left(0\right)-E_1\left(0\right)-E_2\left(0\right)-A\left(0\right)-I_{\mathrm{n}}\left(0\right)-I_{\mathrm{s}}\left(0\right)-J\left(0\right)-R\left(0\right)-D\left(0\right).$ A schematic of the transmissions is provided in Figure 1.

Figure 1. Schematic of the transmission of COVID-19 for one population.

The basic, control, and effective control reproduction numbers $R_0$, $R_0^{\mathrm{c}}$ and $R_{\mathrm{e}}^{\mathrm{c}}$ quantify the epidemic impact of infectious disease in a population and the effectiveness of control measures. The basic reproduction number $R_0$ defines the average number of secondary cases generated by primary infectious individuals during the early transmission period when the population is completely susceptible and without control interventions. The control reproduction number $R_0^{\mathrm{c}}$ defines the average number of new cases generated by infectious individuals when the population susceptible is with some control measure or intervention. The effective control reproduction number $R_{\mathrm{e}}^{\mathrm{c}}$ is the reproduction number that varies with time, and it is defined as partially susceptible population, where if $R_{\mathrm{e}}^{\mathrm{c}}>1$, the disease transmission continues, and if $R_{\mathrm{e}}^{\mathrm{c}}<1$ the disease transmission declines. Therefore, while the basic reproduction number $R_0$ and control reproduction number $R_0^{\mathrm{c}}$ captures the transmission potential of infectious disease during the early epidemic, the effective control reproduction number $R_{\mathrm{e}}^{\mathrm{c}}$, captures changes in the transmission potential over time.

To calculate these numbers, we employ the well-known next-generation matrix approach [59]. This yields the expression $R_0^{\mathrm{c}}=\rho \left(\mathcal{F}{\mathcal{V}}^{-1}\right),$ where ρ denotes the spectral radius and for the present model, the matrices $\mathcal{F}$ and $\mathcal{V}$ are given by $\begin{array}{rl}\mathcal{F}& =\left[\begin{array}{cccccc}0& \text{β}{q}_{\mathrm{e}}& \text{β}{q}_{\mathrm{a}}q& \text{β}q& \text{β}q& \text{β}h\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right],\\ \mathcal{V}& =\left[\begin{array}{cccccc}{\kappa }_1& 0& 0& 0& 0& 0\\ -{\kappa }_1& {\kappa }_2& 0& 0& 0& 0\\ 0& -{\kappa }_2{\rho }_{\mathrm{a}}& {\gamma }_1& 0& 0& 0\\ 0& -{\kappa }_2(1-{\rho }_{\mathrm{a}})(1-{\rho }_{\mathrm{s}})& 0& {\gamma }_1& 0& 0\\ 0& -{\kappa }_2(1-{\rho }_{\mathrm{a}}){\rho }_{\mathrm{s}}& 0& 0& \alpha & 0\\ 0& 0& 0& 0& -\alpha & {\gamma }_2+\delta \end{array}\right].\end{array}$ The result is (3) $R_0^{\mathrm{c}}=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}+\frac{q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha }+\frac{h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta }).$(3) Moreover, the corresponding effective reproduction number is defined as $R_{\mathrm{e}}^{\mathrm{c}}\left(t\right)=R_0^{\mathrm{c}}S\left(t\right)/N,$ where $S\left(t\right)/N$ is the proportion of susceptible individuals in the population at time t, where the total population N is assumed to be constant (see (1(1) $N:=S+Q+E_1+E_2+A+I_{\mathrm{n}}+I_{\mathrm{s}}+J+R+D=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$(1) )). The explicit expression (3(3) $R_0^{\mathrm{c}}=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}+\frac{q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha }+\frac{h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta }).$(3) ) implies the following contributions of the individual compartments: $\begin{array}{rl}{R}^{E_2}& =\frac{\text{β}{q}_{\mathrm{e}}}{{\kappa }_2},R^A=\frac{\text{β}{q}_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1},R^{I_{\mathrm{n}}}=\frac{\text{β}q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1},\\ R^{I_{\mathrm{s}}}& =\frac{\text{β}q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha },R^J=\frac{\text{β}h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta },\end{array}$ where (4) $R_0^{\mathrm{c}}=R^{E_2}+R^A+R^{I_{\mathrm{n}}}+R^{I_{\mathrm{s}}}+R^J.$(4) This parameter $R_0^{\mathrm{c}}$ is a key epidemic parameter to measure the transmissibility of a pathogen; by (3(3) $R_0^{\mathrm{c}}=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}+\frac{q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha }+\frac{h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta }).$(3) ) we see that $R_0^{\mathrm{c}}$ is a function of several parameters of the epidemic model (2(2) $\begin{array}{rl}{S}^{\prime }& =-\frac{\text{β}S}{N}(q_{\mathrm{e}}{E}_2+q_{\mathrm{a}}qA+q{I}_{\mathrm{n}}+q{I}_{\mathrm{s}}+hJ)-pS+\lambda Q,\\ Q^{\prime }& =pS-\lambda Q,\\ E_1^{\prime }& =\frac{\text{β}S}{N}(q_{\mathrm{e}}{E}_2+q_{\mathrm{a}}qA+q{I}_{\mathrm{n}}+q{I}_{\mathrm{s}}+hJ)-{\kappa }_1{E}_1,\\ E_2^{\prime }& ={\kappa }_1{E}_1-{\kappa }_2{E}_2,\\ A^{\prime }& ={\kappa }_2{\rho }_{\mathrm{a}}{E}_2-{\gamma }_{1}A,\\ I_{\mathrm{n}}^{\prime }& ={\kappa }_2(1-{\rho }_{\mathrm{a}})(1-{\rho }_{\mathrm{s}})E_2-{\gamma }_1{I}_{\mathrm{n}},\\ I_{\mathrm{s}}^{\prime }& ={\kappa }_2(1-{\rho }_{\mathrm{a}}){\rho }_{\mathrm{s}}{E}_2-\alpha I_{\mathrm{s}},\\ J^{\prime }& =\alpha I_{\mathrm{s}}-({\gamma }_2+\delta )J,\\ R^{\prime }& ={\gamma }_1(A+I_{\mathrm{n}})+{\gamma }_{2}J,\\ D^{\prime }& =\delta J,\\ C^{\prime }& =\alpha I_{\mathrm{s}}.\end{array}$(2) ), including transmission rates and infectious periods of the epidemiological states that contribute to new contagions. In addition, we can observe that the basic reproduction number $R_0$ corresponds to the control reproduction number $R_0^{\mathrm{c}}$ in the absence of interventions, i.e. when q = 1, h = 0, and ${\rho }_{\mathrm{s}}=0$ (where q is interpreted as the level of no protection and 1−q is the level of protection due to the use of PPE). These assumptions lead to the expression $R_0=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}).$

2.2. Simulated data

Using our model proposed for COVID-19, we simulate various scenarios fixing their parameters and initial conditions. These curves will help analyze the ability of our model to capture such dynamics. Then a parameter estimation process is necessary together with a bootstrapping procedure to examine the identifiability of several of its parameters, as will be explained in the following. For this, we assume the values and references summarized in Table 1 for the parameters, and for the initial conditions the values $\begin{array}{rl}J\left(0\right)& =C\left(0\right)=1,Q\left(0\right)=D\left(0\right)=R\left(0\right)=0,I_{\mathrm{n}}\left(0\right)=1,I_{\mathrm{s}}\left(0\right)=20,\\ E_1\left(0\right)& =2J\left(0\right),E_2\left(0\right)=4J\left(0\right)-E_1\left(0\right),A\left(0\right)=4,\\ S\left(0\right)& =N-E_1\left(0\right)-E_2\left(0\right)-I_{\mathrm{n}}\left(0\right)-I_{\mathrm{s}}\left(0\right)-A\left(0\right)-J\left(0\right)-R\left(0\right)-D\left(0\right)\end{array}$ for a total population N = 500, 000. These values are used to simulate the following scenarios.

Table 1. Values of model parameters selected to generate simulated data.

Parameter	Description	Value	Source
		$1.0676,$	Computed
β	transmission rate	$1.2435,$	using (3(3) $R_0^{\mathrm{c}}=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}+\frac{q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha }+\frac{h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta }).$(3) )
		1.7794
h	relative isolation transmissibility of infected individuals	0	Assumed
$q_{\mathrm{e}}$	relative transmissibility of exposed individuals	0.1	[17]
$q_{\mathrm{a}}$	relative transmissibility of asymptomatic cases	0.4	[17]
q	relative transmissibility of asymptomatic and infected cases that wear PPE	0.4	Assumed
$1/{\kappa }_1$	length of latent period	2.5	[49, 72]
$1/{\kappa }_2$	length of infectiousness prior to symptom onset	2.5	[49, 72]
${\rho }_{\mathrm{a}}$	proportion of exposed individuals who become asymptomatically infected	0.4	[42, 46]
${\rho }_{\mathrm{s}}$	proportion of fully infectious individuals who undergo testing	0.6	Assumed
$1/\alpha$	time from symptom onset to isolation or hospitalization	10	Assumed
$1/{\gamma }_1$	time from illness onset to recovery	7	[49, 71]
$1/{\gamma }_2$	time from diagnosis to recovery	5	[55]
δ	death rate within hospitals	0.021	Computed Chilean data
p	proportion of the susceptible population in quarantine declared	0.01, 0.05	Assumed
$1/\lambda$	duration of quarantine	15	Assumed
$R_0^{\mathrm{c}}$	control reproduction number	3,4,5	Assumed

Note: The parameter $R_0^{\mathrm{c}}$ is fixed, and β is computed for each one using the expression (3(3) $R_0^{\mathrm{c}}=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}+\frac{q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha }+\frac{h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta }).$(3) ).

2.2.1. Scenario 1: initial phase for three values of $R_0^{\mathrm{c}}$ without quarantine

We consider three values of $R_0^{\mathrm{c}}$ that generate three values of β via (3(3) $R_0^{\mathrm{c}}=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}+\frac{q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha }+\frac{h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta }).$(3) ). These values along with other parameters are summarized in Table 1. The corresponding incidence curves $t\mapsto C^{\prime }\left(t\right)$ are generated for $t\in [0,200]$, see Figure 2(a). These incidence curves will be our datasets for the identifiability study. Notice that the curve for $R_0^{\mathrm{c}}=3$ effectively grows more slowly than the others, and its maximum is smaller than the curves with $R_0^{\mathrm{c}}=4$ and $R_0^{\mathrm{c}}=5$.

Figure 2. Incidence curves $C^{\prime }=C^{\prime }\left(t\right)$ for synthetic data based on the parameters of Table 1(a) of Scenario 1, (b) of Scenario 2 with p = 0.01, (c) of Scenario 2 with p = 0.05. Magenta, green, and blue curves correspond to $R_0^{\mathrm{c}}=3$, 4, and 5, respectively.

2.2.2. Scenario 2: initial phase for three values of $R_0^{\mathrm{c}}$ with quarantine.

For this scenario, we consider the same values for $R_0^{\mathrm{c}}$ and corresponding values of β as in Scenario 1 since expression (3(3) $R_0^{\mathrm{c}}=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}+\frac{q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha }+\frac{h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta }).$(3) ) does not depend on the parameters p or λ. The pairs (p = 0.01, $\lambda =1/15$) and (p = 0.05, $\lambda =1/15$), with the other parameters as given in Table 1, produce the datasets shown in Figure 2(b,c), respectively. Again, the curve with $R_0^{\mathrm{c}}=3$ grows more slowly than the others, and its maximum is smaller than the curves with $R_0^{\mathrm{c}}=4$ and $R_0^{\mathrm{c}}=5$. Besides, the maxima of the curves with p = 0.05 are consistently smaller than those of the corresponding curves with p = 0.01, which evidences the effect of the quarantine modelled.

Using the bootstrapping approach, we intend to use these curves defined for each scenario to generate multiple samples from the best-fit model. Before doing this, we need to define the parameter estimation process to explain the bootstrap approach and finally the identifiability study.

2.3. Parameter estimation

To estimate the parameter values of our COVID-19 model, we need to use one optimization process to fit the model to each simulated dataset. Then, the method used for this purpose is outlined as follows, executed in Matlab (Mathworks, Inc).

2.3.1. Mixed simulated annealing method with a least-squares approach (SA-LSQ)

The SA-LSQ approach is defined as a parameter estimation process that combines the simulated annealing and least-squares optimization algorithms to minimize the Euclidean distance between the incidence curve $C^{\prime }(t,\mathrm{\Theta })$, defined by our COVID-19 model (2(2) $\begin{array}{rl}{S}^{\prime }& =-\frac{\text{β}S}{N}(q_{\mathrm{e}}{E}_2+q_{\mathrm{a}}qA+q{I}_{\mathrm{n}}+q{I}_{\mathrm{s}}+hJ)-pS+\lambda Q,\\ Q^{\prime }& =pS-\lambda Q,\\ E_1^{\prime }& =\frac{\text{β}S}{N}(q_{\mathrm{e}}{E}_2+q_{\mathrm{a}}qA+q{I}_{\mathrm{n}}+q{I}_{\mathrm{s}}+hJ)-{\kappa }_1{E}_1,\\ E_2^{\prime }& ={\kappa }_1{E}_1-{\kappa }_2{E}_2,\\ A^{\prime }& ={\kappa }_2{\rho }_{\mathrm{a}}{E}_2-{\gamma }_{1}A,\\ I_{\mathrm{n}}^{\prime }& ={\kappa }_2(1-{\rho }_{\mathrm{a}})(1-{\rho }_{\mathrm{s}})E_2-{\gamma }_1{I}_{\mathrm{n}},\\ I_{\mathrm{s}}^{\prime }& ={\kappa }_2(1-{\rho }_{\mathrm{a}}){\rho }_{\mathrm{s}}{E}_2-\alpha I_{\mathrm{s}},\\ J^{\prime }& =\alpha I_{\mathrm{s}}-({\gamma }_2+\delta )J,\\ R^{\prime }& ={\gamma }_1(A+I_{\mathrm{n}})+{\gamma }_{2}J,\\ D^{\prime }& =\delta J,\\ C^{\prime }& =\alpha I_{\mathrm{s}}.\end{array}$(2) ), and the time series data (${\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}}_t$) that corresponds to time series data with $t\in [0,T]$, where T is the final time (date). Here Θ is the set of parameters to estimate of our COVID-19 model.

In our approach, we obtain the set of estimated parameters $\hat{\mathrm{\Theta }}$ that satisfies $\hat{\mathrm{\Theta }}=\underset{\mathrm{\Theta }\in S_{\mathrm{\Theta }}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}f(\mathrm{\Theta },R),$ with initial guess $S_{\mathrm{\Theta }}$ and the objective function $f(\mathrm{\Theta },R):=\frac{1}{R}\sum _{t_i\le T}{\left|C^{\prime }\right(t_i,\mathrm{\Theta })-{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}}_{t_i}|}^2,$ for a given normalizing coefficient R.

We wish to avoid attaining local minima and to guarantee that the optimization process finds the global minimum. Consequently, due to the sensitivity of the solution with respect to the initial guess utilized for these processes, we decided, before starting the SA fitting, to select the best initial guess running several Latin hypercube sampling (LHS) runs. This is done by the built-in Matlab function lhsdesign(), applied to the bounds defined for each parameter to estimate.

Together with simulated annealing (SA), using the Matlab built-in function simulannealbnd(), and the best initial guess obtained with LHS, we estimated the parameters that are transformed as the initial guess to apply the least squares method (LSQ) with lsqcurvefit() Matlab function to speed up the steepest descent phase of the simulated annealing algorithm.

More details about the procedure are summarized in Algorithm 1.

This method is a powerful optimization tool, particularly for fitting models that describe epidemic phenomena with real data, such as those developed in [7, 8, 36, 68, 70] Additionally, the SA technique was inspired by Metropolis et al. [37] to select the new characterized solutions. Even so, its adaptation allows the application to a frequentist paradigm.

Specifically, to study parameter identifiability, we define various experiments that use the generated curves in Scenarios 1 and 2, and we intend to estimate or recover some parameters of interest from these curves. Therefore parameters of interest are selected from $\{\alpha ,\text{β},{\rho }_{\mathrm{s}},I_{\mathrm{s}}(0),J(0),p\}$, while the parameters not selected are fixed to the respective value of Table 1.

2.3.2. Sets of parameters to be estimated

For the data generated using Scenarios 1 and 2, we define the sets of parameters to estimate given in Table 2. Each combination explored will be mentioned as an experiment. For example, Experiment 101 indicates to data from Scenario 1 to estimate the parameter set 1 and Experiment 201-1 corresponds to data from Scenario 2 to estimate the parameter set 1, with p = 0.01, where Experiment 201-2 will be with p = 0.05 (see Tables 3 and 4). Therefore, for Scenario 1, we have 24 experiments, that is, the number of values assumed of $R_0^{\mathrm{c}}$ (three) times the number of parameter sets (eight). Likewise, for Scenario 2, we have 84 experiments corresponding to the number of values of $R_0^{\mathrm{c}}$ (three) times the number of values of p (two) times the number of parameter sets (fourteen).

Table 2. Set of parameters to be estimated in Scenarios 1 and 2.

	Experiment	Parameter set (Θ)	Experiment	Parameter set (Θ)
Scenario 1	1	β	5	$\text{β},I_{\mathrm{s}}\left(0\right),J\left(0\right)$
	2	$(\alpha ,\text{β})$	6	$\alpha ,\text{β},I_{\mathrm{s}}\left(0\right),J\left(0\right))$
	3	$(\text{β},{\rho }_{\mathrm{s}}$	7	$(\text{β},{\rho }_{\mathrm{s}},I_{\mathrm{s}}(0),J(0\left)\right)$
	4	$(\alpha ,\text{β},{\rho }_{\mathrm{s}})$	8	$(\alpha ,\text{β},{\rho }_{\mathrm{s}},I_{\mathrm{s}}(0),J(0\left)\right)$
Scenario 2	1	β	8	$(\text{β},I_{\mathrm{s}}(0),J(0\left)\right)$
	2	$(\alpha ,\text{β})$	9	$(\alpha ,\text{β},I_{\mathrm{s}}(0),J(0\left)\right)$
	3	$(\text{β},{\rho }_{\mathrm{s}})$	10	$(\text{β},{\rho }_{\mathrm{s}},I_{\mathrm{s}}(0),J(0\left)\right)$
	4	$(\alpha ,\text{β},{\rho }_{\mathrm{s}})$	11	$(\alpha ,\text{β},{\rho }_{\mathrm{s}},I_{\mathrm{s}}(0),J(0\left)\right)$
	5	$(\text{β},p)$	12	$(\text{β},I_{\mathrm{s}}(0),J(0),p)$
	6	$(\alpha ,\text{β},p)$	13	$(\text{β},\alpha ,I_{\mathrm{s}}(0),J(0),p)$
	7	$(\text{β},\alpha ,{\rho }_{\mathrm{s}},p)$	14	$(\alpha ,\text{β},{\rho }_{\mathrm{s}},I_{\mathrm{s}}(0),J(0),p)$

Table 3. Summary of results obtained for Scenario 1.

Note: A green box indicates that the parameter estimated satisfies the parameter identifiability conditions (small confidence interval with centred mean value, small MSE). Parameters that do not meet the identifiability criteria are coloured in red.

Table 4. Summary of results for Scenario 2.

Note: The meaning of green and red boxes is as in Table 3.

2.4. Bootstrapping method

This method uses the parametric bootstrapping approach with Poisson error structure (see [14, 51]) to repeatedly sample observations from the best-fit obtained by the SA-LSQ fit of our COVID-19 model to an time series data. This process can be summarized in the following steps:

1. Obtain the deterministic model solution $C(t,\hat{\mathrm{\Theta }})$ or total daily incidence series from the best-fit obtained using the SA-LSQ estimation.

2. Generate S replicated datasets assuming the Poisson error structure. To generate these simulated data, we incorporate the Poisson error structure using the incidence curve $C^{\prime }(t,\hat{\mathrm{\Theta }})$, where for each time $t_i$ we generate a new incidence value using a Poisson random variable with mean $C^{\prime }(t_i,\hat{\mathrm{\Theta }})$. Therefore, this new dataset represents an incidence curve for the system, assuming the time series follows a Poisson distribution centred on the mean at time points $t_i$.

3. Re-estimate model parameters for each of the S simulated realizations (using the SA-LQS method), which are denoted by ${\hat{\mathrm{\Theta }}}_i$ for $i=1,\dots ,S$.

4. Finally, using the set of ${\hat{\mathrm{\Theta }}}_i$ for $i=1,\dots ,S$, we construct the confidence intervals for each estimated parameter. Also, for each set of estimated parameters, $R_0^{\mathrm{c}}$ is calculated to obtain a distribution of $R_0^{\mathrm{c}}$ values as well.

2.5. Parameter identifiability

2.5.1. Preliminaries

The complexity of our mathematical model involves ten epidemiological states or compartments and an auxiliary variable $C\left(t\right)$ that indicates the cumulative number of symptomatic diagnosed or new cases reported. In addition to the ten system states, our model consists of 15 parameters (see Table 1), where five classes contribute to a new infection $(E_2,A,I_{\mathrm{n}},I_{\mathrm{s}},J)$, which also contribute to the representation of $R_0^{\mathrm{c}}$ in these classes, see Equations (3(3) $R_0^{\mathrm{c}}=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}+\frac{q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha }+\frac{h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta }).$(3) ) and (4(4) $R_0^{\mathrm{c}}=R^{E_2}+R^A+R^{I_{\mathrm{n}}}+R^{I_{\mathrm{s}}}+R^J.$(4) ). Then, the estimated individual parameters are used to compute $R_0^{\mathrm{c}}$. Some initial conditions of the model are unknown and possibly need to be estimated along with the other model parameters. For our case we consider the values $I_{\mathrm{s}}\left(0\right)$ and $J\left(0\right)$ to estimate together to the parameters α, β, ${\rho }_{\mathrm{s}}$, and p. This defines the set from parameters to estimate $\{\alpha ,\text{β},{\rho }_{\mathrm{s}},p,I_{\mathrm{s}}(0),J(0\left)\right\}$, where the remaining parameters and initial conditions are known or assumed.

To estimate the parameter sets, we fit the model (with other parameters and initial values fixed) to each dataset using the auxiliary variable $C^{\prime }$ and the curve of the daily cases generated. The parameter estimation procedure employs the methods detailed in Section 2.3. The initial parameter prediction affects the solution for the model as local minima occur. To guarantee the global minimum obtained in the optimization processes, we generate a multi-start (using LHS in Algorithm 1) repeating these processes for a random amount of initial parameters, whose ranges are defined as follows: $\begin{array}{rl}& 0<\text{β}<5,0.1<\alpha <2,0<J\left(0\right)<C\left(0\right),\\ & 0<I_{\mathrm{s}}\left(0\right)<300,0<{\rho }_{\mathrm{s}}<1,\mathrm{a}\mathrm{n}\mathrm{d}0<p<1.\end{array}$ Here, the particular choice of the wide range for $I_{\mathrm{s}}\left(0\right)$ aims to avoid restrictions during the calibration of the model due to its sensitivity.

We focus the parameter identifiability study for (2(2) $\begin{array}{rl}{S}^{\prime }& =-\frac{\text{β}S}{N}(q_{\mathrm{e}}{E}_2+q_{\mathrm{a}}qA+q{I}_{\mathrm{n}}+q{I}_{\mathrm{s}}+hJ)-pS+\lambda Q,\\ Q^{\prime }& =pS-\lambda Q,\\ E_1^{\prime }& =\frac{\text{β}S}{N}(q_{\mathrm{e}}{E}_2+q_{\mathrm{a}}qA+q{I}_{\mathrm{n}}+q{I}_{\mathrm{s}}+hJ)-{\kappa }_1{E}_1,\\ E_2^{\prime }& ={\kappa }_1{E}_1-{\kappa }_2{E}_2,\\ A^{\prime }& ={\kappa }_2{\rho }_{\mathrm{a}}{E}_2-{\gamma }_{1}A,\\ I_{\mathrm{n}}^{\prime }& ={\kappa }_2(1-{\rho }_{\mathrm{a}})(1-{\rho }_{\mathrm{s}})E_2-{\gamma }_1{I}_{\mathrm{n}},\\ I_{\mathrm{s}}^{\prime }& ={\kappa }_2(1-{\rho }_{\mathrm{a}}){\rho }_{\mathrm{s}}{E}_2-\alpha I_{\mathrm{s}},\\ J^{\prime }& =\alpha I_{\mathrm{s}}-({\gamma }_2+\delta )J,\\ R^{\prime }& ={\gamma }_1(A+I_{\mathrm{n}})+{\gamma }_{2}J,\\ D^{\prime }& =\delta J,\\ C^{\prime }& =\alpha I_{\mathrm{s}}.\end{array}$(2) ) on characterizing the identifiable parameters using uncertainty studies to assess the stability of the parameters estimated. Then, the goal is to evaluate the identifiability using the mean value, the $95\mathrm{\%}$ confidence intervals, and the mean square error (MSE). The step-by-step description of the identifiability analysis is detailed in the next subsection.

2.5.2. Identifiability analysis process

In the study of identifiability, we must follow the following steps that rely on definitions and computations presented in detail in [14, 51]:

1. Assume that $\mathrm{\Theta }=\{\alpha ,\text{β},{\rho }_{\mathrm{s}},p,I_{\mathrm{s}}(0),J(0\left)\right\}$ is the set of parameters to be estimated and ${\mathrm{\Theta }}_0$ is the initial guess.

2. Estimate Θ fitting the model to the dataset generated, using the SA-LSQ procedure with the values of the initial conditions and parameters given. The hat symbol indicates an estimated parameter. Then, we have the set $\hat{\mathrm{\Theta }}=\{\hat{\alpha },\widehat{\text{β}},{\hat{\rho }}_{\mathrm{s}},\hat{p},\widehat{I_{\mathrm{s}}}(0),\hat{J}(0\left)\right\}.$

3. Using the bootstrapping process we obtain the $95\mathrm{\%}$ confidence intervals and their frequency distributions for each parameter estimated.

4. Compute the mean squared error (MSE) $\mathrm{M}\mathrm{S}\mathrm{E}:=\frac{1}{S}\sum _{j=1}^S(\mathrm{\Theta }-{\hat{\mathrm{\Theta }}}_j{)}^2$ for each parameter estimated. Here Θ represents the true parameter value, and ${\hat{\mathrm{\Theta }}}_j$ represents the estimated value of the parameter for the jth bootstrap realization.

5. For each set of estimated parameters, $R_0^{\mathrm{c}}$ is calculated to obtain a distribution of $R_0^{\mathrm{c}}$ values as well.

Following the process described, we say that a model parameter is identifiable if its confidence interval lies in a finite range of values and contains the estimated (and real) parameter value. The preferred situations is that these values have a central position or are close to the mean value. Besides, a small confidence interval indicates that the parameter can be identifiable, while a broader range could indicate a lack of identifiability. When an estimated parameter has a low MSE and a small confidence interval, the we may assert that the parameter is identifiable from the model. On the other hand, larger confidence intervals or larger MSE values may suggest non-identifiability. Here confidence interval is considered to be small if its lower and upper bounds are not smaller or greater, respectively, than 10% of the value estimated. Besides, the smallness of the MSE depends on the nature of the estimated parameter. For instance, if the parameter represents a probability, then we say that MSE is small if it is a value close to zero, but if the parameter represents the quantity of population, one acceptable range can be considered no more than 0.01% of the total population.

Specifically, to analyze structural identifiability, we will fit the entire curves generated for each scenario to recover the values from the parameters assumed in each case. Moreover, to investigate practical identifiability, we will consider different fit times. This is done to address situations when only data of the initial phase (before the peak or maximum value) are available.

3. Results

3.1. Structural identifiability

By applying the identifiability methodology, we obtain results for each one of the steps 1 to 5 outlined in Section 2.5.2. The first type of results corresponds to the parameter estimation process that leads to the best-fit model solution. Secondly, we have 250 replicated datasets applying the bootstrap process to re-estimate model parameters. Figures A1 and A2 (Appendix) show the confidence intervals (represented by vertical lines) and the MSE bar plots for the experiments of Scenario 1, and Figures S1 and S2 (in the supplemental material) display analogous results for the experiments of Scenario 2. Observing the results obtained and registered in these figures, we concentrate our investigation on Tables 3 and 4, in which we classify the parameters estimated between identifiable and non-identifiable. A parameter is considered identifiable with narrow confidence interval, a mean value centred and close to the real value, and the MSE is small. If at least one of these identifiability criteria is not satisfied, then the parameter is considered non-identifiable.

From Table 3 we conclude that all estimated parameters are identifiable for Experiments 1 to 4, i.e. those without initial conditions to fit. For Experiments 5 to 7 in which the initial conditions are also estimated, we note that all the parameters except the initial condition $J\left(0\right)$ are identifiable, maintaining mainly the identifiability of the parameter of interest $R_0^{\mathrm{c}}$. Only Experiment 8 exhibits non-identifiability situations for all parameters except for the dataset generated under the assumption $R_0^{\mathrm{c}}=5$. Here we can suspect that the identifiability improves when the force of infection is greater as in this last case. Nevertheless, we will continue to explore this part later. Another observation is that the estimated $R_0^{\mathrm{c}}$ can be robust to variation or bias in the other estimated parameters.

Analogously, we deduce from Table 4 that the parameters estimated in Experiments 1, 2, 3, and 5 are identifiable. The parameters in Experiments 8, 9, 10, and 12 are also identifiable (except for the initial condition $J\left(0\right)$). For Experiments 4, 7, and 9, the parameter $R_0^{\mathrm{c}}$ is identifiable even though no other parameters are identifiable. Particularly in Experiment 7, all the remaining parameters are non-identifiable. Other situations occur; for example, in Experiment 6, all estimated parameters except those based on data generated with $R_0^{\mathrm{c}}=3$ and p = 0.05 are identifiable (red boxes in column Experiment 6 from Table 4). For Experiment 13, all estimated parameters based on data generated with $R_0^{\mathrm{c}}=4$ and 5 are identifiable. For Experiment 14, the results are not good, but the parameter $R_0^{\mathrm{c}}$ estimated using data generated with p = 0.01 is identifiable. A similar situation occurs with $R_0^{\mathrm{c}}$ in Experiment 11, where this quantity is identifiable only by fitting the data generated with $R_0^{\mathrm{c}}=5$ and p = 0.1.

3.2. Practical identifiability

The results obtained for this part are in Supplemental Material, including coloured tables, namely Tables S2 and S3 that display information analogous to Tables 3 and 4. Tables S2 and S3 allow one to analyze the identifiability. For this study, we investigate the performance of our model to capture the initial phases of an epidemic disease, examining the choices of the fit periods, namely 20, 30, or 40 days. According to Table S2, in Experiment 1 the estimated parameters are identifiable. In Experiments 3 and 5 with fit times 30 and 40 days, the parameters $R_0^{\mathrm{c}}$ are identifiable, even though not all estimated parameters are identifiable. In contrast to this result, in Experiment 2 with a fitting period of 20 days, the parameters that define $R_0^{\mathrm{c}}$ are identifiable, but the computed value $R_0^{\mathrm{c}}$ is not. This situation could be linked with the inverse relation of $R_0$ with parameter α (see Equation (3(3) $R_0^{\mathrm{c}}=\text{β}(\frac{q_{\mathrm{e}}}{{\kappa }_2}+\frac{q_{\mathrm{a}}q{\rho }_{\mathrm{a}}}{{\gamma }_1}+\frac{q(1-{\rho }_{\mathrm{s}})(1-{\rho }_{\mathrm{a}})}{{\gamma }_1}+\frac{q{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{\alpha }+\frac{h{\rho }_{\mathrm{s}}(1-{\rho }_{\mathrm{a}})}{{\gamma }_2+\delta }).$(3) )). Remember that $1/\alpha$ represents the time from symptom onset to isolation or hospitalization. Thus a variation in α of 0.1 to 0.2 implies a change of 5 days from the period of symptoms onset to isolation, a very sensitive parameter. Notably, for this experiment, the confidence intervals and the MSE are plotted in Figure S3, where we can see that these intervals are not very large but lack precision, mainly for the case of a fit period of 20 days, and the epidemic growth is given by $R_0^{\mathrm{c}}=3$ and 4. For this and other experiments, the parameter $R_0^{\mathrm{c}}$ tends to be identifiable when fit periods are greater than 20 days, and $R_0^{\mathrm{c}}\le 4$. The parameter ${\rho }_{\mathrm{s}}$ and the initial condition $I_{\mathrm{s}}\left(0\right)$ and $J\left(0\right)$ for all experiments of Scenario 1 are not identifiable. This suggests a lack of identifiability especially of the parameter ${\rho }_{\mathrm{s}}$ when we have fewer data to fit.

Table S2 also shows that the estimated parameters for Experiment 1 are identifiable. This result differs from Experiments 12, 13, and 14, where all parameters are non-identifiable. For Scenario 2, we can observe that for most experiments with a fit period of 20 days (except for Experiments 1 and 9 with $R_0^{\mathrm{c}}=3$), the parameter $R_0^{\mathrm{c}}$ is non-identifiable. Moreover, all estimated parameters are also non-identifiable, especially for these cases, excluding Experiment 2. These results suggest using more fit data to guarantee identifiability, as it occurs more frequently in experiments with longer fit periods.

Another situation to comment on here is that for most experiments where the quarantine parameter p is estimated, the results mainly for $R_0^{\mathrm{c}}$ result non-identifiable, excluding only Experiment 5 for the fit period 40 days, and $R_0^{\mathrm{c}}$ assumed as 4 and 5, and Experiment 6 again for the fit period 40 days and $R_0^{\mathrm{c}}=3$ and 4. Besides, for some quarantine parameters estimated from Experiments 5 and 6 with fit period 40 days, the result is identifiable, mostly when the fit data are generated assuming p = 0.01. We can think that the lack of identifiability increases when a larger part of the population is quarantined in the initial phase.

4. Application to regional data of Chile

We now apply our methodology to the data collected for some Chilean regions. From the structural and practical results obtained for our model, we apply the study to two cases, one to scan the initial phase (first 30 days) and another for an intermediate phase (first 120 days), considering the data before the change of quarantine criteria applied by the Chilean government. To explain the selection of regions and the phenomena of COVID-19 in Chile, we include the following subsection.

4.1. The early transmission of COVID-19 in Chile

We wish to apply the COVID-19 model of Section 2.1 to some administrative regions of Chile (see Figure 3), where there are regions without quarantine periods, others with more than one quarantine period, a situation that the Chilean government named ‘dynamical quarantines’. This policy was in effect until July 27 2020. In addition, this measure was applied only to municipalities with a greater incidence of positive cases for COVID-19. For this reason, we consider a decoupled analysis to select the regions to analyze until July 27 2020, assuming the day of the first positive case as the first day for the timeline of each region and July 27 2020 is the last day. Since the onset of COVID-19 in Chile, the Ministry of Health has been reporting new cases of infection and death daily for each region and new cases of recovered persons at the national level. The official information (see [24, 41] for details) also includes the daily numbers of polymerase chain reaction (PCR) tests applied, as well as the daily numbers of critical patients and occupancy of intensive care units (ICUs). A timeline of strategies and measures implemented to contain the pandemic is summarized in Table S1. Several times the Chilean Ministry of Health changed the case and death definitions. Table 5 lists the most relevant of these changes along with their dates. Considering the difficulties associated with the changes in case definitions, we decided to include data reported for symptomatic and asymptomatic cases in the COVID-19 model. Total data correspond to the addition of symptomatic and asymptomatic cases reported because an asymptomatic (tested) case may become symptomatic later.

Figure 3. The 16 administrative regions of Chile [63] and their population according to the 2017 census [10]. The Roman numbers are the official administrative numbers of the geographic regions. The greater Santiago area is the Metropolitan region (RM) and is counted as Region 13. Note that numbering is not strictly ordered from north to south; regions 14 to 16 have been created by dividing existing regions. The total population at the end of 2020 is estimated at 19.1 million.

Table 5. Changes in criteria of reported data.

Date	New criterion
29 Apr	Incorporation status ‘asymptomatic,’ i.e. people without symptoms but with PCR positive.
1 Jun	Deaths without a diagnosis but with suspected symptoms of COVID-19 or indeterminate PCR are added to official daily reports.
2 Jun	Report of PCR tests applied for each region is suspended (for approx. a week).
3 Jun	Report of recovered cases is suspended.
7 Jun	653 deaths were added as possible cases of COVID-19 (patients with indeterminate PCR or with symptoms similar to those caused by COVID-19), and 96 deaths of people diagnosed do not have a precise regional distribution in the reports, due to a change in the region of residence
17 Jun	The status ‘without notifying’ is created in the daily reports, summing 31422 cases of people with PCR + who were not reported in the system on diagnosis (reports delayed).

4.2. Experiments applied to regional data

Analyzing the type of information and data published by the Chilean government we identified regions with more than one quarantine period. Then, considering these multiple quarantine periods and the conclusions obtained from our experimental identifiability study, we decided to select only those regions with exactly one and without quarantine periods from the Chilean data. Specifically, we selected those administrative regions that satisfy at least one of these conditions, namely those with

1. exactly one unique quarantine period for the first 30 or 120 days,

2. no quarantine periods in the first 30 or 120 days.

We also prefer data with an increasing number of reported cases, avoiding having consecutive null case counts, i.e. we select robust datasets for capturing epidemic growth. Overall, in light of the experimental results based on simulated data to analyze the (structural and practical) identifiability of the parameters, we know that the best results are achieved for the sets of estimated parameters $(\text{β},R_0^{\mathrm{c}})$ and $(\text{β},J(0),I_{\mathrm{s}}(0),R_0^{\mathrm{c}})$. However, to observe the quarantine parameter, we also consider the parameter sets $(\text{β},p,R_0^{\mathrm{c}})$ and $(\text{β},p,J(0),I_{\mathrm{s}}(0),R_0^{\mathrm{c}})$, for the case that the region has a single quarantine period, even knowing that these parameters are not always identifiable from the experimental results in the initial phase. Table 6 provides more details about parameter identifiability results.

Table 6. Summary of the identifiable parameters concluded from the structural and practical identifiability study used to Chilean data.

	Without Quarantine		With Quarantine
Parameter set	Structural	Practical	Structural	Practical
$(\text{β},R_0^{\mathrm{c}})$	Identifiable	Identifiable	Identifiable	Identifiable
$(\text{β},J(0),I_{\mathrm{s}}(0),R_0^{\mathrm{c}})$	Identifiable	Identifiable	Identifiable	Identifiable
		$T_{\mathrm{f}\mathrm{i}\mathrm{t}}>20$ days		$T_{\mathrm{f}\mathrm{i}\mathrm{t}}>20$ days
$(\text{β},p,R_0^{\mathrm{c}})$			Identifiable	Identifiable $T_{\mathrm{f}\mathrm{i}\mathrm{t}}=40$
				$R_0^{\mathrm{c}}=5$ and p = 0.01
$(\text{β},p,J(0),I_{\mathrm{s}}(0),R_0^{\mathrm{c}})$			Identifiable	Nonidentifiable

Subsequently, we filter the regional data for two periods, one for the first 30 days and the other for 120 days, and we select some of those that satisfy the above mentioned conditions. We have established the following regions for our application.

• Regions 1 (Tarapacá) and 8 (Bío-Bío), without quarantine periods in the initial phase and one quarantine in the final phase.

• Regions 3 (Atacama) and 4 (Coquimbo) were without quarantine during the 120 days.

• Region 10 (Los Lagos) had one quarantine in the initial phase.

Then, with this regional data, we apply our identifiability study to two events,

1. For regions without quarantine, we estimate the parameter sets $(\text{β},R_0^{\mathrm{c}})$ (Experiment 1) and $(\text{β},J(0),I_{\mathrm{s}}(0),R_0^{\mathrm{c}})$ (Experiment 2). We select a fit period of 30 or 120 days depending on the phase analyzed.

2. For regions with quarantine, we estimate the parameter sets $(\text{β},p,R_0^{\mathrm{c}})$ (Experiment 1) and $(\text{β},p,J(0),I_{\mathrm{s}}(0),R_0^{\mathrm{c}})$ (Experiment 2).

We are preserving the same consideration for the time fit.

The left plots in Figures S4 and S5 display the regional data distributions, including some dates of quarantine measures and important holidays.

Establishing the experiments and the Chilean data for the analysis, we need to fix the values of those parameters that are not estimated. In Table 7, we summarize the values assumed for the parameters. For example, we consider the parameters for the time from symptom onset to isolation or hospitalization $1/\alpha$, and the proportions of fully infectious individuals who undergo testing ${\rho }_s$ are averaging. In addition, we assume the quarantine parameters $\lambda =0$ and p = 0 when the experiment does not involve quarantine periods. Alternatively, when there is quarantine, the quarantine parameter p is estimated depending on the region, being the parameter $1/\lambda$ the quarantine duration applied to it (see Table 7).

Table 7. Fixed parameters for running experiments and duration of quarantine periods.

Parameter	Selected for simulations	Source
h	0	Assumed
$q_{\mathrm{e}}$	0.1	[17]
$q_{\mathrm{a}}$	0.4	[17]
q	0.4	Assumed
$1/{\kappa }_1$	2.5	[49, 72]
$1/{\kappa }_2$	2.5	[34, 67]
${\rho }_{\mathrm{a}}$	0.4	[42, 46]
${\rho }_{\mathrm{s}}$	0.5	Assumed
$1/\alpha$	5	Assumed
$1/{\gamma }_1$	7	[49, 71]
$1/{\gamma }_2$	5	[55]
δ	0.021	National death case data
Duration of quarantine periods $1/\lambda$
Region	Initial phase	Final phase
1	0	67
8	0	11
10	21	–

Establishing the experiments and the Chilean data for the analysis, we need to fix the values of those parameters that are not estimated. In Table 7, we summarize the values assumed for the parameters. For example, we consider the parameters for the time from symptom onset to isolation or hospitalization $1/\alpha$, and the proportions of fully infectious individuals who undergo testing ${\rho }_s$ are averaging. In addition, we assume the quarantine parameters $\lambda =0$ and p = 0 when the experiment does not involve quarantine periods. Alternatively, when there is quarantine, the quarantine parameter p is estimated depending on the region, being the parameter $1/\lambda$ the quarantine duration applied to it (see Table 7).

The initial conditions satisfy the considerations assumed in the identifiability study, i.e. $\begin{array}{rl}& R\left(0\right)=Q\left(0\right)=D\left(0\right)=0,E_1\left(0\right)=2J\left(0\right),\\ & E_2\left(0\right)=4J\left(0\right)-E_1\left(0\right),I_{\mathrm{n}}\left(0\right)=min\left\{J\right(0),1\},\\ & S\left(0\right)=N-E_1\left(0\right)-E_2\left(0\right)-I_{\mathrm{n}}\left(0\right)-I_{\mathrm{s}}\left(0\right)-A\left(0\right)-J\left(0\right),\end{array}$ where for the case when $J\left(0\right)$ and $I_{\mathrm{s}}\left(0\right)$ are estimated, these quantities are assumed to satisfy $0<J\left(0\right)<C\left(0\right)$ and $0<I_{\mathrm{s}}\left(0\right)<300$. In those cases where $J\left(0\right)$ and $I_{\mathrm{s}}\left(0\right)$ are assumed, the values $J\left(0\right)=1$ and $I_{\mathrm{s}}\left(0\right)=10$ are utilized.

Table 8. Parameter estimation results for the application with Chilean data.

Exp.	$T_{\mathrm{f}\mathrm{i}\mathrm{t}}$	region	β	p	$I_{\mathrm{s}}\left(0\right)$	$J\left(0\right)$	$R_0^{\mathrm{c}}$ comp	presnorm
without quarantine
1	30	1	1.10459				2.30859	1.40E+02
	30	8	1.55171				3.24307	3.17E+03
	120	3	0.69428				1.45104	5.28E+03
	120	4	0.73344				1.53288	8.88E+04
2	30	1	1.08156		11.23401	1.00000	2.26045	1.40E+02
	30	8	1.12597		52.89334	1.00000	2.35327	2.37E+03
	120	3	0.68112		14.52313	0.99994	1.42353	5.23E+03
	120	4	0.61346		148.95176	2.00000	1.28214	5.63E+04
with quarantine
1	30	10	2.80707	0.08965			5.86678	9.20E+02
	120	1	2.54383	0.05503			5.31660	1.21E+05
	120	8	5.00000	0.54340			10.45000	1.44E+05
2	30	10	3.70074	0.12446	5.43419	0.99998	7.73455	9.11E+02
	120	1	1.50472	0.03014	53.68752	0.61428	3.14486	1.05E+05
	120	8	3.53811	0.41259	62.03312	0.98066	7.39466	1.10E+05

4.3. Application results

Finally, applying our methodology to Chilean regions data selected for the experiments defined, we obtained the parameter estimations (Table 8), the fits plotted in Figures S4 and S5 and Table 9 which we constructed following the same procedure as for Tables 3 and 4. We observe the 95% confidence intervals obtained by the 250 replicates of the bootstrapping process and the MSE (see Figures S6 to S14 in Supplemental Material). According to the coloured table, the parameters β and $R_0^{\mathrm{c}}$ (computed) are usually identifiable, except for Region 10, which has a quarantine period in the first 30 days. Quarantine leads some points to fit, so perhaps our model is unsuccessful in this case, and the estimates generated are non-identifiable. As expected, given the summary of results in Table 6, Region 8 has a short quarantine period in the middle of the final phase, and it turns out difficult to fit data from this region during 120 days to estimate the sets of identifiable parameters. On the other hand, for Region 1, with a prolonged quarantine period in the final phase, conceivably, the initial conditions incorporate more noise in the estimations because it is difficult to estimate the parameters in Experiment 2.

Table 9. Parameter identifiability verification. A green box indicates that the parameter estimated satisfies the parameter identifiability conditions (confidence intervals obtained via bootstrapping are narrow; their mean values are centered; the MSE is small). A red box means the parameter that does not meet the identifiability criteria. Cases not studied are painted in gray.

In short of these two last regions, we can see that the duration of quarantine and the phase analyzed can be decisive when estimating parameters; these examples suggest that we need to explore more combinations of quarantine periods in our COVID-19 model. In addition, the parameters estimated from experiments applied to Regions 3 and 4 are identifiable, as expected, given the structural and practical identifiability demonstrated for our model using our methodology. In the same way, we verify the identifiability obtained for Experiments 1 and 2 applied to Region 8, and Region 1 has identifiable parameters only when we use Experiment 2. However, their confidence intervals are not distant when we apply Experiment 1. The noise from initial conditions is possibly more prevalent in Region 1 than Region 8.

5. Discussion and conclusions

We present and apply a computational approach developed by Roosa and Chowell [51] for simple epidemic models. However, we propose adapting it to a complex compartmental model inspired by the first wave of the COVID-19 outbreak in Chile, where quarantines were declared involving different Chilean regions in early phases. Specifically, we propose a methodology to analyze the structural and practical parameter identifiability for a COVID-19 model. We explore the parameter uncertainty computationally through a parametric bootstrapping approach to achieve that goal. In that exploration, we involve synthetic data generated with the same COVID-19 model to measure the capability of our model to recover the parameters assumed to satisfy the identifiability structural and practical. Those criteria involve observing the $95\mathrm{\%}$ confidence intervals and the MSE (mean square error), constructed from empirical distributions resulting from the bootstrapping process. In this first attempt, we selected some experiments that would allow us to validate the structural and practical parameter identifiability. We consider different variables, for example, the number of parameters to fit, the time to fit, and the quarantine periods. Then, in our methodology, we construct coloured tables (e.g. Table 3) to visualize these qualities. Overall in the results obtained, we observe that the parameter β, and in consequence $R_0^{\mathrm{c}}$, are identifiable (when the rest of the parameters are known). The following best parameter sets in terms of structural identifiability are the sets without initial conditions to estimate, especially in the case without quarantines involved. Wherein the cases with quarantine are more sensitive as a percentage of the population in quarantine increases and the epidemic growth is $R_0<4$, evidence of a lack of identifiability in the model when the quarantine dynamics are considered. Then, we would take quarantine periods studied in the models with caution. Expanding the analysis to the initial phases, we find a significant loss of parameter identifiability. It compared with the structural results. Possibly verifying the conclusion shown in [52] indicates the non-identifiability in pre-peak moments. However, we are interested in analyzing the initial period to recover the control reproduction number. Therefore, we consider applying our computational approach to various initial epidemic periods to evaluate our COVID-19 model. In our experiments, we confirm that the parameter β is recovered, and in the cases when the initial conditions are also estimated and fit time is $T_{\mathrm{f}\mathrm{i}\mathrm{t}}>20$ days (see Tables S2 and S3 in the supplemental material) identifiability is guaranteed. Showing a characteristic dependence between the fit times and the growth rates, besides quarantine periods, which also affect identifiability, but maintain the parameter β identifiable.

Finally, collecting the conclusions from the experiment with synthetic data, we select the parameter sets (β), and ($\text{β},J\left(0\right),I_{\mathrm{s}}\left(0\right)$) for periods without quarantine, and ($\text{β},p$), and ($\text{β},p,J\left(0\right),I_{\mathrm{s}}\left(0\right)$) for periods with quarantine, to fit with Chilean data for some regions. The model parameters are not identifiable in the Chilean regions with a quarantine period in the initial phase. On the other hand, regions without quarantine periods have better results. In particular, when more data ($T_{\mathrm{f}\mathrm{i}\mathrm{t}}>20$) is fit to estimate initial conditions is guaranteed practical identifiability. (see Table 6). For future work, we see the necessity of selecting the fit times for each Chilean region because each one has a particular quarantine declared; besides, the quality data also is variable. The study developed by Freire-Flores et al. in [21] also analyzes the quarantine periods in some Chilean regions using computations of reproduction numbers on different periods characterized by extraordinary government measures. Then we can use these times of fitting with our approach and compare with these results, where, in particular, the results for the parameter β in Region 8 are close.

Our findings further underscore the importance of investigating the identifiability of model parameters through structural and practical identifiability analysis to interpret the results better and avoid misleading model inferences during an epidemic outbreak [22, 51, 56]. On the one hand, structural identifiability analysis should be conducted before the model is applied to real data since it can help guide the modeller about any parameters in the model that are not identifiable even in the scenario that the entire epidemic trajectory has been observed. Parameters deemed not structurally identifiable should be flagged by the modeller and determine whether the model can still be used to infer other quantities of interest. Analyzing a parameter's practical identifiability is unnecessary if it is not structurally identifiable. On the other hand, our results indicate that modellers should be especially careful when calibrating complex epidemic models using limited data of the early epidemic phases, as one or more parameters are often difficult to identify from limited data depending on the model's complexity level. Parameter identifiability occurs more frequently in cases involving more complex epidemic dynamics where parameters change over time due to the implementation and relaxation of social distancing interventions. Unsurprisingly, a parameter is unidentifiable when its function is to model dynamics for an event (e.g. quarantine measures) that has not occurred. As more data on an unfolding epidemic becomes available, model parameters tend to become identifiable, especially when the epidemic has already reached its peak. This highlights the challenges modellers face when using models to generate real-time epidemic predictions as the epidemic evolves.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Availability of data and materials

The data were recorded by the Chilean Government [24], which can be accessed via a GitHub repository supported by Ministery of Science, Technology, Knowledge, and Innovation of Chile [41].

Additional information

Funding

RB is supported by project MATH-Amsud 22-MATH-05 ‘NOTION: NOn-local conservaTION laws for engineering, biological and epidemiological applications: theoretical and numerical’ and from ANID (Chile) through Fondecyt project 1210610; Anillo project ANID/ACT210030; Centro de Modelamiento Matemático (CMM), project FB210005 of BASAL funds for Centers of Excellence; and Centro de Recursos Hídricos para la Agricultura y la Minería (CRHIAM), project ANID/FONDAP/15130015. GC is partially supported by National Science Foundation (NSF) grants #2026797, #2034003, and National Institutes of Health (NIH) R01 GM 130900, and by project ANID/PCI/MEC80170119. IK would like to thank the German Research Foundation (DFG) for financial support of the project within the Cluster of Excellence ‘Data-Integrated Simulation Science’ (EXC 2075), DFG Project 432343452. LYLD is supported by Agencia Nacional de Investigación y Desarrollo (ANID) scholarship ANID-PCHA/Doctorado Nacional/2019-21190640 and Centro de Modelamiento Matemático (CMM), project FB210005 of BASAL funds for Centers of Excellence.

Appendix. Tables and Figures

Figure A1. Summary of results using conditions of experiment of scenario 1 (Parameters (α, β, ${\rho }_{\mathrm{s}}$)). Here and in Figure A2 in per row, the left plot corresponds to the 95% confidence intervals (vertical lines) for the distributions of each estimated parameter obtained through 250 realizations of the synthetic data generated for Scenario 1. Each red point denotes the mean estimated parameter value. The light-blue dashed horizontal line represents the true value (or assumed) for each parameter estimated; the experiments applied to each parameter are fixed on the horizontal axis. The right plot corresponds to the mean squared error (MSE) distribution of parameter estimates considering each experiment and dataset. Finally, Yellow, blue, and green lines or bars are grouped for each experiment and represent each $R_0^{\mathrm{c}}$ assumed for creating the synthetic data, i.e. $R_0^{\mathrm{c}}=3$, 4, and 5, respectively.

Figure A2. Summary of results using conditions of scenario 1 (Parameters ($I_{\mathrm{s}}\left(0\right)$, $J\left(0\right)$, $R_0^{\mathrm{c}}$)).