From the hospital scale to nationwide: observability and identification of models for the COVID-19 epidemic waves.

Two mathematical models of the COVID-19 dynamics are considered as the health system in some country consists in a network of regional hospital centers. The first macroscopic model for the virus dynamics at the level of the general population of the country is derived from a standard SIR model. The second local model refers to a single node of the health system network, i.e. it models the flows of patients with a smaller granularity at the level of a regional hospital care center for COVID-19 infected patients. Daily (low cost) data are easily collected at this level, and are worked out for a fast evaluation of the local health status thanks to control systems methods. Precisely, the identifiability of the parameters of the hospital model is proven and thanks to the availability of clinical data, essential characteristics of the local health status are identified. Those parameters are meaningful not only to alert on some increase of the infection, but also to assess the efficiency of the therapy and health policy.

Introduction

COVID-19 is an aerial virus which strickens humans through a respiratory infection, see Zhu et al., 2020. It belongs to the Coronavirus family which was discovered in the 60′s and has already infected humans through SARS and MERS. SARS is an atypical pneumonia which appeared for the first time in 2002 in China as described in Ksiazek et al. (2003). MERS which appeared in 2012 in China too (see Zaki, van Boheemen, Bestebroer, Osterhaus, & Fouchier, 2012), is very similar to SARS but with a higher mortality. COVID-19 was declared to WHO (the World Health Organization) at the end of 2019 from cases in Wuhan, China, see Zhu et al., 2020. On the 20th of February 2020 the declared positive cases were 76,000 with almost 2500 deaths, mainly in China. At the same time, there were about 50 declared positive cases in Northern Italy, while in France there were less than 30 declared positive cases. The health status was far from being homogeneous in those or other countries. The WHO declared on March 11th, 2020, that COVID-19 was a pandemic on the website (World Health Organization, 2019). On June 1st, 2020, the declared positive cases around the world were 6.15 million with 375,000 deaths, with almost 233,000 positive cases and 33,000 deaths in Italy and 152,091 positive cases and 28,833 deaths in France. On September 21st 2020, the confirmed cases globally amount to 30,909,405 with 958,754 deaths communicated to WHO by national authorities. The total number of cases in France on September 21st is 420,855 with 31,109 deaths, while in Italy the total number of cases is 296,569 with 35,692 deaths (website World Health Organization, 2020).

Mathematical models are fundamental to understand and to predict the mechanisms of the spread of an epidemic. The most popular and widely used is the SIR model in Kermack & McKendrick, 1927 for human-to-human transmission. For the Covid-19 pandemic, many models have been built to explore the epidemic at the scale of a country like the SIDARTHE model which is an upgrade of the SIR model and is found in Giordano et al., 2020. Recall that SIDARTHE stands for a Susceptible, Infected, Diagnosed, Ailing, Recognized, Threatened, Healed and Extinct model. The SIDARTHE model has highlighted that restrictive social-distancing measures need to be combined with widespread testing and contact tracing in Giordano et al., 2020. Gevertz, Greene, Hixahuary Sanchez Tapia, and Sontag (2020) incorporate explicit social distancing via separate compartments for susceptible and asymptomatic individuals in the SIR model, while in Di Giamberardino, Iacoviello, Papa, & Sinisgalli, 2020 infected people are split into infected with low viral load, undiagnosed ones, diagnosed ones and in quarantine ones. Alternatively the SARS-CoV-2 dynamics is also described at a within-host level in Hernandez-Vargas and Velasco-Hernandez (2020).

In this paper, a new continuous-time macroscopic model, valid at the scale of a country, is introduced and it is argued that suitable delays are mandatory to reflect the dynamics of the infection while maintaining a certain simplicity in the modelling. An observability analysis is processed for further insight of this model.

Since the health status is far from being homogeneous over a full country and since the health system consists in a network of major regional hospitals, it is worth to take advantage of the availability and agility of local hospital data rather than just merging them in some centralized information system.

