Recursive state and parameter estimation of COVID-19 circulating variants dynamics | Scientific Reports – Nature.com

Predicting the dynamics of an epidemiological evolution is of utmost importance to control its spread in a population. Modeling by a SIR-based model is conventional in control applications because of its simplicity and real-time applicability. The SARS-CoV-2 virus, however, presents high mutability, which affects model parameters over time. In addition, a relevant percentage of the population has been getting vaccinated in 2021, which also affects those parameters. Both uncertainty sources were irrelevant in the early stages of the COVID-19 pandemic, but their systematic increase make long-term forecasts unreliable. Hence, an accurate system estimate over an extended analysis period requires state feedback. In this work, the state feedback was done by simultaneous state and parameter estimation using an augmented vector.

In this section, a compartmental model is adapted to improve estimation performance considering data availability in Brazil. The model in Equation (1) was adapted from the SIDARTHE model proposed by Giordano et al.30 Hence, it assumes homogeneous states without age structure or the effect of vaccination coverage. The estimated parameters \(\alpha _0\), \(x_c\), and \(x_m\) related to transmissibility, severity of the disease, and lethality, respectively, are described later in this section. These parameters correlate with the dynamics analyzed in COVID-19 vaccine effectiveness studies31,32,33,34. We have selected a federative unit per Brazilian region to evaluate epidemic progression countrywide, but the southeast region, the most populated one, is an exception with two units. Amazonas (AM) was chosen for the north region; Mato Grosso do Sul (MS) for the central-west; Rio Grande do Norte (RN) for the northeast; Rio Grande do Sul (RS) for the south; Rio de Janeiro (RJ) and São Paulo (SP) for the southeast. The total population \(N_i\) from each federative unit i consists of the following compartments:

  • Susceptible (S): individuals prone to infection;

  • Exposed (E): individuals infected in the incubation period, while they are not infectious;

  • Infected (I): undetected asymptomatic individuals;

  • Quarantined (Q): detected asymptomatic individuals who self-quarantine after detecting the disease;

  • Ailed (A): undetected symptomatic individuals;

  • Recognized (R): detected symptomatic individuals who self-quarantine after detecting the disease;

  • Threatened (T): individuals hospitalized in nursery or intensive care units (ICU);

  • Healed detected (\(H_d\)): detected individuals cured without treatment;

  • Deceased (D): individuals deceased due to the disease;

  • Healed with treatment (\(H_t\)): individuals cured after a hospitalization period;

  • Healed undetected (\(H_u\)): individuals cured of the disease without being detected.

$$\begin{aligned} \frac{d S}{dt}&= – \nu S \end{aligned}$$

(1a)

$$\begin{aligned} \frac{d E}{dt}&= \nu S – \rho E \end{aligned}$$

(1b)

$$\begin{aligned} \frac{d I}{dt}&= p \rho E – (\lambda +\varepsilon ) I \end{aligned}$$

(1c)

$$\begin{aligned} \frac{d Q}{dt}&= \varepsilon I – \lambda _{d} Q \end{aligned}$$

(1d)

$$\begin{aligned} \frac{d A}{dt}&= (1-p) \rho E – (\theta +\mu +\kappa ) A \end{aligned}$$

(1e)

$$\begin{aligned} \frac{d R}{dt}&= \theta A – (\mu _{d}+\kappa _{d}) R \end{aligned}$$

(1f)

$$\begin{aligned} \frac{d T}{dt}&= \mu A +\mu _{d} R – (\sigma +\tau ) T \end{aligned}$$

(1g)

$$\begin{aligned} \frac{d H_d}{dt}&= \lambda _{d} Q + \kappa _{d} R \end{aligned}$$

(1h)

$$\begin{aligned} \frac{d D}{dt}&= \tau T \end{aligned}$$

(1i)

$$\begin{aligned} \frac{d H_{t}}{dt}&= \sigma T \end{aligned}$$

