Monitoring and Tracking the Evolution of a Viral Epidemic Through Nonlinear Kalman Filtering: Application to the COVID-19 Case

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 26, NO. 4, APRIL 2022 1441 Monitoring and Tracking the Evolution of a Viral Epidemic Through Nonlinear Kalman Filtering: Application to the COVID-19 Case Antonio Gómez-Expósito , Fellow, IEEE, Jose A. Rosendo-Macías , Senior Member, IEEE, and Miguel A. González-Cagigal Abstract—This work presents a novel methodology for systematically processing the time series that report the number of positive, recovered and deceased cases from a viral epidemic, such as Covid-19. The main objective is to unveil the evolution of the number of real infected people, and consequently to predict the peak of the epidemic and subsequent evolution. For this purpose, an original nonlinear model relating the raw data with the time-varying geometric ratio of infected people is elaborated, and a Kalman Filter is used to estimate the involved state variables. A hypothetical simulated case is used to show the adequacy and limitations of the proposed method. Then, several countries, including China, South Korea, Italy, Spain, U.K. and the USA, are tested to illustrate its behavior when reallife data are processed. The results obtained clearly show the beneficial effect of the severe lockdowns imposed by many countries worldwide, but also that the softer social distancing measures adopted afterwards have been almost always insufficient to prevent the subsequent virus waves. Index Terms—Nonlinear Kalman filtering, parameter estimation, Covid-19, geometric series. I. INTRODUCTION DESPITE the spectacular medical advances of the 20th century, and the practical eradication of viral diseases that in the past caused great mortality (e.g., smallpox), modern societies are still very vulnerable to the sudden appearance of new viruses, such as the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), cause of the coronavirus disease 2019 (Covid-19), for which there is still no vaccine. In addition, once a viral outbreak originates in a region of a country (in the case of the Covid-19, the Chinese region of Hubei, where the first reported case was dated on December 2019), the globalization of the economy and mass tourism spread it almost inevitably and quickly to the rest of the world. Manuscript received May 11, 2020; revised October 15, 2020 and November 5, 2020; accepted February 24, 2021. Date of publication March 3, 2021; date of current version April 13, 2022. The Miguel A. González-Cagigal thanks the financial support of the Spanish Ministry of Education and Professional Training under Grant FPU17/06380. (Corresponding author: Antonio Gómez-Expósito.) The authors are with the Department of Electrical Engineering, Laboratory of Engineering for Energy and Environmental Sustainability, University of Seville, 41004 Sevilla, Spain (e-mail: [email protected]; [email protected]; [email protected]). Digital Object Identifier 10.1109/JBHI.2021.3063106 In the absence of effective treatments, once a certain threshold has been passed, the main and almost sole remedy against the spread of the disease to the entire population is social distancing, the objective of which is to minimize the contact between people, and therefore morbidity [1]. In extreme cases, when the speed of propagation of the outbreak is very high, massive lockdowns of entire countries may be needed, which cannot last indefinitely owing to their drastic impact on the economic activity. For this reason, all the agents involved (governments, international organizations, institutions, companies and individuals) have the greatest interest in knowing how the number of affected and deceased people will evolve over time, with a view, on the one hand, to verifying the beneficial effects of social distancing, and on the other to scheduling the already saturated health resources and taking the economic measures intended to mitigate as far as possible the devastating effects of an epidemic like that of Covid-19. Scientists, engineers, economists, etc. are acquainted with several mathematical and statistical toolkits (recently renamed collectively as “data analytics”) for the treatment and filtering of time series, with a view to extracting useful information from the available data, uncertain by definition, such as trends, patterns, average values, expected variances, etc. In the specific case of a viral epidemic, such as that of Covid-19, there are basically two categories of models for processing the information: 1. Models that try to characterize the “physical” reality explaining the observed data. In the case of a viral epidemic, these models [2], [3] consider, for example, what fraction of people are at work, in teaching or travelling, how long it takes for an infected person to manifest symptoms, what is the mortality rate according to age groups, etc. This type of modeling is widely used in engineering, because the dynamics of the underlying systems or devices are generally well characterized, through mathematical relationships obtained from the physical laws that govern them (such is the case, for example, of electrical networks or an artificial satellite). 2. Models that try to determine explanatory parameters or variables from a purely mathematical point of view (“black box” approach), without going into the causes or interactions between components that explain the resulting data. Given uncertain data, which enter the system regularly (in our case, every day), the aim is to characterize © IEEE 2021. This article is free to access and download, along with rights for full text and data mining, re-use and analysis. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 29,2022 at 09:49:10 UTC from IEEE Xplore. Restrictions apply. 1442 IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 26, NO. 4, APRIL 2022 its temporal evolution by adjusting the parameters of a mathematical model, so that the differences between what is observed and what is estimated are minimized. Two variants can be considered in this category: a) Mathematical models that do not assume a priori what the shape of the temporal evolution of the involved magnitudes will be, but rather use a state transition equation, which tries to capture the dynamics of the problem in question by relating the variables in an instant of time to the variables in the previous instant. In this case, it is a matter of determining how the coefficients that define this equation evolve over time. In the case of epidemics, among the most popular models are those derived from the SIR (Susceptible-Infected- Recovered) model, [4], such as the one used for example in [5] to analyze the evolution of Covid-19 in Italy. This model is also considered in [6], where the evolution of the epidemic is forecasted using a novel state filtering algorithm. b) Mathematical models based on the assumption that the evolution of infected people, deceased, etc. obeys a predetermined curve (based on the experience of previous epidemics), whose coefficients are estimated based on the time series of reported data. For example, the evolution of the accumulated number of infected people can be satisfactorily approximated by means of a sigmoid curve, as assumed in [7], where the curve proposed by Gompertz [8] is used. The methodology proposed in this work belongs to the second category. As explained in the next section, we depart from the basic SIR model, by considering that the number of susceptible people, being large enough and changing relatively slowly, does not have to be explicitly considered in the model, but can be rather embedded in other equally significant parameters, such as the time-varying ratio of the geometric series characterizing the progression of affected people. Moreover, the proposed model explicitly distinguishes between people who have proved positive in a test, and people actually infectious, who are many more and for whom there is no reliable information available. In this work, a Kalman filter (KF) is used to process both the assumed dynamic model and the information available throughout the outbreak. The KF, proposed for linear dynamic systems in the early 1960s, is considered one of the fundamental tools that allowed men to walk on the moon, as it was successfully used in guiding the Apollo program space missions [9]. This filter, which constitutes a generalization of the technique known as “recursive least squares”, estimates the maximum likelihood evolution (that is, the most statistically probable, according to the assumed uncertainties and the observed samples) of the state of a dynamic system, and can be generalized to the non-linear case, including situations where the model parameters are also to be estimated. Reference [10] applies the KF for the estimation of the evolution of AIDS, while several recent studies related to the Covid-19 have arisen. In [11] the KF is used to deal with the estimation of the reproductive number of the virus. A shortterm prediction model is proposed in [12], where the time update equations of the estimator are used for future forecasts of the pandemic spread. ARIMA models are combined with a KF in [13] to track the evolution of the Covid-19 in Pakistan. Unlike in those references, where the parameters involved in the state estimation process are supposed to be known, in this work such assumptions are not required. This is the major distinguishing feature of the proposed methodology, compared to the state of the art, and the main contribution of the paper. II. PROPOSED MODEL We start from the well-known and simple SIR model [4], mathematically described by: ˙ Sc(t)=−β·Sc(t)I(t)/N ˙ I(t)=β·Sc(t)I(t)/N −γ·I(t) ˙ Rc(t)=γ·I(t) where Sc(t)and Rc(t)are, respectively, the cumulative or total susceptible and recovered people, I(t)represents the active infectious (not to be confused with cumulative infectious), β and γare the transmission and recovery rates, and Nis the total population of the studied region, satisfying N=Sc(t)+ I(t)+Rc(t). Note that, in this compact model, the deceased cases are paradoxically included in Rc(t)(alternatively, they could be subtracted from N). For practical purposes, the discrete counterparts obtained by numerical integration (forward Euler) are rather of interest. Moreover, as dead people are separately reported, they can be explicitly modeled, leading to a discrete-time SIRD (Susceptible-Infected-Recovered-Deceased) model, as used in [14]–[16]: Sc(n+1)=Sc(n)−β·Sc(n)I(n)/N (1) I(n+1)=I(n)+β·Sc(n)I(n)/N −(γ+μ)I(n)(2) Rc(n+1)=Rc(n)+γ·I(n)(3) Dc(n+1)=Dc(n)+μ·I(n)(4) where nis the elapsed time (in days) from a given origin, Dc(n) is the cumulative dead and μis a mortality ratio. The data publicly reported (available in references such as [17], [18]), typically comprise the following three items: rFraction of infectious people who, subject to a test, yield a positive outcome. This considers the fact that there may be many more infected than those reported positives, as happens with a large number of asymptomatic people. The cumulative positives will be denoted Pc(n). rFraction of recovered people who have been previously identified as positive. For simplicity of notation, the same symbol as in the basic SIR model, Rc(n), will be used in the sequel, even though we are referring here to a subset of Rc. rCumulative number of deceased, Dc(n), which is assumed to be the same as in the SIRD model, even though Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 29,2022 at 09:49:10 UTC from IEEE Xplore. Restrictions apply. GÓMEZ-EXPÓSITO et al.: MONITORING AND TRACKING THE EVOLUTION OF A VIRAL EPIDEMIC 1443 the actual number of dead by the virus may differ from the reported figures. Some sources directly provide the active positive cases, P(n), defined as: P(n)=Pc(n)−Rc(n)−Dc(n).Note that both I(n)and P(n)tend to zero as nincreases sufficiently (end of the viral outbreak), while the remaining cumulative magnitudes asymptotically reach a maximum or steady-state value. Epidemiologists use the so-called basic reproductive number, R0(average number of people infected by a single infectious person during the infective period at the onset of the outbreak) to characterize whether and how fast an epidemic spreads at the very beginning. If R0>1, then the epidemic will progress exponentially. As time elapses, though, the number of susceptible people decreases, either by the virus evolution itself or as a consequence of social distancing measures, and R0is replaced by the effective reproductive number, Rt. In terms of SIRD coefficients, R0and Rtare given by: R0=β/(γ+μ); Rt(n)=R0Sc(n)/N However, as thoroughly discussed in [19], the basic reproductive number R0is not free from ambiguity and controversy. For instance, it is stated in [19] that “using R0as a threshold parameter for a population-level model could produce misleading estimates of the infectiousness of the pathogen, the severity of an outbreak, and the strength of the medical and/or behavioral interventions necessary for control.” Moreover, if R0is estimated from time series of reported data, as in [11], then there is no way, at least for a new virus such as Covid-19, to subsequently check or contrast the accuracy of the estimates. This probably explains the wide confidence intervals so far reported for R0values [20]. Similar arguments apply to Rt. For this reason, instead of or in addition to R0, we postulate in this work the use of a more intuitive and measurable index, related with the growth rate of the infected class, to duly and unambiguously characterize a viral epidemic. Let the daily evolution of the active infectious be expressed as a geometric time series: I(n+1)=r(n)·I(n)(5) where r(n) is the time-varying ratio of the series. Then the daily growth rate is obtained from: Growth rate (p.u.)=ΔI(n)/I(n)=r(n)−1 Clearly, as long as r(n)>1, the viral outbreak will continue its expansion, whereas the disease extinguishes when r(n)<1. There is no ambiguity in using r(n) as a threshold, when referred to a whole population. Note however that, if (5) were expressed in terms of cumulative magnitudes, rather than daily or active cases, then r(n) would tend asymptotically to 1. By direct comparison of (5) with (2), the following relation is obtained, r(n)=1+β·Sc(n)/N −(γ+μ) or, in terms of Rt: r(n)=1+(γ+μ)[Rt(n)−1] Given r(n),one still would have to guess the values of the parameters involved in the SIRD model (1)–(4), to obtain R0. We contend that there is no need to worry in the short term about R0,asr(n)suffices to duly track the epidemic evolution on a daily basis. This work is aimed at estimating, from the daily reported data, the evolution of r(n)and, as a consequence, the growth rate of the infectious people. Note that, if r(n)can be somehow estimated, then (1) becomes unnecessary. In our approach, the impact of susceptible people, a factor which varies smoothly, is also embedded into r(n). In order to take advantage of the reported numbers of positive, deceased, and recovered cases, the following relationships are considered, taking into account (3)–(4): P(n)=t(n)·I(n)(6) D(n)=μ(n)·I(n)(7) R(n)=t(n)·γ(n)·I(n)=γt(n)·I(n)(8) where t(n)is a testing or reporting ratio that models the fraction of those infectious who are subject to tests and yield positive, D(n)=Dc(n)−Dc(n−1) is the daily increase in the number of deaths, and R(n)=Rc(n)−Rc(n−1) is the daily variation in the number of recovered cases. From (7) at two consecutive instants, keeping (5) in mind: D(n+1)=μ(n+1)·r(n)·I(n) D(n)=μ(n)·I(n) and dividing: rD(n)= r(n)·Kμ(n)(9) where rD(n)=D(n+1)/D(n)is the ratio of consecutive daily deaths and Kμ(n)=μ(n+1)/μ(n)is in turn a ratio of consecutive mortality ratios. Similarly, from (6) and (5): P(n+1)=t(n+1)·r(n)·I(n) P(n)=t(n)·I(n) and dividing: rP(n)= r(n)·Kt(n)(10) where rP(n)=P(n+1)/P(n)is the ratio of consecutive daily positives and Kt(n)=t(n+1)/t(n)is a ratio of consecutive t(n)coefficients. Finally, from (8) and (5): R(n+1)=γt(n+1)·r(n)·I(n) R(n)=γt(n)·I(n) and dividing: rR(n)= r(n)·Kγ(n)(11) where rR(n)=R(n+1)/R(n)is the ratio of daily recovered cases and Kγ(n)=γt(n+1)/γt(n)is a ratio of consecutive γt(n)coefficients. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 29,2022 at 09:49:10 UTC from IEEE Xplore. Restrictions apply. 1444 IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 26, NO. 4, APRIL 2022 III. KALMAN FILTER APPLICATION AND IMPLEMENTATION From the above simple model, given by (9)–(11), a system of state equations can be formulated allowing the sequence r(n)to be estimated by means of a non-linear KF, such as the Extended KF (EKF), Unscented KF (UKF) or Ensemble KF (EnKF) [9], [21]. This type of filter is capable of estimating the dynamic evolution of both parameters and state variables. Even though its maximum likelihood has only been proven for linear problems, it is applied successfully in nonlinear problems, such as the one in hand. In our proposal, the state vector is composed of the infectious geometric ratio, r(n), not directly measured, and the gains Kμ(n),Kt(n)and Kγ(n),all of them initially assumed to evolve with a random walk. Therefore, the state vector is x=KμKtKγrt and the state transition equation: x(n)=x(n−1) + w(n) where w=[ wμwtwγwr]tis a Gaussian noise vector with a covariance matrix Q(n), which accounts for model errors such as possible time lags not duly considered in the model. The state is to be estimated with the help of a measurement vector, composed of rD(n),rP(n)and rR(n), along with the pseudo-measurements of Kμ(n),K t(n)and Kγ(n).Thus, assuming that the coefficients μ(n),t(n)and γt(n)change slowly, their variation ratios can be considered to lie around the unity. Therefore, the measurement vector is, z(n)=rD(n)rP(n)rR(n)111 t and the measurement equation, z(n)=h[x(n)] + v(n) where v=[ vDvPvRvμvtvγ]tis a Gaussian noise vector that models the measurement error with a covariance matrix R(n), and the nonlinear measurement function is h(x)=r·Kμr·Ktr·KγKμKtKγt With this formulation, the EKF is able to deal properly with the non-linearity arising in the measurement equation, performing at each iteration a linear prediction step, followed by the nonlinear correction step. So, the KF provides the sequence of states estimates ˆx(n), whose last component, the estimate of r(n), will be denoted as rK(n). The sequence rK(n)incorporates, in a statistically optimal fashion, the noisy information provided by raw data ratios, such as rD(n),rP(n)or rR(n), which results in a more reliable estimation than the raw data ratios themselves. After several initial trials, an enhanced model has been finally implemented, considering two improvements: 1) The random walk in r(n)can be advantageously replaced, while keeping the linearity of the state equation, by a linear prediction based on the recent past history, e.g., Based on r(n−1),r(n−2) and r(n−3).Inthis Case, x=KμKtKγrr n−1rn−2rn−3t and x(n)=A·x(n−1) + w(n)(12) where: A= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 1 1 4 3 −1 3 1 1 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,w= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ wμ wt wγ wr 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ 1) The information available at instant ncan be used to improve previous estimations. For this purpose, the Rauch- Tung-Striebel (RTS) smoother [9], is implemented in two steps: rForward pass performed with EKF. rBackward recursion smoother based on the linear state transition equation. Regarding the tuning of the filter, the diagonal elements of R(n)related to vD,vPand vRare self-tuned according to the respective sample variances of the last available days of rD(n), rP(n)and rR(n). The terms related to the pseudo-measurements are set to 10−4to increase their weights. The diagonal elements of Q(n)are set to constant values Qi,i =10 −3, according to the observed variance of rK(n)−rK(n−1), except in the case of history variables, where Qi,i =0. The initial values of x(0) and P(0), required by the filter, are not relevant as their effect quickly vanishes. The proposed values are: rK(0) and history variables initialized to the mean of the raw ratios, Kigains initialized to 1, and P(0) initialized to an identity matrix. Additionally, the results of the KF allow the estimation or prediction of two magnitudes of great interest: 1) The day in the future when the infection peak will be reached, np, which will occur when rk(np)=1.For this purpose, the sequence rK(n)can be fitted to a predetermined evolution in order to predict its future behavior from the past history. According to what can be observed empirically, we have fitted the sequence rK(n) to a decreasing exponential, characteristic of first-order systems. It is worth stressing that, at this point, we are talking about the peak of infected, the peak of positives being a proxy for it. 2) The estimated number of active infectious people, IK(n). This magnitude can be computed from those infected at a given reference day, n0, as follows: IK(n)= n−1  i=n0 rk(i)·I(n0)(13) This means that, besides the errors of the sequence rK(n), the uncertainty of IK(n)inherits that of the initial guess I(n0). Therefore, as happens also with the SIR model, and in fact with any other model based on differential equations, precisely estimating IK(n)based on the above expression requires that an accurate number of infectious Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 29,2022 at 09:49:10 UTC from IEEE Xplore. Restrictions apply. GÓMEZ-EXPÓSITO et al.: MONITORING AND TRACKING THE EVOLUTION OF A VIRAL EPIDEMIC 1445 TABLE I VALUE OF THE PARAMETERS IN THE SIMULATION people be known on a given day, which can be very challenging. For our results, we will obtain the initial guess from IK(n0)=P(n0)/t(n0), for an assumed value of t(n0).For instance, t(n0)=0.2if we believe the first day there were 5 infectious people per reported positive. IV. SIMULATION RESULTS In this section, the performance of the proposed implementation of the KF is tested on a set of simulated scenarios, where the SIR model described in Section II is considered for the propagation of a virtual virus. While the total simulation time is 90 days, a lockdown is assumed to take place at day 15. This restrictive policy is aimed to completely stop non-essential public mobility, resulting in a quick reduction of the transmission rate, β(which is assumed to evolve exponentially, from an initial value 0.5 to 0.09) and, therefore, also of the ratio r(n). The remaining simulation parameters are given in Table I. In all simulated cases presented in the sequel, artificial Gaussian noise has been added to the measurements used by the KF algorithm in order to represent a more realistic scenario where the reported information presents inaccuracies. A. Base Case In the base case, the testing ratio, t(n), is assumed to remain constant (except for the noise). The time evolution of the estimated geometric ratio, r(n), is represented in Fig. 1 along with the raw noisy measurements provided by the simulation and the actual value of r(n). Note that the estimation provided by the KF is very close to the simulated value, giving evidence of the good performance of the proposed method. Fig. 2 shows the benefit attained from the incorporation of the smoothing filter mentioned above to the basic EKF algorithm, in terms of more damped oscillations. To compare the proposed KF implementation with other methods customarily employed as filters in these cases, Fig. 3 shows the estimated value of r(n)along with the results provided by three moving-average filters respectively applied to each of the noisy measurements, rP,rDand rR. As can be seen, the KF more closely tracks the evolution of r(n). From the estimated r(n)sequence, and the initial testing ratio, t(0), an estimation is obtained for the evolution of infectious people, which is compared in Fig. 4 with the simulated value of I(n)and the reported positives. The maximum estimation error (around 4.5%) takes place, as expected, at the peak of the epidemic. Fig. 1. Estimation of r(n) in the base case. Fig. 2. Estimation of r(n) with and without the smoother. B. Error Assessment Once the proposed estimation technique is validated in the base case, where only Gaussian noise is considered, the effect of different error sources is studied in the following scenarios. rStep in t(n) An abrupt change is simulated in the testing ratio from t(n)=0.2 to 0.3 at day 25, representing an increase in the availability of the tests (this has been observed in practice in several countries). Fig. 5 shows the estimation of r(n), along with the simulated value and the measurements of the geometric ratios rP,rDand rR. Note that the step in t(n)is observed as an impulse in the ratio rP, which is quite effectively filtered out by the proposed KF implementation. The estimation of I(n)is represented in Fig. 6, where the actual simulated value is again very close to the estimated one, Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 29,2022 at 09:49:10 UTC from IEEE Xplore. Restrictions apply. 1446 IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 26, NO. 4, APRIL 2022 Fig. 3. Comparison of the proposed method with moving average filters. Fig. 4. Estimation of the infectious people in the base case. giving evidence of the good performance of the method in the presence of a step in the testing ratio. rDeviations in t(0) Finally, the last scenario considered in this section shows how errors in the initial guess of the testing ratio, t(0), with respect to the assumed value t(0) = 0.2, affect the results. Given that the errors in this factor only affect the estimation of the infectious people, I(n), the representation of the estimated r(n)is not repeated. Fig. 7 shows the estimation of I(n)for t(0) = 0.15 and t(0) = 0.25 (±25% error). The results in Fig. 7 clearly show the importance of having an accurate guess of the number of infectious at the onset of the outbreak, as this initial error propagates proportionally up to the peak. Note, however, that any epidemiological model, such as SIR, faces the same challenge. Fig. 5. Estimation of r(n) with a step on the testing ratio. Fig. 6. Estimation of the infectious people with a step on the testing ratio. V. CASE STUDIES In this section, the proposed KF-based estimation technique is applied to the real data reported by different countries. For convenience, the results have been divided in two subsections: 1) the first period, when massive and aggressive lockdowns occurred in most countries, denoted in the media as the “first wave” of the pandemic [22], and 2) the subsequent period, once the first lockdown was relaxed, characterized by one or more additional waves, usually interleaved with several de-escalation phases. A. First Wave A total of four countries have been considered in this period: China, South Korea, Spain and the U.K.. At the early stage of Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 29,2022 at 09:49:10 UTC from IEEE Xplore. Restrictions apply. GÓMEZ-EXPÓSITO et al.: MONITORING AND TRACKING THE EVOLUTION OF A VIRAL EPIDEMIC 1447 Fig. 7. Estimation of the infectious people for a range of t(0) values. Fig. 8. Estimation of r(n) in China in the first period considered. the pandemic, the information provided by these countries was sufficient to allow the application of the proposed methodology. Figs. 8–15 represent the estimated sequence of the geometric ratio, r(n), and the number of infected people, I(n),for the four countries. The KF implementation is tuned as described in Section III for the covariance matrices Qand R, and the initial values of the vector x(0) and the covariance matrix P(0). In order to estimate the number of infectious people, I(n), according to (13), an initial value for the parameter t(0) is needed. In absence of a better clue, t(0) = 0.2is considered in all cases, expect for Spain (see the discussion of this particular case below). The points for which r(n)=1(peak of the epidemic) are highlighted with a dot. The following remarks can be made from those results: rA different evolution of r(n)can be observed for the Asian countries (China and South Korea), where the effects of Fig. 9. Estimation of infectious people in China in the first period considered. Fig. 10. Estimation of r(n) in South-Korea in the first period considered. Covid-19 started earlier. Once the geometric ratio r(n)< 1,the trend for South Korea is to remain roughly constant throughout the considered period, whereas for China a certain rebound can be noticed after March 10. rThe estimation results obtained for Spain show an asymptotic trend towards r(n)= 0.95. A slight increase is observed in r(n)between April 10 and 15, probably influenced by a sudden increase in the number of tests. Regarding the parameter t(0), in the Spanish case we have taken into account the results of a massive seroprevalence test performed by the government in the first half of May [23], from which it was concluded that the total number of infected people was around 5.2% of the population (approximately 2.3 million Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 29,2022 at 09:49:10 UTC from IEEE Xplore. Restrictions apply. 1448 IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 26, NO. 4, APRIL 2022 Fig. 11. Estimation of infectious people in South-Korea in the first period. Fig. 12. Estimation of r(n) in Spain in the first period considered. people). In view of this valuable information, the initial value t(0) has been adjusted so that the cumulative number of infected people matches the result of the survey on the date it was released (May 13), leading to t(0) = 0.12. This provides the estimation of I(n)shown in Fig. 13, where a maximum value of the active infectious people of around 1.3 million can be noticed by mid-April. Fig. 16 represents the estimation of the cumulative infectious people for the Spanish territory, where the total number of infected people matches the results of the survey. rWith the available information in mid-May, some countries had already left behind the peak of the epidemic (i.e., r(n)<1). For those cases (China, South-Korea and Spain), a rather accurate early forecasting of the epidemic evolution can be made, around 10 to 14 days before the peak, by fitting a decreasing exponential to a window of Fig. 13. Estimation of infectious people in Spain in the first period considered. Fig. 14. Estimation of r(n) in U.K. in the first period considered. past estimated data. This prediction is shown with green dotted lines in Figs. 8–13. rRegarding the U.K., where the peak of the number of infectious people had not been reached in the period considered (i.e., r(n)>1), an exponential fitting (made between around mid-April and early May) and the corresponding extrapolation is considered for this case. According to such fitting, the first peak should have taken place in the second half of May, provided the social distance measures were not relaxed. Finally, Fig. 17 shows the evolution of the above-mentioned fitted exponential curves for different countries (Italy has been included in the representation in order to establish a more complete comparison), all of them represented from a common Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 29,2022 at 09:49:10 UTC from IEEE Xplore. Restrictions apply. GÓMEZ-EXPÓSITO et al.: MONITORING AND TRACKING THE EVOLUTION OF A VIRAL EPIDEMIC 1449 Fig. 15. Estimation of infectious people in U.K. in the first period considered. Fig. 16. Estimation of cumulative infectious people in Spain. threshold r(n)=1.2,so that the corresponding time constants can be easily compared. In light of this representation, it can be noticed that the reduction of the geometric ratio is faster in China (just 13 days from r(n)=1.2to r(n)=1), possibly as a consequence of a more severe lockdown, followed by Spain and South Korea (between 25 and 27 days to reach r(n)=1), showing similar trends, and finally Italy (40 days to reach r(n)=1). B. Subsequent Waves As the pandemic evolves, it becomes more difficult to properly report on a regular basis all the information involved in the estimation of active positives. Many countries (notably Spain) stopped reporting the number of recovered people, probably Fig. 17. Fitted geometric ratios from common threshold. Fig. 18. Estimation of r(n) in USA in the second period considered. owing to the remarkable increase in the number of asymptomatic positive cases, which never entered a hospital and hence never counted as recovered or dead. For this reason, it is not possible to accurately update the estimations of the geometric ratios of active positives, r(n),for some of the countries considered in the early stages. Instead, Figs. 18–25 represent the estimated geometric ratio, r(n),and the number of active infectious people, I(n),for four countries (USA, Italy, India and Brazil), all of them specially affected by the pandemic and still reporting the information required by the proposed estimation technique. Similar assumptions as in the previous section are made regarding the KF tuning. The following remarks can be made from the results presented in Figs. 18–25: rThe number of infectious people in the USA was on the verge of reaching a peak by the end of May (r(n)≈1). However, Figs. 18 and 19 show that, afterwards, r(n) Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 29,2022 at 09:49:10 UTC from IEEE Xplore. Restrictions apply.