At this point, the mathematical modelling of the local population flows at the level of one single hospital appears to be definitely relevant. Such models may be interconnected to include the transfer of COVID-19 patients from one hospital to another to distribute the pressure at some peaks of the infection, as done in Italy with transfers from the north to the south, or in France with transfers from the east to the west and even to Germany.

The identification of the parameters of the second local model is agile, fast and relevant not only to alert on some increase of the infection, but also to evaluate the flows of severely infected patients, and thus to assess the efficiency of the therapy and the local health policy.

Data from Nantes University Hospital were collected on a daily basis for about 6 months from March 16th to September 17th 2020, and are used here to identify essential characteristics of the local health conditions. These data include the daily values of patients in conventional care, intensive care, patients who died from COVID-19 and patients which were discharged from the hospital because completely recovered or partially recovered but in the condition of continuing their treatment at home. The originality of the data from Nantes university hospital lies in the fact that the lockdown started early, the emergency services were not saturated and they focus on a population of patients with mild or severe Covid-19 cases. The macroscopic and local models are illustrated in Fig. 1 .

Fig. 1

Epidemiologic models for two different scales. (In dark grey : the map of France, in red : Nantes University Hospital).

Control systems theoretic tools have shown their efficiency to give a new insight for various biomedical systems including for instance HIV infection in Chang, Moog, and Astolfi (2014). It is argued that new solutions to cope with the COVID-19 pandemic may also benefit of those engineering science tools.

The outline of the paper is as follows. A new infection model involving delays is derived from the standard SIR model in Section 2. To fit to the mainstream of the current literature, this model is described in continuous time. The observability analysis of this continuous time-delay model is processed in Section 3. A subsystem, consisting in a node of the global health system network, is extracted and further detailed in Section 4 as it models the dynamics inside a hospital center. Since such a model is obtained from the daily data available from a University hospital, this second local model is directly designed in discrete-time. In Section 5 a discussion on the identification of the parameters characterizing the local model and their interpretation is carried out. Conclusions are pointed out in Section 6.

The dynamics in the general population

The evolution of the illness is described through the following dynamical model which represents a modification of the well-known SIR model introduced in 1927 in Kermack & McKendrick, 1927 to describe the diffusion of an epidemic disease by considering the evolution of three classes of people (compartments) the Susceptible individuals, those who can become infected, the Infected, those who spread the disease around, and the Recovered, those who recovered from the disease. With respect to the SIR model, hereafter we split the Infected people into two compartments due to the specificity of the disease, and we consider also the compartment of dying people. In particular, apart the high transmission rate, other two aspects were immediately pointed out by the physicians which did strongly influence the diffusion of the disease and the medical resources: first it was estimated that a large delay of time (10 to 14 days) is present between the moment in which a person becomes infected and can infect, and the instant in which symptoms become evident and the person is isolated and sent to quarantine. Secondly it took a long time for many patients to recover from the disease. Some of them were in hospital and in intensive care for a big amount of time (many more than one month) thus maintaining the resources unavailable to help other patients. To highlight these two aspects we then refer in the following to a time delay system. As it will be shown in the next section the delay introduced will have an important role on the observability properties of the dynamics and thus cannot be neglected. Such a dynamics can then be described by the following delay-differential equations:

I are the infected patients aware of their disease and who are thus in quarantine. They include all hospitalized patients, but not only.

R is the sub-population which has recovered from the infection.

D denotes the cumulative number of patients who deceased from the infection and were already identified so they are part of the I population.

S is the amount of naive individuals among a given population which are susceptible to become infected. The S population is infected by the I infected individuals which are not in quarantine.

I are asymptomatic infected people which represent the main source of infection and who spread out the infection among the general population.

βSI is the amount of newly infected individuals per time unit. This term splits into two parts, a smaller part of very sensitive people which affects the dynamics of I and the majority of newly infected individuals will increase the number I.

α denotes the death rate due to the infection and mainly affects the I population.

is the amount of patients in quarantine who recover from the disease.

