A stochastic inventory model of COVID-19 and robust, real-time identification of carriers at large and infection rate via asymptotic laws.

A critical operations management problem in the ongoing COVID-19 pandemic is cognizance of (a) the number of all carriers at large (CaL) conveying the SARS-CoV-2, including asymptomatic ones and (b) the infection rate (IR). Both are random and unobservable, affecting the spread of the disease, patient arrivals to health care units (HCUs) and the number of deaths. A novel, inventory perspective of COVID-19 is proposed, with random inflow, random losses and retrials (recurrent cases) and delayed/distributed exit, with randomly varying fractions of the exit distribution. A minimal construal, it enables representation of COVID-19 evolution with close fit of national incidence profiles, including single and multiple pattern outbreaks, oscillatory, periodic or non-periodic evolution, followed by retraction, leveling off, or strong resurgence. Furthermore, based on asymptotic laws, the minimum number of variables that must be monitored for identifying CaL and IR is determined and a real-time identification method is presented. The method is data-driven, utilizing the entry rate to HCUs and scaled, or dimensionless variables, including the mean residence time of symptomatic carriers in CaL and the mean residence time in CaL of patients entering HCUs. As manifested by several robust case studies of national COVID-19 incidence profiles, it provides efficient identification in real-time under unbiased monitoring error, without relying on any model. The propagation factor, a stochastic process, is reconstructed from the identified trajectories of CaL and IR, enabling evaluation of control measures. The results are useful towards the design of policies restricting COVID-19 and encumbrance to HCUs and mitigating economic contraction.

infection rate or rate or new carriers contracting the disease in period t, a = steady state level

= U carriers at large at the end of time period t. See also U.

the pandemic due to SARS-CoV-2

exit distribution from the group of CaL to HCUs or directly to ICUs

outflow or exit rate from carriers at large (CaL) entering health care units in the beginning of period t. analytic expression in terms of IR in eq. (T12), in Table A, Appendix A (SM).

steady state value of E

exit probability distribution (g=1,2,…, ν, = probability to exit to health care units in period t = t* + T-μ + i, i = 1,2,..,ν, of a carrier at large infected in time period t*)

complementary cumulative distribution of g

health care unit.

intensive care unit.

infection rate = number of new carriers per time period

propagation factor, infections at time period t per item in CaL at the end of time period t -1

mean absolute percentage error

Markov chain - Monte Carlo

= τ = mean residence time in the CaL group of all carriers (early exodus and entries to HCUs)

operational research

severe acute respiratory syndrome corona virus related to the COVID-19 pandemic

susceptibles, exposed, infected, recovered, deceased (state variables of SEIRD models)

supplementary material in the web

probability of abandonment of the group of CaL in period t as early exodus = Ω(U + Ω) (independent from the exit to HCUs / ICUs)

center axis of the exit distribution (ED) to HCUs from infection time (time periods)

CaL (inventory) level = number of carriers at large at the end of period t. It includes infected, symptomatic, asymptomatic, carriers that will eventually be self-cured, carriers that will eventually perish from the disease and carriers that will be taken to HCUs for treatment. Given analytically in terms of the infection rate by eq. (T15), in Table A, Appendix A (SM).

steady state value of U

1-s = probability of not leaving CaL as early exodus throughout period t.

dimensionless rate of entry to health care units = E/a

fraction of deceased among early leaves, ,

steady state value of ε

mean residence time in CaL of alive carriers (population average) at the end of period t. It can be monitored from samples of symptomatic carriers stand alone, from ongoing tests in the general population.

steady state value of η

mean residence time in CaL of carriers that enter HCUs (population average) in the beginning of time period t = (carriers to health care units: E). Monitored from samples of entries to HCUs, E.

steady state value of θ

half spread of the exit distribution to HCUs and ICUs, time periods

= 2μ+1=spread of the exit distribution to HCUs, time periods

mean residence time in the CaL group, or mean residence time as CaL of all carriers, including symptomatic, asymptomatic, self-cured, entering HCUs and deceased (population average)

early leaves or early exodus from the CaL group = exiting carriers that do not enter HCUs = rate of disinfected, self-cured carriers and of deceased carriers not entering HCUs.

equal by definition

time period, current time

complementary cumulative distribution

Introduction

With more than 40 million cases and 1 million deaths reported by October 15, 2020, figures that grew to about 170 million cases and 3.55 million deaths by May 31, 2021 worldwide, COVID-19 has threatened human existence. The SARS-CoV-2 is highly contagious, effectuating explosive growth of the pandemic at a global scale (Fig. 1) (Yuan et al., 2020a). In parallel to the self-sacrificial, heroic struggle of the personnel in health care units (HCUs) an unprecedented scientific effort has been undertaken in the medical and pharmaceutical industry to develop efficient preventive and treatment means to our defence (Deftereos et al., 2020). Authorities were coerced into painful decisions and policy measures, including curfews that have stifled economic activity (Nicola et al., 2020; Tsiodras, 2020). Quantitative representation of epidemic evolution may assist towards organization and planning intervention and use of resources, especially in the early stages (Oliva, 2019). Up-to-date work has provided valuable analytics, hinging on the endogenous characteristics of contagion and propagation of epidemics, for tracking the disease and devising control measures. The ineffable mark on society of COVID-19 compels further research and advancement of quantitative tools.

Fig. 1

COVID-19: Construal as a time-varying, uncertain inventory system with stochastic inflow, randomly varying early exodus and delayed/distributed outflow with shifting delay and spread and time-varying fractions of the exit distribution.

An important problem in this vein and target of this research effort is cognizance of: (a) the number of the carriers at large (CaL) of SARS-CoV-2 and (b) the number of new infected individuals, that is, the infection rate (IR) or rate of new carriers. To this end a stochastic model is to be developed representing actual COVID-19 evolution and a method for identifying CaL and IR in real time. Denoted herein by U the number of CaL is crucial for evaluating the state of the epidemic and devising control action, since it reflects the extent of the spread of SARS-CoV-2 in the general population. The IR is also critical, since new, infected carriers enter CaL and contribute to its explosive increase and repercussions. Significant uncertainty and under-reporting is present in both figures (Krantz and Rao, 2020). CaL encompass any hosts that have not entered HCUs for treatment. They include infected, exposed, symptomatic, positive or negative testing and asymptomatic, recurrent, restricted, confined, under quarantine, self-cured, deceased carriers not corroborated by tests and patients not in HCUs. Asymptomatic carriers may reach 50-60% of all carriers according to field studies in Iceland and on the USS TH. ROOSEVELT (CVN-71) reported on April 17, 2020. CaL is the pool from which all incidences emanate, including highly infectious ones, self-cured (carriers of the Ig-G and Ig-M antibodies) severe cases and victims of COVID-19. CaL is viewed here as the host “inventory’ ('inventory' in the SARS-CoV-2 perspective) which is unknown and unobservable. The IR is unknown and unobservable as well. Viewed as the inflow to inventory and denoted by a, the IR consists the main target of restrictive measures (e.g., social distancing) aiming at controlling the IR and indirectly, at controlling CaL. A third variable, important to the health care system is the outflow from CaL (E in our notation) that is, the number of patients who enter HCUs, or, less often, who enter directly intensive care units (ICUs). Traditionally viewed as ‘customer’ arrivals, patient flow is associated with waiting times in hospitals or emergency departments. It has been investigated using queuing theory (Armony et al., 2015). In our perspective, E is the throughput of the ‘inventory system’, an additional concern of policy and proactive measures against COVID-19: increasing values of E, combined with long treatment periods, threaten any health care system with immediate collapse, as in the on-going pandemic.