(1j)

$$\begin{aligned} \frac{d H_{u}}{dt}&= \lambda I + \kappa A \end{aligned}$$

(1k)

where all states are fractions of a total population \(N_i\), informed by the Ministry of Health of Brazil2. \(\nu\) is the infection rate, \(\rho\) is the incubation rate, and p is the fraction of infected individuals who remain asymptomatic. \(\varepsilon\) and \(\theta\) are the detection rates of I and A, respectively. \(\lambda\), \(\lambda _{d}\), \(\kappa\), \(\kappa _{d}\), and \(\sigma\) are the recovery rates of I, Q, A, R, and T, respectively. \(\tau\) is the mortality rate, whereas \(\mu\) and \(\mu _{d}\) are severe illness rates of A and R, respectively. Fig. 1 shows a scheme of the state transitions.

Figure 1
figure 1

Schematic diagram of the proposed compartmental model.

The analysis of several cases studies within a larger region provides spatial dynamics concerning virus spread. The chosen federative units for the study are known to be heterogeneous among each other35,36,37 since Brazil is a large country where there were several different outbreaks dates, local government policies, and population behavior. Brazilian spatial epidemic progression studies35,36 indicated multiple initial outbreaks spread progressively to neighboring territories. The numerous government policies and population behavior are pointed out by the higher variance of the first wave duration of Brazilian states when compared to the United States and India variances37.

The healed compartment from the Giordano et al. model30 is subdivided into three compartments: \(H_{d}\), \(H_{u}\), and \(H_{t}\). \(H_{t}\) is a measurable state by the Brazilian severe acute respiratory syndrome (SARS) database38,39. \(H_{d}\) is an unmeasured state because the Brazilian SARS database accounts only for the hospitalized individuals, and the Ministry of Health of Brazil only provides recovered estimate countrywide. However, the cumulative confirmed cases provided by the latter are composed mainly of \(H_{d}\) for any analysis post the first wave. \(H_{u}\) is an unmeasured state containing most post-infection individuals for all studied federative units. The closed system assumption of the compartmental model leads to the following constraint: \(S + E + I + Q + A + R + T + H_d + D + H_t + H_u = 1\); thus, we substituted Equation (1k) by Equation (2).

$$\begin{aligned} H_{u} = 1 – S – E – I – Q – A – R – T – H_d – D – H_t \end{aligned}$$

(2)

The additional state E corresponds to a natural time delay of the system, which is usual in control-oriented models8,11 and forecasts to a lesser extent20,40,41. Presymptomatic infected individuals are within this state as \((1-p)E\), but their infection rate is assumed insignificant to simplify parameter estimation.

Reinfections play a significant role in the resurgence of COVID-19 infection waves since new lineage might evade immunity from previous infections42. Gamma43 and Delta44 mutations allow them to infect individuals recovered from other variants. Hence, reinfections from natural immunity decrease are assumed negligible to the emergence of another variant. The latter, however, can not be forecasted as they happen in occasional events. Hence, S in the model Equation (1) is unconnected with healed compartments \(H_u\), \(H_t\), and \(H_d\), and the reinfection dynamics are assumed to be comprised in the state estimation.

The infection rate is simplified into a single parameter to guarantee observability. Hence, infections caused by presymptomatic and detected infected are assumed to be insignificant compared to infections caused by undetected infected. In addition, the same infection rate is applied to symptomatic and asymptomatic, although the first is acknowledged as more infectious45.

$$\begin{aligned} \nu = \alpha (I + A) \end{aligned}$$

(3)

where \(\alpha\) is the contagion rate, consisting of the probability that a susceptible individual contracts the disease from possible contact with an infectious individual. It is a function of non-pharmaceutical interventions (NPI), vaccination coverage, and circulating variants. NPI and vaccination mitigate virus spread in the short-term, while virus mutations might affect its transmissibility, as happened for the Gamma43 variant.

