Recursive state and parameter estimation of COVID-19 circulating variants dynamics.

COVID-19 pandemic response with non-pharmaceutical interventions is an intrinsic control problem. Governments weigh social distancing policies to avoid overload in the health system without significant economic impact. The mutability of the SARS-CoV-2 virus, vaccination coverage, and mobility restriction measures change epidemic dynamics over time. A model-based control strategy requires reliable predictions to be efficient on a long-term basis. In this paper, a SEIR-based model is proposed considering dynamic feedback estimation. State and parameter estimations are performed on state estimators using augmented states. Three methods were implemented: constrained extended Kalman filter (CEKF), CEKF and smoother (CEKF & S), and moving horizon estimator (MHE). The parameters estimation was based on vaccine efficacy studies regarding transmissibility, severity of the disease, and lethality. Social distancing was assumed as a measured disturbance calculated using Google mobility data. Data from six federative units from Brazil were used to evaluate the proposed strategy. State and parameter estimations were performed from 1 October 2020 to 1 July 2021, during which Zeta and Gamma variants emerged. Simulation results showed that lethality increased between 11 and 30% for Zeta mutations and between 44 and 107% for Gamma mutations. In addition, transmissibility increased between 10 and 37% for the Zeta variant and between 43 and 119% for the Gamma variant. Furthermore, parameter estimation indicated temporal underreporting changes in hospitalized and deceased individuals. Overall, the estimation strategy showed to be suitable for dynamic feedback as simulation results presented an efficient detection and dynamic characterization of circulating variants.

Introduction

The first official cases of COVID-19 were dated in December 2019 in Wuhan, China. Its spread worldwide in later months resulted in a pandemic classification from the World Health Organization (WHO) on 11 March 2020[1]. The first official case in Brazil was reported two weeks earlier, on 26 February 2020, from a man returning to São Paulo from Italy[2]. Social distancing presents effective mitigation over virus spread[3]; however, it generates negative impacts on the economy and on the mental health of the population[4]. Vaccination is a control action that progressively reduces virus transmissibility aiming for disease elimination defined by zero community infections[5]. Nonetheless, vaccination coverage is delimited by the vaccine acceptance rate, which makes its goal unfeasible even if vaccination provides 100% efficacy against transmission.

The SARS-CoV-2 virus is highly mutable, with thousands of variants documented since its origin in December 2019[6]. Mutations might change system dynamics; thus, model updating is required for reliable predictions. Genomic surveillance of SARS-CoV-2 virus in Brazil indicated four predominant circulating variants from February 2020 to July 2021. B.1.1.33 and B.1.1.28 were predominant from the pandemic beginning to September 2020; the Zeta variant (P.2), which originated in Rio de Janeiro, was predominant from October 2020 to February 2021; and the Gamma variant (P.1), which originated in Amazonas, was predominant from mid-February 2021 to July 2021[7]. In the pandemic modeling, the dynamics from each variant correspond to a set of parameters that must be estimated to ensure an accurate prediction over an extended period of analysis.

Modeling epidemiological evolution by a compartmental model is standard for control-oriented models since its simplicity suits real-time applications[8-12]. Optimality is usually defined to mitigate virus spread within health system capacity, while an input or manipulated variable is correlated to contagion rate. The input variable is discrete for a definition based on previously implemented government restrictive measures[8,10,11]; and continuous for a definition based on mobility data[12] or possible government measures (e.g., complete lockdown and no countermeasures)[9]. Nonetheless, model parameters are not constrained to functions of the manipulated variables. Olivier et al.[8] defined several compartmental model parameters as time-varying functions. Köhler et al.[9] described hospitalized parameters as a function of the state variables. Morato et al.[10] used a three-step parameter estimation of contagion, recovery, and mortality rates.

Mobility data is regarded generally in models focused on the forecast[13-19]. There are many available database, among which there are data related to local infection probability and restrictive government measures. For instance, SafeGraph details node-related information fit for a network model[13]. Facebook details geographic movement metrics suitable for spatial models[14,15]. Apple and Google detail location mobility trends correlated to government measures which are used coupled[16,17] or standalone[18,19]. In addition, mobility dynamics might be identified and described on time-varying functions from previous government measures[20,21].

