A new SEIR type model including quarantine effects and its application to analysis of Covid-19 pandemia in Poland in March-April 2020
AA new SEIR type model including quarantine effects and its application to analysis ofCovid-19 pandemia in Poland in March-April 2020
Tomasz Piasecki , Piotr B. Mucha , Magdalena Rosi´nska
1. Institute of Applied Mathematics and Mechanics, University of Warsaw,ul. Banacha 2, 02-097 Warszawa, Poland.E-mails: [email protected], [email protected]. Department of infectious Disease epidemiology and Surveillance,National Institute of Public Health - National Institute of Hygiene,ul. Chocimska 24, Warsaw, Poland.E-mail: [email protected]
Abstract
Contact tracing and quarantine are well established non-pharmaceutical epidemic controltools. The paper aims to clarify the impact of these measures in COVID-19 epidemic. A newdeterministic model is introduced (SEIRQ: susceptible, exposed, infectious, removed, quaran-tined) with Q compartment capturing individuals and releasing them with delay. We obtaina simple rule defining the reproduction number R in terms of quarantine parameters, ratio ofdiagnosed cases and transmission parameters. The model is applied to the epidemic in Polandin March - April 2020, when social distancing measures were in place. We investigate 3 scenarioscorresponding to different ratios of diagnosed cases. Our results show that depending on thescenario contact tracing could have prevented from 50% to over 90% of cases. The effects ofquarantine are limited by fraction of undiagnosed cases. Taking into account the transmissionintensity in Poland prior to introduction of social restrictions it is unlikely that the control ofthe epidemic could be achieved without any social distancing measures. MSC Classification: 92D30, 34D20Keywords: Covid-19, SEIR model, quarantine, epidemiology, ordinary differential equations, stabil-ity analysis
The epidemic of SARS-CoV-2 infection triggered an unprecedented public health response. Giventhe lack of effective vaccine and treatment in 2020, this response included a variety of travel re-strictions and social distancing measures [2]. While these measures help to slow down the epidemicthey come at significant economical and societal cost [1]. As an alternative an approach focusing onrapid diagnosis is increasingly recommended [21] and prior to lifting social distancing measures large-scale community testing should be in place [36]. Testing efforts are complemented by identifyingand quarantining contacts of the diagnosed cases. Of note, by isolating the asymptomatic contactsfrom their social networks, this strategy takes into account the pre-symptomatic and asymptomaticspread of the infection [9, 10], believed to be one of the key drivers of fast spread of COVID-19. Asan example, wide spread testing in general population followed by isolation of the infected helpedto reduce COVID-19 incidence by 90% in an Italian village of Vo’Euganeo [14]. A modelling studyin France offers similar conclusions arguing that relaxing social lock-down will be only feasible incase of extensive testing [22]. While there is already a number of studies estimating the effects ofgeneral social distancing measures [25, 2, 22, 26], less is known about the impact of quarantine.Hellewell et al. [27] investigated the potential of rapid isolation of cases and contact tracing tocontrol the epidemic, finding that prohibitively high levels of timely contact tracing are necessaryto achieve control. However, new technologies may offer sufficiently fast alternative to traditionalcontact tracing, in which case the epidemic could be still controlled by contact tracing [28].Our aim is to develop a SEIR-type model which incorporates the effects of quarantine andvalidate it in a setting in which measures to reduce contacts are in place. We apply it to investigate1 a r X i v : . [ q - b i o . P E ] M a y he role of quarantine in Poland. The first case of COVID-19 in Poland was diagnosed on March4th. Social distancing measures were rapidly introduced during the week of 9 - 13th March includingclosure of schools and universities, cancellation of mass events and closure of recreation facilitiessuch as bars, restaurants, gyms etc. as well as shopping malls. Religious gatherings were limited.Finally, borders were closed for non-citizens [24]. These measures were fully in place on March16th. Further, beginning at March 25th restrictions on movement and travel were introduced (lock-down). Wearing face covers became obligatory on April 14th. The restrictions were gradually liftedbeginning at April 20th. We focus on modelling the time period when the social distancing measureswere in place and then consider different scenarios of relaxation of the restrictions with possibleimprovement of testing and contact tracing. We note that the procedures for quarantine were inplace even before the social distancing measures. They initially focused on individuals arriving fromCOVID-19 affected areas in China. When the epidemic started spreading in European countriespeople who came back to Poland from these countries were advised to immediately seek medicalattention if they experienced any symptoms consistent with COVID-19. However, adherence tothese recommendations was not evaluated. As soon as the first case was diagnosed in Polandquarantine for close contacts was also implemented.This paper aims to define a deterministic population model describing the epidemic in classicalterms of susceptible, exposed, infectious, removed. In our model the quarantine becomes a separatestate that removes individuals from susceptible and exposed states. We show that the reproductivenumber in our model is given by a simple formula referring to the parameters of transmissionand transition, but also to parameters describing the quarantine. We demonstrate that in a reallife scenario (case study of Poland) the quarantine effectively reduces the growth of infectiouscompartment. Increasing the efficiency of contact tracing and testing may may to some extentcompensate lifting up the social distancing restrictions. We introduce a modification of the classical SEIR model including effects of quarantine. To underlineimportance of that extension we call it SEIRQ. Formally the model is described by a system ofordinary differential equations with delay dedicated to the quarantine. 𝑆 𝑄 𝐸 𝐼 𝑑 𝐼 𝑢 𝑅 𝑑 𝑅 𝑢 𝛽 𝑑 𝐼 𝑑 + 𝛽 𝑢 𝐼 𝑢 (1 − 𝜅)𝜎 𝜅𝜎 𝛾 𝑢 𝛾 𝑑 (1 − 𝜃)𝐾 𝜃𝐾 (1 − 𝜃)𝐾 𝑇 𝜃𝐾 𝑇 Figure 2.1:
Schematic representation of the states included in the model . The solid linesrepresent the transition parameters and the dashed line indicate that the specific quantity is added.The following states are included in the model: S ( t ) – susceptible E ( t ) – exposed (infected, not infectious) I d ( t ) – infectious who will be diagnosed I u ( t ) – infectious who will not be diagnosed R d ( t ) – diagnosed and isolated 2 u ( t ) – spontaneously recovered without being diagnosed Q ( t ) – quarantinedThe figure 2.1 presents the schematic representation of the model. A susceptible individual(state S ), when becoming infected first moves to the state E , to model the initial period, when theinfected individual is not yet infectious. Next the cases progress to one of the infectious states I d , I u at the rates κσ and (1 − κ ) σ , respectively. In general, moving through the I d pathway concernsthese individuals who (independently of quarantine) would meet the testing criteria, as relevant tothe local testing policy, e.g. testing of people with noticeable symptoms. We shall emphasize thatfrom the point of view of analysis of spread of infection, the quantity I u shall be regarded rather asnot recognized infections, not necessarily asymptomatic or mild. With this interpretation the valueof κ can be influenced by intensity of testing.The creation of state E is via I d and I u with transmission rates β d and β u , respectively, nor-malized to the total population size N = S + E + I d + I u + R d + R u + Q , which is assumed to beconstant in time, births and deaths are neglected.The transition parameter σ is assumed identical for both groups, relating to the time betweeninfection and becoming infectious. The infectious individuals then move to the state R d , whichis the state of being diagnosed and isolated (and later recovered or deceased), with the rate γ d corresponding to the observed time between onset and diagnosis. On the other hand R u containspeople who spontaneously recovered with rate γ u . Our model includes an additional state of beingquarantined ( Q ). To mimic the situation of contact tracing, individuals can be put in quarantinefrom the state S (uninfected contacts) or the state E (infected contacts). These individuals stayin the quarantine for a predefined time period T . We assume that the number of people whowill be quarantined depend on the number of individuals who are diagnosed. An average numberof individuals quarantined per each diagnosed person is denoted as α . However, as the epidemicprogresses some of the contacts could be identified among people who were already infected, butwere not previously diagnosed, i.e. the state R u . We note that moving individuals between thestates Q and R u has no effect on the epidemic dynamics, therefore we assume that only individualsfrom S and E are quarantined and we reduce the average number of people put on quarantine bythe factor S ( t ) S ( t )+ R u ( t ) . Further, to acknowledge the capacity limits of the public health system toperform the contact tracing, we introduce a quantity K max , describing the maximum number ofpeople who can be put in quarantine during one time step.We also assume that among the quarantined a proportion θ is infected. After the quarantine,the infected part θK ( t − T ) goes to R d and the rest (1 − θ ) K ( t − T ) returns to S .Taking all of the above into account, the model is described with the following SEIRQ system:˙ S ( t ) = − S ( t )( β d I d ( t ) + β u I u ( t )) − (1 − θ ) K ( t ) + (1 − θ ) K ( t − T ) , ˙ E ( t ) = S ( t )( β d I d ( t ) + β u I u ( t )) − σE ( t ) − θK ( t )˙ I d ( t ) = κσE ( t ) − γ d I d ( t ) , ˙ I u ( t ) = (1 − κ ) σE ( t ) − γ u I u ( t ) , ˙ R d ( t ) = γ d I d ( t ) + θK ( t − T ) , ˙ R u ( t ) = γ u I u ( t ) , ˙ Q ( t ) = K ( t ) − K ( t − T ) , where K ( t ) = min { S ( t ) S ( t )+ R u ( t ) αγ d I d ( t ) , K max } , and α, β d , β u , γ d , γ u , θ, T ≥ . (2.1)We assume that the parameters α, β d , β u , θ , γ u and γ d depend on the country and time-specificpublic health interventions and may therefore change in time periods. Due to proper interpretationof the equation on E we require that β d ≥ θαγ d to ensure positiveness of E .3 .2 Basic reproductive number, critical transmission parameter β ∗ . Based on the general theory of SEIR type models [23], we introduce the reproductive number R = κ (cid:18) β d γ d − θα (cid:19) + (1 − κ ) β u γ u . (2.2)It determines the stability of the system as R < R > κ ) and quarantine(in terms of α ), but also gives an indication on levels of optimal testing and contact tracing. Weunderline that this formula works for the case when the capacity of the contact tracing has not beenexceeded ( K ( t ) < K max ). The details of derivation of (2.2) are provided in the Appendix, sectionA.4. We shall emphasize the formal mathematical derivation holds for the case when I and E aresmall comparing to S . Therefore the complete dynamics of the nonlinear system (2.1) is not fullydetermined by (2.2). However in the regime of epidemic suppression, which is the case of COVID-19epidemic in Poland, I and E are small compared to S and so the formula (2.2) reasonably prescribesspreading of infection in the population.The critical value R = 1 defines the level of transmission which is admissible, taking into accountthe existing quarantine policy, in order to control epidemic. As the level of transmission depends onthe level of contacts, this provides information on the necessary level of social distancing measures.The formula (2.2) indicates that improving the contact tracing may compensate relaxation of contactrestrictions. The key quantity is θα . Indeed the system with the quarantine has the same stabilityproperties as one without K , but with the new transmission rate β newd = β d − θαγ d . In order toguarantee the positiveness of E , β newd must be nonnegative. It generates the constraint θαγ d ≤ β d . (2.3)The above condition also implies the theoretical maximal admissible level of quarantine. We defineit by improving the targeting of the quarantine, i.e. by the highest possible level of θ : θ max = β d γ d α . (2.4)As long as the K max threshold is not exceeded the effect of the increase in θ or in α play the samerole at the level of linearization (small I, E ). However, in general it is not the case and for thepurpose of our analysis we fix α .For our analysis we assume β d = β u = β . The reason is that, both I d and I u contain a mixture ofasymptomatic and symptomatic cases and although there might be a difference we lack informationto quantify this difference. Then using formula (2.2) we compute critical values β ∗ ( κ, θ, α ) definedas R ( β ∗ ) = 1 , namely β ∗ ( κ, θ, α, γ d , γ u ) = (1 + θακ ) γ d γ u γ u κ + γ d (1 − κ ) . (2.5)It shows the upper bound on transmission rate β which still guarantees the suppression of pandemic.We shall omit the dependence on γ d , γ u as these are fixed in our case, and denote briefly β ∗ ( κ, θ, α ).In the case of maximal admissible quarantine (2.4) we obtain β ∗ ( θ max , κ ) = γ u − κ , (2.6)which can be regarded as theoretical upper bound for β if we assume ”optimal admissible” quaran-tine for fixed κ , for which the epidemic could be still controlled. It must be kept in mind thoughthat the condition (2.3) means that we are able to efficiently isolate all persons infected by everydiagnosed, therefore is unrealistic. The resulting β ∗ ( θ max , κ ) should be therefore considered as atheoretical limit for transmission rate. 4 .3 Fitting procedure All simulations are performed using GNU Octave ( ). Theunderlying tool for all computations is a direct finite difference solver with a 1 day time step.
Basic assumptions for data fitting.
We estimate the transmission rates β by fitting themodel predictions to the data on the cumulative number of confirmed cases. Since people withconfirmed diagnosis are efficiently isolated, they are immediately included into R d . Therefore, thequantity fitted to the data is R d ( t ).The crucial assumption behind our approach is that the parameter β changes twice duringthe period of analysis. The reason is that we can distinguish two important time points in thedevelopment of epidemic in Poland. The first is initial restrictions including school closure effectiveMarch 12, which was accompanied with restrictions on other social activities. As we do not takemigration into account in our model, we assume that the effect of border closing is reflected in β .The second turning point was a lockdown announced on March 25.For simplicity we comprise the effect of above measures in two jump changes in β in t ∈ { t , t } and choose t = 14 , t = 28. With t = 1 corresponding to March 3 it means small delay withrespect to the above dates which can be justified by the fact that new cases are reported with adelay of approximately 2 days. Choice of fixed parameters (Tab. 2.1).
We assume that the parameters σ, γ u represent thenatural course of infection and their values could be based on the existing literature. The parameter σ describes the rate of transition from non-infectious incubation state E into the infectious states I d or I u . The value of σ takes into account the incubation period and presymptomatic infectivityperiod. γ u relates to the period of infectivity, which we select based on the research regarding mildercases, assuming that serious cases are likely diagnosed. Further, κ is a parameter related both to theproportion of asymptomatic infection and the local testing strategies. Since the literature findingsprovide different possible figures, for κ we examine three different scenarios.Parameters γ d , θ and α are fixed in our model for the purpose of data fitting, but informedby available data. One of the scenarios of future dynamics of the epidemic (section 3.3) considerspossible increase of θ . Parameter γ d was estimated basing on time from onset to diagnosis fordiagnosed cases, and θ as rate of diagnosed among quarantined. Furthermore we fix the parameter α by comparing the number of quarantined people obtained in simulations with actual data. Thecapacity level of public health services is set in terms of possible number of quarantined per day K max , as double the level observed so far. Detailed justification of the values of fixed parameterscollected in the following table, is given in the Appendix, section A.2.Parameter Value Source σ . Literature: incubation time [3, 4, 5]+ presymptomatic spread [7, 9, 10] γ d . Observed data: appendix section A.2 γ u Literature: [11], WHO mission report from China κ { } Literature: proportion asymptomaticor undocumented [6, 13, 14, 15] θ α
75 Observed data: appendix section A.2 K max
50 000 2 × the maximum level observed so far (arbitrary decision)Table 2.1: Fixed parameters used in the model Optimization algorithm.
In order to fit the values β , β , β we use a standard gradient5escent algorithm. The error function is defined as mean square difference between the cumulativenumber of diagnoses and the R d ( t ) predicted from the model.For the initial values the error function is optimized only for a limited number of possibleconditions, as these mostly impact β , which is less relevant for future predictions. To estimateconfidence intervals we use a method of parametric bootstrap. The optimisation procedures aredescribed in the Appendix, section A.1, where we also show precise errors of data fitting. Dataset.
The data series contains cumulative number of confirmed cases of COVID-19 in Polandfrom March 3 (first confirmed case in Poland) till April 26, which amounts to 54 observations. Thedata are taken from official communications of the Ministry of Health. As explained in table 2.1and appendix (section A.2 additional data sources were used for choosing θ , α and γ d . In Table 3.1 we show estimated values of β i , where i = 1 , , R for the third time interval. Given the social distancingmeasures in place early April 2020, as well as the quarantine levels, the reproductive number wasbelow 1, independently of the value of κ , which relates to testing effectiveness. The figure 3.1 showsthe fit of the models assuming different levels of κ . Good fit is found for all three models althoughpredictions start to differ in the middle-term prognosis. κ = 0 . κ = 0 . κ = 0 . β β β R ( β , . ,
75) 0.817 0.802 0.772(0.651 , 0.977) (0.648 , 0.915) (0.569 , 0.874)Table 3.1: Estimated values of β i and values of R corresponding to the latest estimation periodwith 95% confidence intervalsWe proceed with predictions assuming that the restrictions are continued, i.e.keeping β = β (note that the estimated β is different for each κ ). We calculate the epidemic duration ( t max ), thepeak number of infected ( I maxd , I maxu ) and the final size of the epidemic ( R d ( t max ) , R u ( t max )). Inorder to show the influence of quarantine we compare the situation with quarantine, keeping thesame θ, α , to the situation without quarantine, setting αθ = 0. The results of the development ofthe epidemic during the first 120 days are shown on Figure 3.2.For κ = 0 . κ = 0 . κ = 0 . t = 40 leading to huge difference in the total time of epidemic and total number of cases. Thesevalues are summarized in the table 3.2. We note that given the epidemic state in the first halfof April 2020 for all values of κ the model predicts epidemic extinction both with quarantine andwithout quarantine. However, since the epidemic is very near to the endemic state, the predictedduration is very long, especially if no quarantine is applied.6 Figure 3.1: Results of model fit to cumulative diagnosed cases ( R d ) for κ = 0 . , . , . R u (panel B). Colouredshades correspond to 95% confidence intervals for the respective colour line. κ quarantine factors R d ( t max ) R u ( t max ) I maxd I maxu t max θ = 0 . , α = 75 31 85 1.9 10.8 450 θ, α = 0 44 175 2.1 12.5 8300.5 θ = 0 . , α = 75 29 20 1.9 2.7 330 θ, α = 0 1078 1078 5.1 9.2 32000.8 θ = 0 . , α = 75 24 4 1.9 0.7 230 θ, α = 0 6317 1579 10.6 47.6 1280Table 3.2: Duration of epidemic ( t max ) in days, the final values of R d and R u , in thousands( R d ( t max ) , R u ( t max )) and peak values of I d and I u , in thousands ( I maxd , I maxu ) according to quaran-tine and testing scenarios. 7 Figure 3.2: Predicted values of R d , R u (panels A – C) and I d , I u (panels D – E), as depending onthe value of κ and whether or not the quarantine is implemented. For t > β = β estimated foreach κ , with the same quarantine parameters or without quarantine at all. β ∗ for the current situation Using the formula (2.5) we can compute critical values β ∗ . In Table 3.3 we show the values of β ∗ ( κ, . ,
75) and for convenience recall also estimated values of β and R , listed already inTable 3.1. Moreover we compute β ∗ ( κ, , R for ourestimated values of β and the same γ d , γ u but without quarantine. Comparing the estimated valuesof β (table 3.3) for all cases of κ are only slightly below β ∗ . κ β β ∗ ( κ, . , R ( β , κ, . , β ∗ ( κ, , R ( β , κ, , β − θαγ d β ∗ and R ( β ) with quarantine ( i = 0 . , α = 75) and without quarantineEliminating the quarantine, for the estimated values of β , we have different situations dependingon the actual value of κ . In case κ = 0 .
2, so assuming that currently only 20% of infections arediagnosed, the low values of R are due to low β rather than the effect of quarantine (controlling8pidemic by social contact restrictions). In effect even if we remove the quarantine we have still R <
1, but very close to 1. On the other hand if κ = 0 . κ = 0 . β ,which corresponds to the situation of controlling the epidemic by extensive testing and quarantine.For these cases, if we remove the quarantine, we end up with R >
1. The quantity β − θαγ d represents effective transmission rate due to diagnosed cases. In particular it shows by how muchthe transmission could be reduced by improved contact tracing ( θα ) and faster diagnosis ( γ d ).These results confirm that the higher is the ratio of undiagnosed infections, the weaker is influ-ence of quarantine. In the next section we verify these results numerically. Our second goal is to simulate loosening of restrictions. In particular we want to verify numericallythe critical thresholds β ∗ listed in table 3.3. For this purpose we assume that at t = 60 we change β . For each value of κ we consider 3 scenarios: (a) Current level quarantine: i.e. quarantine parameters θ = 0 . , α = 75 are maintained; (b) No quarantine is applied starting from t = 60; (c) The maximal admissible quarantine is applied, meaning that θ max = βαγ d (see (2.2)). In thiscase α = 75. As long as the limit K max is not reached there is no difference whether we increase α or θ , the decisive parameter is αθ . Increasing α would lead to reaching K = K max earlier andhence worse outcomes.Figures 3.3-3.5 show the final values of R = R d + R u and time till the end of epidemic dependingon the value of β for t ≥
60. for above 3 scenarios and different values of κ . The theoretical valuesof β ∗ are shown by black lines.The results confirm that around β ∗ a rapid increase in the total number of infected occurs,coinciding with the peak total epidemic duration. Thus the numerical computations confirm thatthe critical β ∗ calculated for the linear approximation in the section 2.2 are adequate, with a smallbias towards lower values. Figure 3.3: Duration of epidemic and the final epidemic size as dependent on β , for κ = 0 . κ = 0 . Figure 3.4: Duration of epidemic and the final epidemic size as dependent on β , for κ = 0 . Figure 3.5: Duration of epidemic and the final epidemic size as dependent on β , for κ = 0 . β observed for κ = 0 . κ = 0 .
8, both in case θ = 0 .
006 and θ = θ max . The values of R d and R u before and after these qualitative changes are summarized in Table 3.4.A closer investigation for these values of β shows that in all 4 cases the jump occurs for thefirst value of β for which the limit number of quarantined, K max = 50000, is achieved. Notice thatimmediately after passing the threshold the values become very close to those without quarantine.10 R d ( t max ) R u ( t max ) t max κ = 0 . , θ = 0 . κ = 0 . , θ = θ max ( β ), 0.211 1423 666 48000.212 11458 10971 840 κ = 0 . , θ = 0 .
006 0.218 1137 236 47400.219 13706 3365 1060 κ = 0 . , θ = θ max ( β ) 0.566 1762 108 28500.567 27602 6729 570Table 3.4: Critical values of β obtained in simulations and corresponding final numbers of diag-nosed/undiagnosed (in thousands) and total time of epidemic.Therefore the effect of quarantine is immediately and almost completely cancelled after passing thecritical value of β . The transition is milder in the case κ = 0 . β .Results of our simulations confirm the theoretical prediction that the margin in relaxation ofrestrictions is very narrow if we want to avoid a blow up of the number infections. Strengtheningof quarantine allows to remain in a stable regime while increasing β . We propose a simple SEIR-type model (SEIRQ), which includes the effects of testing and contacttracing. The model formulation allows to calculate an interpretable formula for the reproductivenumber R (2.2). As typical for this class of models, R depends on transmission parameters β .Increasing β corresponding in e.g. to higher frequency of social contacts increases R . Decreasing β , for example in consequence of widespread use of face masks, has the opposite effect. On theother hand γ d reflects the time to diagnosis and the formula indicates that more rapid diagnosisis associated with lower R . In addition, our model offers a clear interpretation of the quarantineeffect. The transmission rate due to diagnosed cases, β d , is decreased by the factor θαγ d indicatingthat both the number of quarantined per diagnosed individual ( α ) and proper targeting of thequarantine (the infection rate among the quarantined θ ) equally contribute to this factor. Also theparameter related to testing: the delay in diagnosis, γ − d , plays similar role. This quantifies thepotential of a wide range of interventions to improve testing and contact tracing, as outlined in e.g.in ECDC recommendations [31]. In particular, as the number of people put in quarantine per eachcase and the infection rate among the quarantined impact R in similar fashion, our results supportthe recommendations to focus on the high risk contacts when the resources do not allow to followall contacts.Our model takes into consideration only the effective contact tracing, i.e. the situation when theinfected contacts are identified and put in quarantine before they become infectious. People whoare identified later would be modelled as passing through one of the I states to the R states. Thismeans that the number of quarantined in our model can be also increased by faster contact tracing.The timely identification of contacts may be a significant challenge in the quarantine approachgiven that the incubation time can be as short as 2 days in 25% of cases [3]. As mentioned byother Authors [28], the delays in manual contact tracing are usually at least 3 days and under suchcircumstances the contact tracing and quarantine alone may be insufficient to control the epidemic.This could be improved with digital contact tracing. Notably, mixed contact tracing strategiesimplemented in South Korea indeed helped to control the epidemic at the early stages [32]. Theuse of ”smart contact tracing” with mobile phone location data and administrative databases werealso key to rapid identification and self-quarantine of contacts in Taiwan [33] and implementation of11uch strategy helped Singapore to control the epidemic without major disruptions of social activities[34].We note that the quarantine effect relates only to transmission due to diagnosed cases. Asexpected, in order to control the epidemic the transmission due to undiagnosed cases has to benegligible. This can be controlled by general measures such as lockdown , which universally decreasethe frequency of social contacts and are therefore likely to reduce β u . In our model the part of R representing transmission due to undiagnosed cases is scaled by (1 − κ ), the parameter relating tothe efficiency of the testing system. Again, the examples of Singapore as well as the Italian villageof Vo’Euganeo show that the widespread testing complementing the efficient contact tracing wasessential to suppress epidemic. Testing unrelated to epidemiological links decreases (1 − κ ) factor,thus making the factors impacting transmission due to diagnosed cases, such as quarantine, morepowerful to decrease R .Further, our model allows to study the effect of the situation, in which the contact tracingcapacities are exceeded. In this situation the epidemic is likely to quickly develop to the levelsobserved without quarantine. It is therefore quite crucial to implement the aggressive contacttracing system, when the epidemic is still at low levels and it is possible to bring the epidemic tosuppression phase.We demonstrate the high impact of contact tracing and quarantine on the observed numbers ofcases in Poland. This effect was coupled with substantial reduction in the transmission parameter β resulting from social contact restrictions. Depending on the scenario, β decreased by 76% to 84%,bringing R below 1. The estimated effect of the quarantine in Poland would depend on which ofthe considered scenarios regarding testing efficiency was the most relevant to our situation. In ourmodel the quarantine is estimated to be the most effective for the scenario in which most of the casesare diagnosed ( κ = 0 . κ of the order of 0.8. The Polish clinical recommendations specifically mention only testingall individuals with severe infections [35]. In addition testing is provided to health care workers.The severe course corresponds to approximately 18% of all infections [3]. Therefore, the κ = 0 . κ in our countrylies between 0.2 and 0.5. For these scenarios the model shows that the control of the epidemic islargely achieved through suppression of β . In case of relaxation of social contact restrictions, theefforts should be focused on increasing the level of testing in order to decrease the proportion ofundiagnosed cases as well as maintaining or increasing the effectiveness of quarantine. For smaller κ , even substantially increasing the effectiveness of quarantine does not allow to go back to the levelof social contacts from before the epidemic ( β ).Finally, the contact tracing effort was manageable in Poland due to relatively small number ofcases. Should the case load increase substantially longer delays in contact tracing would occur,which can substantially decrease the effects of quarantine [27, 28]. Limitations and future directions of research.
We do not consider the likely reducedtransmission from undiagnosed cases who are more likely to be asymptomatic or paucisymtopmaticcases ( β u < β d ). The reduction factor for infectiousness of asymptomatic is still under investigation.One study found a 60-fold lower viral loads in asymptomatic cases [29], but another estimated thetransmissibility reduction by 50% [6]. Moreover, we did not have sufficient data to include thisadditional parameter. We calibrated our model only to diagnosed cases, which were officially avail-able. Calibration to mortality data is another approach successfully implemented in e.g. [25] thatpotentially removes bias due to different testing policies. As there were relatively fewer fatalities inPoland and little data on clinical progression we decided on simplified model without explicit mod-elling of the outcomes. Furthermore, we did not consider the sub-optimal adherence to quarantine.It is likely that some individuals would not fully comply to strict quarantine rules. However, onlyanecdotal evidence for such phenomenon is available at this time. In our model it would decreasethe effective αθ , which was chosen to fit to observed number of people put in quarantine. Finally,the analysis of R is suitable for small size of epidemic, when S ≈ N . For other cases the results are12till useful, but the approximation may be biased, as we have shown for β ∗ . Due to little availabledata and policy changes we did not have sufficient data to determine which κ scenario is the mostappropriate.In conclusion we present a simple model, which allows to understand the effects of testing,contact tracing and quarantining of the contacts. We apply the model to the data in Poland andwe show that despite a substantial impact of contact tracing and quarantine, it is unlikely that thecontrol of the epidemic could be achieved without any reduction of social contacts. Acknowledgments.
This work was partially supported by the Polish National Science Centre’sgrant No2018/30/M/ST1/00340 (HARMONIA).
References
A Appendix
A.1 Optimization algorithm.
In order to fit the values β , β , β we use a standard gradient descent algorithm. Namely, we defineerror function as E ( κ ) = (cid:34) (cid:80) t =1 [ R d ( κ, t ) − data ( t )] (cid:80) t =1 | data ( t ) | (cid:35) / , (A.1)where { R d ( κ, t ) } t =1 is the vector of computed values of R d and { data ( t ) } t =1 the vector of data(cumulative number of confirmed cases).At each step we approximate the gradient of the error function with respect to β , β , β bydifferential quotients and move in the direction opposite to the gradient. The algorithm revealsa good performance provided we start sufficiently close to the minimum, which is not difficult toensure in our case.It remains to choose the initial data. A closer look on results of simulations shows that thechoice of initial data mostly influence the fitting in the beginning of period under considerationand hence the value of β , while for analysis of future scenarios β is the most important. Takingall this into account we do not struggle for sharp optimization of data fitting with respect toinitial data and restrict to the following heuristic choice. It is natural to assume I u (0) = − κκ I d (0).Concerning the choice of E (0) we assume it in a form E (0) = m ( I d (0) + I u (0)). We set initialvalues I d (0) ∈ { , , } and for each value we set I u (0) according to the above formula andthree values of E (0) corresponding to m ∈ { , , } . For each of these 9 combinations we run theoptimization algorithm looking for the best fit of β i , i = 1 . . .
3. We have repeated this approach for κ ∈ { . , . , . } . It turns in that for all values of κ the best fit was obtained for I d (0) = 20 and m = 2. More careful analysis around I d (0) = 20 did not improve the quality of fitting, therefore: I d (0) = 20 , I u (0) = 1 − κκ I d (0) , E (0) = 2( I d (0) + I u (0))is our final choice. We obtain the following fitting error defined by (A.1): κ E ( κ ) 0.0077 0.0079 0.0084Table A.1: Errors of data fit. A.2 Choice of fixed parameters
1. The parameter σ describes the rate of transition from non-infectious incubation state E into theinfectious states I d or I u . The median incubation time from exposure till the onset of symptoms wasestimated at 4 to 5 days [3, 4, 5]. However, there exists evidence that typically infectivity preceedssymptoms, by 1 to 3 days [7, 9, 10].A modelling study identified the rate of transition between the non-infectious and infectiousstates at . [6], which corresponds to an average time lag of 3.69 days untill the case becomesinfectious.2. The parameter γ u represents the period of infectivity during the natural course of disease.We discuss the period of infectivity, especially as applied to mild cases. The median duration of16iral shedding was estimated among 113 Chinese hospitalized patients. Overall it was 17 days, butit was shorter among cases with milder clinical course [16]. A study among 23 patients in HongKong confirmed viral shedding longer than 20 days among a third of patients, although the peaklevel of shedding was noted during the first week of infection [17]. In the mission report from ChinaWHO reports viral shedding in mild and moderate cases to last 7 - 12 days from symptom onset.Among younger and asymptomatic or mild cases the shedding may be shorter: in a study among24 initially asymptomatic youngsters the median duration was 9.5 days [11].3. The value of κ generally depends on the testing policy. However, recommended testingpolicies often rely on the presence of respiratory symptoms. This is also the case in Poland. Itwas observed that some infected people never develop symptoms, although the precised rate of suchtruly asymptomatic infections is still under investigation. Some studies may be biased by a too shortfollow-up time. A small study among residents of a long-term care skilled nursing facility found thateven though more than half of individuals with confirmed infection were asymptomatic at the timeof test, majority of them subsequently developed symptoms. The proportion of people who remainasymptomatic may be higher among younger individuals [Hu]. A study among Japanese nationalsrepatriated from Wuhan suggests the proportion of asymptomatic infections is about 30% [13]. Ananalysis among the passengers of Diamond Princess ship, where a COVID-19 outbreak occurred,taking into account this delayed onset of symptoms estimated the proportion of asymptomaticinfections to be about 18%, even though almost 50% were asymptomatic on initial test. In addition,large scale screening implemented in Italian village Vo’Euganeo indicated that 50% to 75% ofinfected individuals did not report symptoms [14]. Similarly, in population screening in Iceland50% were asymptomatic at the time of screening [15]. It may be stipulated that some of the peoplediagnosed through screening developed symptoms latter, consistently with the findings from theDiamond Princess study.On the other hand a sizable proportion of infected people, especially at younger ages, experienceonly mild symptoms, for which they may not seek medical attention. In the study of Li [6], theproportion of undocumented cases was estimated as 86%.4. The parameter γ d was estimated basing on a sample of case-based data available in routinesurveillance, by fitting gamma distribution to the time from onset to diagnosis, for cases whowere not in quarantine before diagnosis. Time from onset to diagnosis was estimated based onsurveillance data available in the Epidemiological Reports Registration System for COVID-19, as of28.04.2020.The system collects epidemiological data on cases diagnosed in Poland and is operatedby local public health departments. All cases eventually are entered into the database. However,substantial reporting delays are noted. There were altogether 4976 cases registered in the system,including 1995 (40.1%), who did not have symptoms at the time of diagnosis. Plausible onsetdate and plausible diagnosis date were available for 2884 cases ( 96.7% of 2981 cases that were notasymptomatic)Gamma distribution was fitted by maximum likelihood to cases who were not diagnosed inquarantine. The observed and fitted distributions are shown below (figure A.1).We next fitted gamma-regression model with week of diagnosis as an explanatory variable. Wefound no significant trend in time. We therefore adopted the average time from onset to diagnosisto be 4.6 days, and taking into account the probability of asymptomatic spread we assumed theparameter γ d to be 1/5.5.5. Next we base θ on available data. We calculate prevalence of infection among the quarantinedindividuals, according to data published by the Chief Sanitary Inspectorate on the number of casesdiagnosed among quarantined people and the total number of quarantined. We used a series ofdata 8.04 – 20.04 to estimate a likely value of θ . We chose this time period due to data availability.Data are shown on the figure below. During this time period there was an increasing trend in theproportion of diagnosed from 0.5% to 0.8% A.2. We presume that this parameter could change withchanging procedures of contact tracing and testing. However, since no detailed data were available,for the modelling purposes we chose a simplifying assumption that θ is stable (i.e. we always takea similar group of contacts under quarantine) selecting an average value of 0.6%.This proportion could be also viewed as attack rate among the contacts of cases. The proportionin Poland is in line with what was observed in Korea, where an estimated attack rate was 0 . . . . . D en s i t y Fitted gamma distribution and observed distribution of the time between onset and diagnosisFitted parameters: Shape = 5.06 Scale = .914Excluding cases diagnosed in quarantine
Figure A.1: Observed and fitter distribution of time between onset and diagnosis among symp-tomatic cases, who where diagnosed outside of quarantine. The fitted gamma distribution has thefollowing parameters: shape = 0.91, scale = 5.06 P E R C E N T D I A G N O S E D W I T H C O V I D - C U M U L A T I V E N U M B E R O F Q U A R A N T I N E D Number ever under quarantine Proportion diagnosed with COVID-19 among quarantined
Figure A.2: Data on the population undergoing quarantine and the percent diagnosed with COVID-19 in this populationoverall and 7 .
56% among household contacts [32], although household attack rate was higher ( > α . Here we make another simplification assuming this18arameter to be constant. The main difficulty is a lack of precise data concerning the number ofnewly quarantined people per day, distinguishing between reasons of quarantine (travel related orcontact tracing relate). At the beginning of epidemic in Poland the average amount of quarantinedfollowing one diagnosed case was definitely higher. Moreover, people coming back from abroadwere subject to obligatory quarantine starting from March 16 and constituted a considerable partof quarantined in the second half of March and beginning of April. In particular, around 54 000Polish citizens staying abroad came back within a special program of charter flights operated byPolish Airlines which ended on April 5. We can assume that after this date the ratio of peoplecoming from abroad among all people subject to quarantine was negligible. As our model does nottake migration into account, we have to take into account only quarantine from contact. For abovereasons, for fitting α we restrict our analysis only to a period of 2 weeks of April. Assuming already θ = 0 .
006 we then choose α minimizing the square error between the number of quarantined fromthe data and computed K(t). This way we obtain α = 75.Taking these into consideration we set the values of parameters collected in table 2.1. A.3 Confidence intervals
Bootstrap.
To estimate confidence intervals we use a method of parametric bootstrap. We generate M = 200 sequences of perturbed data assuming that for each time t ∈ { , } the increment of R (i.e. daily number of new diagnoses) is a random number from Poisson distribution with mean valueequal to increment of observed data.Model parameters are also perturbed, see below. For each series of perturbed data we estimatethe values of β i and take estimated confidence intervals as appropriate quantiles of obtained sets.In order to estimate confidence intervals for R d ( t ), shown on panel A of figure 3.1, we proceedas follows. For each sequence of perturbed data we compute fitted R d ( t ). This way we obtain a setof curves { R ( k ) d ( t ) } t =1 ,..., k =1 ,..., . Then for each time instant t ∈ { , } we set lower and upper bounds of the confidence intervalof R d ( t ) as appropriate quantiles of the set { R ( k ) d ( t ) } k =1 . Analogously we compute the confidenceintervals for R u ( t ) shown on panel B of figure 3.1. Distribution of parameters.
Following other Authors [30] as well as experimental data, forthe uncertainty analysis we used the following distributions of the parameters.1. 1 /γ d ∼ Gamma ( a , b ), where the shape parameter, a = 1 .