State estimation of the time-space propagation of COVID-19 using a distributed parameter observer based on a SEIR-type model.

The real-time prediction and estimation of the spread of diseases, such as COVID-19 is of paramount importance as evidenced by the recent pandemic. This work is concerned with the distributed parameter estimation of the time-space propagation of such diseases using a diffusion-reaction epidemiological model of the susceptible-exposed-infected-recovered (SEIR) type. State estimation is based on continuous measurements of the number of infections and deaths per unit of time and of the host spatial domain. The observer design method is based on positive definite matrices to parameterize a class of Lyapunov functionals, in order to stabilize the estimation error dynamics. Thus, the stability conditions can be expressed as a set of matrix inequality constraints which can be solved numerically using sum of squares (SOS) and standard semi-definite programming (SDP) tools. The observer performance is analyzed based on a simplified case study corresponding to the situation in France in March 2020 and shows promising results.

Introduction

Lately, mathematical models in epidemiology have attracted considerable attention, and extensive research has been carried out in this field due to their important role in human life [1], [2], [3], [4]. Social networks provide practical methodologies to simulate epidemics spreading by developing mathematical models to predict and monitor critical features of outbreaks.

In late 2019, a disease outbreak emerged in Wuhan, China. The culprit was a certain strain called Coronavirus Disease 2019 or COVID-19 in brief. This virus has been identified to cause fever, cough, shortness of breath, muscle ache, confusion, headache, sore throat, rhinorrhoea, chest pain, diarrhea, nausea and vomiting [5]. Although only seven coronaviruses are known to cause disease in humans, three of these, COVID-19 included, can cause a much severe infection, and sometimes fatal to humans.

The rapid spread of COVID-19 all over the world has become a matter of grave concern and has hugely altered the lifestyle and social behavior of the populations from the beginning of . Indeed, it poses considerable economic, environmental, and political challenges all over the world. Consequently, considerable research effort has been made to investigate precise mathematical models for the outbreak of this newborn virus and rapid estimation of its future transmission and mortality rates. For the time being and in spite of the massive vaccination, COVID-19 infection is still active in many countries at the time of writing this article. Governments intensify the vaccination campaigns to combat the disease whereas research institutions try to find out the effectiveness of the vaccines regarding the new variants of the virus circulating in many countries.

Several mathematical models have been proposed from various epidemiological groups [4], [6], [7]. These models help governments as an early warning device about the size of the outbreak, to assess how quickly it will spread, and how effective control measures may be. However, due to the limited emerging understanding of the new variants and their transmission mechanisms, results may only reflect specific scenarios.

Susceptible–Infectious–Recovered (SIR) models, with compartments representing susceptible, infectious, and recovered portions of a population are often considered as variations of the well-known Kermack–McKendrick model [1]. While SIR-type models are commonly used in governing the forward dynamics of epidemiological systems, a variety of state estimation techniques have been used in the literature for control or monitoring purposes [8], [9], [10]. Particularly for the COVID-19 pandemic, one can mention, among others, the recent agent-based modeling studies [11], [12] and those reported in [13], [14], [15] that perform state estimation of time propagation SIR models. On the other hand, time–space propagation applications have not been explored much despite their importance for monitoring the evolution of compartment distribution in a determined host region.

In this work, employing a new generalized epidemiological (susceptible–exposed–infected–recovered) model of COVID-19 spread [7], a model-based state estimator design is developed for the epidemic spread dynamics in a prescribed host population. The presence of time-variant parameters and incomplete measurements are tackled through a robust state estimator based on the Lyapunov theory. In this context, the contribution of this work is a computationally tractable method for the nonlinear Luenberger-like observer design of the SEIR model of COVID-19 spread [7]. This method provides a guaranteed estimation error decay rate and is formulated in terms of LMI constraints and a polynomial parametrization in such a way that the corresponding solution provides a convergent estimator. The observer allows monitoring the susceptible, exposed and infected distribution in a rectangular spatial region with Neumann boundary conditions. Stability analysis in 2D domains sets a challenging problem with regard to the derivation of possibly not too conservative conditions to be solved by semidefinite programming (SDP). The proposed methodology uses integration by parts as a fundamental tool as well as the Poincaré–Wirtinger inequality and the S-procedure for local stability analysis based on sector conditions.