The problem of concurrent identification in real-time of CaL, U and IR, a is addressed, using concepts of inventory theory (Axsaeter, 2000; Porteus, 2002; Zipkin, 2000) together with asymptotic relations between the key variables, namely, inflow, inventory level and throughput. Up-to-date work in epidemic evolution relies on mechanistic models of lumped, finite-dimensional systems of differential equations, with presumed structure and fixed parameters. No structure, fixed parameters, or explicit quantitative relations of endogenous dynamics are presumed in our phenomenological view. The inventory (CaL) is considered time-varying, distributed, stochastic and unobservable. The inflow (IR) is stochastic, unknown and unobservable. The throughput, or exit from CaL and entry to HCUs, E, is registered. The inventory system is decoupled from the group of patients in HCUs/ICUs, the dynamics of which are fundamentally different, featuring low, or no data uncertainty. The quest is: (a) for a quantitative representation of cases and national incidence curves via the uncertain inventory perspective and (b) for a method, which, based on real-time, low uncertainty data could provide reliable estimates of CaL and IR independently of any model. Besides its own merit, the stochastic inventory model will be employed to assess the efficacy of the novel identification method for CaL and IR, based on asymptotic laws. Driven by real-time data, the identification method is to be validated by actual data and model fit, yet it will not rely on the details of any model to provide reliable estimates of CaL and IR. Section 3 describes the problem and the main differences between our perspective and up-to-date models. Section 4 expounds the characteristics of the quantitative description and presents empirical cases, modeled by our inventory construal. They encompass single, e.g., China and multiple-surge patterns, e.g., US, in distinct empirical cases. Section 5 unveils the fundamental results, potentially useful in our endeavor for real-time identification and presents the expressions to be employed to this end. Section 6 presents the identification method and validation, followed by discussion in Section 7.

Literature review

Regarding existing work, a well-developed body of literature exists on dynamics of contagion and disease propagation (May and Anderson, 1991; Capasso, 1993). Inspired by the work of Daniel Bernoulli (1766) (Bernoulli and Blower, 2004) on smallpox and by the malaria models of Ross (1911, 1915, 1916) the legendary SIR model (Kermack and McKendrick, 1927) and further extensions (Hethcote, 2000; Altizer and Nunn, 2006; Brauer, 2008; Vynnycky and White, 2010) assume a compartmental view of susceptible, exposed, infectious, recovered and deceased (compartments, S, E, I, R, D, respectively). The view is similar to chemical reaction kinetics, where transition from one state to another is associated with a specific rate. The essential parameters of compartmental models include the (assumed constant) rates of contagion between various state variables (S, E, I, R, D) (Diekmann et al., 1995; Pell et al., 2018; Giordano et al., 2020) from which the basic reproduction ratio, R, is found (van den Driessche and Watmough, 2002; Chowell et al., 2007; Brauer 2017; Zhao et al., 2020) expressing the number of infected individuals per carrier, when all subjects are equally susceptible. Such models lead to single pattern peak of contagious disease outbreaks (Chowell et al., 2019). However, multiple non-periodic surge patterns and/or stable incidence patterns with periodic behavior may occur (Hoen et al., 2015; Chowell et al., 2015) followed by strong resurgence outbreaks. Short-term sequential forecasts for devising successive control policies should allow for multiple surge patterns, as in actual COVID-19. Several additional issues, including unobservable heterogeneity (Yan, 2018) delay in incidence reporting, spatiotemporal distribution and data uncertainty spurred scepticism regarding mechanistic finite-dimensional models and forecasting of transmission of contagious, lethal diseases (Allen, 2017; Funk et al., 2018). The need for further investigation and analytics is pronounced in COVID-19 (Wang and Jia, 2019; Sen and Sen, 2021) due to the high incidence and fatality rates, despite restrictions and economic sacrifices.

The present work will adopt an alternative perspective to structured models, without attempting to substitute for the detailed interpretation of internal contagion dynamics. A novel, non-compartmental, quantitative representation will be proposed, which does not rely on fixed-structure internal transmission between state variables. Instead, it views CaL as an overall uncertain inventory, an autonomous system, augmented by the (randomly varying) newly infected carriers and depleted by the random early exodus and by the non-stationary, distributed-lag exit to HCUs. CaL entails the myriad of interrelated, randomly evolving, unobserved factors that influence the propagation of SARS-CoV-2 at different spatial and temporal scales (Yan, 2018; Chowell et al., 2019). Such a view considers the infection tree(s) (Piontti et al., 2014) a ‘black box” without seeking a detailed branching process model (Jagers, 1975) or an underlying compartmental, or networked epidemic model (Keeling and Eames, 2005; Nadini et al., 2020).

In summary, up-to-date work is based on lumped, homogeneous, fixed-structure, linear or bilinear models with no dead-time, to represent the evolution of the dominant variable of contagious diseases (incidence rate) and of the resulting recovery rate, death-rate and discharge rate from HCUs. The models feature several compartments of endogenous dynamics and networked state variables, with constant kinetic rates and result in single-pattern evolution and peak of incidence rates. They identify presumably fixed kinetic rates, based on past data of the states (infected, recovered, deceased). Such data feature significant uncertainty and underreporting, regarding CΟVID-19, due to the presence of a large number of asymptomatic carriers and may be potentially obsolete, due to the recurrent, evasive nature and continuous metamorphosis of SARS-CoV-2. Our endeavour is to decouple CaL (the high uncertainty part) from patients treated in HCUs/ICUs and to investigate a stochastic, non-stationary inventory representation of CaL and IR. The quest is to capture the multiple pattern evolution of COVID-19, as effectuated by temporal and spatial variability of contagion, environmental conditions, viral mutations, control measures and vaccination. The identification method will be scoping directly on the values of CaL and IR in real-time, based on low uncertainty data.

Problem definition

The epidemic models of Kermack and McKendrick (1927) and advancements, including SIRD and SEIR models have contributed significantly in understanding and representing quantitatively the spread of several epidemics, e.g., malaria, cholera, dengue fever (Zirka virus) plague, ebola, tuberculosis and recently of SARS and SARS-CoV-2. Several shortcomings have been pointed out, mainly due to the complex, time-varying and random nature of the actual SARS transmission process, spatial heterogeneity, age differentiation, individual immunity factors, the recurrent nature of the disease and mutations of SARS-CoV-2. Additional deficiencies include environmental and demographic factors and the effects of interventions and concomitant behavioral changes. The large number of asymptomatic cases and the large numbers of self-cured and deceased in all age groups, who did not enter HCUs result in significant data uncertainty. The latter hinders parametrization of such models in real time for purposes of identification, prediction and control. A tentative list of possible issues and directions for improvement is as follows.

A mechanistic structure, similar to chemical reaction kinetics is assumed. The latter hinge on the premises of Brownian (random) motion of molecules in a homogeneous mixture, an assumption hardly met in societies. For instance, outbreaks in household communities spread more slowly and culminate at lower peaks than anticipated by homogenous mixing models such as SIR.