NPI dynamics are inserted into a compartmental model by time-varying functions20, independent variables9,11,12, or time-varying parameters estimated over time24,25. We focused on these last two as they are better suited to a control-oriented model. First, we separated social distancing from other NPI by defining \(\alpha\) according to Equation (4).

$$\begin{aligned} \alpha = \alpha _{0} (1 – u) \end{aligned}$$

(4)

where \(u \in [0,1]\) is the manipulated variable related to social distancing and \(\alpha _{0}\) is the estimated contagion rate. The linearity applied over \(\alpha\) and u in Equation (4) gets the correct direction between the contagion rate and social distancing. NPI unrelated to social distancing (e.g., mass gathering restrictions and mask requirements) are comprised in \(\alpha _{0}\).

Social distancing is measurable by Google mobility data as percentage changes concerning a baseline defined from data sets before the COVID-19 outbreak46. Google mobility data are divided into six categories: recreation, essentials, parks, transit, workplace, and resident. A linear combination among the two most independent categories is used to define u. The similarity was measured by a zero-lag cross-correlation matrix through data from all federative units studied between February 2020 and July 2021. The normalized cross-correlation, whose results are presented in Supplementary Table S1, was calculated using the xcorr function from MATLAB. The absolute difference from zero characterizes the similarity between two signals, where independence is defined. The essentials signal had a cross-correlation closer to zero for all categories except itself; however, it is a monthly periodic signal while the others are weekly reported. Hence, the cross-correlation closest to zero, disregarding essentials, is related to parks and workplace; thus, they were selected to define u. Additionally, u was limited in the range [0,1], assuming each mobility category has lower and upper bounds on -100% and 100%, respectively. A weighted sum to assimilate location-dependent correlations concerning each mobility category was used to evaluate u according to Equation (5).

$$\begin{aligned} \begin{aligned} u(t)&=\dfrac{-w_u\left\langle Parks(t)\right\rangle -(1-w_u)\left\langle Workplace(t)\right\rangle -(w_u Parks_{min}+(1-w_u) Workplace_{min})}{w_u (Parks_{max}-Parks_{min})+(1-w_u) (Workplace_{max}-Workplace_{min})}\\&= \dfrac{-w_u \left\langle Parks(t)\right\rangle -(1-w_u)\left\langle Workplace(t)\right\rangle -(w_u (-100)-(1-w_u)(-100))}{w_u(100-(-100))+(1-w_u)(100-(-100))}\\&= \dfrac{-w_u \left\langle Parks(t)\right\rangle -(1-w_u) \left\langle Workplace(t)\right\rangle +100}{200} \end{aligned} \end{aligned}$$

(5)

where \(\left\langle Parks(t)\right\rangle\) and \(\left\langle Workplace(t)\right\rangle\) are weekly moving averages of parks and workplace, respectively, and \(w_u\) is the relative weight concerning parks mobility, which is an additional parameter to be estimated in the model identification.

The fraction of individuals who do not experience symptoms is defined as \(p \in [0.15,0.7]\) according to the U.S. Centers for Disease Control and Prevention (CDC)47. Virus mutations, testing policies, and different age distribution explain the broad range. The parameter p is not estimated over time because it is not observable from available data since there is no classification of symptomatic and asymptomatic. Hence, we defined \(p=0.5\) as an intermediate value whose error is mitigated by the state estimator with estimations of I and A. The vaccine efficacy against infection is correlated to both \(\alpha _{0}\) and p since it only measures symptomatic cases, according to U.S. Food and Drug Administration (FDA).

Following the vaccine efficacy against severe and mild diseases, a parameter \(x_c\) is defined as the fraction of symptomatic individuals who develop severe or mild symptoms. Assuming that severe and mild illnesses imply hospitalization, thus \(x_c\) is the fraction of individuals moving from A and R to T. Summing up Equations (1e) and (1f):