Virus mutations and vaccination coverage affect model dynamics, adding uncertainties to the system. There are examples of recursive state and parameter estimation applications in the COVID-19 pandemic. Sun et al.[22] estimated parameters at each discrete time with a grid search, while Menda et al.[23] estimated them with a neural network. Liao et al.[24] and Morato et al.[10] estimated parameters with moving horizon estimation based on a least-square method followed by a regressive method; moreover, the latter[10] proposed an additional moving average on the estimation structure. Tsay et al.[25] estimated unmeasured states with an unscented Kalman filter. Zhu et al.[26] estimated states and parameters into an augmented state with an extended Kalman filter (EKF). Song et al.[27] estimated states with an EKF and parameters with a proposed strategy based on maximum likelihood. State and parameter estimations in the literature focused on overall system dynamics. The authors used estimation strategies to estimate unmeasured states, capture reinfection dynamics or adjust model parameters for more accurate estimations. Hence, virus mutations and vaccination dynamics have not been study objects with similar estimation strategies.

Transmissibility, severity of the disease, and lethality are three properties of interest for study in an epidemiological model. They are defined by the probability of an infected individual moving from one given compartment to another. Marziano et al.[28] proposed an age-structured model to analyze the Italia epidemic evolution during its first wave for possible outcomes from easing restrictive measures. The transmissibility was a function of google mobility data, the probability of developing severe disease was a fitting parameter per age group, and the lethality was defined as a function of the latter and hospitalized data. Kemp et al.[21] proposed a compartmental model with fitting parameters for each probability of split in the model configuration to analyze herd immunity in Austria, Luxemburg, and Sweden. The transmissibility was a function of mobility fitted for each previous government measure, while other parameters were fitted as constants for each wave of COVID-19 infection.

In this work, we propose a comprehensive compartmental model for detecting epidemiological dynamics in terms of transmission, severity of the disease, and lethality equivalent to vaccine efficacy studies. Classical vaccination coverage modeling through a SIR-based model supposes the vaccinated state is 100% immune to reinfection; however, recent literature contradicted this assumption[29]. The modeling through correlated vaccination parameters is an alternative formulation to comprehend the vaccine dynamics. Hence, it is suitable for analyzing vaccine efficacy or intervention measures. The proposed model considers a recursive estimation approach in which simplifying assumptions focuses on detecting the aforementioned dynamics with parameter estimation. Model accuracy is improved using temporal prevalence distributions from the seroprevalence survey EPICOVID19-BR in the first wave of the pandemic. The proposed estimation strategy identifies COVID-19 variant emergence and characterizes its dynamics on epidemiological evolution based on dynamic feedback, which is suitable for online applications.

First, we describe the proposed compartmental model assumptions and parameter identification of the first wave of the pandemic. Then, we describe the implemented state and parameter estimation strategies. Next, we show numerical results from simulations on several Brazilian federative units and analyze the estimated parameters. Finally, we make our conclusions and discuss possible future works.