Existing phenomenological growth models only support single-peak outbreak dynamics, whereas real epidemics often display more complex transmission.

Effects of demographic heterogeneity (Liu and Stechlinski, 2012; Davila-Pena et al., 2021) and age stratification (Cantó et al., 2017).

Effects of shifting community trust (or mistrust) to mitigation measures and vaccination, and resulting population behavior and compliance (Yan, 2018; Ni and de Brujin, 2020).

Time delays and differing heterogeneous spatial and social factors render infectious disease transmission a truly distributed system. For instance, the death rate in the ICU in a large city in Western Greece was reported over 97% (02 June 2021) much higher than the national average (50%). The overall death rate among confirmed cases in the same region was reported around 30%; national statistics report a fatality rate of 3% since the outbreak.

The models assume that S + I + R + D = constant (it may not be valid in an evolving disease).

Separate effects on contagion dynamics of quarantined and confined, including self-recovered.

Deaths, not registered as COVID-19 deaths. Although this type of uncertainty has little effect on the structure of the SIRD model, since D is an outcome of the last compartment, it has a strong effect on the estimated parameters (rates) found from data of S, R and D.

Deaths of non-hospitalized. A number of carriers, mainly of older age, with parallel diseases (e.g., heart, liver or kidney malfunction) perish without being hospitalized. They reached 70% in Greece (2nd wave, Nov.–Dec. 2020, Nat. Statistics Bureau, 02 April 2021).

Infection is assumed immediate. The dynamic responses are shaped by the time constants (rates) which may represent initial periods with slow dynamics, but they do not include time-delays, i.e., the essential dead-times of incubation, measurement and action. Apart from the latent period after which the symptoms of the disease manifest, there are incident reporting delays, measurement delays (until tests detect the disease) and action delays.

Constant rates do not reflect varying conditions, public health interventions, restrictive measures, regulations and resulting complex systems, producing multidimensional outputs (Ghaderi, 2021).

The rate constants are determined from statistics (Cantó et al., 2017; Anastassopoulou et al., 2020; Lin et al., 2020; Li et al., 2020, Stewart et al., 2020, Stidham, 1974). On-going factors are not included. Reinitialization methods spurred by growing discrepancy and error accrual, include ad-hoc interventions.

Constant values of the rates are determined ex-tempora, by statistically averaging all past data S, R and D via minimization of least squares error. Such methods seeking time-varying parameters, lead to the true parameter values if the parameters converge to constant values. If the true parameters are time-varying (Calafiore et al., 2020a) they may lead to erroneous parameter values, trapped in a linear manifold, though minimizing the sum of error squares.

Real-time estimation of the states, S, E, I, R and D faces difficulties. Well known, procedures such as non-linear Kalman filtering algorithms for state estimation, in which updated estimates depend on current data and previous estimates, are applicable when the exact model is known. Due to structural uncertainty in the model and temporal variability in the rates, application of nonlinear state estimation to SIRD COVID-19 models may lead to accrual of estimation error.

The decisive effect of varying environmental factors (air temperature, humidity, wind) and seasonality (Lin et al., 2006; Groseth et al., 2007; Wu et al., 2016). Such conditions are crucial, particularly for zoonotic infectious diseases (e.g., SARS-CoV-2).

Effect of vaccination strategies (Sinha et al., 2021) and incidents of infected after vaccination.

No inclusion of recurrent cases (affecting susceptibles); they were included in the SEAIR model.

The widely differing time constants of self-cured (a high fraction of infected) symptomatic and asymptomatic, who recover in one to four weeks without being hospitalized.

The different transmission routings regarding patients in HCUs/ICUs than CaL and the effect of hospital infection and nosocomial propagation (Deftereos et al., 2020; Lai et al., 2020 ).

Effect of mutations and genetic evolution of SARS-CoV-2, rendering the rates time-varying.

Linear, first order dynamics, leading to single growth patterns are assumed in all state equations apart from the equations in S and E, which are assumed bilinear (The SIR and SIRD model are bilinear in S and the SEIR model is bilinear in S and E. The modified SIRD model is bilinear in S, E and A). Contagion dynamics are infinite-dimensional, spatiotemporally distributed and nonlinear, with unknown nonlinearities (Keeling et al., 2001; Korobeinikov, 2006) structural instabilities (van den Driessche and Watmough, 2000) and diffusion (Willis et al., 2020).

Several improvements have been suggested to the above. The modified SEIR model (Yang et al., 2020; Mwalili et al., 2020) attempts to improve on demographic homogeneity by considering geographical migration index (In and Out groups of the population to and from Susceptibles and Exposed, respectively). It assumes that the rate of transmission of E → S is five-fold that of I → S. The semi-mechanistic model of Ebola dynamics (Chowell et al., 2015) is a modified SEIR model, accounting for delays in reporting of infected and time-varying transmission. Tsay et al. (2020) presented a minimal cost optimization method for controlling the spread of COVID-19, taking into account the cost of restrictive measures, with a-priori set relative weighting factor. The procedure is based on the SEAIR model with ex-tempora parameter estimation via minimization of the least squares of the estimation error (LSE). The SEAIR model allows for asymptomatic and recurrent and notes the time- variability of several rates. Recovered from asymptomatic are considered and estimation of hidden states (A and E) is proposed by re-initialization based on the (uncertain) data. However, besides relying on the structure and assumptions of SIR, it is well known that if used for real-time estimation, LSE methods provide moving, time averages of the time-varying parameters they attempt to estimate, but they fail to track peaks and valleys. The latter are of crucial importance in COVID-19 (Ahumada et al., 2020). Calafiore et al. (2020b) mention significant problems with the identification of parameters of dynamic models, including the high sensitivity on the uncertain initial conditions, the non-convexity of the optimization problem and the poor performance of the resulting local solutions. They elicit the non-valid assumption of known infected cases. They proposed a SIRD model, in which a constant fraction of infected is known from observations and tests and recast the parameter estimation to a convex optimization problem. Chowell et al., (2019) proposed superposition of multiple sub-epidemics described by a structured generalized logistic growth model, with exponential and polynomial growth patterns that represent overlapping sub-epidemics. Rubio-Herrero and Wang (2021) proposed rolling regression for time-varying SIRD models. The effect of vaccination was included via additional compartments and states in Dashtbali and Mirzaie (2021). Based on the models above, feedback control theory and methods including stochastic optimal control (Tsay et al., 2020; Lesniewski, 2020); Choi and Shim, 2021) and model predictive control (Stewart et al., 2020; Alleman et al., 2020; Morato et al., 2020; Kohler et al., 2021; Carli et al., 2020; Sen and Sen, 2021) were used to derive step-wise optimal distancing policies. However, there are no guaranteed stability margins (stabilizing effect) of fixed-structure, linear quadratic optimal control measures (Doyle, 1978) in the face of model uncertainty. Time-varying rates, delays, limited or biased measurements due to asymptomatic cases, render such policies even less effective, to optimally balance public health and economic contraction.

Debarring from finite-dimensional model predictions, Pastor-Satorras et al. (2015) elicited the potential of coevolving, time-varying networks to represent propagation. Allen (2017) proposed a stochastic differential model using continuous time Markov processes to account for demographic variability and environmental factors, related to terrestrial or aquatic settings. It allows for random behavior of the parameters, yet its basic structure is the SIRD model. Other stochastic epidemic models (Greenwood and Gordillo, 2009; Funk et al., 2018) assume the structure of basic mechanistic models, with inclusion of latency periods (Karako et al., 2020) asymptomatic cases (Bardina et al., 2020) and government control measures (Chang et al., 2020; He et al., 2020). Zheng et al. (2021) suggested abruptly changing, piecewise constant infection rates. Taylor and Taylor (2021) suggested combination of probabilistic forecasts based on reported data in order to predict hospitalization and mortality rates. Dynamical density functions were proposed in te Vrugt et al. (2020) to model the effect of social distancing.

Modeling and formulation: a stochastic inventory model for COVID-19

Figs. S1 and S2 in the Supplementary Material (SM) in the Web present COVID-19 national incidence profiles and SEIRD and SEAIR epidemic models. A discrete time approach is adopted in our stochastic inventory view, (Tsiliyannis, 2016; 2017; 2018; 2020). The underlying idea is to replace the continuous time with discrete points in time, at which the system state is observed (Schwarz et al., 2016). One transition per interval is possible, which, depending on the system and the length of the discretization interval (e.g., daily) may account for multiple incidences, i.e., new SARS-CoV-2 carriers (arrivals in bulk queues) and/or multiple entries to HCUs (“service completions” in queueing systems). Two novel concepts will be employed towards realistic representation of empirical data of COVID-19 evolution. (i) Inclusion of random inventory losses (Tsiliyannis, 2016; 2017) or abandonments of queue -the phenomenon in queueing systems that customers get impatient and renege from the system, when the waiting time exceeds their tolerance - (Feldman et al., 2008; Hampshire and Massey, 2010; Dietz, 2011; Kim and Ha, 2012; Bertsimas and Doan, 2010). Abandonments of queue with retrial represent customers choosing to go back to the service system after some time - (Kapodistria, 2011; Dai and He 2012; Ata and Tongarlak, 2013). (ii) Delayed and distributed exit to HCUs, with randomly varying fractions of the exit distribution (ED) and possibly shifting delay and spread. Both features are present in COVID-19, strongly affecting CaL and IR: losses represent the random exodus from CaL of no longer active carriers, disinfected, self-cured (symptomatic or asymptomatic) without treatment in HCUs and of deceased carriers that have not entered HCUs. The number of no longer infected, self-cured and deceased due to SARS-CoV-2 who do not enter HCUs is unknown and unobservable. Retrials represent cases when a disinfected, HCU cured, or self-cured carrier, is re-infected or re-contracts the disease, a realistic possibility stressed out for COVID-19. Denoted by Ω, the early exodus, is a nonnegative random variable, not identically distributed with the exit E. The early exodus and the IR are random, in the broader sense of a sequence of random variables, that is, a stochastic process (Gallager, 2014) rather than noise-like variations (e.g., white noise). Fig. 1 presents the system concept and main variables. The second novelty representing COVID-19 evolution regards carriers that enter HCUs. It reflects the fact that carriers infected at the same time may enter HCUs at different time periods in the future, according to individual condition and priorities. This empirical fact is modeled via an ED of carriers entering HCUs, which spreads in a finite number of time periods, say, ν days in the future. The delayed ED is centered at a future time, say T time periods after contagion and spreads in μ time periods prior and after T, i.e., from period up to period . Thus, any new carrier entering CaL in period t may exit either (i) as part of the early exodus in the time periods from t+1 up to period , or (ii) as a patient entering an HCU or directly an ICU for treatment (exit, E.) in the beginning of any of the periods . The delayed ED encompasses delays (dead times) and allows overtaking. With these basic considerations, the following issues will be explored and actions taken thereto. (a) Establish the quantitative framework and variables describing epidemic evolution. (b) Provide a rigorous quantitative representation in dynamic time of COVID-19, matching national epidemic curves. (c) Set-up the identification problem of CaL and IR in real time and ascertain the key variables/parameters and the minimum number thereof that must be monitored to this end. (d) Present an efficient identification method and test it against empirical incidence cases, closely represented by the stochastic inventory model.

Modeling of random early exodus or losses of the SARS-CoV-2 inventory

The random early exodus from CaL, , accounts for disinfected, self-cured and deceased, not entering HCUs. It is represented via the early exodus ratio, or the probability of leave, , which equals the ratio of the early exodus flow in period t, , divided by the pool from which this exodus flow emanates during period t, i.e., by , Eq. (T13) in Table A, Appendix A (SM). In general, the ratio may take random values in [0,1] but in practice . The early exodus events are assumed to be independent events and thus the probability of not leaving the group of CaL as part of the early exodus in the time span from t-j+1, to t, is given by . For instance, if x follows a uniform distribution in [0.94, 1] i.e., , as given by the Markov chain Monte Carlo (MCMC) procedure, then, on the average, after ten days a fraction 0.9710 = 0.74 (or 74%) of carriers, who were infected ten days ago and have not entered HCUs or ICUs, remain in CaL. The rest are no longer carriers of the virus, i.e., they are no longer infected, or have self-recovered (the great majority) or have been deceased.

Modeling of the delayed, distributed and non-stationary entry rate to HCUs

The ED itself does not have to be described by any of the well known (or any known) constant probability distributions, or probability mass functions. In fact, it is unknown and time-varying, switching, possibly, from left skew to symmetric, or to right skew, featuring varying kurtosis. This is taken into account in order to faithfully represent empirical cases of non-stationary exit distributions such as in COVID-19. Each fraction of the ED, g=1,2,…,v, 0 ≤ g ,  ≤ 1 is a sequence of random variables generated via MCMC simulation. Details of the representation and a simple example are given in Appendix B in SM. Both the center axis, T, and the spread, ν, of the ED may be shifting with time, thus embodying further structural uncertainty as well. The time variability, as expressed by the random early exodus, by the non-stationary ED of ‘served’ items and by structural variation, allows to accurately represent stochastic routings of disease transmission and uncertainty associated with COVID-19. A direct consequence of the distributed nature of the outflow to HCUs at time t, E is that it consists of carriers who have been infected in various time periods in the past. More specifically, it consists of carriers possibly infected in any of the time periods . For instance, with T = 10 days, μ = 5 days, E includes patients that may have contracted the disease in any of the time periods t-5, t-6,…, t-15, i.e., five to fifteen days prior to entering HCUs.

Modeling the evolution of CaL and IR