$$\begin{aligned} \frac{d(A+R)}{dt}= (1-p) \rho E – (\mu +\kappa ) A – (\mu _{d}+\kappa _{d}) R \end{aligned}$$

which is simplified by assuming \(\mu \approx \mu _{d}\) and \(\kappa \approx \kappa _{d}\) to:

$$\begin{aligned} \frac{d(A+R)}{dt}= (1-p) \rho E – (\mu +\kappa ) (A + R) \end{aligned}$$

Thus:

$$\begin{aligned} x_{c} = \frac{\mu }{\mu + \kappa } \end{aligned}$$

Let us rewrite \(\mu\) and \(\kappa\) as a probability function of the symptomatic individual to follow their ways, then:

$$\begin{aligned} x_{c} = \dfrac{(1-x_k-x_{\theta }) {\tilde{\mu }}}{(1-x_k-x_{\theta }) {\tilde{\mu }} + x_k {\tilde{\kappa }}} \end{aligned}$$

(6)

where \(x_k\) and \(x_{\theta }\) are the probabilities of an individual in A to recover or to get detected, respectively. The average rates of severe illness \({\tilde{\mu }}\) and symptomatic recovery \({\tilde{\kappa }}\) correspond to properties studied in the literature. In this work, we defined \({\tilde{\mu }} = 1/5\ \text {d}^{-1}\) and \(\rho = 1/5.2\ \text {d}^{-1}\) from CDC47, and \({\tilde{\kappa }}\) was based on a study of the detection window and test sensitivity of IgG/IgM tests48. The testing rate is a local and time-dependent property that affects both the probabilities \(x_k\) and \(x_{\theta }\). Let us define the correlated parameter \(x_s\) as the fraction of recovered undetected individuals. We have from Equation (1e):

$$\begin{aligned} x_{s} = \dfrac{\kappa }{\mu + \theta + \kappa } = \dfrac{ x_k {\tilde{\kappa }}}{(1-x_k-x_{\theta }) {\tilde{\mu }} + x_{\theta } {\tilde{\theta }} + x_k {\tilde{\kappa }}} \end{aligned}$$

(7)

Rewriting Equation (7) for \(x_k\):

$$\begin{aligned} x_{k} = \dfrac{x_s x_{\theta } {\tilde{\theta }} + (1-x_{\theta }) x_s {\tilde{\mu }} }{(1-x_s) {\tilde{\kappa }} + x_s {\tilde{\mu }}} \end{aligned}$$

and substituting it in Equation (6) rewritten for \(x_{\theta }\):

$$\begin{aligned} x_{\theta } = \dfrac{(1-x_c-x_s) {\tilde{\kappa }} {\tilde{\mu }}}{(1-x_c-x_s) {\tilde{\kappa }} {\tilde{\mu }} + (1-x_c) x_s {\tilde{\theta }} {\tilde{\mu }} + x_c x_s {\tilde{\theta }} {\tilde{\kappa }}} \end{aligned}$$

(8)

Considering that \(x_{c}\) and \(x_{s}\) represent fractions of the symptomatic infected, then \(x_{s}+x_c \in [0,1]\). Locations with a steadier testing policy could estimate \(x_s\) as a constant. However, rapid tests and RT-PCR were not available in public health services in the early stages of the COVID-19 pandemic in Brazil. Defining \(x_s\) as a logistic equation in the function of time according to:

$$\begin{aligned} x_{s} = a_{x_{s}} \left( 1-\frac{a_{\zeta }}{1+\exp {\left( -b_{\zeta } (t – c_{\zeta })\right) }}\right) \end{aligned}$$

(9)

where \(a_{x_{s}}\), \(a_{\zeta }\), \(b_{\zeta }\), and \(c_{\zeta }\) are identified model parameters. The definition of Equation (9) is based on heuristics that \(x_s\) is initially high and decreases progressively to a steady state following test availability to the population. These model parameters also comprises uncertainties regarding test policy.

