key: cord-0858581-3pawdjok authors: d’Onofrio, Alberto; Manfredi, Piero; Iannelli, Mimmo title: Dynamics of partially mitigated multi-phasic epidemics at low susceptible depletion: Phases of COVID-19 control in Italy as case study date: 2021-07-21 journal: Math Biosci DOI: 10.1016/j.mbs.2021.108671 sha: a90dc793887bcec3e26800f29efc25be5346c291 doc_id: 858581 cord_uid: 3pawdjok To mitigate the harmful effects of the COVID-19 pandemic, world countries have resorted - though with different timing and intensities - to a range of interventions. These interventions and their relaxation have shaped the epidemic into a multi-phase form, namely an early invasion phase often followed by a lockdown phase, whose unlocking triggered a second epidemic wave, and so on. In this article, we provide a kinematic description of an epidemic whose time course is subdivided by mitigation interventions into a sequence of phases, on the assumption that interventions are effective enough to prevent the susceptible proportion to largely depart from 100% (or from any other relevant level). By applying this hypothesis to a general SIR epidemic model with age-since-infection and piece-wise constant contact and recovery rates, we supply a unified treatment of this multi-phase epidemic showing how the different phases unfold over time. Subsequently, by exploiting a wide class of infectiousness and recovery kernels allowing reducibility (either to ordinary or delayed differential equations), we investigate in depth a low-dimensional case allowing a non-trivial full analytical treatment also of the transient dynamics connecting the different phases of the epidemic. Finally, we illustrate our theoretical results by a fit to the overall Italian COVID-19 epidemic since March 2020 till February 2021 i.e., before the mass vaccination campaign. This show the abilities of the proposed model in effectively describing the entire course of an observed multi-phasic epidemic with a minimal set of data and parameters, and in providing useful insight on a number of aspects including e.g., the inertial phenomena surrounding the switch between different phases. Since its detection in the city of Wuhan, China, by end December 2019, the novel SARS-CoV-2 coronavirus has rapidly spread worldwide with a dramatic direct disease burden (WHO situation reports 1-209 and Weekly epidemiological updates, [1] ), and a dramatic impact on the economy, health and the society 5 as a whole. During the first COVID-19 pandemic wave, an increasing number of countries worldwide have opted, following China, for massive social distancing measures -what is now universally identified as lockdown, integrated with a range of supplementary measures such as tracing and isolation of confirmed cases. The aim of lockdown is essentially that of abruptly bringing the basic 10 reproduction number of the infection, R 0 , below threshold, therefore stopping sustained transmission. If this actually occurs, and the measures are not prematurely interrupted, the epidemic will be brought on a suppression path (Ferguson et al. [2] ). The China initial experience, especially the one of the Hubei province, has shown that a combination of severity and appropriate duration 15 of intervention can achieve high degrees of suppression even for epidemics that have reached a non-negligible scale. A somewhat different story occurred in Western Europe, where in many cases the lockdown was declared with substantial delays, with the aggravating circumstance that its implementation was done through several steps instead than abruptly, therefore yielding much larger epi- 20 demics than in China (see [1] ). For example in Italy, whose Northern regions (especially Lombardia) experienced one of the most dramatic local epidemic worldwide, full lockdown was achieved by a long number of subsequent decrees of the Prime Minister covering a span of more than twenty days. These delays not only caused the epidemic to overwhelm available public health resources, 25 with huge human costs, but also made eventually unavoidable to unlock long before the achievement of adequate suppression levels, due to the dramatically growing economic loss. Commom to several European countries, the subsequent unlocking phase have shown prolonged honey-moon epochs of low epidemic activity, possibly resulting from a combination of reasons such as e.g., an enduring were partly successful i.e., the epidemic curve was just brought to stationarity with the current reproduction index R t in the region of one at the national level, until the new variants of concern initiated their replacement of the original ones [1] , giving rise to a third wave [3] , and forcing new restrictions, including regional lockdowns. A number of key questions were continuously debated by the 40 public during each pandemic wave about the effects of the enacted interventions such as e.g., (i) when the effects of social distancing measures will be detectable in the data, or when the epidemic peak of that particular phase will be achieved, (ii) which the further epidemic growth will be, in terms of hospitalizations, and deaths, after the measures implementation, etc. 45 Previous questions have represented major concerns throughout the COVID-19 pandemic and were responsible of widespread anxiety in the public opinion. They therefore call for general answers, even more so in the current moment when a third wave has started and it is clear that the benefits of vaccinationat current rates of administration -will take time to manifest. And, obviously, 50 in view of future pandemics as they are not anymore classifiable as rare events. The COVID-19 pandemic has called for an unprecedented interest of modelers for epidemic dynamics both per-se, to investigate the direct epidemiological consequences of the pandemic and related mitigation measures, and also to investigate its broader implications on society. This has led to an endless number 55 of papers and preprints of which we quote here only a few focusing on COVID-19 epidemiological modeling [2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] . Most of the cited investigations have resorted to simulation. However, even under fairly general conditions, many of these phenomena can be investigated in a fully analytical manner, therefore providing a range of further 60 robust insights about epidemic control. Therefore, motivated by the main feature of the COVID-19 epidemic i.e., being characterised by a sequence of phases driven by the intensity and duration of the undertaken mitigation measures, in this paper we attempt to analyse, in an analytical manner, the dynamics of a generic multi-phasic epidemic, much in the 65 same way as we observed for COVID-19 prior the current mass vaccination campaign started to deploy its main effects. These phases include an early invasion phase, the subsequent mitigation interventions (e.g., generalised lockdown), the unlocking/releasing phase, the second wave phase, the further mitigating interventions etc. To make the problem tractable, we assumed that the interventions 70 considered are sufficiently long-lasting and effective so to prevent that the susceptible proportion of the population excessively departed from its initial value (presumably around 100% for COVID- 19) or from any other reference value. This hypothesis can be always made valid by considering suitable observation windows. Besides the rapid control achieved in China already by February 2020 of confirmed infections was in the range of 5%. This is strongly suggestive of the fact that the potential impact on the course of the epidemics, due to the depletion of the susceptible population, would always be secondary with respect to the impact due to the undertaken mitigation measures. More in general we 85 expect this to be largely true in all countries and subnational settings where social distancing arrived at a sufficiently earlier stage of the epidemic so to halt it before it depleted the susceptible proportion. Note additionally that this approach is valid for any combination of epidemic phases where interventions are able to approximately "freeze" the susceptible fraction around some given level. In these circumstances, the analysis of a multi-phasic epidemic yields to a linear problem which can be analysed by the tools of age-structured populations mathematics ( [23, 24, 25] and references therein). Our principal aim here will be to provide a unified framework to multi-phase epidemics while keeping our results as more analytical as possible. Our approach will be kinematic in the 95 sense that will look at how the different epidemic phases unfold over time reshaping the epidemic, without minding at whether, for example, a certain value of the current reproduction number R t achieved during e.g., a lockdown phase, was the mere outcome of government-enforced restrictions or was also reflecting the beneficial (or harmful) effect of virtuous (resistant) behavior spread in 100 the population, as would instead be the approach of behavioral epidemiologists ( [26, 27] ). Consistently, we will simply consider a SIR epidemic model including class age i.e., age since infection, in the infective compartment and piece-wise constant transmission and recovery rates over time reflecting different intervention epochs. The reason why we use such a simple model (e.g., without including 105 more compartments as symptomatic cases, hospitalization, deaths etc as a faithful description of COVID-19 would require), is that this simple model provides the backbone description of a general multi-phasic epidemic. This backbone can be easily refined by adding further compartments representing delayed events following infection [28] . Despite its simplicity, our SIR model is general in that we consider both a latency delay and general forms for the infectiousness and recovery processes. By the proposed model, we first provide a general characterization of a multi-phasic epidemic. Second, by adopting a wide class of generation time and recovery kernels allowing reducibility to ordinary or delayed differential equations ( [29] ), 115 we were able to provide a detailed investigation of a non-trivial low-dimensional case allowing a full analytical treatment of the transient dynamics of a multiphase epidemic. Finally, we complemented our theoretical results by a number of illustrations, using realistic parameter constellations drawn from available estimates of COVID parameters. In particular, we fitted the model to the overall 120 Italian COVID-19 epidemic since March 2020 until late February 2021. This illustrates the abilities of the proposed model in effectively explaining the entire course of the epidemic, including the effects of interventions, using a minimal set of data and parameters, and in providing useful insight on a number of aspects including e.g., the demographic inertia occurring during the switch between 125 different epidemic phases. The article is organised as follows. Section 2 develops our model for a multi-phasic epidemic. Section 3 reports the main theoretical results. Section 4 presents the proposed class of reducible kernels and reports a number of sub-cases and illustrations. Section 5 reports the application to Italian data. Concluding remarks follow. The model we present here is a SIR epidemic model structured by age since infection or class age. The evolution of the epidemics over time is subdivided into phases, each one characterized by a certain intervention intensity. 135 Our model is based on the state variable where x denotes class age, representing the (improper) age-density of the number of infected individuals. In particular, x + denotes the maximal age since infection and may be finite or infinite. Our main assumption is that interventions are 140 early and effective enough to prevent a large depletion of the susceptible population during the overall epidemic course i.e., the susceptible fraction does not significantly departs from 1. Note that this hypothesis can actually be applied at any level S of the susceptible proportion, provided subsequent interventions are effective enough to "freeze" the susceptible population therafter. In this manner, the model equations read Where • γ(x, t) = the removal rate of infected individuals aged x at time t, that is, the overall rate at which infected individuals aged x are removed by any cause as, e.g., recovery or death (during an uncontrolled epidemic) 150 or by screening, tracing, isolation, hospitalization etc (in the presence of interventions). Note that, thanks to the stated hypotheses, the incidence of new infections depends linearly on the infected density Y (x, t). Both rates β(t, x), γ(x, t) are taken as functions of time as well, to reflect the possibility of interventions (e.g., social distancing, testing, isolating), or the removal of such measures. To mimic the various phases of the epidemic experienced by most European countries we 160 will represent these rates as piecewise function of time. Just for illustration, to represent a 3-phase epidemic, as it might be the case for an epidemic with an initial outbreak phase, a lockdown phase, and a subsequent release phase, we will take and where I [a,b) (t) represents the characteristic function of the interval [a, b), while t L and t R , 0 < t L < t R , are the lockdown and unlocking times respectively. Of course, it is possible to consider any number of epidemic phases. Specification (2) models interventions aiming to affect the age-density of removal in the broad sense defined above. In (3), factor c(t) denotes the number of adequate contacts 170 per person and per unit time, while β 0 (x) tunes the (average) intrinsic infectiousness of an infected individual aged x (i.e, the probability that an infected individual aged x infects a susceptible during an adequate contact), possibly related to her/his viral load. Our formulation therefore assumes that intervention strategies may act either on the number of adequate contacts (as is the case for 175 social distancing) or the scale of transmission (e.g., by using masks), but not on the shape of the distribution of infectiousness by age since infection. Note moreover that the previous form describes abrupt switches between phases, reflecting an ideal world where compliant agents suddenly adjust their behavior upon governmental intervention. However, extension to smooth transitions over time is 180 possible. Overall, the proposed model represents a non-autonomous, infinitedimensional, dynamical system with piecewise constant switching parameters (see for instance [30] , [31] ). Stepwise solution of the model Within the previous framework, the model can be solved in a stepwise man-185 ner, phase by phase. For the sake of simplicity we limit the discussion to three "representative" phases, namely an invasion (or outbreak) phase, a lockdown phase and a release one (respectively labeled by i = 1, i = 2, i = 3), but the analysis can continue to include successive phases. We first note that under our hypotheses each phase of the epidemics (i = 190 1, 2, 3) is characterized by a phase-specific reproduction number given by is the phase-specific survival-to-removal probabilities, defining the probability that a newly infected individual is still infective by age x conditional on the interventions undertaken during any given phase. Concerning the reproduc-195 tion numbers, we will make the assumptions (consistent with the experience of countries which resorted to lockdown) R 1 0 > 1, R 2 0 < 1, R 3 0 > 1 (but other hypotheses can obviously be made). For each phase (i = 1, 2, 3) we have the Lotka-Von Foerster problem where the initial age distribution of the next phase is the terminal distribution , and the solution to (1) reads By integration along characteristic lines (see for instance [24] ), each system (4) can be transformed into an integral equation in the phase-specific incidence Namely we obtain the phase-specific renewal equation where is a normalised kernel representing the generation time density, and On the further assumption that the temporal duration of each phase is sufficient to allow the emergence of long-term behavior, the ultimate behavior of the incidence during the i-th phase is therefore exponential at a Lotka's intrinsic rate α * i and scale constant U * i , also termed the stable 235 equivalent in demographic jargon (see [23] ). The requirement that intervention phases are sufficiently long is by no means unrealistic: during invasion phases convergence to the exponential epidemic path is usually very fast, while for e.g., a lockdown phase aiming to bring reproduction below threshold, it would purely be irrational to stop the intervention before the epidemic has entered its 240 exponentially declining path. The asymptotic behavior of incidence, as stated in (10), in turn implies the emergence of the stable age distribution (i.e., a stable age profile) of the infected population: where, if x + < +∞, in the denominator we intend that the function Γ i is 245 extended as 0 for x > x + (indeed, the stable distribution arises if the function Γ i (x) is Laplace transformable at α * i ). Actually, based on (10), the solution of (4) [see (6) ] will approach the asymptotic form . Precisely, we have [see (6) ] where the convergence occurs pointwise in x. We briefly comment about the components of the stable equivalent U * 1 of COVID-19 dynamics during the exponential phase. The denominator in (9) reads is the average age at which infected individual produce their secondary cases, at the stable age distribution prevailing during invasion. Therefore, the stable equivalent can be represented as Notably, at the numerator of the expression we have the product where represents, still resorting to demographic jargon, Fisher's reproductive value of an infected aged x ( [23] ). Therefore, the numerator in (13) , which scales the magnitude of the epidemic, also has a noteworthy demographic interpretation, i.e., it represents the total reproductive value i.e., the infection potential embed-265 ded in the initial age distribution of infected individuals. This section summarises the main model results in distinct subsections focusing on the distinct phases (i = 1, 2, 3) of the epidemic using equations (5) . During each phase, the infected distribution evolves toward its asymptotic stable 270 age distribution and provides the new initial age distribution for the subsequent phase. Thus each phase is ultimately represented by the dominant real root α * i of the characteristic equation (7) and by the stable equivalent U * i . We discuss in greater detail the invasion phase (i=1), where no intervention is in place, and the subsequent social distancing phase (i=2). Finally we briefly 275 extend the results to the release phase (i=3). We first consider the initial phase of the epidemic, that is the solution of (1) in the interval [0, t L ] where no control measures are in place yet. This is governed by the solution of (4) for i = 1 with R 1 0 > 1, so that the dominant real 280 root α * 1 is strictly positive. Thus, whatever be the initial distribution Y 0 (x), the infected population will reach its stable age distribution [(11)] Note that, since α * 1 > 0, the function Γ 1 (x), being bounded, is Laplace transformable at α * 1 . As known from classical mathematical demography, the stable age distri-285 bution (SAD), which here represents the relative distribution of the infected J o u r n a l P r e -p r o o f Journal Pre-proof by (class) age, is time invariant due to the fact that the absolute numbers of infective individuals in each age group grow at the same rate α * 1 over time. In particular, the normalising term of the SAD, given by the reciprocal of the Laplace transform of the survival to removal function Γ 1 (α * 1 ) has a clear epi-290 demiological meaning. It indeed represents the long-term per-capita birth rate of new infective individuals -generically defined as the ratio between the rate of new births, i.e., the incidence U i (t), and the overall size of the infected population I(t) -in the asymptotic regime of exponential growth of incidence defined in (10). In case the age profile of the initial datum corresponds to the SAD (11), where I(0) is the number of infected individuals at the beginning of the epidemic, then From (17), the corresponding solution to (4) reads i.e., the age-distribution is stable at all times. The latter developments state that if the infective population lies on its SAD already at the initial phase of the epidemic, then the dynamics of incidence is simply given by the exponential propagation at the intrinsic rate α * 1 of the total initial infective population I(0) 305 through its stable (per-capita) birth rate 1/ Γ 1 (α * 1 ). Concerning the behaviour of the total number of infected during the outbreak phase (6) and (12) we see that it splits into two terms Thus, the first term is a transient one accounting for the initial group of infected people, while the second term accounts for the epidemic dynamics resulting 11 J o u r n a l P r e -p r o o f from new cases generated for t > 0 which for large times converges to its stable form. Precisely, the initial group I 0 1 (t) decays at a rate related to the removal probability Γ 1 (x). In fact, if x + < +∞ we simply have Note that the ultimate behavior of the total number of infected during the first phase differs from incidence only in the scale constant, which is now given by U * 1 Γ 1 (α * 1 ). The latter represents the stable equivalent of the total infected population, which differs from that of incidence by a multiplication for the reciprocal of the asymptotic birth rate of new infectives, Γ 1 (α * 1 ). Finally, as regards the lapse of time necessary for observing the emergence of the asymptotic behavior of the solution, we can estimate the approximate Let us now consider the implications of a sufficiently long lasting mitigation phase starting at the time t L , when the reproduction number is abruptly set at R 2 0 < 1. Since in this case we have α * 2 < 0, if social distancing is maintained for 330 a sufficiently long time, the epidemic will eventually set on an exponentially declining path that we will term a suppression path. The suppression path would eventually bring to epidemic extinction if the lockdown phase would continue rather than being halted by unlocking. However, before setting on the suppression path, the lockdown dynamics will be somewhat articulated as it will have to 335 connect the pattern of fast exponential epidemic growth inherited from the first phase with its regime phase of exponential decline. This suggests a transient phase resulting from the abrupt switch between dynamic regimes, which will be dominated by the inertia inherited from the fast growth of the first phase, and might result in an epidemic peak, after which epidemic decline will start. During this phase the age distribution of infective individuals will continuously adjust, and will eventually converge to its new stable limiting form Y * 2 (x) that will promote the exponential decline of the suppression phase. In what follows, we discuss the second-phase dynamics by distinguishing between the limit case of epidemic annihilation (R 2 0 = 0) from the standard case of (eventual) elimi-345 nation 0 < R 2 0 < 1. J o u r n a l P r e -p r o o f In this case, setting R 2 0 = 0, we have U 2 (t) ≡ 0 and Thus, for the prevalence of infected we have only the transient component I 0 that either is identically 0 for t > x + (if x + < +∞) or goes extinct at any rate ω such that 0 < ω < −σ Γ2 [see the analogous (19)]. pointwise in x. Note that, in the case x + < +∞, we certainly have but, if x + = +∞, since α * 2 < 0, the limit on the right hand side of (22) is integrable if and only if α * 2 > σ Γ2 . In this case (23) holds true [compare with (20) ]. However if α * 2 < σ Γ2 we instead have for any ω such that 0 < ω < −α * 2 < −σ Γ2 , i.e., I + 2 (t) decays at the rate −ω. Much information on epidemic trend during the lockdown phase is provided by the two parameters governing the resulting long-term declining exponential path, namely the intrinsic rate α * 2 and the stable equivalent, where the latter tunes the scale of the epidemic curve over the suppression path. We assume that at the lockdown time t L , the epidemic has reached its stable form (12), so 365 that the initial datum at t L is implying a total number of infected individuals Here U * 1 is in general given by (13) while, if the initial population of infectious individuals at the beginning of the epidemic was distributed according to the stable distribution of the first phase, then [see (17) ] . . The stable equivalent of the incidence function for the lockdown phase is [see (13) and (14)] In the particular case of Γ 1 (x) ≡ Γ 2 (x), i.e. if the lockdown measures only involve the lowering of the contact rate from its free-epidemic level c 1 to c 2 < c 1 , 375 the latter quantity reads The previous quantity encapsulates the inertial effects arising from the transition between two regimes of stable growth, namely an initial regime of exponential epidemic growth and the subsequent suppression regime characterised by asymptotic exponential decline. For simplicity we assume that the release of social distancing measures starts when the stable age distribution of the suppression phase was fully established so that we have to consider the incidence equation The analysis proceeds in the same way as the first two phases ending in a longterm regime of exponential growth (R 3 0 > 1) or decay (R 3 0 < 1). In particular the corresponding stable equivalent will embed the inertial effects due to the ageing experienced by the age distribution of infection inherited from the suppression dynamics of the lockdown 390 phase, on the unlocking dynamics. 14 J o u r n a l P r e -p r o o f To better appreciate the general discussion above, it is worth to consider specific parametrizations allowing a more detailed description of the epidemic phases as well as comparisons with data. Moreover, by a parametrized descrip-395 tion of the process, we can exploit Erlang kernels allowing the reduction of infinite-dimension systems to systems of ordinary differential equations (ODEs) or to mixed systems of delayed and ordinary differential equations (DDE-ODE since now on). This is a classical reduction procedure for integral equations of Volterra type, known as the linear chain trick ( [29] ), that can be directly 400 applied to our systems (4). We recall that the compartmental approximation through sequences of ODEs with constant transition rates represent the easiest way to reliably simulate the PDE systems studied here. Our general model depends on two key age-specific epidemiological functions 405 i.e., the (phase-specific) survival-to-infection function Γ i (x), which tunes the age distribution of removal, and the infection reproduction kernel K i (x), which is the normalised product of the infectivity kernel β 0 (x) and Γ i (x), tuning the generation time (age-) distribution. A wide class of reducible problems appears if both distributions are translated Erlang densities, that is Gamma densities 410 with integer index, i.e., i) the infectivity kernel β 0 (x) is a (non-proper) translated Erlang density of order n and rate ϕ: ii) the removal rate γ i (x) is taken as corresponding to the survival function of a translated Erlang density of 415 order m and rate γ i , that is: Under previous specifications, the phase-specific reproduction numbers are given by Formulations (24) and (25) assume a latently infective phase having fixed duration τ, (τ ≥ 0) during which infected individuals are not yet infective and cannot be removed. In relation to COVID-19 data, we remark that although the Gamma kernels best fitting observed data of serial intervals (taken as surrogates 425 of generation times) were obviously characterised by non-integer indices ( [2, 32] ), therefore implying non-reducible kernels, nonetheless integer-indexed Gamma (i.e., Erlang-type) distributions can usefully bound observed distributions. As we now show by means of examples, if duration τ is strictly positive, the linear chain trick will lead to a mixed DDE-ODE system which, in the limit 430 case τ = 0, will collapse into a pure ODE system. To illustrate how the chosen kernels allow reduction (for each phase), we discuss the case n = 2, m = 1 yielding J o u r n a l P r e -p r o o f By considering the basic variables where y i (x, t) is the scaled variable we obtain (see details in Appendix B) the following DDE-ODE system having as state variables the epidemic incidence U i (t) and the auxiliary variable J i (t): Note that the variable J i (t), that appears upon differentiation of U i (t), contributes to determine the behaviour of U i (t) that in turn contains all the rele- Once U i (t) is known, using Y i (x, t) from formula (6), other relevant epidemiological quantities can be computed such as prevalence of infected individuals: For these variables we can also state additional equations that can be added to system (30) . Indeed, proceeding as before we obtain J o u r n a l P r e -p r o o f The previous equations explicitly show the underyling SEIR structure im-450 plicit in the present reduced model, which can indeed be described as a SEIR model with fixed-duration latency time τ . If in particular τ = 0, the delayed system (30) reduces to the 2-dimensional ODE system and I i (t) (note that, since τ = 0, it holds I i (t) = I # i (t), landing on an SIR structure) satisfies the additional equation In this case, the model solution can be calculated explicitly through the roots of the characteristic equation [compare with (7)] and, for each phase From these, we finally calculate I i (t) through (32) . Note that if the initial distribution is assumed to be the stable one (16), then for the first phase we have Next, for the second phase (the "'lockdown"') we have (note this holds for any subsequent phase) I 2 (t) = I 0 2 (t) + I + 2 (t) and, in particular, The first term in (33) is due to the prevalence of infected individuals at the beginning of the second phase, which decays at a rate γ 2 (compare with (18)). The second term accounts for transmission events occurred during the lockdown phase. In this case, since α * 2 < 0, the asymptotic behaviour of I + 2 (t) depends on α * 2 + γ 2 . Indeed we have Some interesting demographic remarks can be made about the dynamics of the infective population I 2 (t) during the lockdown phase (see (33) ) depending on the sign of α * 2 + γ 2 . On the one hand, for α * 2 + γ 2 > 0 it is immediate to see that I 2 (t) asymptotically behaves as U 2 (t)/(α * 2 + γ 2 ) where (α * 2 + γ 2 ) has a clear demographic meaning i.e., it represents the instantaneous per-capita 480 birth rate of the infective population. That is to say, the infective population asymptotically evolves in a stable regime of exponential decline (given that α * 2 < 0 as resulting from the hypothesis R 2 0 < 1) where each infective individual produces on average (α * 2 +γ 2 ) new births (i.e., secondary cases) per unit of time. On the other hand, for α * 2 +γ 2 < 0, I 2 (t) asymptotically behaves as e −γ2t (I 0 2 − 485 U 0 2 /(α * 2 + γ 2 )). The interpretation for this case is that for α * 2 + γ 2 < 0 the infection reproductivity R 2 0 is brought by the lockdown to such a low level that -essentially -the infective population generated at any time t > t L does not differ significantly from the sum of the initial infective population present at the beginning of the lockdown phase I 0 2 , plus their new infective cases (given 490 by −U 0 2 /(α * 2 + γ 2 ) as in this case the per-capita birth rate is exactly given by −(α * 2 + γ 2 )) that decay over time at the removal/mortality rate γ 2 . In other words, in this case subsequent new infections generated at any time t > t L are so few that they are hardly distinguishable from the case of a population lacking reproductive power (i.e., having R 0 = 0), despite we specifically assumed 495 R 2 0 > 0. The separation between these two cases occur for α * 2 + γ 2 = 0, to which 19 J o u r n a l P r e -p r o o f it corresponds a reproduction number given by The latter expression is informative about the role played by the two functions shaping the generation time kernel during the second phase, K 2 (x). If ϕ = 0 i.e., if the transmission rate β 0 (x) is just an increasing power function, then 500 R 2, * 0 = 0. On the other hand, if ϕ > 0 then R 2, * 0 > 0 i.e., it becomes a non trivial threshold. This peculiar situation is due to the fact that in our formulation (31) , together with the prevalence function I(t); note that the incidence is discontinuous at t L . Right panel: plot of the prevalence I(t); for the second phase I 2 (t) is splitten into its two components I 0 2 (t) and I + 2 (t) as in (33); the parameter choice implies α * 2 + γ 2 > 0 so that I + 2 (t) eventually behaves as B + e α * 2 (t−t L ) ; the two components I 0 2 (t) and I + 2 (t) are plotted together with I 2 (t) and B + e α * 2 (t−t L ) . declining infectivity (occurring at the rate ϕ > 0) and removal (occurring at the rate γ 2 > 0) act as additive and independent causes of infection removal. 505 Figure 2 illustrates, for the ODE system (31), the behavior of the main variables U (t), J(t), I(t) when the change from phase 1 (invasion) to phase 2 (lockdown) occurs. The left panel shows the discontinuity of the incidence function U (t) at the switching time t L , resulting from the massive abrupt change in infection reproduction ratios (from R 1 0 = 2.0 to R 2 0 = 0.7). On the right panel 510 we analyze the prevalence I(t) and the two separated terms I 0 2 (t) and I + 2 (t) for the lockdown phase. We reported only a case with α * 2 + γ 2 > 0, the other case with α * 2 + γ 2 < 0 being very similar. The left panel of Fig. 3 illustrates the transition pattern of the prevalence function I(t) for the DDE system (30) considering different values of the re- . intrinsic demographic inertia of the infective population. We termed this type of inertia as "demographic inertia" because it exactly represents the type of inertia that arises in an age-structured population as a consequence of the switch between two stable dynamic regimes caracterized by their SAD's. This notion was 520 popularised in classical mathematical demography by Nathan Keyfitz's through the parameter he called population momentum [23] . In our model this inertia is marked at relatively high values of R 2 0 (say, in the range 0.85 < R 2 0 < 1.0), where prevalence may continue to grow for weeks after the lockdown implementation. This demographic inertia represents the key component of COVID-19 525 inertia following mitigation interventions, which can be further increased by the characteristic delays arising for COVID-19 (e.g., those due to diagnosis and cases confirmation). At lower values of R 2 0 this inertia is obviously mitigated. The right panel of the same figure compares the stable age distributions Y * 1 (x) emerging during the invasion phase to those relative to the second phase (still 530 considering different values of R 2 0 ). The difference between these stable age distributions are the ultimate responsible of the transient phase leading from the stable asymptotic behaviour during phase 1 to that of phase 2. Notably, during the first phase the stable age distribution is a piecewise-declining exponential function reflecting only epidemic growth (at rate α * 1 > 0) at ages x < τ (where 535 no removal is in place) and the additional effect of removal at ages larger than τ , when also removal starts occurring. On the other hand, during the second phase, the stable distribution is increasing at ages x < τ , reflecting the ageing 21 J o u r n a l P r e -p r o o f of the infective population, which is asymptotically declining at the rate α * 2 < 0. . Concerning the dependence of the DDE system (30) on the latency delay τ , the left panel of Figure 4 compares the model solution for infected prevalence I(t) for τ = 0 (the unlagged case) with those resulting from different positive values of τ > 0. In the latter case, the characteristic equation is transcendental and reads Thus, we may have infinitely many roots, though we are mainly interested in the leading root which is real. In the right panel of Figure 4 we plot the first roots of (35) as a function of τ . Finally, Fig. 5 reports the halving time of the various epidemic curves in the stable regime of the lockdown phase versus the reproduction number R 2 0 , 550 for different combinations of values of the latency delay τ and of the removal rate γ 2 . Note that, given the value of R 2 0 , an increase of τ produces an increase of the halving time while an increase of the removal rate γ 2 has the opposite effect. Interestingly, the shape of the relationship in Fig. 5 indicates that the greatest part of the decline in the halving time of the suppression path -which is 555 the most straighforward measure of the speed at which the community controls the epidemic -is achieved for R 2 0 values up to say 0.7. Beyond that point, further substantial achievements in decreasing the halving time would require further substantial decreases in R 2 0 . This would in turn require a dramatic strenghtening of the social distancing conditions, possibly causing an explosion . in the epidemic social and economic costs. This effect should primarily be a consequence of the specific steepness of the involved curve whose shape depends on the generation time kernel adopted. Indeed, in the case τ = 0 the halving time is given by From the latter standpoint, Fig. 5 , even if not including social costs, is illus-565 trative of the main societal tradeoff of a generalised lockdown. On one hand, brute-force epidemic suppression requires to bring R 2 0 strongly below one (perhaps below the threshold R 2, * 0 stated in (34) ). This would require a dramatic haltening of socio-economic activity, with large temporary costs, but at the advantage of getting rid of the epidemic much faster. On the other hand, the 570 attempt to keep the maximum feasible of socio-economic activities opened i.e., bringing R 2 0 just slightly below (or about) 1, allows the epidemic to continue for long time, with a persistently large and slowly declining incidence and prevalence. Though our model was mainly theoretical and aiming to investigate the basic mechanisms of a multi-phasic epidemic, it is nonetheless interesting to try to fit the model to COVID-19 data in order to highlight the potential and limits of such a simple tool with real data. Before passing to results we report a short summary of the Italian epidemic as observed so far. Protezione civile. Since then, data have been updated daily on a dedicated site (https://github.com/pcm-dpc/COVID-19) providing information on some main stock/flows (confirmed cases, hospitalised cases, ICU occupancy, deaths 585 and tests) at a regional level. However, more detailed data such as those on contact tracing and symptom onset used to compute the weekly figures of the reproduction number R t by more sophisticated methods were not made available until very late. This makes it of interest an analysis as the present one also for comparing with Governmental estimates. After a number of purely local measures, such as the creation of hotspots in most affected Northern Italy areas since February 20, the full national lockdown was completed on March 25, when only essential activities and services were allowed to be continued. This was preceded by a number of earlier Government measures including the national closure of all schools grades and universities 595 since March 4, 2020, and continued on March 11 by the first nation-wide closure of non-fundamental economic activities. Pairwise, though the lockdown was declared officially over by May 5th, some measures were continued till July 3 in the light of the still critical situation of hospitals, thereby maintaining some degrees of social distancing and individual protection. However, by the start of holidays (August 2020), the Government publicly encouraged people to avoid international travels and to spend vacations in Italy, contributing to spread optimism and confidence that the epidemic was over, possibly bringing a decline in individuals' attentions. This decline, documented by the press in an endless list of episodes and partly captured by google mobility 605 data, led to a growing pool of new infections clusters that were the premise for the second wave once workplaces and schools re-opened since the beginning of September. The subsequent collapse of the tracing system and the ensuing acceleration of the epidemic, brought to further social distancing measures since the second half of October 2020. These measures were, as a rule, aimed to 610 avoid further generalised global lockdowns by local-level targeted interventions informed by local figures such as the current reproduction number R t , the incidence of confirmed cases, and the hospital pressure. Such interventions were able to downturn the second wave and to bring the R t in the region of 1 during January and February 2021, but proved insufficient when the alpha "UK" 615 variant became predominant initiating a third epidemic wave [3] . In order to fit reported data with the simplest versions of our model, namely those analysed in Section 4.2 using parametrization (27) , we made a few simplifying assumptions. First, we considered data at the National level. This is not inappropriate for the phases after the first one where epidemic circulation 620 involved the entire country, possibly it is not for the first phase where the epidemic was localised in a few provinces of a few regions. Moreover, being our model highly stylised with a minimal parametrization, it does not include further epidemiological compartments (e.g., isolation, hospitalizations, ICU, and deaths) which are necessary for realistic epidemic modeling. Therefore, we sim- . observed every week-end (Fig.6) . We deliberately did not detrend data by using the time series of tests actually carried out. Assuming, quite crudely, that 630 newly confirmed cases are isolated and therefore removed from the active infective population, newly reported cases are considered as removed individuals R i (t), satisfying the equation We therefore can use the formula: where γ i is the removal rate of the corresponding phase, getting a sequence I n 635 to be fitted against the theoretical prevalence I # (t) provided by the model. The generation time kernel K i (x) (27) was chosen based on available estimates of the generation time distribution for COVID-19 in Italy [32] . In that study, data on observed serial intervals (i.e., the distance between symptoms in primary infectors and symptoms in their secondary cases), were fitted by a Guzzetta et al best fit to italian data best least square Erlang(2) fit best mean−preserving Erlang(2) fit Figure 7 : Graphic comparison between the generation Gamma kernel estimated for Italy in [32] vs two possible Erlang approximation with index 2 namely (i) the mean preserving Erlang distribution used in our application in this section, (ii) the best fitting Erlang distribution computed by a least squares criterion. approximation of it (see Fig. 7 ), by simply taking (an analogous approach was used in [33] ) To fit our model throughout the various epidemic phases (Fig.6, left panel) , we checked in the data for evidence of the presence of the key model feature, namely the emergence of epochs of stable exponential behavior eventually resulting from the implementation (or the relaxation) -for a sufficiently long period of time -of mitigation interventions. This is clear from the log-scale diagram 650 in the right panel of Fig. 6 where we could well identify at least four phases characterised by an initial transient dynamics which eventually approaches an exponential trend. These four phases are represented by (i) the initial phase of epidemic growth (roughly days 0-19), (ii) the declining phase after the peak induced buy the lockdown (roughly days 80-100), (iii) the phase of exponential 655 increase during the second wave (days 210-250), (iv) the phase of the exponential decline resulting from the measures adopted to mitigate the second wave (days 275-290). This shape in turn suggests a 2-stage fitting procedure, to be repeated for each phase transition, where: (i) in the first stage the characteristic root α * i 660 for the particular phase considered is estimated, (ii) in the second stage the remaining parameters (τ , γ i , t i .) are adjusted conditional on α * i (actually τ is fitted only once because it is a biological constant pertaining to the infection and it will not plausibly change through the different phases). In what follows, we detail the adopted approach for the first two phases and briefly summarise As discussed in Section 4.2, we assume that the epidemic was initiated by an initial cohort of infective individuals distributed according to the corresponding stable age distribution of the invasion phase. Consequently, the initial growth 670 of the epidemic curve is expected to be exponential since initial time with rate α * 1 . This parameter was estimated by a (least squares) exponential fit chosing the data window yielding the highest value of the determination coefficient R 2 . This yielded I 0 = 181 α * 1 = 0.15 day −1 , over the period from day n = 0 to n = 19, with R 2 ≈ 0.99. The fact that 675 including further data points beyond n = 19 day always worsened the fit, may be interpreted as due to the onset of social distancing gradually showing its effect soon after time n = 19 day, that is since March 15. At a crude glance, the lockdown phase might seem to display its effects on This in turn allowed to also estimate the reproduction number R 1 0 of the first phase as well as that of the second phase: The overall fit of the model ( Concerning the subsequent phases, namely (i) the second wave phase which started due to the relaxation of most measures adopted during the lockdown, (ii) the subsequent mitigation phase with the so-called "multicolor strategy" [3] 705 based on locally targeted (typically: at the regional level) measures in view of the documented seriousness of the local situation, up to (iii) the current emergence of more transmissible variants of the virus, we used the same procedure detailed above. The results are reported in Table 1 and Figure 9 . The estimates of the phase-specific reproduction numbers R i 0 are systemat- Fig. 9 , provides a summary view of the explanatory power of the model (30), by comparing its solution with the data over the entire course of the epidemic observed in Italy so far. Note that we deliberately avoided to include 715 generalised smooth transitions between phase specific reproduction numbers just for purposed of obtaining better-looking fits. Note also that the poorer fit over the last two phases is the consequence of their short duration compared to phases I-IV, preventing the emergence of any meaningful stable exponential phase. 720 To sum up, the overall fit in Figure 9 was very good almost everywhere but in the most recent phases, despite the little information used, namely the data from the sole exponential windows in each phase. Essentially, prior to the recent period made complicated by the emergence of the new variants, there are only 725 two periods where some lack of fit emerges, namely a portion of the declining trend during the lockdown phase and the first epoch of the re-opening phase. This mostly lies in (i) the minimalistic structure of our model which does not explicitly include any compartment directly comparable with real data, (ii) the stylised hypothesis by which switches from one phase to the successive one occur 730 abruptly. Refined data explanations would require to add a number of further compartments to account for the delayed processes of the disease pathways, testing, isolation, hospitalization etc, and to amend our homogenous mixing framework to include heterogeneities in transmission and disease (over age and space) and testing. This said, the fitted model is the simplest setup illustrating the naturally embedded demographic inertia of COVID-19 dynamics despite the lack of any delayed dimensions (such as ICU and deaths). Here, by "demographic inertia" we mean the natural tendency of a multi-phasic epidemic to display its true phase-specific trends with a temporal delay compared to the moment when that . intervention. This demographic inertia is to be distinguished from other components of COVID-19 observed inertia that are due to the fact that its main public health outcomes e.g., cases confirmations (not to say of hospitalizations and deaths) arise with substantial time-delays with respect to the time when 745 transmission occurred. This natural demographic inertia can be macroscopic or negligible, other parameters being equal, depending on the relative magnitude of the reproduction numbers R i 0 , R i+1 0 that are characteristic of two successive phases. This was illustrated in the left panel of Fig. 3 . Therefore, looking at Fig. 9 we can say that: (i) the continued growth in 750 prevalence after the first lockdown declaration is largely attributable to delayed patterns rather than to demographic inertia (as can be understood from the model lack of fit at the switch between the first two phases). The reason is that the lockdown brought the reproduction number representative of the second phase down to a quite low level. On the other hand, the continued decline in 755 prevalence for a very long time beyond the lockdown official end-date at May 5 (about time n = 80 day) can be partly explained by demographic inertia exactly due to the very low level of R i 0 eventually achieved during the lockdown phase. Borrowing from a language typical of vaccination programs, the latter prolonged honey-moon period likely surely contributed to the generation of an 760 optimism wave in the public opinion that eventually was responsible for the second pandemic wave. J o u r n a l P r e -p r o o f Based on a simple hypothesis, namely that the depletion of the susceptible population remains contained over time, in this manuscript we have proposed 765 a minimal SIR model for the overall dynamics of an epidemic evolving in a multi-phasic form due to a sequence of epochs of mitigation intervention and their relaxations, as has been the case for the COVID-19 epidemic. Our model provides, first of all, a unified representation of such multi-phasic epidemics, describing how the different epidemic phases unfold over time as a consequence 770 of the interplay between the transmission process and intervention responses. Second, by using a wide class of generation time kernels allowing reducibility (either to ordinary or delayed differential equations), we investigated in depth a low-dimensional case allowing a non-trivial full analytical treatment also of the transient dynamics connecting the different epidemic phases. The latter 775 model provided an excellent fit to the entire course of the epidemic observed in Italy since February 2020, despite its parsimonious parametric structure, not including any disease compartments and using minimal data. Last, the proposed model represents the simplest setup illustrating the intrinsic demographic component of the inertia of COVID-19 dynamics. This demographic inertia adds 780 to the one due to the various delayed phenomena characteristic of COVID-19 to determine the overall COVID-19 inertia. The limitations of the proposed model lie in its parsimonious structure and in its kinematic nature. As for the first drawback, this can be easily overcome by simply adding disease-related and other relevant compartments, as well as by introducing relevant heterogeneities, 785 as those due to chronological age such as, e.g., in the onset of symptoms and serious disease. This however represents only a part of the story of COVID-19 which is clearly more complicated than described by simple models. For example, we described the incoming story of the second COVID-19 wave, debuted in Italy since the start of September 2020, in a short report (Iannelli et al. [34] ) 790 where we correctly predicted that the observed epidemic was dramatically accelerating since the beginning of October, with R t shifting upward, and would continue to do so for a while, due to the overwhelming of the testing-tracing system by epidemic growth. We also suggested that the subsequent epidemic increase would have continued up to a maximum corresponding to the reproduc-795 tion number generated in the preceding weeks by the behavior of risky groups, and would have halted or slowed down once the more prudent population groups would be reached. This is suggestive of the fact that, at a fine scale, the epidemic description requires to account for a number of further dimensions. These include saturation effects, as the finiteness of tracing resources, capable -if alone 800 -to slow-down the epidemic but only for a while, as well as to go beyond the kinematic representation to include a number of behavioral aspects. Behavioral effects are numerous but -just for illustrative purposes -one of these is surely understandable as a consequence of the inertia phenomenon that we termed the post-lockdown honey-moon period. This prolonged period of low epidemic 805 activity surely contributed to the generation of the optimism wave in the Italian public that contributed to the drop of prudent behaviors, and eventually 31 J o u r n a l P r e -p r o o f accelerated the second pandemic wave. This calls for behavioral epidemiology explanations ( [26, 27] ). Learning from such aspects is of dramatic importance also in prospective terms, namely the control of future pandemic events. and nomenclature coming from the demographic theory of the stable population. Under rather general assumption such as • F (t) ≥ 0 (t ∈ [0, ∞)) is continuous and absolutely Laplace transformable, • K(x) ≥ 0 (x ∈ [0, ∞)) is absolutely Laplace transformable, equation (A.1) has a unique continuous non-negative solution u(t) (t ∈ [0, ∞)) whose asymptotic behaviour has been extensively discussed by Feller [35] and is related (under some general assumption that we suppose fulfilled by the equation coming from our model) to the roots of the characteristic equation 820 K(λ) = 1. where we use symbol f to denote the Laplace transform of the function f . Following Feller [35] , it is known that this equation has one and only one leading root α * , which is real, simple and strictly positive (negative) if the condition K(0) > 1 (respectively K(0) < 1) is satisfied. Precisely we have the sequence α * , α j j=∞ j=1 such that 825 α j+1 ≤ α j < α * (j = 1, 2, . . .) and the solution has the asymptotic expansion u(t) = u * e α * t + is an estimate of the minimum time required to get the limit (A.3) reached within a tolerance . Reduction of the main system (4) to a cascade of ODE-DDE equation is 840 allowed by the special form of the basic parameters (24), (25) . Here we illustrate the method through the special case (27) . The starting point are the variables (28) and (29) . From (4) we see that the latter satisfy the reduced system (B.1) WHO, Covid-10 situation reports 1-209 and weekly epidemiological updates Impact of non-pharmaceutical interventions (npis) to reduce covid-19 mortality and healthcare demand. imperial college covid-19 865 response team, Imperial College COVID-19 Response Team Covid-19 epidemic in italy. national and regional upgrade (in italian) Mod-875 elling the impact of testing, contact tracing and household quarantine on second waves of covid-19 Early insights from statistical and mathematical modeling of key epidemiologic pa-880 rameters of covid-19 The effect of travel restrictions on the spread of the 2019 novel coronavirus Covid-19 and sars-cov-2. modeling the present, looking at the future Quantifying SARS-CoV-2 transmission 890 suggests epidemic control with digital contact tracing Estimating the effects of non-pharmaceutical interventions on covid-19 in europe Spread and dynamics of the covid-19 epidemic in italy: Effects of emergency containment measures Modelling the covid-19 epidemic and implementation of population-wide interventions in italy The effect of human mobility and control measures on the covid-19 epidemic in china Early dynamics of 910 transmission and control of covid-19: a mathematical modelling study Possible effects of mixed prevention strategy for covid-19 epidemic: massive testing, quarantine and social distancing A covid-19 epidemic model with latency period Predicting the number of reported and unreported cases for the covid-19 epidemics in china, south korea, italy, france, germany and united kingdom Applications of predictive modelling early in the COVID-19 epidemic Modelling covid-19 The impact of covid-19 and strategies for mitigation and suppression in low-and middle-income countries Changes in contact patterns shape the dynamics of the covid-19 outbreak in china First results of the national serosurvey on sars The basic approach to age-structured population dynamics Age-structured population dynamics in demography and epi-950 demiology Modeling the interplay between human behavior and the spread of infectious diseases Statistical physics of vaccination Estimation in emerging epidemics: Biases and remedies Biological Delay Systems: Linear Stability Theory Pulse and constant control schemes for epidemic models with seasonality Infectious Disease Modeling The impact of a nation-wide lockdown on covid-19 transmissibility in italy Effects of latency and age structure on the 970 dynamics and containment of covid-19 Se l'epidemia accelera cosa resta da fare? On the integral equation of renewal theory The authors would like to thank the anonymous referees whose suggestions improved this manuscript.