key: cord-1050510-vtrl82k2
authors: Tsiligianni, Christiana; Tsiligiannis, Aristeides; Tsiliyannis, Christos
title: A stochastic inventory model of COVID-19 and robust, real-time identification of carriers at large and infection rate via asymptotic laws
date: 2022-01-08
journal: Eur J Oper Res
DOI: 10.1016/j.ejor.2021.12.037
sha: d33ad7353e63c1b3d53a501cb34f355b33f2bf5a
doc_id: 1050510
cord_uid: vtrl82k2

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 victims. 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 design of policies restricting COVID-19 and encumbrance to HCUs and mitigating economic contraction.

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) . In parallel to the selfsacrificial, 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.

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 t , 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 t , 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 t 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 t is the throughput of the 'inventory system', an additional concern of policy and proactive measures against COVID-19: increasing values of E t , 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 t and IR, a t 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 t , 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.

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 Ronald Ross (1911 Ross ( , 1915 ) the legendary SIR model (Kermack and McKendrick, 1927) and further extensions (Hethcote, 2000; Atizer and Nunn, 2006; Brauer, 2008; Vynnycky and White, 2010) 5 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) (Diekman et al., 1995; Pell et al., 2015; Giordano et al., 2020) from which the basic reproduction ratio, R 0 , is found (van der 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 (Chowel 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 finitedimensional 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, noncompartmental, 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; Chowell et al., 2019) . Such a view considers the infection tree(s) (Pionti et al., 2014) a 'black box" without seeking a detailed branching process model (Jagers, 1975) or an underlying compartmental, or networked epidemic model (Keeling et al., 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 6 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 realtime, based on low uncertainty data.

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 Steclinskhi, 2012; Davila-Pena et al., 2021; ) and

age stratification (Canto et al., 2020) .

 Effects of shifting community trust (or mistrust) to mitigation measures and vaccination, and resulting population behavior and compliance (Yan et al., 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.06 2021) much higher than the national average 7 (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 (2 nd 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 timedelays, 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) . 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 and Rohani, 2001; Korobeinikov, 2006) structural instabilities (van der Driessche and Watmough, 2000) and diffusion (Willis et al., 2020) .

Several improvements have been suggested to the above. The modified SEIR model 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, 9 Kohler et al., 2020; Carli et al., 2020; Sereno et al., 2021) were used to derive step-wise optimal distancing policies. However, there are no guaranteed stability margins (stabilizing effect) of fixedstructure, 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) . 

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 (Tsiliyannis, -2019 . The underlying idea is to replace the continuous time with discrete points in time, at which the system state is observed (Schwartz 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; 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, 2015) . 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 Nomenclature a t infection rate or rate or new carriers contracting the disease in period t, a = steady state level CaL = U t , carriers at large at the end of time period t. See also U t . COVID-19 the pandemic due to SARS-CoV-2 ED exit distribution from the group of CaL to HCUs or directly to ICUs E t 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 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 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 Ω t , the early exodus, is a nonnegative random variable, not identically distributed with the exit E t . 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 

. 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 t .) 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 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.

The random early exodus from CaL, t Recorded daily (known).

Ω t ≥ 0, early exodus rate, or early leave rate of no longer hosts or infectious, of self-cured, or deceased carriers, who did not enter HCUs/ICUs: random abandonments of queue or losses of inventory (unknown and unobservable). 

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 i,t , i=1,2,…,v, 0 ≤ g i,t ≤ 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 t , 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 t 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.

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 t-1 (end of period t-1) to level U t (end of period t) by the net amount of new carriers,