The reminder of this paper is organized as follows. In Section 2, the epidemiological model is presented along with some preliminary assumptions. Section 3 introduces the mathematical framework for the proposed observer synthesis based on the abstract formulation of the error dynamics and the set-up of the sector condition related to the nonlinearity embedded into its dynamics. The Lyapunov convergence analysis and the main theoretical contribution of this paper are given in Section 4. Section 5 is devoted to numerical tests of the proposed observer design scheme in a case study dedicated to monitoring of the spread of COVID-19, based on the model presented in Section 2. Section 6 draws concluding remarks and points out some possible research lines.

Notation. We denote the set of natural, real, positive and non-negative real numbers by , ,  and , respectively. The vector space of  dimensional real vectors and -by- real matrices are respectively denoted by  and .  represents the subspace of symmetric matrices, where the multiplicative and additive identities are denoted by  and , respectively. The superscript “” denotes matrix transposition, “” the Kronecker product,  the column vector by vertically stacking the columns of the matrix ,  a block-diagonal matrix,  the Hermitian operator applied to matrix ,  a connected open subset of ,  the boundary of , ,  and  the time, first and second order spatial derivatives of the function , with respect to , respectively.  denotes the space of square Lebesgue integrable -dimensional vector valued functions, that is, functions with finite norm  also  is the Sobolev space, often denoted by  and defined as  where , is a bidimensional multi-index of order  and where the following notation is used:  Explicitly in this work, the Laplacian of  is the sum of all the unmixed second partial derivatives in the Cartesian coordinates :  denotes the vector of monomial bases of degree  or less, i.e,

Instrumental Tools. The following statements are instrumental for deriving the main theoretical results of this paper.

Let . We define by  the convex hull of the set , i.e.

Differential Mean Value Theorem [16]

Let , and be a differentiable function with respect to . Then, there is an element , such that: where .

Poincaré–Wirtinger Inequality [17]

Let be an open bounded Lipschitz connected subset in . Then there exists a constant , depending only on , such that for every function where is the mean value of over , with standing for the Lebesgue measure of the domain .

The Poincaré constant  depends on the geometry of the domain . Particularly, if  is a bounded, convex, Lipschitz domain with diameter , then the Poincaré constant is at most  [18].

Spatiotemporal model of COVID-19 infection spread

We consider the SIR-like model presented in [19] and further developed in [7] where it is assumed that a host population of individuals is divided into compartments corresponding to disease status, modeling the movement in space and time of the subpopulation in each compartment. Specifically, these compartments are the densities of Susceptible population , Exposed population , Symptomatic Infected population , Asymptomatic Infected population , Under treatment population , Removed population  and the deceased population . Note that  refers only to deaths due to COVID-19. We denote the living host population as . Due to the names of the compartments used, this model is often referred as the susceptible–exposed–infected–recovered (or, in short, SEIR) model. We therefore formulate the problem in terms of the state vector  containing the different compartments. This model assumes that

the spatial mobility is governed by diffusion coefficients according to the mobility restrictions of the host population;

only susceptible, exposed and asymptomatic individuals are moving;

there is a latency period between exposure and the development of symptoms.

Description of the infection flow

The susceptible class contains individuals who do not have temporary immunity to the virus, then might become infected if exposed. The exposed class contains individuals who have been infected but do not have symptoms. The period that starts when the person becomes infected, until the person becomes symptomatic or asymptomatic is the latent period  days. The under treatment class contains individuals who are currently infected and cannot transmit the infection because of adequate isolation. The recovered class contains individuals who returned to a normal state of health after having been infected during the latent period  days. The number of deaths depends only on the death rate as the number of recovered depends only on the recovery rate. Finally, the cumulative number of infected depends only on the exposed and the incubation period. The diffusion parameters are included in the model to spread the disease spatially. Fig. 1 depicts the infection flow according to the explanation above.

Fig. 1

Compartmental representation of the -model.