Analogous to \(x_s\), we define the fraction of recovered undetected asymptomatic individuals \(x_{a}\) from Equation (1c) as:

$$\begin{aligned} x_{a} = \frac{\lambda }{\lambda +\varepsilon } = \frac{x_{id} {\tilde{\lambda }}}{x_{id} {\tilde{\lambda }} + (1- x_{id}) {\tilde{\varepsilon }}} \leftrightarrow x_{id} = \frac{x_{a} {\tilde{\varepsilon }}}{{\tilde{\lambda }} + x_{a} ({\tilde{\varepsilon }} – {\tilde{\lambda }})} \end{aligned}$$

where \(x_{id}\) is the probability of detecting the disease in an asymptomatic individual. Defining \(x_a\) similarly to \(x_s\), we have:

$$\begin{aligned} x_{a} = a_{x_{a}} \left( 1-\frac{a_{\zeta }}{1+\exp {\left( -b_{\zeta } (t – c_{\zeta })\right) }}\right) \end{aligned}$$

(10)

where \(a_{x_{a}}\) is an additional identified model parameter. Equations (9) and (10) have linear dependence between \(x_{s}\) and \(x_{a}\) to avoid overfitting of an excessive number of model parameters. The ratio \(a_{x_{a}}/a_{x_{s}}\) comprises the effect of the viral load on the test sensitivity and the test probability between symptomatic and asymptomatic infected individuals. Related uncertainties are assumed to be mitigated by state estimation among the states I, Q, A, and R.

Finally, we define a parameter \(x_m\) analogous to vaccine efficacy against lethality as the fraction of threatened individuals who decease. We define it from Equation (1g) as:

$$x_{m} = \frac{{\mathbf{\tau }}}{{\sigma + {\mathbf{\tau }}}}$$

(11)

Rewriting Equation (11) as a function of a death probability \(x_e\):

$$\begin{aligned} x_{m} = \frac{x_e {\tilde{\mathbf{\tau } }}}{(1-x_e){\tilde{\mathbf{\sigma} }} + x_e {\tilde{\mathbf{\tau } }}} \end{aligned}$$

and isolating \(x_e\) give us:

$$\begin{aligned} x_{e} = \frac{x_m {\tilde{\mathbf{\sigma} }}}{(1-x_m){\tilde{\mathbf{\sigma} }} + x_m {\tilde{\sigma }}} \end{aligned}$$

(12)

where \({\tilde{\sigma }}\) is the average recovery rate from hospitalization and \({\tilde{\tau }}\) is the average mortality rate. These parameters depend on healthcare demand, medical resources, notification delay, virus mutations, vaccine coverage, and testing policy. Nonetheless, they are simplified as constants to allow future estimations since, by assumption, uncertainties are mitigated by the state and parameter estimation.

The definition of parameters equivalent to vaccine efficacy against transmissibility, severity of the disease, and lethality as functions of state transition rates give comprehensive information about the virus spreading dynamics. The model uncertainties are outweighed by better parameter estimations by considering a fewer number of estimated parameters. The definition of \(x_c\) and \(x_m\) yields additional flexibility in the model formulation. Minor changes applied over \(\alpha _0\), \(x_m\), and \(x_c\) can express specific vaccine dynamics on the model. Hence, their definition comprehends an alternative implementation of vaccination in compartmental modeling.