is the amount of other infected individuals among the general population who recover from the disease.

is the amount of unaware infected individuals who become aware of their infection and go into quarantine. This phenomenon occurs with some delay τ which is typically evaluated from 2 to 3 weeks.

As already highlighted, D represents the COVID deaths from the infected and in quarantine compartment. There may be of course unknown COVID deaths which can be described by an additional dynamics of the formso that the total number of deaths due to COVID is actually given by .

D cannot be measured, and had an important role essentially at the beginning of the pandemic when the number of detected COVID cases was much lower than the real one. Due to the high mortality caused by COVID in that period, a rough estimation could be obtained by comparing the month death rate of a single region with the corresponding one in the current year.

As already underlined, the dynamics (1) is affected by two different kind of delays: the first one characterizes the amount of time that passes between the moment a person is infected and can infect, and the moment the person becomes aware of the illness and is put to quarantine. This delay is intrinsic of the disease and can only be marginally affected. The second delay which characterizes the dynamics (1) is instead linked to the time infected people recover from the disease. This delay instead has drastically changed from the starting of the pandemic and is mainly due to a better knowledge of the disease and the therapies needed by the patients, which allows now to recover faster from the disease.

In Figure 2 a representation of the propagation of the disease is given. Susceptible people (S) remain susceptible in time or become infected. Once infected, they will move to the compartment of people in quarantine (I) or infected asymptomatic people (I), depending whether they have symptoms or not. From the compartment I people will become symptomatic and move to I or recover and move to R. Those in quarantine will recover and move to R or die due to the disease (D). Finally, as it has been recently put in evidence, people who recover may be infected again, so they should be counted in the susceptible people. The model described by Eq. (1) has from this point of view a short term validity. It is possible to take into account the possibility of reinfection of recovered people. In this case, the dynamics (1) would be modified by settingThe coefficient ν seems however to be very small at this stage, and is neglected in this study.

Fig. 2

Flow Diagram associated to the evolution of the disease in a population described by Eq. (1).

In the next Section, the dynamics (1) is further analyzed with respect to its observability properties since this kind of study allows to have an estimation of the state variables. The subsystem (2) consisting by I and D is then further discussed in Section 4: a group of people who are aware of their infection define the flow of admissions in a local hospital and are split into two populations, the patients admitted in conventional hospitalization and the patients admitted in intensive care.

Analysis of observability

The health status of a given population has to be assessed from its number I of attested COVID-19 cases. It is important to guess the real number of infected individuals, and this is the purpose of observability from the measurement of I.

Considering then as output of the dynamics (1), we will see hereafter that successive time differentiations of the measurement y will involve I and S. Thus, the dynamics (1) with the only output is not fully observable and at most I and S can be estimated in addition to the measurement I. On the other hand also the number of deaths due to COVID-19 in the I population can be measured, while the number of recovered people cannot clearly be estimated. In the present section we will thus investigate the observability properties of the given system starting from the measurement of and we will show that given the fact that D is measured and R cannot be estimated, the problem reduces to the study of the observability properties of a subsystem of order 3.

Now, due to the presence of delays, which in this context are considered constant, we can take the differential representation of the dynamics in order to study its behaviour, using the approach introduced in Califano and Moog (2020). To take into account the link between the delayed variables, the backward shift operator δ has to be considered, see Xia, Marquez-Martinez, Zagalak, and Moog (2002). Let us denote by the field of causal meromorphic functions with s ∈ IN. Given a function denotes the function shifted by τ, that is . Let dx(t) denote the differential of x. Then, thanks to the back shift operator δ, . Accordingly, given the functionits differential form can be written in concise form as

Given :Finally is the (left) ring of non commutative polynomials in δ with coefficients in . A general module spanned by the differentials of functions in is then defined over the ring as in Xia et al. (2002). Then, in this framework setting and assuming without loss of generalities ℓ1 and ℓ2 integers, our system is characterized by the differential representation