The dynamics is governed by a system of three partial differential equations (PDE) and three ordinary differential equations (ODE) as follows  for  as spatial and time domains respectively. The total living population is  and the deaths are . No new recruit is added and, if we assume that the region of interest is isolated, we prescribe the following homogeneous Neumann boundary conditions,  where  is the outward normal vector to . This selection of boundary conditions represents the situation where the spatial region under consideration is closed to any in- or out- flow of populations, so that the epidemic spread is only due to local infections, which was the case during lock-down conditions in the first and second wave of the pandemic. If the population is traveling in or out of the considered spatial region, then the boundary conditions need to be formulated in another way using more general boundary conditions or mixed boundary conditions. However the computational procedure proposed in the following is not suited to these conditions, and the approach should be revisited. Specifically, when other kinds of boundary conditions are considered, the additional terms in the integration by parts procedure within the Lyapunov analysis may yield terms that cannot be expressed in an affine way and hence could not be approached using SDP tools directly.

The system parameter definitions are listed in Table 1.

Table 1

Parameter definitions.

Parameter	Definition
wβ_e	Transmission rate from S to E from contact with E (days−1)
wβ_s	Transmission rate from S to E from contact with I_s (days−1)
wβ_a	Transmission rate from S to E from contact with I_a (days−1)
δ_l	Latency rate (days−1)
p_s	Probability of being symptomatic (days−1)
1−p_s	Probability of being asymptomatic (days−1)
γ_r	Recovery rate (days−1)
μ_d	Death rate (days−1)
ν_t	Under treatment rate (days−1)

Latency period and infection period have been estimated as  days and  days respectively [20], and thus  days−1,  days−1. To account for the lockdown and unlockdown, the average number of contacts is updated as follows [21] while the diffusion coefficient is set up to

Here bol stands for beginning of lockdown and eol for end of lockdown. Unlockdown is assumed to be faster than lockdown. The parameter  is a varying coefficient translating respect for distancing. The value of  is fixed according to the average daily commute related to the host population. Then, the numerical values of six parameters  remain to be determined from the measurements  and  for  days considering a least square cost function in such a way that the optimization problem may be solved using Approximate Bayesian Computation combined with a quasi-Newton method [7].

This work was initially conceived in October 2021 when massive vaccination campaigns were not still available and the main control strategy was lockdown by means of the parameters  and . Currently with vaccination as the main retaining strategy of the disease spread, the PDE model in (5) may be modified to attend this counteraction of the pandemic as addressed in [13], thus, considering  as the vaccination rate  which does not impose any constraint in the applicability of the methodology presented in this paper as it may be verified in the mathematical developments.

Basic reproduction numbers

The basic viral reproduction number  of the infection is the expected number of cases directly generated by one case in a population where all individuals are susceptible to infection. In [7], a condition on parameters which defines the basic reproduction number related to the model described by (5)–(6) such that the disease has an exponential initial growth is given and to reflect the spatio-temporal dynamics of the disease, the effective reproduction number is considered:  and its mean with respect to the domain

The density population of any host population, at any time , would remain constant based on the epidemiological dynamics (8) :

Also, it is worth mentioning that all compartmental variables have positive values (populations). Thus, it can be concluded that all of these variables remain bounded during the outbreak management time due to the boundedness of the total density population .

State observer design

In this section, a distributed parameter observer is developed to estimate the compartmental states of the COVID-19 spread dynamics described by (5)–(6). The reason why state estimation is required is two-fold: (a) it is impossible to accurately measure the numbers of susceptible, exposed, and infected asymptomatic people; and (b) other state measurements are typically global values throughout the analyzed region, so it might not be convenient to use them in real-time control design as proposed in [13]. In order to derive the LMI conditions for the observer design arising from the Lyapunov stability analysis and make the solution computationally tractable by means of SDP tools, we address the state estimation problem in a rectangular spatial domain defined by . This choice along with the Neumann-type boundary conditions related to the model, set a feasible analysis for the epidemic propagation in a general host population whenever the rectangle circumscribes the entire population distribution.

Preliminary results

Assuming that the measured output is given by the total symptomatic and death individuals (which corresponds to the data usually provided by health authorities), we consider

Then, (5)–(6) take the form of  by considering  with output measurement  where  since .

Nonlinear Luenberger-type state observer