The Ministry of Health of Brazil2 provides accumulated data on confirmed cases, deceased, and their respective incidences for each federative unit and county. The Brazilian SARS database38,39 provides clinical data from patients with a severe acute respiratory syndrome which comprehend confirmed and suspected cases of COVID-19 and other diseases. It notifies the period of hospitalization, evolution date, the confirmation status of COVID-19, among other information. Summing up all confirmed COVID-19 patients per each federative unit i gives observability on \(T_i\) and \(H_{t,i}\). In addition, overall means of hospitalized evolution between April 2021 and July 2021 were used to define \({\tilde{\sigma }}\) and \({\tilde{\tau }}\). Both databases are daily measured; hence sampling time \(T_s = 1\) d. Average testing rates \({\tilde{\varepsilon }}\) and \({\tilde{\theta }}\) are location-dependent; however, we assumed that correlated uncertainties are comprehended in \(a_{\zeta ,i}\), \(b_{\zeta ,i}\) and \(c_{\zeta ,i}\). Hence, we defined \({\tilde{\varepsilon }} = {\tilde{\theta }}\) = \(1\ \text {d}^{-1}\) to suit sampling time. In summary, the monitored variable \(\mathbf {y_i}\) is defined as:

$$\begin{aligned} \mathbf {y_i}(k) = \mathbf {h}(\mathbf {x}_{\mathbf {i}}) = \begin{bmatrix} Q_i(k)+R_i(k)+T_i(k)+H_{d,i}(k) + D_i(k) + H_{t,i}(k)\\ D_i(k)\\ T_i(k)\\ H_{t,i} (k)\end{bmatrix} \end{aligned}$$

(13)

where \(\mathbf {x_i} = \left[ S_i\ E_i\ I_i\ Q_i\ A_i\ R_i\ T_i\ H_{d,i}\ D_i\ H_{t,i} \right] ^T\).

EPICOVID19-BR provides additional data over temporal distributions in Brazil. It surveyed COVID-19 prevalence in cities from all regions on different timelines49,50. Let us consider the prevalence estimations from federative units given by Marra and Quartin51 based on three phases of EPICOVID19-BR. Furthermore, if we assume \({\tilde{\lambda }}={\tilde{\kappa }}=1/15\ \text {d}^{-1}\), then we can correlate states \(H_{d,i}\) and \(H_{u,i}\) with test sensitivity. EPICOVID19-BR did not test hospitalized patients49 and used an IgM and IgG antibody test more sensitive 15 days after the appearance of symptoms48. Hence, the state transition model \(\mathbf {f}(\mathbf {x_i},u_i)\) for each federative unit i is defined in Equation (14).