Compute the observability matrix which is given bywithand Clearly for any j ≥ 0 does not depend on dx 2 and dx 3 which proves that the whole system cannot be weakly, regularly or strongly observable. We may however be interested in studying the reduced system defined by the variables and given by

Using the previous computations, one gets that the associated observability matrix isThis subsystem will then be weakly observable if has full rank over . It will be strongly observable if it is also unimodular. If it is weakly observable but not strongly we will have to check if it is regularly observable, see Califano and Moog (2020), that is if it is possible to reconstruct the state of the system by using also derivatives of higher order of the output function1

We may essentially distinguish three cases based on the values of the parameters ϵ and η which are discussed hereafter and are represented in Table 1 .

Table 1

Observability properties of subsystem (2).

	Weakly observable	Regularly observable	Strongly observable
ϵ=0	for ηβx_5(−1)≠0	×	×
ϵ ≠ 0, η=0	for βx4x5 ≠ 0	for βx4x5 ≠ 0	for βx4x5 ≠ 0
ϵ ≠ 0, η ≠ 0	for βx5 ≠ 0	×	×

First Case .

This case corresponds to the situation in which there is an important delay between the time people get infected and the time they become aware and get into quarantine. As a consequence, necessarily η ≠ 0 and the observability matrix associated to the third order subsystem (2) will beThe subsystem (2) will be then weakly observable for η ≠ 0, β ≠ 0, .

It is easily verified that it will neither be strongly nor regularly observable.

Second Case ϵ ≠ 0 and .

In this case there is no structural delay between the two classes I and I. People get infected and after a negligible time are moved to quarantine. The delay characterizes only the large amount of time that ill people need to get recovered. In this situation, for the subsystem (2) one getsNo delay affects the observability matrix. The subsystem (2) will be strongly observable if and only if the matrix has full rank, which happens if and only ifthat is whenever β ≠ 0, x 4 ≠ 0, x 5 ≠ 0. The last condition in factis always satisfied since x 5 ≥ 0.

Third Case ϵ ≠ 0 and η ≠ 0.

To check it let us use Smith decomposition. We have thatwherewithThe subsystem (2) is then weakly observable for βx 5 ≠ 0. It cannot be strongly observable. After some tedious computations, it appears that it is not regularly observable either.

Data from a regional university hospital

The originality of this section is to consider a smaller granularity with a subsystem of the previous dynamics: a local hospital center in charge of Covid-19 patients. Its input is the number A(k) of admissions in the hospital on day k. It is then in principle less than or at most equal to I in model (1), and it is further split into the number I(k) of patients admitted in intensive care and the number C(k) of patients admitted in conventional hospitalization on day k. The model describes the evolution of the dynamics defined by the hospitalized patients I(k) and C(k) as well as the number of deaths D(k) in hospital and the number R(k) of patients who have recovered from the disease (completely or partially and continue their treatment outside the hospital).

In Fig. 3 a representation of the management of the COVID–19 patients in the hospital is highlighted. Each day k a set of (A) new patients arrive at the hospital admission centre. Depending on their health status they can be forwarded to the conventional care section (C) or the intensive care one (I). Patients move from intensive care to conventional care and then are released from the hospital (R) once their health allows. Some of them instead do not survive (D) once their condition worsen.

Fig. 3

Flow Diagram associated to the management of the COVID–19 patients in the hospital described by Eq. (3).

As already underlined data from Nantes University Hospital were collected for about 6 months from March 16th to September 17th, 2020. These data include the daily values of C(k), I(k), R(k) and D(k) valid at day k.

Figs. 4 and 5 depict the number D(k) of deaths and of recovered patients R(k) at day k, respectively. The value k = 0 corresponds to March 16, 2020.

Fig. 4

Cumulative deaths over time from Nantes University Hospital Dataset.

Fig. 5

Cumulative number of recovered patients from Nantes University Hospital Dataset.