Within this framework and notation as defined, the system evolves in time according to Eq. (T1) in Table 1 . A fundamental relation of inventory dynamics, Eq. (T1) states that the ‘system’ of CaL at the end of period t has moved from level U 1 (end of period t-1) to level U (end of period t) by the net amount of new carriers, , where a is the IR during period t, E is the exit from CaL in the beginning of period t entering HCUs and Ω is the early exodus from CaL during period t. Subject to different dynamics, patients in HCUs and ICUs are not included. As stated, the IR depends on the number of CaL. This can be expressed by a time-varying, linear relation (Eq. (T3), Table 1) which, for positive values of the propagation factor, k, (k represents infections in day t per carrier at the end of day t-1) induces a positive feedback mechanism, forcing the system to higher and higher values of . This is seen by substituting Eq. (T3) into Eq. (T1) to obtain a time-varying linear system with a time-varying pole at 1+k, (Kamen, 1987; 2010). CaL exhibits exponential growth for constant , hyper-exponential growth under rising and slower, polynomial growth for shrinking. Restrictive measures or vaccination lower k and for values tending to zero the pole tends to one (no growth). Then, the (stochastic) effects of early exodus and of the exit to HCUs become dominant, lowering the number of CaL and the IR. Allowing such arbitrarily random temporal-variability of a single parameter, namely of the propagation factor, k (Eq. (T3)) is a novelty in modeling actual epidemics, compared to SEIRD models, for it allows stepwise fine-tuning of k at the present time (t) given the history , for mimicking multiple-pattern evolution. Tuning accounts for shifting transmission mechanisms, behavioral changes in the community, extent, adoption and impact of control measures, genetic evolution of the virus, environmental factors, seasonality, etc.

Table 1

COVID-19: Evolution of CaL and IR under random early exodus.

U_t: CaL at the end of time period t.	U_t=U_t−1+a_t−Ω_t−E_t(T1)
Ω_t: random exodus, or random leaves (losses of Ut) during time period t; it includes no longer infected, self-cured and deceased carriers, who have not entered HCUs by the end of time step t.E_t: exit to HCUs (recorded).	0≤Ω_t≤U_t−1+a_t−E_tΩ_t=φr(U_t−1+a_t−E_t)(T2)0<φ<1, r = random [0, 1]φ=0.1, COVID-19
a_t: infection rate (new carriers infected in period t (inflow to CaL).k_t: propagation factor: definition.	a_t: a multiple-incidence discrete stochastic processa_t=:k_tU_t−1(T3)

The assumptions of the stochastic inventory construal for COVID-19 are summarized as follows. (1) Finite duration of carrying the virus (at most Τ + μ days) with possible exit to HCUs/ICUs from day Τ − μ up to T + μ; T and μ can be slowly varying (2) being a sequence of independent events, the early exodus from CaL, Ω, depends only on time t and is independent from the outflow to HCUs, E (3) the IR evolves according to Eq. (T3), Table 1.

Data and model validation: empirical cases and matching realizations

Several cases of COVID-19 national incidence evolution are presented, featuring distinct profiles. It is assumed that the entry time to HCUs, E, ranges from one to about four weeks. Then T − μ = 7 and T + μ = 27 days, from which the ED delay and spread range around T = 17 days and μ = 10 days. The early exodus probability and the randomly varying fractions of the ED enable to analytically express the inventory, U and the exit to HCUs, E, in dynamic time, in terms of the IR (expressions T12 and T15 in Table A, Appendix A in SM). The IR is varying according to Eq. (T3), Table 1. Alternative MCMC realizations of are generated by allowing random variability of x and k, some of which faithfully describe the spread of the virus, as expressed by the mean absolute percentage error (MAPE) metric. The efficacy of a realization to represent a real evolution profile is assessed by matching the cumulative number of cases by the integral of CaL, throughout the evolution period. Matching of the (daily) entry rate to HCUs, E, corroborates the efficacy of the realization. A matching realization provides the benchmark against which the identification method is tested (Section 6).

Greece

Fig. 2b presents a model realization matching COVID-19 reported data in Greece, from February 26, 2020 to May 10, 2021 (Fig. 2a) where the disease was initially effectively controlled. On day 24 (March 20 2020) tame preventive measures already taken, have been severely tightened with general curfews. As a consequence, the spread was harnessed (Figs. S3 and S4 in SM) with visible effects on IR starting about two weeks later (days 30–35). The measures were relaxed on May 04,

Fig. 2

COVID-19, Greece (a) official data. (b) A stochastic model realization of carriers at large (CaL) U, infection rate (IR) a, early leaves Ω and entry rate to HCUs, E, resulting in close representation of the reported entry rate to HCUs and of the cumulative number of SARS-CoV-2 cases (see fig. 3) in the period 26.02.2020 to 10.05.2021.

Fig. 3

COVID-19, Greece: Matching of reported data by the realization in fig. 2b in the period 26.02.2020 to 10.05.2020. Left: Matching of the reported daily rate of entries to HCUs. The MAPE between the reported rate and the model output is 9%. Right: Matching of the reported cases by the model in the realization of fig. 2b. The MAPE between the reported cumulative number of cases (dashed line) and the cumulative 0.1CaL of the model realization (dark blue line) is 3%.

2020. In regard to the model MCMC representation, the early exodus probability varied randomly in (0.00, 0.12] that is the probability of not leaving as part of early exodus at any day varied in the range [0.88, 1). The fractions of the ED were taken as random normal deviates around a quantized Gaussian ED, g1, 2,...,21, with standard deviations 100% of the respective mean g. The model output for the entry rate to HCUs (Fig. 3 a) matches actual daily data of E (Fig. 3a). Fig. 3b depicts the representation of reported cases up to day 440, i.e., from 26.02.202 to 10.05.2021, in order to compare with the reported data (Fig. 3b, dashed line). The model output for 10% of the cumulative number of cases, i.e., the cumulative 0.1CaL (time integral, or summation in discrete time, 0.1) matches the cumulative number of reported cases at any day t for t = 1,2,…,440: the profiles in Fig. 3b feature a mean absolute percentage deviation of ∼3% (MAPE). It follows from Figs. 2b and 3a,b that the realization closely matches the reported data of entry rate to HCUs and of the cumulative number of CaL in days 1–440.

Regarding the evolution of the epidemic, eight distinct phases are evident in the profiles of CaL and IR: (i) the outbreak, hypergeometric growth phase, which lasted for less than a month (up to mid March 2020) in which the IR rose to ∼16% of CaL, i.e., a = 0.16 U 1 (Fig. 4 a) (k = 0.16 in Eq. (T3), Table 1) resulting in fast spread of the disease to ∼1000 CaL and IR to ∼120/day (ii) a descending phase, lasting for about three weeks (March–mid April 2020) in which the propagation factor k utilized in the model was sharply lower (k∼0.1) and even lower (k∼0.065) (iii) the decline phase (April - mid May 2020) with k ∼0.08, in which the IR descended to about 20/day and the number of CaL dropped to ∼50, (iv) the leveling phase (late May–June 2020 up to day 120) corresponding to relaxed measures, during which a higher k was necessary (k ∼0.1) resulting in (iv) a slowly resurging phase (late July–September 2020 until day 210). This was followed by an explosive phase, up to day 270, due to imported cases from visitors, repatriates, tourists, vacationers and summer activities, in which CaL rose to 9000 and IR to 1300/day (k rose to ∼0.14). (v) A tame period, up to day 330 (k ∼0.10) in which CaL and IR fell to 7000 and 800/day, respectively, followed by (vi) a strong resurging, oscillatory phase, after the winter holidays, which was weakly controlled via severe restrictive measures throughout winter and spring and in which CaL peaked at 14000 and IR at 1600 (k varied randomly around k ∼0.10). (vii) A multiple wave pattern in April 2021 during Lent and Easter holidays (CaL peaking at ∼19000, IR∼1700) and persisting in May 2021 despite 25% vaccination, was slowly subdued (phase (viii): June 2021).

Fig. 4

National profiles: propagation factor, k (infected in day t per carrier in CaL at the end of day t-1).

COVID-19 evolution in China

Fig. 5 depicts a model realization (top) of carriers at large (CaL) infection rate (IR) early exodus and distributed exit/entry to HCUs for China, in the period January 01, 2020 to November 27, 2020. The early exodus probability and the fractions of the ED in the model were randomly varying as in 4.4.1 above. The realization output for the cumulative 0.1CaL in Fig. 5 (0.1) matches closely the cumulative number of reported cases for t = 1,2,…,330, i.e., from 01.01.2020 to 27.11.2020: the profiles in days 1-330 feature an average deviation of ∼4%. The factor k utilized in the model (Eq. (T3), Table 1) rose initially to 0.53 (Fig. 4b) reflecting the outbreak. Severe control measures imposed after two weeks, tamed the spread, with visible effects starting about two weeks later (days 35–40) Fig. 5. Lower values of k were used to represent the effect of the control action: k was reduced to 0.20 and then, after a short rise, to 0.025 (Fig. 4). The explosive growth of CaL and IR ceased. The IR dropped sharply from 16000/day thereafter and remained negligible, although randomly varying. CaL also fell from 150000 and remained at low levels ever since. Two distinct phases of evolution are evident: (i) explosive growth and (ii) decline due to lock-down. A recurrence phase that emerged around July 2020 was controlled by vaccination. Up-to-date models may reproduce the single pattern evolution of the epidemic up to the recurrent phase.

Fig. 5

COVID-19 China: A model realization (top) of carriers at large (CaL) U, infection rate (IR) a, early leaves Ω and entry rate to HCUs, E, resulting in close representation of the reported cumulative number of SARS-CoV-2 cases (bottom) in the period 01 January through 27 November 2020). The MAPE between the actual reported number of carriers and the corresponding outcome from the MCMC realization (dashed line in embedded fig.) is 4%.

