A time-modulated Hawkes process to model the spread of COVID-19 and the impact of countermeasures
AA time-modulated Hawkes process to model the spreadof COVID-19 and the impact of countermeasures (cid:63)
Michele Garetto a, ∗ , Emilio Leonardi b , Giovanni Luca Torrisi c a Universit`a degli Studi di Torino, C.so Svizzera 185, Torino, Italy b Politecnico di Torino, C.so Duca degli Abruzzi 24, Torino, Italy c IAC-CNR, Via dei Taurini 19, Roma, Italy
Abstract
Motivated by the recent outbreak of coronavirus (COVID-19), we propose astochastic model of epidemic temporal growth and mitigation based on a time-modulated Hawkes process. The model is sufficiently rich to incorporate specificcharacteristics of the novel coronavirus, to capture the impact of undetected,asymptomatic and super-diffusive individuals, and especially to take into ac-count time-varying counter-measures and detection efforts. Yet, it is simpleenough to allow scalable and efficient computation of the temporal evolutionof the epidemic, and exploration of what-if scenarios. Compared to traditionalcompartmental models, our approach allows a more faithful description of virusspecific features, such as distributions for the time spent in stages, which iscrucial when the time-scale of control (e.g., mobility restrictions) is comparableto the lifetime of a single infection. We apply the model to the first and sec-ond wave of COVID-19 in Italy, shedding light onto several effects related tomobility restrictions introduced by the government, and to the effectiveness ofcontact tracing and mass testing performed by the national health service.
Keywords:
COVID-19, Compartmental models, Branching process, Hawkesprocess
1. Introduction and related work
Models of epidemic propagation are especially useful when they provide keyinsights while retaining simplicity and generality. For example, in mathematicalepidemiology the SIR model reveals in very simple terms the fundamental role ofthe basic reproduction number ( R ) which governs the macroscopic, long-term (cid:63) The simulation code and data used in this work are available under GPL v3 on GitHub: https://github.com/michelegaretto/covid.git ∗ Corresponding author
Email addresses: [email protected] (Michele Garetto), [email protected] (Emilio Leonardi), [email protected] (GiovanniLuca Torrisi)
Preprint submitted to Annual Reviews in Control January 5, 2021 a r X i v : . [ q - b i o . P E ] J a n volution of the outbreak in a homogeneous population. The vast majority ofmodels developed for the novel SARS-CoV-2 are extensions to the classic SIR,along the standard approach of introducing additional compartments to describedifferent phases of the infection, the presence of asymptomatic, symptomatic orpauci-symptomatic individuals, the set of quarantined, hospitalized people, andso on. Such models lead to a system of coupled ODE’s with fixed or time-varyingcoefficients to be estimated from traces. A very incomplete list of modelingefforts pursued in this direction, early applied to COVID-19, includes the SEIRmodels in [1, 2, 3, 4, 5], the SIRD models in [6, 7], the SEPIA model in [8], theSIDHARTE model in [9].In this paper we develop a slightly more sophisticated model that allows amore accurate representation of native characteristics of a specific virus, such asactual distributions of the duration of incubation, pre-symptomatic and symp-tomatic phases, for various categories of infected. This level of detail is im-portant when the intensity of applied countermeasures varies significantly overtime-scales comparable to that of an individual infection, and we believe it isessential to address fundamental questions such as: i) when and to what extentcan we expect to see the effect of specific mobility restrictions introduced by anational government at a given point in time? ii) what is the impact of hard vspartial lockdowns enforced for given numbers of days? when can restrictions besafely released to restart economic and social activities while still keeping theepidemic under control?Specifically, we propose and analyse a novel, modulated version of the markedHawkes process, a self-exciting stochastic point process with roots in geophysicsand finance [10]. In a nutshell, in the standard marked Hawkes process eachevent i with mark m i , occurring at time t i , generates new events with stochasticintensity ν ( t − t i , m i ), where ν ( · , m i ) is a generic kernel. The process unfoldsthrough successive generations of events, starting from so-called immigrants(generation zero). In our model events represent individual infections, and amodulating function µ ( t ), which scales the overall intensity of the process atreal time t , allows us to take into account the impact of time-varying mobilityrestrictions and social distance limitations. In addition, we model the transi-tion of infected, undetected individuals to a quarantined state at inhomogeneousrate ρ ( t ), to describe the time-varying effectiveness of contact tracing and masstesting.We mention that branching processes of various kind, including Hawkes pro-cesses, have been proposed in various biological contexts [11, 12, 13, 14]. In [15],a probabilistic extension of (deterministic) discrete-time SEIR models, based onmulti-type branching processes, has been recently applied to COVID-19 to cap-ture the impact of detailed distributions of the time spent in different phases,together with mobility restrictions and contact tracing. In parallel to us, au-thors of [16] have proposed a Hawkes process with spatio-temporal “covariates”for modeling COVID-19 in the US, together with an EM algorithm for parame-ter inference. The work in [16] is somehow orthogonal to us, since they focus onspatial and demographic features, aiming at predicting the trend of confirmedcases and deaths in each county. In contrast, we introduce marks and stages to2atively model the course of an infection for different categories of individuals,with the fundamental distinction between real and detected cases. Moreover,we take into account the impact of time-varying detection efforts and contacttracing. At last, we obtain analytical expressions for the first two moments ofthe number of individuals who have been infected within a given time.We demonstrate the applicability of our model to the novel COVID-19 pan-demic by considering real traces related to the first and second wave of coro-navirus in Italy. Our fitting exercise, though largely preliminary and basedon incomplete information, suggests that our approach has good potential andcan be effectively used both for planning counter-measures and to provide ana-posteriori explanation of observed epidemiological curves.The paper is organized as follows. In Section 2 we provide the mathematicalformulation of the proposed modulated Hawkes process to describe the temporalevolution of the epidemic. In Section 3 we motivate our approach by comparingit to the standard SIR model. Some mathematical results related to the momentgenerating function of our process are presented in Section 4. In Section 5 wedescribe our COVID-19 model based on the proposed approach. We separatelyfit our model to the first and second wave of COVID-19 in Italy in Sections 6and 7, respectively, offering hopefully interesting insights about what happenedin this country (and similarly in other European countries) during the recentpandemic. We conclude in Section 8.
2. Mathematical formulation of modulated Hawkes process
We first briefly recall the classic Hawkes process restricted to the tempo-ral dimension (the spatio-temporal formulation is similar, but we focus in thiswork on the purely temporal version). Events of the process (or points) oc-cur at times T = { T i } i ≥ , which are R -valued random variables. A subset I of these points, called immigrants, are produced by an inhomogeneous Poissonprocess of given intensity σ ( t ). Each immigrant, independently of others, isthe originator of a progeny (or cluster) of other points, dispersed in the futurethrough a self-similar branching structure: a first generation of points is pro-duced with intensity ν ( t − T j ), where T j is the occurrence time of immigrant j ,and ν : R + → R + is a kernel function. Each point of the first generation, in turn,generates new offsprings in a similar fashion, creating the second generation ofpoints, and so on.The above process can be easily extended to account for different typesof points with type-specific kernel functions. Types are denoted by marks M = { M i } i ≥ , which are assumed to be i.i.d. random variables with val-ues on an arbitrary measurable space (M , M ), with a probability distribu-tion Q . Let N ( t, m ) be the counting process associated to the marked points N = { ( T i , M i ) } i ≥ . The (conditional) stochastic intensity λ ( t ) of the overallprocess is then given by: λ ( t ) = σ ( t ) + (cid:90) t ν ( t − s, m ) N ( d s, d m ) = σ ( t ) + (cid:88) T k ∩ (0 ,t ) ν ( t − T k , M k )3here ν : R + × M → R + is a type-dependent kernel function. In the followingwe assume ν ( t ) := E [ ν ( t, M )] to be summable,. i.e. R = (cid:90) ∞ ν ( t )d t < ∞ (1)Note that R is the average number of offsprings generated by each point, whichis usually referred to as basic reproduction number in epidemiology. The processis called subcritical if R <
1, supercritical if R > N isto modulate the instantaneous generation rate of offsprings N (cid:48) = N \ I by apositive, bounded function µ ( t ) : R + → R + , representing the impact of mobilityrestriction countermeasures. By so doing, we obtain the modified stochasticintensity of the process: λ ( t ) = σ ( t ) + µ ( t ) (cid:20)(cid:90) t ν ( t − s, m ) N ( d s, d m ) (cid:21) = σ ( t ) + µ ( t ) (cid:88) T k ∩ (0 ,t ) ν ( t − T k , M k ) (2)Note that, when µ ( t ) = µ is constant, we re-obtain a classic Hawkes process withmodified kernel µ ν (). In general, the obtained process is no longer self-similar.In particular, the average number of offsprings generated by a point becomes afunction of time: R ( t ) = (cid:90) ∞ µ ( t + τ ) ν ( τ ) d τ (3)which provides the infamous real-time reproduction number usually referred toon the media as ‘Rt index’.We emphasize that the process can initially start in the supercritical regime,and then it can become subcritical for effect of a decreasing function µ (), acase of special interest in our application to waves of COVID-19. We mentionthat the great bulk of literature related to the Hawkes process and its applica-tions to geophysics and finance focuses on the subcritical regime, whereas ourformulation applies also to the supercritical regime, which is more germane toepidemics.The conditional intensity (2) can be easily de-conditioned with respect to4 , obtaining the ‘average’ stochastic intensity λ ( t ) := E [ λ ( t )]: λ ( t ) = σ ( t ) + µ ( t ) E (cid:34)(cid:90) (0 ,t ) × M ν ( t − s, m ) N (d s, d m ) (cid:35) = σ ( t ) + µ ( t ) E (cid:34)(cid:90) (0 ,t ) × M ν ( t − s, m ) λ ( s ) d s Q (d m ) (cid:35) = σ ( t ) + µ ( t ) E (cid:20)(cid:90) t ν ( t − s ) λ ( s )d s (cid:21) = σ ( t ) + µ ( t ) (cid:90) t ν ( t − s ) λ ( s )d s. (4)where we recall that ν ( t ) = E [ ν ( t, M )].We observe that (4) is a linear Volterra equation of the second kind, whichcan be efficiently solved numerically. In the special case of constant modulatingfunction, (4) reduces to a convolution equation which can be analyzed and solvedby means of Laplace transform techniques (see Appendix A).At last we introduce the total number of points up to time t , N ( t ) (regardlessof their associated marks), and its average: N ( t ) = (cid:90) t λ ( τ ) d τ (5)
3. Comparison with SIR model
When the kernel function has an exponential shape, i.e., ν ( t ) = Ke − γt , itis possible to establish a simple connection between the Hawkes process andthe classic SIR model [17]. Specifically, consider a stochastic SIR model withan infinite population of susceptible individuals, where each infected generatesnew infections at rate β , and recovers at rate γ (i.e., after an exponentiallydistributed amount of time of mean 1 /γ ). Then the average intensity of theprocess generated by this stochastic SIR, averaging out the times at which nodesrecover, has exactly the (conditional) intensity of a Hawkes process startingwith the same number of initially infected nodes, no further immigrants, and ν ( t ) = βe − γt [17]. This can be intuitively understood by considering that inSIR the average, effective rate with which an infected individual generates newinfections after time t since it became infected equals β times the probabilitythat the node has not yet recovered, which is e − γt .We remark that in [17] authors push this equivalence a bit further, by show-ing that a SIR model with finite population N is equivalent to a modulatedHawkes process similar to (2), where µ ( t ) = 1 − N ( t ) /N . In our model and its For example, the standard trapezoidal rule allows obtaining a discretized version of λ ( t )by a matrix inversion. . Therefore, wedo not need the ‘correction factor’ 1 − N ( t ) /N , and instead use the modulatingfunction µ ( t ) to model the (process-independent) effect of mobility restrictions.A similar effect due to mobility restrictions can be incorporated into the SIRmodel, by applying factor µ ( t ) to the infection rate β .Given the above connection between a (possibly modulated) Hawkes processand the traditional SIR, one might ask what is the benefit of our approach withrespect to SIR-like models. In contrast to SIR, the Hawkes process allows us tochoose an arbitrary kernel ν ( t ), not necessarily exponential. With this freedom,can we observe dynamics significantly different from those that can be obtainedby a properly chosen exponential shape? We answer this question in the positivewith the help of an illustrative scenario.Consider an epidemic starting at time zero with I = 1000 infected indi-viduals. We fix to g = 10 (days) the average generation time, which is themean temporal separation between a new infection belonging to generation i and its parent in generation i −
1. Note that the above constraint implies that (cid:82) ∞ t ν ( t ) d t = g .We also normalize (cid:82) ∞ ν ( t ) d t = 1, since we can use µ ( t ) to scale the infectionrate, in addition to considering time-varying effects due to mobility restriction.Note that in SIR we can only satisfy the above constraints by choosing theexponential kernel ν ( t ) = g e − tg , t >
0. Instead, in our Hawkes model we havemuch more freedom.First, we consider a scenario in which µ ( t ) is constant, equal to either 1.2(supercritical case) or 0.9 (subcritical case). In Fig. 1 we show the temporalevolution of the mean number of infected N ( t ), for three kernel shapes thatallows for an explicit solution of (4),(5) using Laplace transform (Appendix A).The delta shape corresponds to the kernel function ν ( t ) = δ ( t − g ), where δ ( · )is Dirac’s delta function, for which N delta ( t ) = I µ (cid:98) t/g (cid:99) +1 − µ − exp shape corresponds to the exponential kernel ν ( t ) = g e − tg (SIRmodel), for which N exp ( t ) = I (cid:20) µµ − (cid:16) e t ( µ − /g − (cid:17)(cid:21) (7)The hyper shape corresponds to the hyper-exponential kernel: ν ( t ) = p α e − α t + p α e − α t , where 0 ≤ p ≤ p = 1 − p , p /α + p /α = g , This assumption is largely acceptable at the beginning of an epidemic. Note that in our model the number of immigrants has a Poisson distribution with mean (cid:82) ∞ σ ( t ) d t . However to reduce the variance we have preferred to start simulations with thesame deterministic number of infected. Further, note that analytical results for the meantrajectory of infected nodes are not affected by this choice.
0 50 100 150 200 250 300 350 m ea n nu m b e r o f i n f ec t e d time (days)hyper 3-30expdelta µ = 0.9µ = 1.2 Figure 1: Evolution of mean number of infected for different kernel shapes, in the case ofconstant µ = 1 . µ = 0 . I = 1000, g = 10. which also permits obtaining an explicit, though lengthy expression of N ( t ) thatwe omit here. Specifically, the ‘hyper’ curve on Fig. 1 corresponds to the case1 /α = 3, 1 /α = 30.We observe that, in the subcritical case ( µ = 0 . µ = 1 . e ηt )(as t grows large), where, interestingly, η > y axes on Fig. 1). In particular, the hyper-exponential kernel canproduce arbitrarily large η (Appendix A). We conclude that, even when we fixthe average generation time g , different kernels can produce largely different (inorder sense) evolutions of N ( t ). Note that, by introducing compartments, SIR-like models can match higher-order moments of the generation time, but ourresults suggest that N ( t ) depends on all moments of it, i.e., on the precise shapeof the kernel. Moreover, recall that some shapes are difficult to approximateby a phase-type approximation (e.g., the rectangular shape, or more in general,kernels with finite support).The strong impact of the specific kernel shape becomes even more evidentwhen we consider a time-varying µ ( t ), as in our modulated Hawkes process.As an example, consider the COVID-inspired scenario in which the modulating7unction µ ( t ) corresponds to the black curve in Fig. 2: during the first 30 days, µ ( t ) decreases linearly from 3 to 0.3; it stays constant at 0.3 for the next 60days; and it goes back linearly to 3 during the next 30 days, after which it staysconstant at 3. m ea n nu m b e r o f i n f ec t e d µ ( t ) time (days)expuniformdelta µ (t) Figure 2: Evolution of mean number of infected (left y axes) in the COVID-like scenario, fordifferent kernel shapes. Shaded areas denote 95% − level confidence intervals obtained by 100simulation runs. The right y axes refers to modulating function µ ( t ) (black curve). In Fig. 2 we show the mean number of infected estimated by simulation(averaging 100 independent runs), for the three kernel shapes exp , delta , and uniform , where uniform corresponds to the kernel ν ( t ) = g , t ∈ (0 , g ), whilethe shapes exp , delta have been already introduced above. Shaded areas aroundeach curve denote 95% − level confidence intervals.We observe huge discrepancies among the trajectories of N ( t ) obtained underthe three kernels. In particular, after the first ‘wave’ of 30 days, the exp kernelproduces about four times more infections than those produced by the delta kernel. This can be explained by the fact that µ ( t ) varies significantly on a timewindow (30 days) comparable to the average generation time (10 days).It is also interesting to compare what happens on the second wave startingat day 90, after all 3 curves have settled down to an almost constant value:now discrepancies are even more dramatic: although we observe a very fastresurgence of the epidemic in all cases, this happens with significant delays fromone curve to another. This is due to the fact that the number of individuals whoare still infectious after the subcritical period (when µ <
1) is largely different,especially considering that both the uniform and delta kernels have finite support(20 and 10 days, respectively) much smaller than the duration of the subcritical8eriod, whereas the exp shape has infinite support: this implies that under the uniform and delta kernels the virus survives the subcritical period only throughchains of infections belonging to successive generations, whereas under the exp shape in principle the epidemic can restart just thanks to the original immigrantsat time 0, who are still weakly infectious when we re-enter the supercriticalregime (around day 100). Actually, in our experiment, under the delta kernelthe epidemic died out in 68 out of 100 simulation runs, which explains the largeconfidence intervals obtained in this case at end of the observation window.Under the uniform kernel, the epidemic died only in 3 out of 100 runs, while italways survived under the exp kernel.We conclude that, even while fixing the mean generation time, the preciseshape of the kernel function can play an important role in predicting the processdynamics and the impact of countermeasures, especially when the time-scale ofcontrol is comparable to the time-scale of an individual contagion. One mightobtain a good fitting with measured data also by using an exponential shape,and a properly chosen modulating function µ ( t ), but this is undesirable, sincethe required µ ( t ) would no longer reflect the actual evolution of mobility andinterpersonal contact restrictions. In the case of COVID-19, several researchershave actually attempted to incorporate into analytical models detailed informa-tion about people mobility, using for example data provided by cellular networkoperators or smartphone apps [18, 19].
4. Moment generating function
Our modulated Hawkes process is stochastic in nature, hence it is impor-tant to characterize how realizations of the process are concentrated aroundthe mean trajectory derived in Sec. 2. This characterization is instrumental,for example, in designing simulation campaigns with proper number of runs.Moreover, note that it is entirely possible that the epidemic gets extinct at itsearly stages, or in between two successive waves, as we have seen in the sce-nario in Fig. 2. Actually, an epidemic could die out even when starting inthe supercritical regime (consider the case of a single immigrant at time t , whodoes not generate any offspring with probability e − R ( t ) ), something that is notcaptured by deterministic mean-field approaches. Therefore, it is interesting tounderstand the variability of the process at any time.Under mild assumptions an expression for the moment generating functionof the number of points in [0 , t ) can be given, and an iterative procedure can beapplied to compute the moment of any order n in terms of moments of smallerorder, though with increasing combinatorial complexity.In this section we limit ourselves to reporting the main results on this mo-ments’ characterization, without their mathematical proofs, to keep the paperfocused on the application to COVID-19 and its control measures.Hereon, we shall assume K := sup ( t,m ) ∈ (0 , ∞ ) × M ν ( t, m ) ∈ (0 , ∞ ) , K := sup t ∈ (0 , ∞ ) µ ( t ) ∈ (0 , ∞ ) . (8)9lso, we fix a horizon τ ∈ (0 , ∞ ), and, for t ∈ (0 , τ ), we denote with S t,τ thenumber of points, up to time τ , in the cluster generated by an immigrant at t (including the immigrant). We denote with | · | the modulus of a complexnumber.We are interested in N ( τ ), i.e., the total number of points generated up totime τ , irrespective of their mark. Theorem 4.1 (Moment generating function of N ( τ ) ). Assume (1) and (8) .Then there exists θ c > such that, for any z ∈ Θ c , Θ c := { z ∈ C : Re z < θ c } , we have E [e zN ( τ ) ] = exp (cid:18)(cid:90) τ ( G ( t, z ) − σ ( t )d t (cid:19) (9) and sup z ∈ Θ c | E [e zN ( τ ) ] | < ∞ . (10) Here G ( t, z ) := E [e zS t,τ ] , (11) is the solution of the following functional equation. G ( t, z ) = e z E (cid:104) e (cid:82) τ ( G ( v,z ) − µ ( v ) ν ( v − t,M ) d v (cid:105) , ( t, z ) ∈ (0 , τ ) × Θ c . (12)From the moment generating function we can obtain an expression for thefirst and second moment of N ( τ ) (the first moment can also be obtained directlyfrom the average stochastic intensity, as we have done in Sec. 2.). Theorem 4.2 (The first two moments of N ( τ ) ). Assume (1) and (8) . Then E [ N ( τ )] = N ( τ ) = (cid:90) τ λ ( t ) d t = (cid:90) τ E ( t ) σ ( t ) d t (13) and E [ N ( τ ) ] = (cid:90) τ E ( t ) σ ( t ) d t + (cid:18)(cid:90) τ E ( t ) σ ( t ) d t (cid:19) , (14) with E ( t ) = G (cid:48)(cid:48) ( t,
0) = E [ S t,τ ] = 1 + (cid:90) τt E ( v ) µ ( v ) ν ( v − t ) d v ++ 2 (cid:90) τt E ( v ) µ ( v ) ν ( v − t ) d v + (cid:18)(cid:90) τt E ( v ) µ ( v ) ν ( v − t )d v (cid:19) (15) E ( t ) = G (cid:48) ( t,
0) = E [ S t,τ ] = 1 + (cid:90) τt λ ( t ) ( v ) d v (16) and λ ( t ) ( v ) = µ ( v ) ν ( v − t ) + µ ( v ) (cid:90) vt ν ( v − u ) λ ( t ) ( u ) d u. (17)10ote that (15) and (17) are second-type inhomogeneous Volterra equations.In the particular case in which µ () is constant, solutions for (15) and (17) canbe found by applying a standard Laplace transform methodology.
5. COVID-19 Model
Now we describe how we applied the modulated Hawkes process introducedbefore to model the propagation dynamics of COVID-19. The proposed modelcould actually be used to represent the dynamics of other similar viruses as well.First, we take advantage of the fact that we do not need to consider a uniquekernel function for all infected. Indeed, the presence of marks allows us to in-troduce different classes of infectious individuals with specific kernel functions.Specifically, we have considered three classes of infectious: symptomatic , asymp-tomatic and superspreader , denoted by symbols { s, a, h } . We assume that, whena person gets infected, it is assigned a random class C ∈ { s, a, h } with prob-abilities p s , p a , p h , respectively, p s + p a + p h = 1, 0 ≤ p s ≤
1, 0 ≤ p a ≤ ≤ p u ≤ symptomatic people are those who will develop evi-dent symptoms of infection, and we assume that because of that they will beeffectively quarantined at home or hospitalized at the onset of symptoms. Onthe contrary, the asymptomatic mark is given to individuals who will not de-velop strong enough symptoms to be quarantined. Therefore, they will be ableto infect other people for the entire duration of the disease, though at low in-fection rate (unless they get scrutinized by mass testing, as explained later). Atlast, superspreaders are individuals who exert a high infection rate but do notget quarantined due to several possible reasons (unless they get scrutinized bymass testing). This class also includes people with mild symptoms, who becomehighly contagious because of their mobility pattern (e.g., participation to ‘su-perspreading events’). Though the above classification of infectious individualsis a simplified one, a properly chosen mix of the three considered classes canrepresent a wide range of different scenarios.Irrespective of their class, we will assume that all infectious people go throughthe following sequence of stages: first, there is a random incubation time, de-noted by r.v. I , with given cdf F I (). During this time we assume the all infectedexert a low infection rate λ low . Then, there is a crucial pre-symptomatic period,during which infected in classes { s, h } already exert a high infection rate λ high ,while infected in class a still exert low infection rate λ low . For simplicity, wehave assumed the presymptomatic phase to have constant duration w .The following evolution of the infection rate of an individual depends onthe class: since we assume that symptomatic people get effectively quarantined,they no longer infect other individuals after the onset of symptoms. We modelthis fact introducing a quarantined class of people (denoted by q ), and deter-ministically moving all infected in class s to class q after time I + w .People in classes { a, h } continue to be infectious during a disease period ofrandom duration D , with given cdf F D (). The difference between these two11lasses is that infectious in class a ( h ) exert, during the disease period, infectionrate λ low ( λ high ), respectively. At last, we assume that people in classes { a, h } enter a residual period of random duration E , with cdf F E (), during which allof them exert infection rate λ low . We introduced this additional phase becausesome people who recovered from COVID-19 where found to be still contagiousseveral days (even weeks) after the end of the disease period. Note that durations I, D, E are assumed to be independent, and that the complete mark associatedto an infected is the tuple M = ( C, I, D, E ). symptomaticasymptomaticsuperspreader λ high λ low λ low λ high λ low I w tI w D E tI w D E tν ( t, s ) ν ( t, a ) ν ( t, h ) Figure 3: Kernel functions for classes s, a, h (from top to bottom), for given durations
I, D, E of stages.
An illustration of the three class-dependent kernels ν ( t, s ) , ν ( t, a ) , ν ( t, h ),conditioned on the durations ( I, D, E ), is reported in Fig. 3.
Remark 5.1.
We emphasize that our model of COVID-19 is targeted at pre-dicting the process of new infections, rather than the current number of peoplehosting the virus in various conditions. In particular, we do not explicitly de-scribe the dynamics of symptomatic but quarantined people and their exit process(i.e., recovery or, in the worst case, death), since such dynamics have no effecton the spreading process (under the assumption of perfect quarantine). How-ever, if desired, one could describe through appropriate probabilities and timedistributions how quarantined people split between those isolated at home and hose who get hospitalized, the fraction going to intensive care, and those whounfortunately die. Parameter symbol COVID-19 fitted valueincubation period I tri([2,12], mean 6)pre-symptoms period w 2disease period D unif([2,12])residual period E exp(10)low infection rate λ low λ high p s p a p h Table 1: Virus-specific parameters.
The parameters introduced so far, summarized in Table 1, are related tospecific characteristics of the virus. We now model properties of the specificenvironment where the virus spreads, taking into account the impact of coun-termeasures. First, we need to specify the immigration process σ ( t ). To keepthe model as simple as possible, we have assumed that the system starts ata given time with I new infections, i.e., immigrants arrive as a single burstconcentrated at one specific instant.The impact of countermeasures is taken into account in the model in twodifferent ways. First, modulating function µ ( t ) can be used to model the in-stantaneous reduction of the infection rate at time t due to the current mobilityrestrictions, and, more in general, changes in the environment which affect theability of the virus to propagate in the susceptible population (such as seasonaleffects). Typically, µ ( t ) is a decreasing function in the initial part of the epi-demics, in response to new regulations introduced by the government, while itgoes up again when mobility restrictions are progressively released.Second, we assume that any infected not already found to be positive istested (e.g., by means of massive swab campaigns) at individual, instantaneousrate ρ ( t ), which reflects both the amount of resources employed by the healthservice to discover the infected population, and the effectiveness of tracking.Recall that, when an infected is found positive, we assume it transits to thequarantined state q and stops infecting other people.Fig. 4 shows the resulting transitions that can occur for the different classes.Note that class q collects all cases known to the health authorities at a giventime t , and thus coincides with the set of detected cases. Infected in classes s, a, h are still unknown to the health service (but they are detectable). Whenthey stop to be positive, people in classes a, h transit to a class u (undetected),which collects all infected who remain unknown to the health service.Environment-related parameters are summarized in Table 2.13 ew infection after I+wafter I+w+D+Eafter I+w+D+E sa qh u p s p a ρ ( t ) ρ ( t ) ρ ( t ) p h Figure 4: State transitions of an infected individual.
Parameter symbol fitted value for first wave in Italyinitial day − ∆ t -20initial number of infected I µ ( t ) T a = − µ a = 3 . T b = 30, µ b = 0 .
31, (see Fig. 5)detection rate ρ ( t ) ρ ( t ) = 0 . t , (see Fig. 5) Table 2: Environment-specific parameters.
Remark 5.2.
We emphasize that the model described so far does not explicitlymodel spatial effects and sub-populations. As such, it is more suitable to de-scribe a homogeneous scenario, where its environmental parameters (includingmodulating functions µ ( t ) and detection rate ρ ( t ) ) can be reasonably assumedto apply to all individuals of the population irrespective of their spatial location.The natural candidate of such scenario is a single nation or just a sub-regionof it. By so doing, we can assume that common national regulations are ap-plied, as well as common health-care practices. However, our approach based onHawkes processes could be extended to incorporate spatial effects, by consideringspatio-temporal kernel functions ν ( x , t ) , where x is a spatial vector originatedat the point at which a new infection occurs. Another possibility would be toconsider a multivariate temporal Hawkes process, where each component rep-resents a homogeneous region, and the different regional processes can interactwith each other. .1. Computation of the real-time reproduction number R ( t )From our model, it is possible to compute in a native way the average numberof infections caused by an individual who gets infected at time t , i.e., the real-time reproduction number R ( t ). To compute R ( t ), we condition on the duration x of the incubation time, on the duration y of the disease time, on the duration z of the residual time, and on the class assigned to the node getting infected attime t (note that they are all independent of each other): R ( t ) = (cid:90) x (cid:90) y (cid:90) z (cid:32) λ low (cid:90) x µ ( t + τ ) u ( t, τ ) d τ + p a λ low (cid:90) x + w + y + zx µ ( t + τ ) u ( t, τ ) d τ +( p s + p h ) λ high (cid:90) x + wx µ ( t + τ ) u ( t, τ ) d τ + p h λ high (cid:90) x + w + yx + w µ ( t + τ ) u ( t, τ ) d τ + p h λ low (cid:90) x + w + y + zx + w + y µ ( t + τ ) u ( t, τ ) d τ (cid:33) d F E ( z ) d F D ( y ) d F I ( x ) (18)where u ( t, τ ) = e − (cid:82) t + τt ρ ( s ) d s is the probability that a node which gets infected at time t is still undetectedat time t + τ .
6. Model fitting for the first wave of COVID-19 in Italy
We fit the model to real data related to the spread of COVID-19 in Italy,publicly available on GitHub [20]. Italy was the country where the epidemicsfirst spread outside of China into Europe, causing about 34600 deaths at theend of June 2020 during the first wave.Our main goal was to match the evolution of the number of detected cases,represented in the model by individuals in class q . The actual count is providedby the Italian government on a daily basis since February 24th 2020. We takethis date as our reference day zero. However, it is largely believed that theepidemics started well before the end of February. Indeed, it became soon clearthat detected cases were just the top of a much bigger iceberg, as the prevalenceof asymptomatic infection was initially largely unknown, which significantlycomplicated the first modeling efforts to forecast the epidemic evolution. DuringJune-July 2020, a blood-test campaign (aimed at detecting IgG antibodies)was conducted on a representative population of 64660 people to understandthe actual diffusion of the first wave of COVID-19 in Italy [21]. As a mainresult, it has been estimated that 1 482 000 people have been infected. Thisfigure provides a fundamental hint to properly fit our model, and indeed whileexploring the parameter space we decided to impose that the total number ofpredicted cases at the end of June (day 120) is roughly 1 500 000.For the durations I , w , D , E of the different stages of an infection, we havetried to follow estimates in the medical literature. In particular, the duration of15he incubation period is believed to range between 2 and 12 days, with a sort ofbell shape around about 5 days [22, 23, 24], which we have approximated, forsimplicity, by a (asymmetric) triangular distribution with support [2 ,
12] andmean 6. Moreover, we have fixed the duration of the pre-symptomatic phaseto 2 days. For the duration of the disease period we have taken a uniformdistribution on [2 , µ ( t ), we could arbitrarily normalize to 1 the value of λ high , while weset λ low = 0 . − ∆ t (recall that day zero is thefirst day of the trace) at which to start the process with I initially infectedindividuals. While different pairs (∆ t, I ) are essentially equivalent, we decidednot to start the process with too few cases and too much in advance with respectto day 0, to limit the variance of a single simulation run. We ended up setting∆ t = 20, while I = 370 was selected as explained later. ρ ( t ) Ta 0 Tb µµ ba −∆Τ t(t) µ Figure 5: Chosen profiles and parameters of functions µ ( t ) and ρ ( t ) for the first wave ofCOVID-19 in Italy. Mobility restrictions in Italy were progressively enforced by national lawsstarting a few days before day 0, first limited to red zones in Lombardy, andsoon extended to the entire country through a series of increasingly restrictiveregulations, introduced over the next 30 days. Instead of trying the capture thestep-wise nature of such restrictions, we assume µ ( t ) to take the simpler profiledepicted in Figure 5, i.e., a high initial value µ a before time T a , a low finalvalue µ b after time T b , and a linear segment connecting point ( T a , µ a ) to point( T b , µ b ). We decided to set T a = − T b = 30 to reflect the time window inwhich mobility restriction were introduced.In Fig. 5 we also show our choice for the profile of detection rate ρ ( t ), i.e., alinear increase starting from day 0, with coefficient α , ρ ( t ) = max { , αt } . This16rofile is justified by the fact that the first wave caught Italy totally unprepared( ρ = 0 before day zero), while massive swabs were only gradually deployed overtime after day zero.The critical parameters I , µ a , µ b , α were fitted by a minimum mean squareerror (MMSE) estimation technique based on the curve of detected cases in thetime window [0 − nu m b e r o f d e t ec t e d ca s e s d a il y v a r i a ti on time (days)detected cases (model)detected cases (trace)daily variation (model)daily variation (trace) Figure 6: Evolution of the number of detected cases according to model and real data. Cu-mulative number (left y axes) and its daily variation (right y axes).
Fig. 6 shows the final outcome of our fitting, comparing the evolution ofthe number of detected cases according to model and real data. We show boththe cumulative number (left y axes) and its daily increment (right y axes).Analytical predictions for the mean trajectory of detected cases were obtainedby averaging 100 simulation runs. The shaded region around analytical curves(barely visible) shows 95%-level confidence intervals computed for each of the120 days.We emphasize that a similar good match could be possibly obtained by othersets of parameters. Our purpose here was not to compute the best possible fitin the entire parameter space (which would be nearly impossible), but to showthat the model is rich enough to capture the behavior observed on the real traceafter a reasonable choice of most of its parameters, driven by their physicalmeaning.Fig. 7 shows instead the evolution of the real number of cases (both thecumulative number and its daily variation) according to the model only, in theabsence of data. Note however that we have constrained ourselves to obtain atotal number of about 1 500 000 cases on day 120, as suggested by the serological17 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 0 10000 20000 30000 40000 50000 60000 nu m b e r o f r ea l ca s e s d a il y v a r i a ti on time (days)real cases (model)daily variation (model) Figure 7: Evolution of the actual number of cases according to model only. Cumulativenumber (left y axes) and its daily variation (right y axes). test [21].Interestingly, by looking at the values of µ a and µ b computed by our MMSE,it appears that the national lockdown was able to reduce the spreading abilityof the virus within the Italian population by a factor of about 12 (from 3.71 to0.31). We will see later on in Fig. 12 that our fitting of µ ( t ) for the first waveis consistent with mobility trends estimated from the usage of the Apple mapsapplication [19].In Fig. 8 we show the real-time reproduction number computed by (18). Wealso applied the classic method of Wallinga-Teunis [26], implemented as in theR0 package [27], to estimate R ( t ) from the trace of detected cases. To applytheir method, we used a generation time obtained from a Gamma distributionwith mean 6.6 (shape 1.87, scale 0.28), which has been proposed for the firstwave of COVID-19 in Italy [28]. Interestingly, though both estimates of R ( t )exhibit qualitatively the same behavior, the model-based value of R ( t ), whichconsiders also undetected cases, tends to be smaller. Having fitted the model to the available trace, we proceed to exploit themodel to examine interesting what-if scenarios. First, we investigate what wouldhave happened (according to the fitted model) if lockdown restrictions wereshifted in time by an amount of days δ . This means that we keep all parame-ters the same, except that we translate horizontally in time the profile of µ ( t )depicted in Fig. 5. In Figure 9 we show the total number of cases that wouldhave occurred for values of δ ∈ {− , − , , , } .18 R t time (days) from modelfrom trace Figure 8: Real-time reproduction number computed by the model, compared to its estimationbased on the trace of detected case, according to the Wallinga-Teunis method (with 95%confidence interval denoted as shared region).
We observe that a shift of just 3 days corresponds to a factor of roughly 2in the number of cases. This translates, dramatically, into an equivalent impacton the number of deaths, if we assume that the mortality rate would havestayed the same (i.e., 34600 / ∼ . ρ ( t ), bychanging its slope α (recall Fig. 5). We report both the number of real casesand the number of detected cases predicted by the model, when all other fittedparameters are kept the same. We consider what would have occurred with α = 0 (which means that infectious people are never tested), and by doublingthe intensity of the detection rate (a tracking system twice more efficient). Asexpected, with α = 0 only symptomatic cases ( p s = 6%) are eventually detected.This time the effect on the final number of cases (or deaths) is not as dramatic asin the previous what-if scenario. This suggests that the impact of mass testingin Italy during the first wave was marginal, and doubling the efforts would nothave produced significant changes in the final outcome. This is somehow optimistic, since mortality also depends on the saturation level of inten-sive therapy facilities.
20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 0 500000 1x10 nu m b e r o f r ea l ca s e s time (days)lockdown +7lockdown +3actual lockdownlockdown -3lockdown -7 Figure 9: What-if scenario: total number of cases if restrictions were shifted in time by δ days.
7. Model fitting for the second wave of COVID-19 in Italy
The second wave of COVID-19 hit Italy in late summer 2020, as in manyother European countries, mainly as en effect of relaxed mobility in July/Augustand possibly other seasonal effects. It is interesting to compare the second wavewith the first wave, by looking at the daily increments of detected cases anddeaths, see Fig. 11.We observe that the number of daily deaths is similar between the two waves.The daily increment of detected cases, instead, is very different (around 5000at the peak of the first wave, around 35000 at the peak of the second wave, a7-fold increase). This can be explained by the much larger capacity of the healthservice to perform swabs and track down the infected, built on the experiencegained from the first wave. It also suggests that, differently from the first wave,the impact of ρ ( t ) (the individual rate at which an infectious is detected) isexpected to be much more important during the second wave than in the firstwave, requiring a careful treatment of it in the model.We have indeed tried to fit our model to the second wave in Italy, keeping allparameters related to the virus (previously fitted for the first wave) unchanged(Table 1), and adapting only environment-specific parameters (Table 2). Amajor difficulty that we had to face was the unknown (at the time of writing)actual diffusion of the virus in the population during the second wave. Recallthat, for the first wave, we exploited a blood-test campaign to get a referencefor the total number of real cases at the end of the first wave. For the secondwave we employed instead a different approach based on a projection back inthe past of the increment of deaths. A similar idea has been adopted in [29] to20
20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1x10 nu m b e r o f ca s e s time (days) α = 0fitted α fitted α x 2 real casesdetected cases Figure 10: What-if scenario: impact of detection rate ρ ( t ) = αt on the real and detectednumber of cases, varying α . estimate the time-varying reproduction number in different European countryat the onset of the first wave.Specifically, we got from [30] the indication of the median (11 days) and IQR(6-18 days) of the amount of time from symptoms onset to death during Oct-Dec 2020, that we fitted by a Gamma distribution (shape 1.65, scale 8.45). Byconvolving such Gamma distribution with the distribution of incubation time,and the pre-symptoms period, we obtained a distribution of the total time frominfection to death, that we used to estimate the time in the past at whicheach dead person was initially infected. By amplifying the number of infectionsleading to death by the inverse of the mortality rate we eventually obtained anestimate of the daily increment of real cases, as reported in Fig. 11 (green curve,left y -axes).We chose August 1st (day 160 on Fig. 11) as starting date of the newinfection process producing the second wave in Italy. This time, instead ofmanually searching for suitable profiles of µ ( t ) and ρ ( t ) generating the expectedcurves of real and detected cases, we adopted a novel “reverse-engineering”approach: we took the daily increments of real and detected cases as input tothe model, and computed the functions µ ( t ) and ρ ( t ) that would exactly producein the model the given numbers of real and detected cases, on each day . Discretized values of functions µ ( t ) and ρ ( t ) are uniquely determined in the model, oncewe constrain ourselves to produce a given number of real and detected cases at each timeinstant. d a il y i n c r e m e n t o f d e t ec t e d / r ea l ca s e s d a il y i n c r e m e n t o f d ea t h s time (days)daily increment of casesdaily increment of deathsdaily increment of detected Figure 11: Daily increment of the number of deaths (black curve, right y -axes), daily incrementof the number of detected (red curve, left y axes), and estimated daily increment of real casesprojected back from the curve of deaths (green curve, left y -axes). In Fig. 12 we show the obtained ‘reversed’ µ ( t ) (blue curve), starting fromday 0 (August 1st) of the time reference adopted for the second wave. We alsoreport the average mobility measured in Italy by the Apple maps application[19], where we have given equal weight (1/3) to driving, transit and walkingmobility. We observe that the ‘reversed’ µ ( t ) qualitatively follows the samebehavior of mobility measured on the maps application, with a gradual increaseduring the month of August followed by a gradual decrease as people (and thegovernment) started to react to the incipient second wave.For completeness, we have also reported on Fig. 12 the fitted µ ( t ) for thefirst wave (black curve). We observe that, during the first wave, the quick in-troduction of hard lockdown caused an abrupt decay of both measured mobilityand fitted µ ( t ), characterized by a bigger reduction in a shorter time. Thesecond wave, instead, was characterized by a smoother transition, due to thedifferent choice of applying just partial lockdowns and progressive restrictionsmore diluted over time .In Fig. 13 we show instead the ‘reversed’ function ρ ( t ), focusing on thesecond wave. We also report on the same plot the daily increments of realand detected cases, which allow us to better understand the obtained profile of ρ ( t ), highlighting an interesting phenomenon occurred around day 30. Indeed,we observe that, during the first month, the epidemic was closely tracked bythe national health system, but at some point, around day 30, the curve ofdetected cases stops increasing, and stays roughly constant during the entiresecond month, lagging more and more behind the otherwise exploding curve ofreal cases (time window 30-60). As consequence, function ρ ( t ), which describesthe effectiveness of individual tracking, falls down to a minimum reached ataround day 60. This behavior can be interpreted as an effect of the saturation of The fact that, during the second wave, mobility values similar to those of the first wavedo not translate into equally similar values of µ ( t ) can be attributed to increased awarenessof people during the second wave. µ ( t ) m ob ilit y ( A pp l e ) time (days) µ (t) - 1st wavereverse µ (t) - 2nd waveaverage mobility - Apple maps Figure 12: Modulating functions µ ( t ) for the first wave (black, fitted) and second wave (blue,reversed), and average mobility according to Apple maps application. the capacity to perform swabs, resulting in a progressive collapse of the trackingsystem, as actually experienced by many people during those days.Our analysis shows that, in contrast to what might be believed by justlooking at the curve of detected cases, September (days 30-60) was, perhaps, themost critical period for the outbreak of the second wave of COVID-19 in Italy.During this period, the detection capacity of the national health system wassaturated, and could not keep the pace with the rapid growth of real cases, givinginstead the illusion of maintaining the epidemic under control. Our reverse-engineering approach can thus shed some light on what actually happened at theonset of the second wave, and quantitatively assess the collapse of the trackingsystem.
8. Conclusions
We have proposed a time-modulated version of the Hawkes process to de-scribe the temporal evolution of an epidemic within an infinite population ofsusceptible individuals. Our approach allows us to take into account precisedistributions for the time spent in different stages of the infection, which isof paramount importance when the intensity of countermeasures (mobility re-strictions, testing and tracing) varies significantly on time-scales comparable tothat of an individual infection. We have applied the model to the spread of23 r e v e r s e ρ ( t ) d a il y i n c r e m e n t o f d e t ec t e d / r ea l ca s e s time (days)reverse ρ (t)daily increment of detected casesdaily increment of real cases Figure 13: Reverse detection rate ρ ( t ) for the second wave (left y-axes), and daily increase ofreal and detected cases (right y-axes). COVID-19 in Italy, either by a direct fit of its parameter (first wave), or bya novel reverse fit (second wave) which allows us, in retrospect, to understandfrom data the time-varying effectiveness of applied countermeasures. Futurework will extend the model to overcome some of its current limitations, like theimpact of spatial effects and other sources of heterogeneity in the population,such as age groups. We think the proposed approach is promising and could beusefully applied to explain the epidemic progress and forecast/assess the impactof control/mitigation measures.
Appendix A. Explicit solution of λ ( t ) for constant µ ( t ) = µ When µ ( · ) ≡ µ is constant, denoting with ∗ the convolution product, we canrewrite equation (4) as λ ( t ) = σ ( t ) + µ · ( ν ∗ λ )( t ) , t > , which can be easily solved in the transformed domain. For a non-negativefunction f : R + → R + , we denote by (cid:98) f ( s ) := (cid:90) R + e − st f ( t ) d t, s ∈ C the Laplace transform of f . Then we have: (cid:98) λ ( s ) = (cid:98) σ ( s )1 − µ (cid:98) ν ( s ) for Re( s ) > Re( z max ) , ν ( t ) λ ( t ) delta δ ( t − α ) Θ(1) e α log( µ ) t erl- α ) t e − αt Θ(1) e α ( √ µ − t exp α e − αt Θ(1) e α ( µ − t hyper z e − zz +1 αt +2e − z +1 αt ( z +1) Θ(1) e αt (cid:20) ( µ − z − z + √ ( µ z z − − µ )(1+ z )2 (cid:21) hyper z α e − zαt + αe − αz t z ( z +1) Θ(1) e αt ( µ − z − µz + √ µ z z − µ z − µ ( z z +1)+( z +1)2]2 z Table A.3: Dominant term of λ ( t ) for different kernel shapes having the same average gener-ation time g = 1 /α . with z max equal to the zero of 1 − µ (cid:98) ν ( s ) with largest real part.In addition, formally, as long as µ < σ is bounded, we can write λ ( t ) = σ ( t ) + σ ( t ) ∗ ∞ (cid:88) i =1 µ i ν ∗ i ( t )where ν ∗ i ( t ) is the i -th fold convolution of v ( t ).An analytical expression of λ ( t ) can be obtained when σ ( s ) and ν ( s ) are bothrational. Table A.3 reports the dominant term of λ ( t ) (and also of N ( t )) forthe case in which σ ( t ) = β exp( − βt ), and ν ( t ) takes different shapes satisfying: (cid:90) ∞ ν ( t ) d t = 1 ; (cid:90) ∞ t ν ( t ) d t = 1 α Besides the simple deterministic ( delta ) and exponential ( exp ) shapes, weconsider the Erlang-2 ( erl-
2) and two variants of hyper-exponential, whose gen-eral form is ν ( t ) = p α e − α t + p α e − α t . To reduce the degrees of freedomof the general hyper-exponential, we have assumed a particular relationshipbetween p /p and α /α , which allows us to introduce a single parameter z .Specifically, we have set either p /p = α /α = z (denoted by hyper ) or p /p = (cid:112) α /α = z (denoted by hyper ). Note that in both cases the vari-ance of the corresponding hyper-exponential distribution increases with z ≥ z → ∞ .For all kernel shapes, and µ >
1, the growth rate of λ ( t ) is exponential, i.e., λ ( t ) = Θ( e ηt ), and the same holds for N ( t ) = (cid:82) t λ ( τ ) d τ . The rate η of theexponential growth, however, depends on the specific shape.In particular, for the two considered cases of hyper-exponential shape, η isan increasing function of z : while in the hyper case η saturates to 2( µ − α ,as z → ∞ , in the hyper case η is even unbounded, as z → ∞ .We also notice that the rate of the exponential kernel is larger than the rateof the Erlang-2 kernel, which is in turn larger than the rate of the deterministickernel. We recall that a non negative function f ( t ) is said to be Θ(1), iff there exist two constants0 < c ≤ C < ∞ such that c ≤ f ( t ) ≤ C , for sufficiently large t . µ ( t ),the epidemic growth rate depends on the specific shape of the kernel function ν ( t ). References [1] Y. Fang, Y. Nie, M. Penny, Transmission dynamics of the COVID-19 out-break and effectiveness of government interventions: A data-driven analy-sis, Journal of medical virology 92 (6) (2020) 645–659, 32141624[pmid].URL https://pubmed.ncbi.nlm.nih.gov/32141624 [2] A. J. Kucharski, T. W. Russell, C. Diamond, Y. Liu, J. Edmunds,S. Funk, R. M. Eggo, F. Sun, M. Jit, J. D. Munday, N. Davies, A. Gimma,K. van Zandvoort, H. Gibbs, J. Hellewell, C. I. Jarvis, S. Clifford, B. J.Quilty, N. I. Bosse, S. Abbott, P. Klepac, S. Flasche, Early dynamics oftransmission and control of COVID-19: a mathematical modelling study,The Lancet Infectious Diseases 20 (5) (2020) 553 – 558.URL [3] S. Pei, J. Shaman, Initial Simulation of SARS-CoV2 Spread and Interven-tion Effects in the Continental US, medRxiv (2020).URL [4] D. Zou, L. Wang, P. Xu, J. Chen, W. Zhang, Q. Gu, Epidemic ModelGuided Machine Learning for COVID-19 Forecasts in the United States,medRxiv (2020).URL [5] A. C. Miller, N. J. Foti, J. A. Lewnard, N. P. Jewell, C. Guestrin, E. B.Fox, Mobility trends provide a leading indicator of changes in SARS-CoV-2transmission, medRxiv (2020).URL [6] C. Anastassopoulou, L. Russo, A. Tsakris, C. Siettos, Data-based analysis,modelling and forecasting of the COVID-19 outbreak, PLOS ONE 15 (3)(2020) 1–21.URL https://doi.org/10.1371/journal.pone.0230405 [7] G. C. Calafiore, C. Novara, C. Possieri, A time-varying SIRD model forthe COVID-19 contagion in Italy, Annual Reviews in Control 50 (2020)361 – 372.URL [9] G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, A. Di Mat-teo, M. Colaneri, Modelling the COVID-19 epidemic and implementationof population-wide interventions in Italy, Nature Medicine 26 (6) (2020)855–860.URL https://doi.org/10.1038/s41591-020-0883-7 [10] A. G. Hawkes, Spectra of Some Self-Exciting and Mutually Exciting PointProcesses, Biometrika 58 (1) (1971) 83–90.URL [11] A. D. Kimmel M., The Bellman-Harris Process, Branching Processes inBiology, Interdisciplinary Applied Mathematics, Springer, New York, NY19 (2002) 87–102. doi:https://doi.org/10.1007/0-387-21639-1_5 .[12] H. Mei, J. Eisner, The Neural Hawkes Process: A Neurally Self-ModulatingMultivariate Point Process, in: Proceedings of the 31st International Con-ference on Neural Information Processing Systems, NIPS’17, Curran Asso-ciates Inc., Red Hook, NY, USA, 2017, p. 6757–6767.[13] M. Kim, D. Paini, R. Jurdak, Modeling stochastic processes in diseasespread across a heterogeneous social system, Proceedings of the NationalAcademy of Sciences 116 (2) (2019) 401–406.URL [14] H. Xu, H. Zha, A Dirichlet Mixture Model of Hawkes Processes for EventSequence Clustering, in: I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach,R. Fergus, S. Vishwanathan, R. Garnett (Eds.), Advances in NeuralInformation Processing Systems, Vol. 30, Curran Associates, Inc., 2017,pp. 1354–1363.URL https://proceedings.neurips.cc/paper/2017/file/dd8eb9f23fbd362da0e3f4e70b878c16-Paper.pdf [15] M. Akian, L. Ganassali, S. Gaubert, L. Massouli´e, Probabilistic and mean-field model of COVID-19 epidemics with user mobility and contact tracing(2020). arXiv:2009.05304 .[16] W.-H. Chiang, X. Liu, G. Mohler, Hawkes process modeling of COVID-19with mobility leading indicators and spatial covariates, medRxiv (2020).URL [17] M.-A. Rizoiu, S. Mishra, Q. Kong, M. Carman, L. Xie, SIR-Hawkes: Link-ing Epidemic Models and Hawkes Processes to Model Diffusions in FinitePopulations, 2018, pp. 419–428. doi:10.1145/3178876.3186108 .2718] Google LLC, Google COVID-19 Community Mobility Reports, .[19] Apple Inc., Apple COVID-19 Mobility Trends Reports, .[20] Presidenza del Consiglio dei Ministri, Dipartimento della Protezione Civile:Italian survelliance data, https://github.com/pcm-dpc/COVID-19 .[21] Istat, Istituto Nazionale di Statistica, Primi risultati dell’indagine disieroprevalenza sul SARS-CoV-2, .[22] The Incubation Period of Coronavirus Disease 2019 (COVID-19) FromPublicly Reported Confirmed Cases: Estimation and Application, Annalsof Internal Medicine 172 (9) (2020) 577–582, pMID: 32150748.URL https://doi.org/10.7326/M20-0504 [23] N. M. Linton, T. Kobayashi, Y. Yang, K. Hayashi, A. R. Akhmetzhanov, S.-m. Jung, B. Yuan, R. Kinoshita, H. Nishiura, Incubation Period and OtherEpidemiological Characteristics of 2019 Novel Coronavirus Infections withRight Truncation: A Statistical Analysis of Publicly Available Case Data,Journal of Clinical Medicine 9 (2) (2020).URL [24] J. Qin, C. You, Q. Lin, T. Hu, S. Yu, X.-H. Zhou, Estimation of incubationperiod distribution of COVID-19 using disease onset forward time: Anovel cross-sectional and forward follow-up study, Science Advances 6 (33)(2020).URL https://advances.sciencemag.org/content/6/33/eabc1202.full.pdf [25] A. W. Byrne, D. McEvoy, A. B. Collins, K. Hunt, M. Casey, A. Barber,F. Butler, J. Griffin, E. A. Lane, C. McAloon, K. O’Brien, P. Wall, K. A.Walsh, S. J. More, Inferred duration of infectious period of SARS-CoV-2:rapid scoping review and analysis of available evidence for asymptomaticand symptomatic COVID-19 cases, BMJ open 10 (8) (2020) e039856–e039856, 32759252[pmid].URL https://pubmed.ncbi.nlm.nih.gov/32759252 [26] J. Wallinga, P. Teunis, Different Epidemic Curves for Severe Acute Respi-ratory Syndrome Reveal Similar Impacts of Control Measures, AmericanJournal of Epidemiology 160 (6) (2004) 509–516.URL https://academic.oup.com/aje/article-pdf/160/6/509/179728/kwh255.pdf [27] T. Obadia, R. Haneef, P.-Y. Bo¨elle, The R0 package: a toolbox to estimatereproduction numbers for epidemic outbreaks, BMC Medical Informaticsand Decision Making 12 (1) (2012) 147. doi:10.1186/1472-6947-12-147 .2828] D. Cereda, M. Tirani, F. Rovida, V. Demicheli, M. Ajelli, P. Poletti,F. Trentini, G. Guzzetta, V. Marziano, A. Barone, M. Magoni, S. Dean-drea, G. Diurno, M. Lombardo, M. Faccini, A. Pan, R. Bruno, E. Pariani,G. Grasselli, A. Piatti, M. Gramegna, F. Baldanti, A. Melegaro, S. Mer-ler, The early phase of the COVID-19 outbreak in Lombardy, Italy (2020). arXiv:2003.09320 .[29] S. Flaxman, S. Mishra, A. Gandy, H. J. T. Unwin, T. A. Mellan, H. Cou-pland, C. Whittaker, H. Zhu, T. Berah, J. W. Eaton, M. Monod, P. N.Perez-Guzman, N. Schmit, L. Cilloni, K. E. C. Ainslie, M. Baguelin,A. Boonyasiri, O. Boyd, L. Cattarino, L. V. Cooper, Z. Cucunub´a,G. Cuomo-Dannenburg, A. Dighe, B. Djaafara, I. Dorigatti, S. L. van El-sland, R. G. FitzJohn, K. A. M. Gaythorpe, L. Geidelberg, N. C. Grassly,W. D. Green, T. Hallett, A. Hamlet, W. Hinsley, B. Jeffrey, E. Knock, D. J.Laydon, G. Nedjati-Gilani, P. Nouvellet, K. V. Parag, I. Siveroni, H. A.Thompson, R. Verity, E. Volz, C. E. Walters, H. Wang, Y. Wang, O. J.Watson, P. Winskill, X. Xi, P. G. T. Walker, A. C. Ghani, C. A. Donnelly,S. Riley, M. A. C. Vollmer, N. M. Ferguson, L. C. Okell, S. Bhatt, I. C.C.-. R. Team, Estimating the effects of non-pharmaceutical interventionson COVID-19 in Europe, Nature 584 (7820) (2020) 257–261.URL https://doi.org/10.1038/s41586-020-2405-7 [30] Istituto Superiore di Sanit`a, Characteristics of COVID-19 patientsdying in Italy,