Defining the Luenberger-type state observer for  for , subject to  and the initial condition  for . Here the gain  is the output injection gain to be designed. The dynamics of the state estimation error , satisfies

subject to  and the initial condition

The function denoted as  and the estimation error  satisfy a sector condition based on the boundedness of the Jacobian matrix of the nonlinear function . Then, the Differential Mean Value Theorem gives  where  for all . Let  be the constant matrices whose entries are the local lower and upper bounds, respectively, of the Jacobian matrix entries of  and hence the following inequality holds  which implies the following

From the definition of the rate function in (15), its jacobian is given by

From Property 1, , ,  and  , hence one choice for  may be defined as  with .

Abstract formulation

The error dynamics described by (21)–(22) can be rewritten as an abstract first order ordinary differential equation in the Hilbert space  according to  where the operators ,  are defined as

The error dynamics in (29) can be represented as a Lure system, depicted in Fig. 2, where the sector condition for the estimation error  and the deviation function  are expressed through the scalar constraint (26) that can be embedded into the local stability analysis by applying the S-Procedure [22].

Fig. 2

Lure-System representation of the error dynamics.

Lyapunov convergence analysis

The state observer design problem is addressed within a weighted Lyapunov framework, with the weight function as a degree of freedom. The analysis of the corresponding dissipation mechanism leads to an LMI convergence condition, in which depending on the spatial coordinates, the observer gain and the Lyapunov weight function are the decision variables. To this end, let us set the positive-definite weighted candidate Lyapunov functional  as  where  is a strictly positive operator defined by the polynomial matrix  as  for all . The following Lemma shows how two positive semi-definite matrices  and some constant  can be used to define the polynomial matrix  such that the operator  is positive and therefore the functional  is a Lyapunov candidate for the observation error dynamics (29).

Given any positive semi-definite matrices , and where , is a vector of monomials with degree or less and is Kronecker product. Let for all if for some then the functional , defined as is a strictly positive functional over , whenever , and satisfies

The error dynamics in (29) is (locally) exponentially stable with decay rate if there exist

, and real positive scalars ,

block diagonal positive semidefinite matrices , , with such that the polynomial matrix satisfy (36) ,

the th degree polynomial matrix

and the following matrix inequality is feasible: and where with

Consider the linear dissipation expression of the Lyapunov function  then, the substitution of (21) into (43) yields

Applying the first mean value theorem for integration, there exists a scalar  such that  and regarding the constraints on the definitions of  and , . Hence, in order to apply the integration by parts in (44), we take (59) into account as presented in Appendix. Thus, (44) becomes

Notice by the virtue of Poincaré inequality that the following holds  with . Hence, substituting (47) into (46) leads to

We can rewrite (48) as  where . Therefore, in order to ensure the negativity of the right side of (49), it suffices that

Applying the -procedure to (26), (50), we obtain

Then (51) implies that  and from the comparison lemma [23], it follows that  where

Covid spread monitoring

In this section, the distributed parameter observer is evaluated under some realistic scenarios. Particularly, the state estimation of the compartmental variables in a host population defined by a 2D spatial domain corresponding to the French situation in March 2020 is discussed. Based on data available from the National Institute of Statistics and Economic Studies (INSEE), the average commuting implies  and the values corresponding to the model parameters are , , , , , , ,  [7]. The population distribution is depicted in Fig. 3, Fig. 4 shows the initial simulation scenario in which the population distribution corresponds to the infected symptomatic confirmed on 18 March 2020 by the health authorities.

Fig. 3

Population density of France  (peoplekm).

Fig. 4

Initial infection density  (peoplekm) on 18 March 2020.

Firstly, we show the feasibility provided through the solution of the LMIs in (39) for different selections of lower bound decay rates . Then, the convergence features of the proposed observer are assessed through the numerical simulation of the observer system.

As formulated in Section 3, the available online measurements are the global infected  and dead population . Thus, the system detectability is sufficiently set by means of the LMI feasibility in (39) according to the conditions presented in Theorem 1. To obtain the output injection gain  through the application of the sufficient conditions presented in Theorem 1, we make use of SOSTOOLS for Matlab with SeDuMi [24]. The following algorithm for the computation of the output injection gain  is proposed for prescribed values of .