$$\begin{aligned} \begin{aligned} \frac{d S_i}{dt}&= – \nu _i S_i \\ \frac{d E_i}{dt}&= \nu _i S_i – \rho E_i \\ \frac{d I_i}{dt}&= p \rho E_i – (\lambda _i+\varepsilon _i) I_i \\ \frac{d Q_i}{dt}&= \varepsilon _i I_i – \lambda _i Q_i \\ \frac{d A_i}{dt}&= (1-p) \rho E_i – (\theta _i+\mu _i+\kappa _i) A_i \\ \frac{d R_i}{dt}&= \theta _i A_i – (\mu _i+\kappa _i) R_i \\ \frac{d T_i}{dt}&= \mu _i (A_i + R_i) – (\sigma _i+\tau _i) T_i\\ \frac{d H_{d,i}}{dt}&= \lambda _i Q_i + \kappa _i R_i \\ \frac{d D_i}{dt}&= \tau _i T_i\\ \frac{d H_{t,i}}{dt}&= \sigma _i T_i\\ H_{u,i}&= 1 – S_i – E_i – I_i – Q_i – A_i – R_i – T_i – H_{d,i} – D_i – H_{t,i}\\ \nu _i&= \alpha _{0,i} (1 – u_i) (I_i + A_i)\\ x_{a,i}&= a_{x_{a,i}} \left( 1-\frac{a_{\zeta ,i}}{1+\exp {\left( -b_{\zeta ,i} (t – c_{\zeta ,i})\right) }}\right) ,\ x_{s,i} = a_{x_{s,i}} \left( 1-\frac{a_{\zeta ,i}}{1+\exp {\left( -b_{\zeta ,i} (t – c_{\zeta ,i})\right) }}\right) ,\\ x_{id,i}&= \frac{x_{a,i} {\tilde{\varepsilon }}}{{\tilde{\lambda }} + x_{a,i} ({\tilde{\varepsilon }} – {\tilde{\lambda }})}\\ x_{k,i}&= \frac{x_{s,i} x_{\theta ,i} {\tilde{\theta }} + (1-x_{\theta ,i}) x_{s,i} {\tilde{\mu }} }{(1-x_{s,i}) {\tilde{\kappa }} + x_{s,i} {\tilde{\mu }}},\\ x_{\theta ,i}&= \frac{(1-x_{c,i}-x_{s,i}) {\tilde{\kappa }} {\tilde{\mu }}}{(1-x_{c,i}-x_{s,i}) {\tilde{\kappa }} {\tilde{\mu }} + (1-x_{c,i}) x_{s,i} {\tilde{\theta }} {\tilde{\mu }} + x_{c,i} x_{s,i} {\tilde{\theta }} {\tilde{\kappa }}},\ x_{e,i} = \frac{x_{m,i} {\tilde{\sigma }}_i}{(1-x_{m,i}){\tilde{\tau }}_i + x_{m,i} {\tilde{\sigma }}_i}\\ \lambda _i&= (1-x_{id,i}) {\tilde{\lambda }},\ \varepsilon _i = x_{id,i} {\tilde{\varepsilon }},\ \theta _i = x_{\theta ,i} {\tilde{\theta }},\ \kappa _i = x_{k,i} {\tilde{\kappa }},\ \mu _i = (1- x_{k,i} – x_{\theta ,i}) {\tilde{\mu }},\\ \sigma _i&= (1-x_{e,i}) {\tilde{\sigma }}_i,\ \tau _i = x_{e,i} {\tilde{\tau }}_i \end{aligned} \end{aligned}$$

(14)

Each federative unit i under study has prevalence distribution from EPICOVID19-BR formulated as Equation (15) for each phase \(j\in \{1,2,3\}\), assuming a test sensitivity of 100% during \(T_{p}\) days followed by a sudden decay to zero.

$$\begin{aligned} \begin{aligned} Pre_{i,j,min} \le \frac{1}{N_{total,j}} \sum _{k=0}^{N_{total,j}-1} \left( H_{all,i}(N_{ep,j}+k) – H_{all,i}(N_{ep,j}+k-T_{p}))\right) \le Pre_{i,j,max} \end{aligned} \end{aligned}$$

(15)

where prevalence bounds \(Pre_{i,j,min}\) and \(Pre_{i,j,max}\) can be found in Supplementary Table S251, and \(T_{p} = 50\) was the arbitrated value for the detection window. \(t(N_{ep,1})=\text {14 May 2020}\), \(t(N_{ep,2})=\text {4 June 2020}\), \(t(N_{ep,3})=\text {21 June 2020}\) are initial dates from the first, second and third phases of EPICOVID19-BR, respectively, while \(N_{total,1}=8\) and \(N_{total,2}=N_{total,3}=4\) correspond to their respective duration in days, and \(H_{all,i} = H_{u,i}+H_{d,i}+H_{t,i}\). In the early stages of the pandemic outbreak, recovered individuals are approximately null; thus, we defined \(H_{d,i}(N_{ep,j}+k-T_{p})=H_{t,i}(N_{ep,j}+k-T_{p})=H_{u,i}(N_{ep,j}+k-T_{p})=0 \forall \left\{ N_{ep,j}+k<T_{p}| j \in \{1,2,3\} \right\}\).