COVID-19 evolution in the USA

Fig. 6a depicts a multiple pattern (oscillatory with increasing amplitude) model realization for , regarding COVID-19 in the US, from the start of the outbreak, up to the first days of December 2020 (20.01.2020 to 03.12.2020) fitting actual data, as manifested by the embedded Fig. 6a. Fixed structure, time-invariant models can hardly reproduce such a multiple-peak evolution profile. The early exodus probability and the ED were taken as in the previous cases, 4.4.1 and 4.4.2. The realization output for the cumulative number of cases in the embedded Fig. 6a, i.e., the time integral, or summation of 0.1CaL in discrete time, for t=1,2,…,320, closely matches the reported data of the cumulative number of incidences in days 1-320 (20.01.2020 to 03.12.2020): the profiles in the embedded Fig. 6a feature a mean absolute percentage deviation around 4%. The propagation factor k utilized in the model (Fig. 4c) was initially > 0.3, reflecting the explosive growth that led to high CaL (∼300000) and IR (∼50000) by day 90 (April 2020). Lower values (between 0.28 and 0.02) used to match the evolution in May 2020 were increased up to k ∼0.15 to represent the swelling wave in July 2020. Oscillatory trends followed until day 320 (Nov.27.2021) with the peak corresponding to k ∼0.26. As depicted in Fig. 4c, the ranges of k were as follows: 0.02 to 0.1 in the spring-summer of 2020, with a peak at 0.14 around day 200 (early August) and 0.02 to 0.26 in the fall of 2020. The oscillatory evolution of the epidemic led to explosive growth, initiating on day 260 (early October 2020) and culminating around day 300, after the 2020 elections, with high levels of CaL (around 1800000) and IR (around 300000).

Fig. 6

COVID-19: Top (a): USA. A realization of carriers at large (CaL) U, infection rate (IR) a, early exodus Ω and entry rate to HCUs, E, resulting in close representation of the reported cumulative number of SARS-CoV-2 cases (embedded) in the period 20.01.2020 to 03.12.2020. The model output for the cumulative 0.1CaL fits the reported number of cases with a MAPE ∼8% (embedded figure). Bottom (b): A realization of carriers at large, infection rate, early exodus and entrance to HCUs, matching the reported incidence data in Italy in the period 21.02.2020 to 27.11.2020. The model output for the cumulative number of SARS-CoV-2 carriers (cumulative 0.1CaL, dark line in embedded fig.) fits the cumulative incidence rate with a MAPE around 6%.

Italy: COVID-19

Fig. 6b depicts a realization of , of our COVID-19 model for Italy, from 20.01.2020 (active from 21.02.2021) until 27.11.2020, where the disease has taken a dramatic toll in human lives in the first three months of the outbreak. Regarding the model, the early exodus probability and the ED were taken as in the previous cases. The realization fits the reported number of cases: the model output for the cumulative number of cases in the embedded Fig. 6b, for t=1,2,…,310, matches the reported cumulative incidence rate with a MAPE of ∼5%. The profiles of include five phases: (i) the outbreak/exponential growth phase, which lasted for about seventy days, with the first peak at 50000 CaL and IR∼10000 (ii) a descending phase, lasting for about two months (April–May 2020) (iii) a second, lower wave (peaking by May end) (iv) the decline, leveling-off phase in June–August 2020 (assisted by summer temperatures) (v) an explosive growth phase, culminating in late November 2020 (CaL∼450000, IR ∼100000). The factor k utilized in the model (Fig. 4d) was initially around ∼ 0.20 reflecting explosive growth. It was reduced to 0.15 around mid March, then to 0.05, and increased to 0.33 by mid April. Subsequently, it varied randomly between 0.08 and 0.25. The surge of the second wave corresponded to the peak, k = 0.25. The multiple wave profile, faithfully reproduced herein, may not be modeled by constant parameter models (Calafiore et al., 2020a).

Identification of carriers at large and infection rate in real-time

Having presented the stochastic inventory model, reproducing pragmatic data, it is time to present the identification method, which exploits the fact that the entry rate to HCUs, E, is known. The method is based on a recent result (Tsiliyannis, 2016) which relates U, a and E asymptotically under random losses and non-stationary, distributed exit, E:where

θ is the mean time (population average) spent as CaL of the patients in E, that is, of patients who are exiting CaL and entering HCUs in the beginning of time period t (θ = mean duration of sojourn in the group of CaL of all members of E)

η is the mean time (population average) of alive CaL (η = mean duration as CaL of all alive members of U at the end of period t). It can be monitored as the mean time of carrying the virus by symptomatic carriers stand alone, during ongoing tests on the general population.

Eq. (1) can be used directly to identify any of the two variables (CaL or IR) based on the other one (if known). For instance, if in a community hit by SARS-CoV-2, the number of CaL is approximately known from extensive tests, then the IR can be directly determined from Eq. (1). The next and greater challenge however is to concurrently identify both the number of CaL and the IR. To this end, Eq. (1) will be combined with a celebrated law in inventory and queueing systems. Among the basic tenets, used to size inventories, to identify bottlenecks and control queue length is the law of mean residence time (Cobham, 1954; Morse, 1958) or Little's law in queuing theory (Little, 1961, 2011; Jewell, 1967; Eilon, 1969; Stidham, 1972, 1974; Whitt, 1991, Kim and Whitt, 2013). It posits that at steady state and under no losses, the expected number of items in the system U (a time average) and the expected mean time in the system (MTS) τ, a population average of all items in the system, are linearly related viawhere a is the inflow rate (time average of the inflow). Law 2 has been used to identify any of the three variables if the other two are known, e.g., in lead time reduction (de Treville et al., 2014) and, extensively, for sizing manufacturing and processing equipment. Under real conditions and uncertainty in U, more accurate estimates of U can be obtained via Eq. (2), if τ and a are monitored (the customer arrival rate, or rate of order placement) rather than by direct monitoring of U (Glynn and Whitt, 1989). This may be attributed to MTS being a scaled variable.