Algorithm 1.

Define the values of  and .

Declare matrices  and scalar .

Construct polynomial matrices  and .

Construct matrices  and  corresponding to the respective substitution of  by  and  into (40).

Solve the affine constraint (39) for  and  simultaneously.

Compute .

Thus, the solution for different values of  and , along with bisection search, provides  for each case. Some results are presented in Table 2.

Table 2

of the proposed approach for different combinations of  and .

(m,q)	γ_max
(5,3)	0.14
(6,4)	0.26
(7,5)	0.37

For , the observer and system responses are generated via numerical simulation starting from the system initial conditions defined in Fig. 3, Fig. 4 which represent the infection situation on 18 March 2020 whereas the observer initial conditions are set to zero. The simulation scenario considers a first time interval  days,  days which corresponds to the end of lockdown on 11 May 2020 and schedules the time-varying definitions of  and  in (7), (8) respectively.

Fig. 5 shows a few snapshots of the evolution of the actual distribution (on the left) as compared to the estimated distribution (on the right) of the  asymptomatic compartmental variable considering the proposed observer. As depicted in the corresponding figures, the uniform wave front propagates for . This is logical, because  and  increases during this time period. It may also be seen that the state estimation has converged in space with respect to the actual variables by these corresponding time instants.

Fig. 5

Spatial distribution of  and  (peoplekm) at time instants  days.

Fig. 6 shows the time evolution of the total values of the compartmental variables, while Fig. 7 shows the evolution of the estimation error norm.

Fig. 6

Time evolution of the total number of individuals in the several compartments of the generalized epidemiological SEIR model of COVID-19.

Fig. 7

Time evolution of the estimation error norm .

Although the initial estimation profiles were all set to be zero, the estimation error norm converges quickly (approximately after  days), and hence, provides very satisfactory estimates. In the same figure, the result corresponding to the inclusion of measurement white noise (variance ) in the  variable is depicted. It should be observed that the procedure outlined here guarantees that the state estimate robustly converges to the true state.

One of the main virtues of the SEIR distributed model based state estimation is to attain the monitoring of the pandemic throughout the host region which plays a fundamental role as warning system of the current epidemiological situation in each sub-population in contrast with state estimation techniques based on lumped SIR models such as [15]. Thus, the distributed estimation of the effective reproduction number may be helpful and a good strategy for authorities to take measures regarding the geographical propagation of the disease. Hence, applying the state estimation methodology proposed in this work, Fig. 8 displays such variable in the host region of study on 11 May 2020 (end of the lockdown).

Fig. 8

Estimation of the reproduction number  at the end of lockdown on 11 May 2020.

Conclusion

In this work,a nonlinear Luenberger-like state observer was designed to estimate the compartmental variables of the generalized epidemiological SEIR model of COVID-19 spread proposed in [7]. The proposed approach aims at obtaining a mechanism to monitor the COVID-19 epidemic in a prescribed host region with different population distributions by considering the realistic scenario of time-variant parameters of diffusion and transmission and only measurement data provided by health authorities, namely positive tests and deaths. The Lyapunov method has been applied to derive a set of LMI based conditions with SOS parametrization to ensure the local stability of the error dynamics which has been modeled as a Lure type system with a multivariable sector condition, thus the appropriate choice of the observer gains results from a SDP problem solution. In this framework, the presence of bounded time variant parameters and incomplete measurement feedback has been tackled through the polytope analysis of the proposed LMIs which guarantees a robust convergence. Numerical experiments are presented to illustrate the method efficiency. Future work entails the extension of the approach to time and space varying diffusion coefficients relative to the distinct mobility of the compartmental population, as well as the effect of vaccination.

Funding

This research received no external funding.

CRediT authorship contribution statement

Ivan F.Y. Tello: Conceptualization, Methodology, Formal analysis, Writing – original draft, Writing – review & editing. Alain Vande Wouwer: Conceptualization, Methodology, Formal analysis, Writing – review & editing, Supervision. Daniel Coutinho: Conceptualization, Methodology, Formal analysis, Writing – review & editing, Supervision.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.