Mathematical model

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 , , and 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 studies[31-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 from each federative unit i consists of the following compartments: where all states are fractions of a total population , informed by the Ministry of Health of Brazil[2]. is the infection rate, is the incubation rate, and p is the fraction of infected individuals who remain asymptomatic. and are the detection rates of I and A, respectively. , , , , and are the recovery rates of I, Q, A, R, and T, respectively. is the mortality rate, whereas and are severe illness rates of A and R, respectively. Fig. 1 shows a scheme of the state transitions.

Figure 1

Schematic diagram of the proposed compartmental model.

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 (): detected individuals cured without treatment;

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

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

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

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 other[35-37] since Brazil is a large country where there were several different outbreaks dates, local government policies, and population behavior. Brazilian spatial epidemic progression studies[35,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 variances[37].

The healed compartment from the Giordano et al. model[30] is subdivided into three compartments: , , and . is a measurable state by the Brazilian severe acute respiratory syndrome (SARS) database[38,39]. 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 for any analysis post the first wave. 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: ; thus, we substituted Equation (1k) by Equation (2).

The additional state E corresponds to a natural time delay of the system, which is usual in control-oriented models[8,11] and forecasts to a lesser extent[20,40,41]. Presymptomatic infected individuals are within this state as , 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 infections[42]. Gamma[43] and Delta[44] 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 , , and , 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 infectious[45].where 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 Gamma[43] variant.

NPI dynamics are inserted into a compartmental model by time-varying functions[20], independent variables[9,11,12], or time-varying parameters estimated over time[24,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 according to Equation (4).where is the manipulated variable related to social distancing and is the estimated contagion rate. The linearity applied over 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 .

Social distancing is measurable by Google mobility data as percentage changes concerning a baseline defined from data sets before the COVID-19 outbreak[46]. 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).where and are weekly moving averages of parks and workplace, respectively, and 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 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 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 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 is defined as the fraction of symptomatic individuals who develop severe or mild symptoms. Assuming that severe and mild illnesses imply hospitalization, thus is the fraction of individuals moving from A and R to T. Summing up Equations (1e) and (1f):which is simplified by assuming and to:

Thus:

Let us rewrite and as a probability function of the symptomatic individual to follow their ways, then:where and are the probabilities of an individual in A to recover or to get detected, respectively. The average rates of severe illness and symptomatic recovery correspond to properties studied in the literature. In this work, we defined and from CDC[47], and was based on a study of the detection window and test sensitivity of IgG/IgM tests[48]. The testing rate is a local and time-dependent property that affects both the probabilities and . Let us define the correlated parameter as the fraction of recovered undetected individuals. We have from Equation (1e):

Rewriting Equation (7) for :and substituting it in Equation (6) rewritten for :

Considering that and represent fractions of the symptomatic infected, then . Locations with a steadier testing policy could estimate 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 as a logistic equation in the function of time according to:where , , , and are identified model parameters. The definition of Equation (9) is based on heuristics that 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 , we define the fraction of recovered undetected asymptomatic individuals from Equation (1c) as:where is the probability of detecting the disease in an asymptomatic individual. Defining similarly to , we have:where is an additional identified model parameter. Equations (9) and (10) have linear dependence between and to avoid overfitting of an excessive number of model parameters. The ratio 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 analogous to vaccine efficacy against lethality as the fraction of threatened individuals who decease. We define it from Equation (1g) as:

Rewriting Equation (11) as a function of a death probability :and isolating give us:where is the average recovery rate from hospitalization and 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 and yields additional flexibility in the model formulation. Minor changes applied over , , and 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 Brazil[2] provides accumulated data on confirmed cases, deceased, and their respective incidences for each federative unit and county. The Brazilian SARS database[38,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 and . In addition, overall means of hospitalized evolution between April 2021 and July 2021 were used to define and . Both databases are daily measured; hence sampling time d. Average testing rates and are location-dependent; however, we assumed that correlated uncertainties are comprehended in , and . Hence, we defined = to suit sampling time. In summary, the monitored variable is defined as:where .

EPICOVID19-BR provides additional data over temporal distributions in Brazil. It surveyed COVID-19 prevalence in cities from all regions on different timelines[49,50]. Let us consider the prevalence estimations from federative units given by Marra and Quartin[51] based on three phases of EPICOVID19-BR. Furthermore, if we assume , then we can correlate states and with test sensitivity. EPICOVID19-BR did not test hospitalized patients[49] and used an IgM and IgG antibody test more sensitive 15 days after the appearance of symptoms[48]. Hence, the state transition model for each federative unit i is defined in Equation (14).

Each federative unit i under study has prevalence distribution from EPICOVID19-BR formulated as Equation (15) for each phase , assuming a test sensitivity of 100% during days followed by a sudden decay to zero.where prevalence bounds and can be found in Supplementary Table S2[51], and was the arbitrated value for the detection window. , , are initial dates from the first, second and third phases of EPICOVID19-BR, respectively, while and correspond to their respective duration in days, and . In the early stages of the pandemic outbreak, recovered individuals are approximately null; thus, we defined .

Gene sequences reported in GISAID[7] indicate Zeta variant appearance in mid-October 2020. Hence, the identification step is bounded at 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 . Therefore, only and were considered optimization variables for the model identification. Hospitalized individuals were used to define from the solution of a system composed by , , , 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 IPOPT[52] via CasADI/MATLAB[53].where , are the measured variables and 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 CVODES[54] via CasADI/MATLAB. The initial guess was set as . Results and location-dependent parameters are shown in Table 1.

Table 1

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

	AM	MS	RN	RS	RJ	SP
t_0,i	03 April 2020	19 April 2020	19 April 2020	19 April 2020	15 April 2020	27 March 2020
S_i(t_0,i)	0.9449	0.9997	0.9951	0.9992	0.9762	0.9956
E_i(t_0,i)	0.0551	0.0003	0.0049	0.0008	0.0238	0.0044
α_0,i	0.1890	0.3953	0.5217	0.2848	0.2628	0.3223
a_x_a,i	1.0000	0.8999	1.0000	0.9999	1.0000	1.0000
a_x_s,i	0.9730	0.8999	0.9370	0.8963	0.9491	0.8861
x_c,i	0.0270	0.1001	0.0630	0.1037	0.0509	0.1139
x_m,i	0.3796	0.2586	0.4606	0.2833	0.4530	0.2694
a_ζ,i	0.2502	0.6539	0.6649	0.6419	0.2207	0.5591
b_ζ,i	0.2500	0.0900	0.1654	0.0429	0.2499	0.0540
c_ζ,i	39.3836	68.3508	51.2077	71.2191	38.0807	71.0390
σ~_i	0.0782	0.0859	0.0721	0.0822	0.0491	0.0794
τ~_i	0.0672	0.0645	0.0766	0.0614	0.0726	0.0673
w_u	0.2604	1.000	1.0000	0.0000	0.7048	1.0000

State and parameter estimation

State estimation is essential for a model with uncertainties and without measurements from all states. It comprehends estimates of unknown properties based on available measures while filtering them to reduce the noise effects. The proposed model in Equation (14) has , , , and as unmeasurable, , , and as unmeasured, and , , and as measured states. Besides, the sum of the states , and , which represent the confirmed cases, is also a measured variable. Furthermore, the epidemiological model parameters have uncertainties related to time-varying NPI, its acceptance from the population, circulating virus variants, and vaccine coverage. Hence, a state estimator must accurately forecast the epidemiological evolution of COVID-19 on each analyzed federative unit i. In this work, we selected the parameters , , and to estimate over time. Parameter estimation was performed using an augmented state within a state estimator. The state is defined as:where are the parameters to be estimated.

Time-varying dynamics from are unknown; thus, we assumed their differential equations equal to zero, and they are subject to artificial noise. Therefore, the state transition model for the augmented state is:

Measurements in process control usually constraint real-time applicability for state estimators within seconds or minutes. Hence, the sampling time d allows analysis over different state estimation strategies. In this work, we evaluated the same scenario for each analyzed federative unit with a constrained extended Kalman filter (CEKF)[55], a constrained extended Kalman filter and smoother (CEKF & S)[56], and a moving horizon estimator (MHE)[57]. We used constrained observers to satisfy the feasible region .

CEKF is an extension of the Kalman filter for nonlinear models. It uses a first-order Taylor expansion of the system model to estimate the current value based on the latest measurement and estimated state. The COVID-19 data, however, are given in weekly cycles in which weekends have fewer notifications that are updated on working days. The CEKF & S is an intermediate option between a regular CEKF and a MHE regarding computational time and performance. First, it forwards estimates from a moving horizon with a CEKF followed by backward estimation with a smoothing equation. The weekly oscillations are attenuated in the resulting state and in the covariance update. The MHE uses a moving horizon of estimates and measured variables in a nonlinear optimization problem, which is solved at each sampling time. This optimization problem has times the degrees of freedom of the CEKF, where is the horizon size. Therefore, it provides better estimation at the cost of a significantly higher computational time.

The error covariance matrix and the initial estimated state of each federative unit i are defined at . and can be found in Supplementary Table S3, while is shown in Table 1. The matrix was defined as , the covariance matrix of process noise was defined as , and the covariance matrix of observation noise was defined as .

Let us define the model with uncertainties:where and are the process and measurement noises, respectively.

The linearization of Equation (18) into a state-space model yields: where the output matrix and the state transition matrix are defined as: whose analytical expressions for these Jacobian matrices can be found in Supplementary Equation S1. We remark that Equation (19b) is equivalent to as it is a linear function.

The initial conditions for each federative unit i were defined as , and for all state estimators. All simulations with state estimation started at and ended at . The performance of the state estimation was evaluated using the mean absolute percentage error (MAPE) calculated for each output as:

CEKF

For the sake of notation simplicity, the subscript i, denoting each federative unit, was suppressed from the description of the estimators. For the CEKF, the optimization problem in Equation (20) to update at each discrete time k corresponds to a quadratic programming, which was solved at each iteration with qpOASES[58] via CasADI/MATLAB.

The state covariance matrix was updated via the Riccati equation in discrete time as follows:

Thereafter, the discrete-time is advanced to , and Equations (20) and (21) are solved again to update and .

CEKF & S

The CEKF & S was implemented according to the formulation from Salau et al.[56] The state estimation was initially done times with the CEKF from the previous section.

The Rauch-Tung-Stribel (RTS) smooth equations[59] were applied from to the simulation end . Each discrete time started with an additional CEKF iteration to calculate and . Let us define , , , and for estimating backward with the Rauch-Tung-Striebel (RTS) smooth equations[59]. Solving Equation (22) for , yields the solution and , which is the initial conditions and for forward estimation until the current step k through iterations of the CEKF.

State and covariance estimations from each step are used to update their respective values in the vectors and (k|k). Two horizon sizes and were used in the simulations to evaluate the effect of on the estimator performance.

MHE

The MHE was implemented according to the formulation from Rawlings et al.[60] The past horizon at each discrete-time k, where is the given horizon size for the estimator. Hence, we can set the initial condition for the optimization problem in Equation (23) as and . The state is updated at each discrete time k with the solution of Equation (23) through IPOPT[52] via CasADI/MATLAB.where are the estimated states and . The state covariance matrix is updated via the Riccati Eqs. (21).

The MHE gives estimations over a horizon based on the initial conditions and . Current estimations at a discrete time k are and . The advance in discrete time is carried out by solving Equations (23) and (21) based on previous estimations from to . The MHE was implemented with .

Simulation results and discussion

In this section, we present the results from simulations for federative units Amazonas (AM), Mato Grosso do Sul (MS), Rio Grande do Norte (RN), Rio Grande do Sul (RS), Rio de Janeiro (RJ) and São Paulo (SP). States and parameters were estimated from 1 October 2020 to 1 July 2021. Confirmed cases () and deceased () measures were obtained from the Ministry of Health of Brazil[2], whereas hospitalized () and healed with treatment () were obtained from the Brazilian SARS database[38,39]. The input variable , calculated using Google mobility data, is presented in Fig. 2 for each federative unit.

Figure 2

Time evolution of input variable related to social distancing.

All three state estimators drove the estimation toward the measure for each federative unit, as shown in Table 2 with MAPE results smaller than 5%. Hence, COVID-19 dynamic evolution on regional populations was captured despite the model assumptions. The tuning of and based on values resulted in better estimates of and , and higher estimation error on for all studied cases. Using the same tuning formulation for all analyzed federative units implied some suboptimal sets of tuning parameters. Table 2 lets us identify the worst estimation from CEKF followed by CEKF & S and MHE according to expectations. Moreover, the increase in horizon size from 7 to 28 showed loss of estimation accuracy of the CEKF & S. The lack of long-term correlation for estimating state and parameter backward is probably a cause for this result; however, additional studies are required to verify the existence of other sources. The time evolution of output measures from all federative units can be found in Supplementary Figs. S1–S5, except for Amazonas, which is shown in Fig. 3.

Table 2

Mean absolute percentage error for simulation with state estimation for each federative unit i.

	State Estimator	AM	MS	RN	RS	RJ	SP
y_1	CEKF	0.40	0.59	1.05	1.49	0.56	0.61
	CEKF & S (N_p=7)	0.34	0.51	0.90	1.10	0.51	0.52
	CEKF & S (N_p=28)	0.45	0.56	0.98	1.19	0.66	0.64
	MHE (N_p=7)	0.32	0.45	0.85	0.99	0.49	0.48
y_2	CEKF	0.37	0.42	0.31	0.35	0.36	0.29
	CEKF & S (N_p=7)	0.34	0.37	0.28	0.32	0.35	0.28
	CEKF & S (N_p=28)	0.38	0.39	0.30	0.33	0.36	0.28
	MHE (N_p=7)	0.34	0.33	0.27	0.30	0.32	0.25
y_3	CEKF	3.94	3.57	2.05	3.66	1.98	2.35
	CEKF & S (N_p=7)	3.14	2.48	1.55	2.53	1.52	1.60
	CEKF & S (N_p=28)	3.58	2.55	1.45	2.57	1.62	1.55
	MHE (N_p=7)	2.52	1.66	1.26	2.07	1.06	1.25
y_4	CEKF	0.21	0.21	0.21	0.17	0.11	0.12
	CEKF & S (N_p=7)	0.20	0.18	0.16	0.15	0.11	0.11
	CEKF & S (N_p=28)	0.25	0.19	0.16	0.17	0.14	0.13
	MHE (N_p=7)	0.18	0.16	0.15	0.16	0.10	0.11

Figure 3

Time evolution of measures and estimated outputs from Amazonas (AM).

Amazonas had only two coronavirus waves identifiable through , as can be seen in Fig. 3, unlike other analyzed federative units. GISAID data[7] indicated that variants B.1.1.33, B.1.1.28, and a local B.1.378 were significantly circulating from 1 October 2020 to 4 December 2020 when the first Gamma variant sequence was identified. Hence, the predominant circulation of the Zeta variant starting in mid-October 2020, was quickly overlapped by the Gamma variant resulting in a single wave. The estimated parameters, presented in Figs. 4-6, represent this profile specifically with MHE estimations and indicate that CEKF & S with had closer parameter estimations with MHE than CEKF & S with . It is important to emphasize that MHE was applied with , which could explain this similarity.

Figure 4

Time evolution of estimated contagion rate from all analyzed federative units.

Figure 6

Time evolution of estimated fraction of threatened individuals who decease from all analyzed federative units.

All other federative units also had an overlap of the Zeta and Gamma waves; however, there were higher periods with the circulation of the Zeta variant. Figs. 4–6 infer rougher parameter estimation from CEKF & S with , despite having a larger horizon. The moving horizon with estimated states backward and forward might smooth estimations overall, but each estimate still considers a single measure. Further studies on the effects of in CEKF & S are required, but these were not the subject of this work. Since MHE provided better state and parameter estimations, descriptions henceforth are related to these estimates.

Age distribution and local healthcare do not explain the discrepancy between and among the analyzed federative units observed in Figs. 5 and 6. This discrepancy points out the violation of the model assumption considering the complete identification of hospitalized individuals infected by COVID-19. Some case studies presented higher altogether with lower , which do not affect the infection fatality rate (IFR) but indicate underreporting of COVID-19 cases among hospitalized individuals. Amazonas, Rio Grande do Norte, and Rio de Janeiro presented testing policies focused on more severe hospitalizations. The IFR calculus defined in Equation (24) highlights this dynamic.

Figure 5

Time evolution of estimated fraction of symptomatic individuals who develop mild and severe symptoms from all analyzed federative units.

The results in Fig. 7 indicate a guaranteed underreporting of deaths in Amazonas. In addition to Figs. 5 and 6, the decrease shows that underreporting of hospitalized individuals in Rio de Janeiro and Rio Grande do Norte reduced over time. The lethality evaluation of a variant by is unfeasible because it decreased in Rio Grande do Norte and Rio de Janeiro due to testing policy. In addition, its peaks in Mato Grosso do Sul, Rio Grande do Sul, and São Paulo are explained by delayed notifications after an overload of the health system, which temporarily increases mortality. Hence, we used to conclude that the Zeta variant increased lethality from in Rio de Janeiro to in Rio Grande do Norte based on calculated with values from Table 1. In addition, estimation indicates that the Zeta transmissibility increased from in Rio de Janeiro to in Rio Grande do Norte.

Figure 7

Time evolution of IFR considering estimated parameters.

Amazonas faced oxygen shortage on the Gamma variant wave, which implicated in mortality increase beyond virus mutations. Nonetheless, estimates of had an increase of over its initial value, while had an increase of . In addition, its estimations on are smoother, which allowed identifying the increase in the severity of the disease ranging from in Rio Grande do Sul to in São Paulo. The analysis through indicated a lethality increase between in Rio Grande do Norte and in Amazonas for the Gamma variant. Moreover, transmissibility increased between in Rio de Janeiro and in Rio Grande do Sul based on first-wave values of from Table 1. Further investigation on points out variant spread countrywide, being Amazonas its source. The Gamma variant spread initially to the farthest case study from Amazonas: Rio Grande do Sul. Afterward, it spread to the Brazilian economic center, São Paulo, and thereafter to the rest of the federative units at a similar rate. The quicker propagation of the Gamma variant to the farthest location from Amazonas, Rio Grande do Sul, is explained by a less rigid NPI highlighted by the highest estimate. São Paulo had only the fourth highest among the studied cases. However, it was the second to significantly contract the Gamma variant, which enforces the theory that it spread countrywide afterward and corroborates its classification as a super-spreader city par excellence by Nicolelis et al.[35] Overall, more clearly indicated the emergence of circulating variants in the system. Finally, decreased in later times for most federative units until 1 July 2021, indicating that vaccination coverage does reduce mortality in infected individuals.

The circulating variant dynamics assumed the unique circulation of the lineage, whose uncertainty was reduced by a manual definition of the analysis period for each variant. Most studied cases had the Zeta wave overlapped by the Gamma variant; thus, dynamic estimations are expected to be lower or equal to the actual value. The Zeta evaluation period was defined in the last 15 days before Gamma variant emergence. The Gamma variant estimations are expected to be more accurate since it was predominant over some time for all cases studies. The Gamma evaluation period was defined from the first stationary point after variant emergence to the end of the simulation.

The computational time of the three state estimators was evaluated throughout the average simulation time among all analyzed federative units for 273 d, from 1 October 2020 to 1 July 2021. Simulation time was measured by the tic and toc functions in MATLAB. All simulations were carried out on an AMD Ryzen 5 5600X 3.70 GHz in a sequence to mitigate computational noises. The average simulation time was 25.4 s for CEKF, 136.7 s for CEKF & S with , 498 s for CEKF & S with , and 10265 s for MHE. Even the average execution time of 38 s per sampling time for the MHE implies real-time applicability of the state estimators with the selected tuning in the COVID-19 pandemic scenario since all of them have execution times lower than the sampling time d. CEKF & S performance and computational time were between the CEKF and the MHE, which enforce it as an alternative for processes with faster sampling times.

The definition of a compartmental model inherits limitations regarding the closed system and homogeneous compartment assumptions. In addition, all numerical results are dependent on the initial condition, which was determined from a nonlinear optimization in this work. Age distribution was neglected in the model formulation to aim for real-time applicability and fulfill available data of confirmed cases. Model assumptions uncertainties are mitigated by the state and parameter estimation; however, they do not guarantee realistic estimations. For instance, the mitigation of the variants reinfection mostly through compartments and instead of estimated parameters , , is a consequence of the tuning. Hence, a fine-tuning procedure may be required for severe assumption violations to avoid unrealistic estimations. Mitigation of multiple uncertainties in the model formulation is achieved by a conservative tuning concerning small dynamic changes. Hence, smooth vaccination coverage and gene sequence dynamics might be noise to the estimator. Data quality also limits a more aggressive tuning for state estimators, such as sudden data updates of underreporting for cases and deceased (e.g., Rio Grande do Norte on 23 July 2021). Nonetheless, the proposed method allowed the study of overall dynamics in each studied case.

Conclusion

In this work, we proposed a mathematical model able to identify underreported cases of COVID-19 from hospitalized and deceased individuals by comparing the fraction of symptomatic individuals who develop severe or mild symptoms, the fraction of threatened individuals who decease, and the infection fatality rate among analyzed federative units. In addition, the model identified circulating variant dynamics in the aforementioned parameters, and characterize them under some assumptions. We remark that this model is suitable for control strategies, assuming there are available hospitalized data.

The performance among estimators confirmed MHE as a more suitable state estimator for COVID-19 due to daily sampling time. Nonetheless, CEKF & S presented reasonable estimations for comparison, and a significant reduction in computational time, which make it applicable in real-time applications.

Parameter estimations identified a lethality increase ranging from 11 to and a transmissibility increase between 10 and for the Zeta mutation. In addition, we found that the Gamma mutations caused a lethality increase ranging from 44 to and a transmissibility increase between 43 and . The estimation strategy successfully detected and estimated dynamics affected by the emergence of COVID-19 variants, which improves model accuracy for further predictions. Moreover, an initial decrease in lethality due to vaccination was also observed. Hence, the parameter estimation within recursive state estimation can deal with dynamic uncertainties from the COVID-19 pandemic.

Future works account for implementing an economic model predictive control and studies on inserting vaccination into the proposed model. Delta variant has been predominant in Brazil since August 2021. It was disregarded from an initial analysis because its mutation highly increases contagion among vaccinated people, which are measured. Therefore, a model comprising vaccinated individuals should generate better estimations of Delta dynamics.

Supplementary Information.