Due to the format of the daily data, a discrete time model with four parameters is in order as follows.Raw data from the Nantes University Hospital (2020) include, on a daily basis, the number I(k) of patients in intensive care, the number C(k) of patients in conventional hospitalization, the cumulative number of deaths and the cumulative number of patients which have recovered since day 1. The daily Covid-19 admissions are not directly measured but are computed as

Thus, all those state variables are considered to be measured and the input A(k) is computed from the measurements over two days.

The new admissions A(k) are split into θA(k) which increases the number of patients in conventional hospitalization and in which are directly entering intensive care. A value of θ close to 1 is representative of a high level monitoring of the disease by the health care outside hospital. When θ is closer to 0, then a major flow of new admissions enter directly into intensive care and, the other way around, λI(k) is the number of patients leaving intensive care and entering conventional hospitalization.

The number of patients γC(k) who recover from the disease is assumed to be proportional to the number of patients in conventional hospitalization.

μC(k) denotes the number of patients who are moved from conventional hospitalization to intensive care.

The number αI(k) of daily deaths is assumed to be proportional to the number I(k) of patients in intensive care.

The daily number A(k) of new admissions highly depends on the social and medical environment, as well as on some political lockdown regulations or meeting restrictions (with some delay). Thus, A(k) may have a significant variability from one regional hospital center to an other. It is proportional to the number of individuals which are susceptible of becoming infected, in some standard SIR population model.

Note that the local scale model (3) models the dynamics of hospitalized patients. Taking the sum of in (3) over all regional hospitals allows to get the hospitalized patients, nationwide. The latter is a part of the population I in model (1), as not all detected patients are hospitalized. The sum of hospitalized patients and non hospitalized infected people equals in model (1). It will be argued next that some parameters that characterize a single hospital may vary with respect to time. These parameters may also differ consistently from one hospital/region to another, depending on several sociological conditions and local regulations.

While in the model (1) delays were considered to characterize the evolution of the infection among the population, since they play a fundamental role in particular on the observability of the dynamics, in the local model (3) which refers to the hospitalization of COVID-19 patients, delays could be avoided. On one hand, even if data are available daily, the discrete nature of the model allows to work on a large time window (namely 7 days, as it will be discussed in the next section for the identification of the parameters) which is comparable with the delays involved in (1). On the other hand, the eventual presence of a large delay could be handled by splitting the patients in intensive care or in conventional hospitalization in more compartments moving the patients from one compartment to the successive one as the disease evolves, thus avoiding again the use of delays.

Discussion

Note that the infection dynamics described by a standard or modified SIR model, such as the dynamics (1), are external to the hospital and feed the dynamics (3). Thus, the input variable A(k) of (3), is a percentage of I. The parameters α, γ, θ μ and λ of model (3) can be identified over a short or longer period so that their evolution can be tracked. The argument in this section is that these parameters are

easy to derive from standard data from the hospital;

are representative of the regional health status.

Identifiability of model (3) parameters

The death rate α and the recovery rate γ are easily obtained from the measured state variables, since:

From Nantes University Hospital (2020), it is not significant to complete an identification on a daily basis. So, let us identify over an horizon of h days, which is done rewriting the third and fourth equation of (3) for days k, ..., . Standard computations lead toand

Note that parameters α and γ are argued to be time-varying, so that Eqs. (4) and (5) rather compute an average value over h days.

The parameters θ and μ are not simultaneously identifiable in this model for constant sequences of A(k) and C(k). Nevertheless, for some input sequence such that all parameters are identifiable. In other words, the system is generically identifiable, or identifiable for almost all sequence of admissions and for almost all flows of patients.

The identification of θ may be processed as follows. Rewrite the first equation of (3) at three different time instants:From (6), one can easily eliminate the coefficients and λ and after some elementary but tedious computations, θ is computed as follows. Rename the following quantitiesso thatObviously, there exist some singular input sequences A(k) which cancel out the denominator of the right-hand side of (8): for instance, when the sequences A(k) and I(k) are constant, then formula (8) is inapplicable for the estimation of parameter θ.