In the problem addressed with U representing CaL and a representing the IR, the MTS is the mean residence time of all carriers at large, including those in the early exodus (symptomatic and asymptomatic ones, carriers eventually self-cured or eventually deceased) and those entering HCUs. A higher MTS leads to higher a, since SARS-CoV-2 is highly contagious and the longer carriers carry the virus, the higher the IR will be. The MTS depends on the exit and early leave profiles, that is, on the evolving non-stationary ED and on the residence probability in CaL, which is varying with time. If τ is known, then Law 2 could give any of the variables U and a, if the other variable were known. To this end, Law 2 should be proven valid under early exodus, or in general, under losses of inventory. As shown in Appendix C in the SM, being inflow based, Law 2 remains valid under random losses of inventory and thus it carries on as a useful means in our array.

In general, Law 2 stand alone, entailing at least two unknown and unobservable variables and one possibly known variable (τ) cannot provide concurrent estimates for U and a. On the other hand, Eq. (1) provides an asymptotic relation between two unobservable variables, that is, CaL, U and IR, a, in terms of the variables η and E. Eq. (1) is an independent relation, not antithetic and not a pleonasm to Law 2. As with Law 2, a relation between two unknown and unobservable random variables is better than no relations at all. The exit / entry rate to HCUs, E, is recorded (known) and the mean residence time in CaL of entries to HCUs can be readily monitored. More specifically, the variables η and θ, are population averages defined analytically in the conventional statistical way at each point in time. They are distinct and different than the MTS (Tsiliyannis, 2016). The mean time of carrying the virus of patients in E (average sojourn time in CaL of patients entering HCUs) i.e., θ, can be monitored from representative samples of E (samples of patients arriving for treatment in HCUs). The mean time of alive CaL, η, could also be monitored, when SARS-CoV-2 tests are conducted, tracing the Ig-G and Ig-M antibodies, among carriers who are testing positive, from the various vintages in the sample of U. Both η and θ are scaled variables (as is the MTS) and thus they can be monitored from small or decentralized samples with small monitoring error, compared to obtaining direct estimates of CaL and IR. Being a scaled variable, η is about the same for symptomatic and asymptomatic carriers in CaL and thus it can be based only on symptomatic carriers. Thus, Eq. (1) exploits the fact that E is known, θ, can be readily and reliably monitored and η can be monitored from symptomatic carriers. Then Eq. (1) includes two unknowns: U and a: if one of them is known, the other one can be found, i.e., one degree of freedom. If both are unknown, then combination of Eq. (1) with eq. T1, Table 1 and with Law 2 gives the independent relations T4-T8 (Table 2 ) which enable identification. This is the cornerstone of the identification method, expounded in Appendix D (SM) and codified in the sequel for various data sets. A data-driven, machine–learning procedure, rather than based on detailed modeling, the identification method does not rely on any model in order to track CaL and IR. Eqs. (1) and (2) are useful in obtaining conditions for social immunity (Table 3 , Proof in Appendix E, SM).

Table 2

Asymptotic expressions used in the identification method (Eqs. (T4)–(T8)).

U−θE=(η−1)(a−E) (T4)

τ−θε=(η−1)(1−ε) (T5)

E=εa (T6)

τ(1−x)=x(1−ε) (T7)

ε((η−θ)x+(1−η+θ))=ηx+(1−η) (T8)

Table 3

Asymptotic conditions for social immunity.

Condition for social immunity:

δ(1−ε)≈0 and ε≈0, or (T9)

δ(θ−τ)+τ−η+1θ−η+1≈0 (T10)

Or

δ≈0 and θ>>τ (T11).

Implementation of the identification method

A 5-step identification procedure

Unknown: U. Known, recorded E monitored . END. ◊

In each successive time interval, t acquire system data: :

Use Eq. (T5) in Table 2 to obtain, the dimensionless rate (fraction with respect to IR) of patients to HCUs, ε.

Use Eq. (T6), Table 2, to obtain the IR, a.

Use Eq. (T4), Table 2, to obtain CaL, U.

Use Eq. (T3), in Table 1, to update the propagation factor k.

A similar procedure is followed if alternative data sets of scaled or dimensionless variables are acquired: In particular, with the data quadruple, , (or or ) use eq. T6 to obtain the IR, a; then go to step 4. With data set: use Eq. (T8) to obtain the dimensionless rate of new carriers to HCUs, ε; then go to step 3.

Validation of robust identification of CaL and IR

USA

Fig. 7a depicts the robust performance of the identification method for the USA (actual incidence rate and model realization as depicted in Fig. 6a) by monitoring the parameters and applying the procedure in the previous section in real time. Zero-mean error in the monitored values of , uniformly and independently distributed within ±10% of the actual values is included. The multiple pattern, resurging phases in Fig. 6a are reflected in the profiles of CaL and IR estimates. The method returns satisfactory values of CaL and IR in real time (MAPE∼15%) tightly following the multiple waves of the epidemic with a time lag of about one week. The time-lag is due to the fact that an asymptotic method is used to identify the dynamic inventory system.

Fig. 7

COVID-19 Real-time identification of the number of SARS-CoV-2 carriers at large (CaL) and of the infection rate (IR). Top (a). USA: Under the presence of zero-mean error in the monitored values of , uniformly and independently distributed, within ±10% of the actual values. The CaL estimate, U’ and the IR estimate, a’ follow the model outputs, U and a with MAPE ∼15%. (The model outputs, U and a correspond to the realization in Fig. 6a, which matches the actual incidence profile, Fig. 6b). Bottom (b). Italy: the estimates U’ and a’ follow the model outputs U and a (Fig. 6b) with MAPE 10% and 17%, respectively.

Italy