where a t is the IR during period t, E t is the exit from CaL in the beginning of period t entering HCUs and Ω t 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 t , (k t 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 timevarying linear system with a time-varying pole at 1+k t , (Kamen, 1987; 2010 

Restrictive measures or vaccination lower k t 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 t (eq. T3) is a novelty in modeling actual epidemics, compared to SEIRD models, for it allows stepwise fine-tuning of k t at the present time (t) given the history (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, Ω t , depends only on time t and is independent from the outflow to HCUs, E t (3) the IR evolves according to eq. T3, Table 1 .

Several cases of COVID-19 national incidence evolution are presented, featuring distinct profiles.

It is assumed that the entry time to HCUs, E t , 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 t and the exit to HCUs, E t , in dynamic time, in terms of the IR (expressions T12 and T15 in Table A 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 t = 0.16 U t-1 ( fig. 4a ) (k t =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 (Marchmid April 2020) in which the propagation factor k t utilized in the model was sharply lower (k t~0 .1) and even lower (k t~0 .065) (iii) the decline phase (April -mid May 2020) with k t ~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 t was necessary (k t ~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 t rose to ~0.14). (v) A tame period, up to day 330 (k t ~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 t varied randomly around k t ~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). 

~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 t ~0.26. As depicted in fig. 4c , the ranges of k t 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. 6b depicts a realization of 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 t 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 t =0.25. The multiple wave profile, faithfully reproduced herein, may not be modeled by constant parameter models (Calafiore et al., 2020a) .

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 t , 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 t :

where 21 θ t is the mean time (population average) spent as CaL of the patients in E t , that is, of patients who are exiting CaL and entering HCUs in the beginning of time period t (θ t =mean duration of sojourn in the group of CaL of all members of E t ) η t , is the mean time (population average) of alive CaL (η t =mean duration as CaL of all alive members of U t 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; Jewell, 1967; Elion, 1969; Stidham, 1972 Stidham, , 1974 Whitt, 1991, Kim and Whitt, 2013) x E , , ,   use eq. T8 to obtain the dimensionless rate of new carriers to HCUs, ε t ; then go to step 3. 

Condition for social immunity: is due to the fact that an asymptotic method is used to identify the dynamic inventory system. lag. In periods of slow dynamics, the MAPE is less than 7% and overall less than 17%. 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. 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 modelling, decision makers may gain insight by utilizing the model (eqs T2-T8) 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 t , e.g., uniformly in [0.9,1) for COVID-19 (b) selecting the mean ED, g i , i=1,2,…,ν and the distribution of each fraction g i,t , of the ED, with mean g i , as a sequence of random variables (c) selecting stepwise the randomly varying sequence of the propagation factor, k t , given the history, 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 t , 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 t , 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.

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, 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 29 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.

COVID-19 with uncertain phases: Estimation issues with an illustration for Argentina

Covid-19: from model prediction to model predictive control

A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis

Databased analysis, modelling and forecasting of the COVID-19 outbreak

On patient flow in hospitals: A data-based queueing-science perspective

On scheduling a multiclass queue with abandonments under general delay costs, Queueing Systems

Infectious diseases in primates: behavior, ecology and evolution. Oxford Series in Ecology and Evolution

A stochastic epidemic model of COVID-19 disease

An attempt at a new analysis of the mortality caused by smallpox and of the advantages of inoculation to prevent it

Robust and data-driven approaches to call centers

Compartmental models in epidemiology

A final size relation for epidemic models of vector transmitted diseases Infectious Disease Modelling 2

A time-varying SIRD model for the COVID-19 contagion in Italy

Modified SIR Model for the COVID-19 Contagion in Italy

Model predictive control to mitigate the COVID-19 outbreak in a multi-region scenario

Mathematical Structure of Epidemic Systems

Modelling transmission and control of the COVID-19 pandemic in Australia

Optimal strategies for social distancing & testing to control covid-19

Priority Assignment in Waiting Line Problems

Estimation of the reproduction number of dengue fever from spatial epidemic data

The Western Africa ebola virus disease epidemic exhibits both global exponential and local polynomial growth rates

A novel sub-epidemic modeling framework for short-term forecasting epidemic waves

Many-server queues with customer abandonment: A survey of diffusion and fluid approximations

A compartmental model that predicts the effect of social distancing and vaccination on controlling COVID-19

Assessment of the influence of features on a classification problem: an application to COVID-19 patients

The GReek study in the Effects of Colchicine in COvid-19 complications prevention (GRECCO-19 study): rationale and study design

Optimal Sourcing and Lead Time Reduction under Evolutionary Demand Risk

The legacy of Kermack and McKendrick

Guaranteed margins for LGQ regulators

Practical scheduling for call center operations

Staffing of time-varying queues to achieve time-stable performance

Real-time forecasting of infectious disease dynamics with a stochastic semi-mechanistic model

Discrete Stochastic Processes, MIT Lecture Notes

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Public health interventions in the face of pandemics: Network structure, social distancing, and heterogeneity

Indirect estimation via L=λW

Stochastic epidemic modelling

The ecology of Ebola virus

Tutorials in Operations Research, A Tutorial on Dynamic Optimization with Applications to Dynamic Rate Queues: Risk and Optimization in an Uncertain World (Publication for tutorials

The Mathematics of Infectious Diseases

A discrete stochastic model of the COVID-19 outbreak: Forecast and control

Epidemic wave dynamics attributable to urban community structure: a theoretical characterization of disease transmission in a large network

Branching processes with biological applications

A Simple Proof of L=AW

The poles and zeros of a linear time-varying system

Fundamentals of linear time-varying systems

Analysis of COVID-19 infection spread in Japan based on stochastic transition mode

The M/M/1 queue with synchronized abandonments

Seasonally forced disease dynamics explored as switching between attractors

Networks and epidemic models

A Contribution to the Mathematical Theory of Epidemics

Advanced workforce e-management for effective customer services

Robust and optimal predictive control of the covid-19 outbreak

Lyapunov functions and global stability for sir and sirs epidemiological models with non-linear transmission

Level of under-reporting including under-diagnosis before the first peak of COVID-19 in various countries: Preliminary Retrospective Results Based on Wavelets and Deterministic Modeling

Asymptomatic carrier state, acute respiratory disease, and pneumonia due to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2): Facts & myths

Epidemic control via stochastic optimal control

Propagation analysis & prediction of COVID-19

A conceptual model for the coronavirus disease (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action

Environmental factors on the SARS epidemic: air temperature, passage of time & multiplicative effect of hospital infection

Little's Law as Viewed on Its 50th Anniversary" (PDF)

Infectious disease models with time-varying parameters and general nonlinear incidence rate

Infectious diseases of humans: dynamics and control

A parametrized nonlinear predictive control strategy for relaxing COVID-19 social distancing measures in Brazil

Queues, inventories, and maintenance: the analysis of operational system with variable demand and supply

SEIR model for COVID-19 dynamics incorporating the environment and social distancing

Epidemic spreading in temporal and adaptive networks with static backbone

Towards understanding socially influenced vaccination decision making: An integrated model of multiple criteria belief modelling and social network analysis

The socioeconomic implications of the coronavirus pandemic (COVID-19): A review

Intervention as a research strategy

Epidemic processes in complex networks

The infection tree of global epidemics

Using phenomenological models for forecasting the 2015 Ebola challenge

Foundations of Stochastic Inventory Theory

Some a priori pathometric equations

An application of the theory of probabilities to the study of a priori pathometry -II

A flexible rolling regression framework for time-varying SIRD models: Application to COVID-19

Use of a Modified SIRD Model to Analyze COVID-19 Data

2021. Model predictive control for optimal social distancing in a type SIR-switched model

How control theory can help us control COVID-19 IEEE Spectrum

Strategies for Ensuring Required Service Level for COVID-19 Herd Immunity in Indian Vaccine Supply Chain Eur

Performance analysis of time-dependent queuing systems: Survey and Classification

Combining probabilistic forecasts of COVID-19 mortality in the United States

Effects of social distancing & isolation on epidemic spreading modeled via dynamical density functional theory

Modeling, state estimation and optimal control for the US COVID-19 outbreak

A fundamental law relating stock and end-of-life flow in cyclic manufacturing

Mean retention time and end-of-life rate identification in cyclic manufacturing

Markov chain modeling and forecasting of product returns in remanufacturing based on stock mean age

Prognosis of product returns for enhanced remanufacturing

Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission

An Introduction to Infectious Disease Modelling

Stationary distribution of a stochastic SIRD epidemic model of Ebola with double saturated incidence rates and vaccination

Insights into the dynamics and control of COVID-19 infection rates

Impact of climate change on human infectious diseases: Empirical evidence and human adaptation

A frailty model for intervention effectiveness against disease transmission when implemented with unobservable heterogeneity

Modified SEIR and AI prediction of the epide-mics trend of COVID-19 in China under public health interventions

Monitoring transmissibility and mortality of COVID-19 in Europe

Association of radiologic findings with mortality of patients infected with 2019 novel coronavirus in Wuhan, China

Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, 2019-2020: A data-driven analysis in the early phase of the outbreak

Foundations of Inventory Management

where 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 t , more accurate estimates of U t can be obtained via eq. 2, if τ and a t are monitored (the customer arrival rate, or rate of order placement) rather than by direct monitoring of U t (Glynn and Whit, 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 t , 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 t and a t . On the other hand, eq. 1 provides an asymptotic relation between two unobservable variables, that is, CaL, U t and IR, a t , in terms of the variables η t , θ t , and E t . 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 t , is recorded (known) and the mean residence time in CaL of entries to HCUs can be readily monitored. More specifically, the variables η t and θ t , are population averages defined analytically in the conventional statistical way at each point in time. They are distinct and different than the MTS (Tsiliyannis, 2020) . The mean time of carrying the virus of patients in E t (average sojourn time in CaL of patients entering HCUs) i.e., θ t , can be monitored from representative samples of E t (samples of patients arriving for treatment in HCUs). The mean time of alive CaL, η t , 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 t . Both η t and θ t 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 24 variable, η t 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 t is known, θ t , can be readily and reliably monitored and η t 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 T9-T13 (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).

1. In each successive time interval, t acquire system data:: 2. Use eq. T5 in Table 2 to obtain, the dimensionless rate (fraction with respect to IR) of patients to HCUs, ε t .3. Use eq. T6, Table 2 , to obtain the IR, a t . 4. Use eq. T4, Table 2 , to obtain CaL, U t . 5. Use eq. T3, in Table 1 , to update the propagation factor k t .

A similar procedure is followed if alternative data sets of scaled or dimensionless variables are acquired: In particular, with the data quadruple, eq. T6 to obtain the IR, a t ; then go to step 4. With data set: