The epidemic of SARS-CoV-2 infection triggered an unprecedented public health response. Given the lack of effective vaccine and treatment in 2020, this response included a variety of travel restrictions and social distancing measures [2] . While these measures help to slow down the epidemic they come at significant economical and societal cost [1] . As an alternative an approach focusing on rapid diagnosis is increasingly recommended [21] and prior to lifting social distancing measures largescale community testing should be in place [36] . Testing efforts are complemented by identifying and quarantining contacts of the diagnosed cases. Of note, by isolating the asymptomatic contacts from their social networks, this strategy takes into account the pre-symptomatic and asymptomatic spread of the infection [9, 10] , believed to be one of the key drivers of fast spread of COVID-19. As an example, wide spread testing in general population followed by isolation of the infected helped to reduce COVID-19 incidence by 90% in an Italian village of Vo'Euganeo [14] . A modelling study in France offers similar conclusions arguing that relaxing social lock-down will be only feasible in case of extensive testing [22] . While there is already a number of studies estimating the effects of general 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 to control the epidemic, finding that prohibitively high levels of timely contact tracing are necessary to achieve control. However, new technologies may offer sufficiently fast alternative to traditional contact 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 and validate it in a setting in which measures to reduce contacts are in place. We apply it to investigate the role of quarantine in Poland. The first case of COVID-19 in Poland was diagnosed on March 4th. Social distancing measures were rapidly introduced during the week of 9 -13th March including closure of schools and universities, cancellation of mass events and closure of recreation facilities such 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 March 16th. Further, beginning at March 25th restrictions on movement and travel were introduced (lockdown). Wearing face covers became obligatory on April 14th. The restrictions were gradually lifted beginning at April 20th. We focus on modelling the time period when the social distancing measures were in place and then consider different scenarios of relaxation of the restrictions with possible improvement of testing and contact tracing. We note that the procedures for quarantine were in place even before the social distancing measures. They initially focused on individuals arriving from COVID-19 affected areas in China. When the epidemic started spreading in European countries people who came back to Poland from these countries were advised to immediately seek medical attention if they experienced any symptoms consistent with COVID-19. However, adherence to these recommendations was not evaluated. As soon as the first case was diagnosed in Poland quarantine for close contacts was also implemented. This paper aims to define a deterministic population model describing the epidemic in classical terms of susceptible, exposed, infectious, removed. In our model the quarantine becomes a separate state that removes individuals from susceptible and exposed states. We show that the reproductive number in our model is given by a simple formula referring to the parameters of transmission and transition, but also to parameters describing the quarantine. We demonstrate that in a real life scenario (case study of Poland) the quarantine effectively reduces the growth of infectious compartment. Increasing the efficiency of contact tracing and testing may may to some extent compensate lifting up the social distancing restrictions.
We introduce a modification of the classical SEIR model including effects of quarantine. To underline importance of that extension we call it SEIRQ. Formally the model is described by a system of ordinary differential equations with delay dedicated to the quarantine. 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 R u (t) -spontaneously recovered without being diagnosed Q(t) -quarantined
The 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 the infected 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 concerns these individuals who (independently of quarantine) would meet the testing criteria, as relevant to the local testing policy, e.g. testing of people with noticeable symptoms. We shall emphasize that from the point of view of analysis of spread of infection, the quantity I u shall be regarded rather as not recognized infections, not necessarily asymptomatic or mild. With this interpretation the value of κ 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, normalized to the total population size N = S + E + I d + I u + R d + R u + Q, which is assumed to be constant in time, births and deaths are neglected.
The transition parameter σ is assumed identical for both groups, relating to the time between infection and becoming infectious. The infectious individuals then move to the state R d , which is 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 contains people who spontaneously recovered with rate γ u . Our model includes an additional state of being quarantined (Q). To mimic the situation of contact tracing, individuals can be put in quarantine from the state S (uninfected contacts) or the state E (infected contacts). These individuals stay in the quarantine for a predefined time period T . We assume that the number of people who will be quarantined depend on the number of individuals who are diagnosed. An average number of individuals quarantined per each diagnosed person is denoted as α. However, as the epidemic progresses some of the contacts could be identified among people who were already infected, but were not previously diagnosed, i.e. the state R u . We note that moving individuals between the states Q and R u has no effect on the epidemic dynamics, therefore we assume that only individuals from S and E are quarantined and we reduce the average number of people put on quarantine by the factor S(t) S(t)+Ru(t) . Further, to acknowledge the capacity limits of the public health system to perform the contact tracing, we introduce a quantity K max , describing the maximum number of people 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:
where K(t) = min{ S(t) S(t)+Ru(t) αγ d I d (t), K max }, and α, β d , β u , γ d , γ u , θ, T ≥ 0.
(2.1)
We assume that the parameters α, β d , β u , θ, γ u and γ d depend on the country and time-specific public health interventions and may therefore change in time periods. Due to proper interpretation of the equation on E we require that β d ≥ θαγ d to ensure positiveness of E.
2.2 Basic reproductive number, critical transmission parameter β * .
Based on the general theory of SEIR type models [23] , we introduce the reproductive number
It determines the stability of the system as R < 1 and instability for R > 1 (the growth/decrease of pandemia). This quantity not only explains the importance is testing (in terms of κ) and quarantine (in terms of α), but also gives an indication on levels of optimal testing and contact tracing. We underline that this formula works for the case when the capacity of the contact tracing has not been exceeded (K(t) < K max ). The details of derivation of (2.2) are provided in the Appendix, section A.4. We shall emphasize the formal mathematical derivation holds for the case when I and E are small comparing to S. Therefore the complete dynamics of the nonlinear system (2.1) is not fully determined by (2.2) . However in the regime of epidemic suppression, which is the case of COVID-19 epidemic in Poland, I and E are small compared to S and so the formula (2.2) reasonably prescribes spreading of infection in the population. The critical value R = 1 defines the level of transmission which is admissible, taking into account the existing quarantine policy, in order to control epidemic. As the level of transmission depends on the 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 contact restrictions. The key quantity is θα. Indeed the system with the quarantine has the same stability properties as one without K, but with the new transmission rate β new d = β d − θαγ d . In order to guarantee the positiveness of E, β new d must be nonnegative. It generates the constraint
The above condition also implies the theoretical maximal admissible level of quarantine. We define it by improving the targeting of the quarantine, i.e. by the highest possible level of θ:
As long as the K max threshold is not exceeded the effect of the increase in θ or in α play the same role at the level of linearization (small I, E). However, in general it is not the case and for the purpose 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 of asymptomatic and symptomatic cases and although there might be a difference we lack information to quantify this difference. Then using formula (2.2) we compute critical values β * (κ, θ, α) defined as
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
which can be regarded as theoretical upper bound for β if we assume "optimal admissible" quarantine for fixed κ, for which the epidemic could be still controlled. It must be kept in mind though that the condition (2.3) means that we are able to efficiently isolate all persons infected by every diagnosed, therefore is unrealistic. The resulting β * (θ max , κ) should be therefore considered as a theoretical limit for transmission rate.
All simulations are performed using GNU Octave (https://www.gnu.org/software/octave/). The underlying 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 the model predictions to the data on the cumulative number of confirmed cases. Since people with confirmed diagnosis are efficiently isolated, they are immediately included into R d . Therefore, the quantity fitted to the data is R d (t).
The crucial assumption behind our approach is that the parameter β changes twice during the period of analysis. The reason is that we can distinguish two important time points in the development of epidemic in Poland. The first is initial restrictions including school closure effective March 12, which was accompanied with restrictions on other social activities. As we do not take migration 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 1 , t 2 } and choose t 1 = 14, t 2 = 28. With t = 1 corresponding to March 3 it means small delay with respect to the above dates which can be justified by the fact that new cases are reported with a delay of approximately 2 days.
Choice of fixed parameters (Tab. 2.1). We assume that the parameters σ, γ u represent the natural 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 infectivity period. γ u relates to the period of infectivity, which we select based on the research regarding milder cases, assuming that serious cases are likely diagnosed. Further, κ is a parameter related both to the proportion of asymptomatic infection and the local testing strategies. Since the literature findings provide different possible figures, for κ we examine three different scenarios.
Parameters γ d , θ and α are fixed in our model for the purpose of data fitting, but informed by available data. One of the scenarios of future dynamics of the epidemic (section 3.3) considers possible increase of θ. Parameter γ d was estimated basing on time from onset to diagnosis for diagnosed 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. The capacity 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 parameters collected in the following Optimization algorithm. In order to fit the values β 1 , β 2 , β 3 we use a standard gradient descent algorithm. The error function is defined as mean square difference between the cumulative number 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 possible conditions, as these mostly impact β 1 , which is less relevant for future predictions. To estimate confidence intervals we use a method of parametric bootstrap. The optimisation procedures are described 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 Poland from March 3 (first confirmed case in Poland) till April 26, which amounts to 54 observations. The data are taken from official communications of the Ministry of Health. As explained in table 2.1 and 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, 2, 3 correspond to the time intervals when different measures were in place, and the R for the third time interval. Given the social distancing measures in place early April 2020, as well as the quarantine levels, the reproductive number was below 1, independently of the value of κ, which relates to testing effectiveness. The figure 3.1 shows the fit of the models assuming different levels of κ. Good fit is found for all three models although predictions start to differ in the middle-term prognosis. We proceed with predictions assuming that the restrictions are continued, i.e.keeping β = β 3 (note that the estimated β 3 is different for each κ). We calculate the epidemic duration (t max ), the peak number of infected (I max d , I max u ) and the final size of the epidemic (R d (t max ), R u (t max )). In order to show the influence of quarantine we compare the situation with quarantine, keeping the same θ, α, to the situation without quarantine, setting αθ = 0. The results of the development of the epidemic during the first 120 days are shown on For κ = 0.2 the difference between the scenarios with and without quarantine is visible but not striking. However for κ = 0.5 and κ = 0.8 a bifurcation in the number of new cases occurs around t = 40 leading to huge difference in the total time of epidemic and total number of cases. These values are summarized in the table 3.2. We note that given the epidemic state in the first half of April 2020 for all values of κ the model predicts epidemic extinction both with quarantine and without quarantine. However, since the epidemic is very near to the endemic state, the predicted duration is very long, especially if no quarantine is applied. 
Using the formula (2.5) we can compute critical values β * . In Table 3 .3 we show the values of β * (κ, 0.006, 75) and for convenience recall also estimated values of β 3 and R, listed already in Table 3 .1. Moreover we compute β * (κ, 0, 0), i.e. without quarantine and show values of R for our estimated values of β 3 and the same γ d , γ u but without quarantine. Comparing the estimated values of β 3 (table 3. 3) for all cases of κ are only slightly below β * . Eliminating the quarantine, for the estimated values of β 3 , we have different situations depending on the actual value of κ. In case κ = 0.2, so assuming that currently only 20% of infections are diagnosed, the low values of R are due to low β 3 rather than the effect of quarantine (controlling epidemic 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.5 or κ = 0.8 we estimate higher β 3 , 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 β 3 − θαγ d represents effective transmission rate due to diagnosed cases. In particular it shows by how much the 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 influence 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 numerically the 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.006, α = 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 this case α = 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 and hence worse outcomes. 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 that the critical β * calculated for the linear approximation in the section 2.2 are adequate, with a small bias towards lower values. The case κ = 0.2 shows that the influence of quarantine is not high, even for the maximal admissible case, when we are able to efficiently isolate all persons infected by every diagnosed. A striking feature in the behaviour of total number of infected are jumps for certain critical value of β observed for κ = 0.5 and κ = 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 the first value of β for which the limit number of quarantined, K max = 50000, is achieved. Notice that immediately after passing the threshold the values become very close to those without quarantine. 
We propose a simple SEIR-type model (SEIRQ), which includes the effects of testing and contact tracing. The model formulation allows to calculate an interpretable formula for the reproductive number 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 the other hand γ d reflects the time to diagnosis and the formula indicates that more rapid diagnosis is associated with lower R. In addition, our model offers a clear interpretation of the quarantine effect. The transmission rate due to diagnosed cases, β d , is decreased by the factor θαγ d indicating that both the number of quarantined per diagnosed individual (α) and proper targeting of the quarantine (the infection rate among the quarantined θ) equally contribute to this factor. Also the parameter related to testing: the delay in diagnosis, γ −1 d , plays similar role. This quantifies the potential 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 each case and the infection rate among the quarantined impact R in similar fashion, our results support the recommendations to focus on the high risk contacts when the resources do not allow to follow all contacts.
Our model takes into consideration only the effective contact tracing, i.e. the situation when the infected contacts are identified and put in quarantine before they become infectious. People who are identified later would be modelled as passing through one of the I states to the R states. This means 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 approach given that the incubation time can be as short as 2 days in 25% of cases [3] . As mentioned by other Authors [28] , the delays in manual contact tracing are usually at least 3 days and under such circumstances 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 strategies implemented in South Korea indeed helped to control the epidemic at the early stages [32] . The use of "smart contact tracing" with mobile phone location data and administrative databases were also key to rapid identification and self-quarantine of contacts in Taiwan [33] and implementation of such 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. As expected, in order to control the epidemic the transmission due to undiagnosed cases has to be negligible. This can be controlled by general measures such as lockdown, which universally decrease the 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 to the efficiency of the testing system. Again, the examples of Singapore as well as the Italian village of Vo'Euganeo show that the widespread testing complementing the efficient contact tracing was essential to suppress epidemic. Testing unrelated to epidemiological links decreases (1 − κ) factor, thus making the factors impacting transmission due to diagnosed cases, such as quarantine, more powerful to decrease R. Further, our model allows to study the effect of the situation, in which the contact tracing capacities are exceeded. In this situation the epidemic is likely to quickly develop to the levels observed without quarantine. It is therefore quite crucial to implement the aggressive contact tracing system, when the epidemic is still at low levels and it is possible to bring the epidemic to suppression phase.
We demonstrate the high impact of contact tracing and quarantine on the observed numbers of cases 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 of the considered scenarios regarding testing efficiency was the most relevant to our situation. In our model the quarantine is estimated to be the most effective for the scenario in which most of the cases are diagnosed (κ = 0.8). Testing strategies that comprise testing of all individuals with symptoms of respiratory illness could theoretically identify up to 82% of infected, assuming they would all present to medical care. This could be coupled with random screening of high risk individuals, in e.g. health care workers, or -in case of high incidence -even random screening of entire community to achieve the κ of the order of 0.8. The Polish clinical recommendations specifically mention only testing all 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.8 scenario is unlikely to be realistic in Poland. We believe that the plausible current κ in our country lies between 0.2 and 0.5. For these scenarios the model shows that the control of the epidemic is largely achieved through suppression of β. In case of relaxation of social contact restrictions, the efforts should be focused on increasing the level of testing in order to decrease the proportion of undiagnosed 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 level of social contacts from before the epidemic (β 1 ).
Finally, the contact tracing effort was manageable in Poland due to relatively small number of cases. 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 reduced transmission from undiagnosed cases who are more likely to be asymptomatic or paucisymtopmatic cases (β 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 the transmissibility reduction by 50% [6] . Moreover, we did not have sufficient data to include this additional parameter. We calibrated our model only to diagnosed cases, which were officially available. Calibration to mortality data is another approach successfully implemented in e.g. [25] that potentially removes bias due to different testing policies. As there were relatively fewer fatalities in Poland and little data on clinical progression we decided on simplified model without explicit modelling 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, only anecdotal evidence for such phenomenon is available at this time. In our model it would decrease the 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 are still useful, but the approximation may be biased, as we have shown for β * . Due to little available data and policy changes we did not have sufficient data to determine which κ scenario is the most appropriate.
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 and we show that despite a substantial impact of contact tracing and quarantine, it is unlikely that the control of the epidemic could be achieved without any reduction of social contacts.
[ 
A.1 Optimization algorithm.
In order to fit the values β 1 , β 2 , β 3 we use a standard gradient descent algorithm. Namely, we define error function as
where {R d (κ, t)} 54 t=1 is the vector of computed values of R d and {data(t)} 54 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 β 1 , β 2 , β 3 by differential quotients and move in the direction opposite to the gradient. The algorithm reveals a good performance provided we start sufficiently close to the minimum, which is not difficult to ensure in our case.
It remains to choose the initial data. A closer look on results of simulations shows that the choice of initial data mostly influence the fitting in the beginning of period under consideration and hence the value of β 1 , while for analysis of future scenarios β 3 is the most important. Taking all this into account we do not struggle for sharp optimization of data fitting with respect to initial data and restrict to the following heuristic choice. It is natural to assume I u (0) = 1−κ κ 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 initial values I d (0) ∈ {10, 20, 30} and for each value we set I u (0) according to the above formula and three values of E(0) corresponding to m ∈ {2, 3, 4}. For each of these 9 combinations we run the optimization algorithm looking for the best fit of β i , i = 1 . . . 3. We have repeated this approach for κ ∈ {0.2, 0.5, 0.8}. 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: A.2 Choice of fixed parameters 1. The parameter σ describes the rate of transition from non-infectious incubation state E into the infectious states I d or I u . The median incubation time from exposure till the onset of symptoms was estimated at 4 to 5 days [3, 4, 5] . However, there exists evidence that typically infectivity preceeds symptoms, by 1 to 3 days [7, 9, 10] . A modelling study identified the rate of transition between the non-infectious and infectious states at 1 3.69 [6] , which corresponds to an average time lag of 3.69 days untill the case becomes infectious.
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 of viral shedding was estimated among 113 Chinese hospitalized patients. Overall it was 17 days, but it was shorter among cases with milder clinical course [16] . A study among 23 patients in Hong Kong confirmed viral shedding longer than 20 days among a third of patients, although the peak level of shedding was noted during the first week of infection [17] . In the mission report from China WHO 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 among 24 initially asymptomatic youngsters the median duration was 9.5 days [11] .
3. The value of κ generally depends on the testing policy. However, recommended testing policies often rely on the presence of respiratory symptoms. This is also the case in Poland. It was observed that some infected people never develop symptoms, although the precised rate of such truly asymptomatic infections is still under investigation. Some studies may be biased by a too short follow-up time. A small study among residents of a long-term care skilled nursing facility found that even though more than half of individuals with confirmed infection were asymptomatic at the time of test, majority of them subsequently developed symptoms. The proportion of people who remain asymptomatic may be higher among younger individuals [Hu] . A study among Japanese nationals repatriated from Wuhan suggests the proportion of asymptomatic infections is about 30% [13] . An analysis 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 asymptomatic infections 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% of infected individuals did not report symptoms [14] . Similarly, in population screening in Iceland 50% were asymptomatic at the time of screening [15] . It may be stipulated that some of the people diagnosed through screening developed symptoms latter, consistently with the findings from the Diamond Princess study.
On the other hand a sizable proportion of infected people, especially at younger ages, experience only mild symptoms, for which they may not seek medical attention. In the study of Li [6] , the proportion of undocumented cases was estimated as 86%.
4. The parameter γ d was estimated basing on a sample of case-based data available in routine surveillance, by fitting gamma distribution to the time from onset to diagnosis, for cases who were not in quarantine before diagnosis. Time from onset to diagnosis was estimated based on surveillance data available in the Epidemiological Reports Registration System for COVID-19, as of 28.04.2020.The system collects epidemiological data on cases diagnosed in Poland and is operated by 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 onset date and plausible diagnosis date were available for 2884 cases ( 96.7% of 2981 cases that were not asymptomatic) Gamma distribution was fitted by maximum likelihood to cases who were not diagnosed in quarantine. 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. We found no significant trend in time. We therefore adopted the average time from onset to diagnosis to be 4.6 days, and taking into account the probability of asymptomatic spread we assumed the parameter γ d to be 1/5.5.
5. Next we base θ on available data. We calculate prevalence of infection among the quarantined individuals, according to data published by the Chief Sanitary Inspectorate on the number of cases diagnosed among quarantined people and the total number of quarantined. We used a series of data 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 the proportion of diagnosed from 0.5% to 0.8% A.2. We presume that this parameter could change with changing 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 take a 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 proportion in Poland is in line with what was observed in Korea, where an estimated attack rate was 0.55% [32] , although household attack rate was higher (> 19%) in other studies [20] .
6. Furthermore we fix the parameter α. Here we make another simplification assuming this parameter to be constant. The main difficulty is a lack of precise data concerning the number of newly quarantined people per day, distinguishing between reasons of quarantine (travel related or contact tracing relate). At the beginning of epidemic in Poland the average amount of quarantined following one diagnosed case was definitely higher. Moreover, people coming back from abroad were subject to obligatory quarantine starting from March 16 and constituted a considerable part of quarantined in the second half of March and beginning of April. In particular, around 54 000 Polish citizens staying abroad came back within a special program of charter flights operated by Polish Airlines which ended on April 5. We can assume that after this date the ratio of people coming from abroad among all people subject to quarantine was negligible. As our model does not take migration into account, we have to take into account only quarantine from contact. For above reasons, 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 from the 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.
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 ∈ {1, 54} the increment of R (i.e. daily number of new diagnoses) is a random number from Poisson distribution with mean value equal to increment of observed data. Model parameters are also perturbed, see below. For each series of perturbed data we estimate the 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 proceed as follows. For each sequence of perturbed data we compute fitted R d (t). This way we obtain a set of curves {R Distribution of parameters. Following other Authors [30] as well as experimental data, for the uncertainty analysis we used the following distributions of the parameters.
1. 1/γ d ∼ Gamma(a 1 , b 1 ), where the shape parameter, a 1 = 1.05 and scale parameter b 1 = 5.23 2. 1/γ u ∼ Gamma(a 2 , b 2 ), where the shape parameter, a 2 = 2 and scale parameter b 2 = 5 3. 1/σ ∼ Gamma(a 3 , b 3 ), where the shape parameter, a 2 = 2 and scale parameter b 2 = 1.75 4. α ∼ P oisson(α 0 ), where we assume constant α 0 = 75; 5. θ -is not sampled for the uncertainty analysis. As the results depend on the quantity αθ, we rely on the distribution of α.
We take N = 200 (approximate average of daily number of diagnosed cases from the data). We approximate the mean value of N samples from Gamma distribution using Central Limit Theorem. Namely, we generate
Based on the classical approach to epidemiological models we address the basic question concerning the propagation of the disease. Namely, how many persons are infected by one infectious individual, a quantity which is usually called reproductive number, R. In order to compute this quantity we use the approach from [23] . We look at the system assuming S ∼ N and E, I d , I u are close to zero, then we consider the following linearizatioṅ Then one deduces (see [23] ) that if we define R = max{eigenvalues of T Σ −1 } then the system is stable for R < 1 and it is unstable for R > 1.
Stability of system (A.5) means that the whole vector (E, I d , I u ) is going to zero, it follows that the main system (2.1) also tends to the zero solution for (E, I d , I u ). Instability implies that for "almost all" small data, the vector (E, I d , I u ) is growing in time (exponentially fast), causing the nonlinear system also evolves far away from the trivial state, i.e. E, I d , I u rapidly grow. By (A.4) we have
Hence the stability of our system is determined by the following factor:
To make a final comment, let us note that in case of spread of pandemia, as R d , R u grow, the above analysis become less reliable. Recall that β is normalized by N , so as S/N is not close to one and the analysis of stability becomes more complex. This behavior is illustrated by figures 3.3, 3.4 and 3.5, where we observe deviation from predictions based on (2.2),