Fig. 7b presents the results of identification for Italy (actual incidence rate and model realization depicted in Fig. 6b under the presence of error as above (zero mean error in the monitored values of uniformly and independently distributed within ±10% of the true values). The method returns satisfactory estimates of CaL and IR in all phases of COVID-19 evolution with 7-day time-lag. In periods of slow dynamics, the MAPE is less than 7% and overall less than 17%.

China

Fig. S5 in the Supplementary material (SM) in the Web depicts the performance of the identification method for China (actual incidence rate and model realization as in Fig. 5 by monitoring the parameters in real time and applying the procedure in the previous section. It shows that, in the explosive period, the procedure returns satisfactory values of CaL and IR with a delay of about one week, despite intense and random variation.

Fig. S6 in SM presents the identification of CaL and IR for Greece (actual incidence rate and model realization depicted in Figs. 2 and 3 under zero mean error in the monitored parameters, , uniformly and independently distributed within ±10% of the actual values. It manifests tight closure (MAPE ∼20%) in the identified values of CaL and IR.

Managerial insight and practical use

Regarding industrial practice, the beneficial parts of the present work include: (a) The time-varying, stochastic inventory model Eqs. (T1)–(T3) and Eqs. (T12)–(T16) in Appendix A in SM which debars from up-to-date mechanistic epidemic models.

There are no difficulties in applying the model: the analytic convolution expression T12 in SM includes 21 terms and the convolution expression T15 in SM includes 26 terms (T∼17 days and ν∼21 days) readily programmable in widely available software, e.g., EXCEL. of fixed structure (b) the real-time, data-driven identification method giving the number of CaL and the IR. In regard to modeling, decision makers may gain insight by utilizing the model (Eqs. (T1)–(T3) and Eqs. (T12)–(T16)) to obtain alternative realizations of epidemic evolution via Markov-chain Monte Carlo simulation. The model can be tuned to match multi-pattern national incidence data, under any, or no control measures. An advantage with respect to existing work is that distinct phases of evolution and multiple pattern profiles, featuring exponential, hyper-exponential, or polynomial growth, oscillations and resurgences may be followed. This is accomplished by: (a) setting the range of random variability of the early loss probability, x, e.g., uniformly in [0.9,1) for COVID-19 (b) selecting the mean ED, g = 1,2,…,ν and the distribution of each fraction g, of the ED, with mean g, as a sequence of random variables (c) selecting stepwise the randomly varying sequence of the propagation factor, k, given the history, , to match actual data up to time t, by the resulting realizations, within a desired MAPE. The system evolves according to Eqs. (T1)–(T3) and Eqs. (T12)–(T16).

A new feature, regarding the practical use of the identification method giving the number of infectious cases and the infection rate, is that policy makers know in advance the minimum number of variables that must be monitored (four). Several options of such quadruples were proposed that included sets of scaled and dimensionless variables, besides the rate of entry to HCUs, which is known. Valuable insight (cognizance of the actual number of CaL and IR) may be gained, regarding the trajectories of CaL and IR without resorting to dynamic models, since the identification method, is based on asymptotic laws and does not require the inventory model, or any underlying model to run and give the trajectories of CaL and IR. Policy action may benefit from combining the model and the identification method, since the latter leads to the reconstruction of the randomly varying sequence of the propagation factor, k, based on the identified trajectories of CaL and IR. The reconstructed propagation factor can be exploited in tuning the stochastic inventory model to match real data, in the presence of structural variability, noise, uncertainty and the myriad of time-varying factors affecting propagation. Mitigating action and policies such as intensity of social distancing and mobility measures may be based on associating the monitored values of with the ranges of the propagation factor, k, and with past CaL and IR levels, control measures, level of immunity (as reflected by vaccination) economic repercussions and public opinion. Model predictive–adaptive control policies could hinge on the above principles.

Conclusion

Epidemic models are fixed structured, compartmental, finite-dimensional, lumped-parameter dynamic models, with constant parameters (kinetic rates). They are inspired by chemical reaction kinetic models that assume homogeneity (Brownian motion in homogeneous mixtures) and linear, or bilinear kinetics. Model parameters are determined via regression on past data of the states, which entail significant uncertainty. Epidemic models have contributed significantly in understanding outbreak, surge and fading out of contagious diseases. However, such a postulate leads to single pattern outbreaks, while future evolution is predetermined starting at any point in time (initial conditions). This is not the case with COVID-19: infinite routings of evolution of the disease, initiating at each point in time are possible. In particular, existing models may not faithfully represent the complex process of the spread of SARS-CoV-2 under actual conditions, which is infinite dimensional, nonlinear with unknown nonlinearities, delayed, distributed, stochastic and time-varying. The present work contributed towards a representative modeling and identification of COVID-19. The main findings and contributions are the following. (a) A ‘black-box’ inventory model is apposite, in which the inventory is paralleled with the number of carriers at large CaL (i.e., the inventory as ‘viewed’ by the SARS-CoV-2) and the infection rate (IR) is paralleled with supply in classical inventory theory, or arrival rate in queueing systems. (b) CaL and IR are key unknowns of the pandemic. (c) CaL includes all infectious carriers conveying the virus (exposed, transmitting, infected, sick, contagious, symptomatic, asymptomatic, testing positive or negative) under any (or no) restrictive/protective measures, confined, quarantined, vaccinated and recurrent. (d) Random losses of inventory, possibly with retrials, represent no longer infected, self-cured (potentially recurrent) and deceased that do not enter HCUs. (e) The outflow from inventory (demand or job completion rate, or order delivery in classical inventory) is the (accurately recorded) entry of patients to HCUs. It can be modelled as a delayed/distributed exit, with randomly varying (stochastic process-like) fractions of the exit distribution and, possibly shifting delay and spread. (f) Such an infinite-dimensional view unleashes the potential of using the powerful OR array of inventory control and identification methods to track and identify the spread of the disease towards efficient use of health-care resources and control. (g) The key monitored variables in the proposed identification method herein are scaled variables, i.e., characteristic mean times: the mean residence time in CaL of alive carriers (a population average –- monitored from samples of symptomatic carriers via the ongoing tests in the general population) and the mean residence time in CaL of carriers entering HCUs, monitored from samples of entries to HCUs. (h) Based on asymptotic laws, the identification method performs satisfactorily in periods of tame evolution, while it follows closely the trend of CaL and IR under explosive, or controlled growth, providing satisfactory estimates with a lag of one week. (i) The method provides robust identification within 30% of actual values, under zero-mean, uniformly distributed error in the monitored parameters.

Limitations of the analytic model may regard successive drastic structural changes of the ED, e.g., drastic changes in the delay and spread of the ED, due to viral mutations, or abrupt change of conditions. Nevertheless, the fundamental inventory model and the identification method, which is based on asymptotic laws, are hardly affected. Extension of the model in series to include patients treated in HCUs and ICUs (CaL→ HCUs → ICUs, i.e., three successive inventory modules) is a next step. It may lead to improved prediction and reduction of fatality rates related to local conditions, e.g., atmospheric pollution, or ICU treatment protocols. A challenging avenue of future work is to better track explosive outbursts of contagious diseases. To this purpose, work may focus on transient laws of inventory, which could enhance the speed and timeliness of the identification (i.e., with no lag) while restraining overshoot at high peaks, during the explosive and decline periods, under the same inventory perspective of COVID-19. The intriguing question then will be whether robustness under monitoring error is maintained, in the presence of dynamic accrual of estimation error. Another emerging issue is how local inventory models can be combined to yield an overall stochastic inventory model of the same structure. Local models may be differentiated, e.g., if tuned to different areas to represent the respective social and environmental profiles of recurrent surges and controls, or if they correspond to SARS-CoV-2 variants. The quest then is whether asymptotic law 1, featuring the characteristic times of the corresponding areas, could be unified to a single law of the same simple form, engulfing both areas and whether the scaled variables (characteristic times) in the combined law can be taken as the population-weighted mean of the areas. The same questions regard Law 2. This could enable implementation of the same real-time identification method for the overall system. Finally, the results of the present work may provide guidance for advancing structured mechanistic compartmental models to reflect heterogeneous, spatial, social, environmental and genetic transmission factors that give rise to multiple pattern evolution. Possible venues to this end are via inferential monitoring of states, based on the characteristic times, or via inferential update of temporally-varying compartmental rates, using the propagation factor.

ncited reference

Yuan et al., 2020b.