Following a similar elimination process, it is possible to eliminate parameters λ and α among (6). Thus μ is identifiable for almost all input sequences, though the explicit expression becomes a bit more involved. The generic identifiability of λ is proven following similar lines.

Computation of the death rate α and the recovery rate γ from the hospital dataset

Data from Nantes University Hospital were collected from March 16th to September 17th, 2020 and are available in Nantes University Hospital (2020). These data include the daily values of C(k), I(k), R(k) and D(k). In France the lockdown took place from March 17, 2020 to May 11, 2020. When some data were missing or inaccurate then the corresponding periods were skipped.

Since α is computed from the available data over 7 days, its evaluation displayed in Fig. 6 starts on the second week. The empty information, instead, from May 20th to 28th corresponds to errors which are due to inaccurate data, such as the number of admissions which were sometimes incorrectly estimated in the hospital data, thus not allowing a correct estimation of α.

Fig. 6

Time varying death rate α in model (3).

In Fig. 6 it is shown that the death rate increased during the early days after the lockdown. Two weeks later it was significantly decreasing and then it stabilized for over two months.

Surprisingly, the death rate α seems to reach unprecedented levels weeks after the end of the lockdown. A peak of the death rate α is noticed in early June, about one month after the end of the lockdown. This paradoxical result is easily interpreted as the general younger population left (intensive care) hospitalization and only elderly patients were remaining in intensive care, sometimes after several months of hospitalization with weak probability of health improvement.

From mid-June to mid-August, no patient was in intensive care which yields a zero death rate. This is due to a very low number of new Covid-19 patients admissions. Since the end of July 2020, there are new conventional hospitalizations, followed by intensive care hospitalizations at the beginning of August 2020. This indicates an increase in the circulation of the virus and eventually the beginning of a second epidemic wave. Parameter α requires a period of 7 days to be estimated, so that the increase in mortality is confirmed on mid-August.

In Fig. 7 a higher level of the recovery rate is noticed during the first month of the lockdown, i.e. before the peak of the infection was reached. Once the spread out of the infection was under control, the recovery rate was somehow stabilized.

Fig. 7

Time varying recovery rate γ in model (3).

Population levels were low in July 2020: no patients in intensive care and a decrease in the number of patients in conventional hospitalization from 23 patients at the beginning of the month to 2 patients at the end of the month. Therefore the recovery of one or two patients will have a significant impact on the recovery rate estimate.

Similarly in August 2020, there is an increase in the recovery rate due to a low headcount. These results point to the suspicion of a better patient management by the physicians in connection with an earlier treatment and to a massive screening. Summarizing, a higher recovery rate in August/September is due to a better knowledge of the disease by the physicians.

Conclusion

The main message in this paper is certainly that fast, agile, decentralized and relevant actions can be taken on some decentralized local level both to predict a new wave of the COVID-19 pandemic and to assess the real-time efficiency of treatments or of the health policy.

A local mathematical model of flows of patients will not stand in competition with a standard or upgraded centralized SIR type model but it will give additional insights to the dynamics to help decisions at the local level of a regional hospital.

The information provided by local data can be easily collected from any hospital. This information is rich enough to capture essential indicators about the local health status of the regional population. Not only this information is easy to obtain and does not require the involvement of any national health agency, but it is also more precise in the sense that it gives a smaller granularity picture of the health situation. It is well-known that this situation can be dramatically different from one region to the other; for instance, the Bergamo region in Italy had a very different status when compared with the south of Italy as shown in Alicandro, Remuzzi, & La Vecchia, 2020. Similarly, the Mulhouse region in France was one of the major clusters, see Kuteifan, Pasquier, Meyer, Escarment, & Theissen, 2020, with a much higher infection rate than western french regions along the Atlantic coast.

Data were collected during and after the lockdown period from the Medical Information Department of Nantes University Hospital Center. Such a low cost information is able to provide an early alert about a second wave of the infection. In this sense, it is an agile and precise tool able to complement official figures published by central health agencies.

Declaration of Competing Interest

We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.