Gene sequences reported in GISAID7 indicate Zeta variant appearance in mid-October 2020. Hence, the identification step is bounded at \(t(N_f) = \text {1 October 2020}\) to guarantee the steady circulation of variants B.1.1.28 and B.1.1.33. The lower bound aims at an imported infection neglectful in the system when \(\left[ I_i(t_{0,i})\ Q_i(t_{0,i})\ A_i(t_{0,i})\ R_i(t_{0,i})\ T_i(t_{0,i})\ H_{d,i}(t_{0,i})\ D_i(t_{0,i})\ H_{t,i}(t_{0,i})\ H_u(t_{0,i}) \right] ^T \approx {{\textbf {0}}}\). Therefore, only \(S_i (t_{0,i})\) and \(E_i (t_{0,i})\) were considered optimization variables for the model identification. Hospitalized individuals \(T_i\) were used to define \(\{N_{0,i}|t_{0,i} = t(N_{0,i})\}\) from the solution of a system composed by \(T_i(N_{0,i})>0.00003\), \(T_p-N_{ep,3}-N_{0,i} \ge 0\), \(N_{0,i} \in {\mathbb {N}}\), for each federative unit i. The model identification is evaluated through an integral time-squared error performance criteria for reducing the contribution of the initial error of imported infections.

The nonlinear optimization problem in Equation (16) was solved for each federative unit i with IPOPT52 via CasADI/MATLAB53.

$$\begin{aligned} \begin{aligned} \min _{\mathbf {ident_i}}&\quad \sum _{k=N_{0,i}}^{N_f} (k-N_{0,i}+1) \left\Vert \mathbf {y_i}(k)-\mathbf {z_i}(k) \right\Vert _{\mathbf {Q_{id,i}}}^2\\ \text {Subject }&\text {to Equation}~(15)\text { and:}\\&\quad \mathbf {x_i}(k+1) = \mathbf {x_i}(k) + \int _{k}^{k+1} \mathbf {f}(\mathbf {x_i}(t),u_i(t)) \, d{t}\\&\quad \mathbf {y_i}(k) = \mathbf {h}(\mathbf {x_i}(k))\\&\quad S_i(t_{0,i}) \in [0.9,1],\ E_i(t_{0,i}) \in [0,0.1],\ \alpha _{0,i} \ge 0,\ a_{x_{a,i}} \in [0,0.1],\ a_{x_{s,i}} \in [0,0.1],\ x_{c,i} \in [0,1],\ x_{m,i} \in [0,1]\\&\quad a_{\zeta ,i}\ \in [0,1],\ b_{\zeta ,i} \in [0,0.25],\ c_{\zeta ,i} \ge 0,\ x_{a,i} \in [0,1],\ x_{s,i} \in [0,1],\ x_{s,i} + x_{c,i} \in [0,1],\ x_{a,i} \ge x_{s,i} \end{aligned} \end{aligned}$$

(16)

where \(\mathbf {ident_i} = \left[ S_i(t_{0,i})\ E_i(t_{0,i})\ \alpha _{0,i}\ a_{x_{a,i}}\ a_{x_{s,i}}\ x_{c,i}\ x_{m,i}\ a_{\zeta ,i}\ b_{\zeta ,i}\ c_{\zeta ,i}\ w_{u,i}\right] ^T\), \(\mathbf {z_i}\) are the measured variables and \(\mathbf {Q_{id,i}} \in {\mathbb {R}}^{4\times 4}\) is a weight matrix calculated to normalize measurements from the early stages of the pandemic to July 2021 according to Equation (17). All numerical integration in the state estimators were solved with CVODES54 via CasADI/MATLAB. The initial guess was set as \(\mathbf {ident_{0,i}} = [0.95\ 0.05\ 0.1\ 0.9\ 0.9\ 0.02\ 0.1\ 0.1\ 0.1\ 20\ 0.5]^T\). Results and location-dependent parameters are shown in Table 1.

$$\begin{aligned} \mathbf {Q_{{id,i}_{j,j}}} = \dfrac{1}{\left( y_{j,max}-y_{j,min}\right) ^2}, j \in \{1,2,3,4\} \end{aligned}$$

(17)

Table 1 Initial states and parameters of the proposed model for each federative unit i.

Source