Modelling Excess Mortality in Covid-19-Like Epidemics

We develop an agent-based model to assess the cumulative number of deaths during hypothetical Covid-19-like epidemics for various non-pharmaceutical intervention strategies. The model simulates three interrelated stochastic processes: epidemic spreading, availability of respiratory ventilators and changes in death statistics. We consider local and non-local modes of disease transmission. The first simulates transmission through social contacts in the vicinity of the place of residence while the second through social contacts in public places: schools, hospitals, airports, etc., where many people meet, who live in remote geographic locations. Epidemic spreading is modelled as a discrete-time stochastic process on random geometric networks. We use the Monte–Carlo method in the simulations. The following assumptions are made. The basic reproduction number is R0=2.5 and the infectious period lasts approximately ten days. Infections lead to severe acute respiratory syndrome in about one percent of cases, which are likely to lead to respiratory default and death, unless the patient receives an appropriate medical treatment. The healthcare system capacity is simulated by the availability of respiratory ventilators or intensive care beds. Some parameters of the model, like mortality rates or the number of respiratory ventilators per 100,000 inhabitants, are chosen to simulate the real values for the USA and Poland. In the simulations we compare ‘do-nothing’ strategy with mitigation strategies based on social distancing and reducing social mixing. We study epidemics in the pre-vacine era, where immunity is obtained only by infection. The model applies only to epidemics for which reinfections are rare and can be neglected. The results of the simulations show that strategies that slow the development of an epidemic too much in the early stages do not significantly reduce the overall number of deaths in the long term, but increase the duration of the epidemic. In particular, a hybrid strategy where lockdown is held for some time and is then completely released, is inefficient.

In the last decades, the classical epidemic models have been reformulated in the framework of complex networks science [13]. Complex networks [14][15][16][17] are very well-suited to encoding heterogeneity of spatial distribution [18] and mobility of population [19][20][21]. New techniques, which go beyond the classical mean-field approach, have been developed and successfully applied to modelling of epidemic spreading in heterogeneous systems such as degree-based mean-field theory [22,23], models of clustering [24], spatial and mobility networks [19][20][21] and meta-population approach [25,26] where one can superimpose hierarchical transportation network on the population distribution in 1. In the absence of a vaccine, immunity can only be obtained through infection.

2.
People who get infected become infectious for about ten days. 3.
People who recover are immune to reinfection. 4.
About one percent of all infections lead to SARS. 5.
The occurrence and course of SARS is correlated with the health conditions and age of the infected person. 6.
SARS is likely to lead to respiratory failure and death unless the person receives appropriate medical attention. 7.
Respiratory ventilation decreases the probability of death. 8.
The death probability is correlated with general health conditions of the patient. 9.
The healthcare system has a limited capacity. Especially, the number of doctors and the number of trained medical personnel, and the number of intensive care beds and mechanical ventilators is limited. 10. The mortality rate from non-Covid causes, like cancer, cardiovascular diseases, or other chronic diseases, increases during epidemic because of epidemic restrictions in hospitals and health clinics. 11. The epidemic may spread in two distinct modes: via local transmission or non-local (global) transmission. The local transmission mode corresponds to geographic epidemic spreading through person-to-person contacts near the place of residence. The non-local transmission mode, in turn, corresponds to epidemic spreading through contacts in public places like: hospitals, cinemas, sport arenas, schools, universities, churches, airports, means of communication, workplaces and many others, where people, who live in different geographic locations, meet.
We are not attempting to develop a realistic model of Covid-19. Such a model would have to take into account many detailed medical and demographic factors, as well as detailed information about geographical population distribution, migration, social mixing, etc. Instead, we use Occam's razor to develop a model that is as simple as possible and that can be used to qualitatively estimate mortality for a variety of strategies used to restrict epidemic spreading. The idea is to examine dominant factors shaping the death statistics during epidemic whose spread is inhibited by large-scale restrictions on social contacts. We are mainly interested in statistical effects. Let us discuss how the above assumptions are implemented in the model. The numbering below refers to the assumptions.

1.
We are interested in outbreaks like Covid-19, for which there is initially no immunity or vaccine, and vaccine development and validation takes several years. We model an epidemic in the pre-vaccine era over a time span of 1000 days, when the only protection mechanism is herd immunity.

2.
The duration of the infection varies from person to person. The model assumes that it is a random variable defined by the geometric law with the mean τ = 10 days. The value τ = 10 should not be understood literally as ten but rather as the order of magnitude. The mean incubation period for Covid-19 was estimated to be 5.2 ± 1.8 days [34]. Based on existing literature, the incubation period is 2 to 14 days [35]. We do not distinguish the incubation, latent and infectious periods.
This simplification does not significantly influence the epidemic dynamics in large scale, which is simulated as a variant of MSIR dynamics [3]. The rate of epidemic is controlled by the basic reproduction number R 0 . According to early estimates [34,36,37], the basic reproduction number for Covid-19 ranged from 2.2 to 2.7, while according to an analysis of scientific literature on Covid-19 [38], the mean value of R 0 was found to range from 1.90 to 6.49. We set R 0 = 2.5 as the default value in the simulations. 3.
The model assumes that a recovered person is immune to the disease or, alternatively, that reinfections are rare or do not lead to SARS. For some diseases, reinfections are marginal and may be neglected in the description of epidemic spreading. There is currently an ongoing debate as to whether this is the case with Covid-19 [39]. The model does not apply to epidemics in which immunity fades over time and reinfections are likely to result in SARS.

4.
The estimated SARS fatality rate for Covid-19 ranges from 0.9% to 2.1% [40]. In the model, the SARS case frequency rate is 1%. This value should be treated as an order of magnitude, not a specific number.

5.
It is known that the occurence and course of SARS caused by Covid-19 is correlated with the age and co-existing diseases of the patient. We want to introduce such correlations to the model in a minimalist way. For this purpose the agent population is divided into a part where SARS cases are less frequent and a part where they occur more frequently. The first part can be thought of as healthy or young people, and the second as chronically ill or the elderly. We label the parts by H and C. The parts differ not only by the frequency of SARS occurence, but also by the SARS mortality rates. The split into the H and C groups may look artificial at first sight but it seems to be the simplest way of implementing the observed variation for Covid-19 of the probability of occurrence and the course of SARS for different age groups. In principle, the model is well suited to implement many age groups having different patterns for the occurence of Covid-19 SARS and other diseases, but this would make the model more complicated and it would not dramatically change the qualitative picture. So we stick to the simplest solution. We denote the part sizes p H and p C , respectively. To fix attention we choose p H = 75% of the total population and p C = 25%. This split roughly overlaps with the split of the population of a European country (like Poland), into people younger than 60 years (H), and older than 60 years (C). The implementation of mortality rates for SARS and other diseases as well as the frequency of SARS occurrence for H and C subpopulations will be presented in detail in the next section where we discuss the stochastic processes that describe the dynamics of the epidemic. 6.
The occurence of SARS cases is also simulated as a stochastic process. An infected person may develop SARS with a certain probability during the infectious period. In the model, this probability depends on whether the person is mechanically ventilated or not. The respiratory ventilation decreases the death probability, so it is important to support by ventilation as many SARS patients as possible. 7.
The number of ventilators, or more generally the healthcare system capacity is limited. Not everyone who needs a support may obtain it in time, when there are too many SARS cases at once. The limited capacity of the healthcare system is simulated in the model by a simple stochastic process of distributing ventilators between SARS patients. Once a patient receives a ventilator he or she will continue to use it until he or she recovers or dies. The ventilator is then transferred to a new person who is randomly selected from all SARS patients who are in need of one. 8.
The death probability is correlated with the general health conditions of the patient. In the model, it is simulated by differentiating the death probability according to the H and C groups to which the patient belongs. The details are given in the next section. 9.
The healthcare system capacity is simulated by the number of respiratory ventilators available (or the number intensive care beds). The effect is mainly expected in developed countries where it is related to delayed diagnoses and late admissions of patients with cancer [46,47] and coronary heart diseases [48]. For example, it was estimated that the diagnosis delays caused by one year of epidemic conditions would lead within 5 years in the UK to an increase of the number of deaths for colorectal cancer by 15.3-16.6%, for breast cancer by 7.9-9.6%, lung cancer by 4.8-5.3% and esophageal cancer by 5.8-6.0% [46,47]. Cancer Research UK has estimated that 2000 fewer cancers were being diagnosed per week in April 2020 as compared to three years earlier [46]. In Poland, in April 2020 only 50% of patients, compared to April 2019, used the system of rapid therapy, which had been introduced some years ago to speed up the treatment of oncological patients. Every third visit to an oncology doctor was canceled, the number of diagnoses using MRI, computed tomography, PET-CT decreased by 30% [49].
The death risk for people with cardiovascular diseases significantly increases during an epidemic. The number of patients ST-elevation myocardial infarction dropped during the lockdown. More than 40% of patients with a heart attack were admitted beyond the optimal time window [48].
These two examples show that protracted epidemic conditions in the healthcare system may have a significant impact on statistics of non-Covid deaths. Some effects will be seen with a time-lag. Cardiovascular diseases and cancer account for the largest share of death statistics. For example there were 647,457 deaths from heart disease, 599,108 from cancer out of a total of 2,813,503 deaths in the USA in 2017 [35]. This is roughly 44.3%. Thus, an increase in the number of deaths from these causes by a few percent may be significant for the entire population. People with other chronic conditions will be statistically more exposed to the death risk due to restricted access to healthcare resources during a long-lasting epidemic. Elderly people and people with chronic diseases are fearful of exposure to the virus so they avoid public places including hospitals and health clinics. In effect, they are more exposed to health risks. These phenomena are difficult to model, since they depend on many factors, which cannot be easily quantified, like the organization of the healthcare system, redeployment of resources during epidemic, quarantine procedures in hospitals etc. Instead of seeking a complicated model with many parameters which would describe all these factors we propose to investigate what happens when the rate of deaths, due to causes other than those related to the virus, increases on average by a factor x during the epidemic, where x is just an input parameter of the model. In particular, we study an increase of the daily mortality from non-Covid causes by x = 1%, 2%, . . . , 5%. 11. Geographic distribution of the population is simulated in the model by geometric 2d random networks [50,51], see Section 2.1 for details. Compared to classic random graphs [52], growing networks [14] or other classes of random networks which are constructed in a non-geometric way, such networks are much better at mimicking the distribution of social distances between people in a situation when social contacts in public places and non-local transmissions are limited. In the model, one can distinguish local and non-local modes of disease transmission. It is a simplified version of the meta-population dynamics [25,26]. Local transmissions are modelled by infections of neighbouring nodes of the network, while non-local ones by infections of randomly selected nodes, independently of their position in the network. The non-local mode of disease transmission simulates intense social contacts in public places where many people meet, who then move to distant places. The effect leads to outbreaks in remote places and therefore significantly accelerates the spread of the epidemic.

Methods
In this section we provide a detailed mathematical description of the model.

Random Geometric Networks
Random geometric networks are constructed by the proximity rule [50,51]. Two nodes are connected by an edge if they lie within the given distance from each other. The simplest example is a network constructed by connecting randomly distributed points in a d-dimensional Euclidean space. We are using this construction here for d = 2 to mimic geographical distribution of the population which defines a network of everyday social contacts. For sake of simplicity we assume that the points are uniformly distributed on a two-dimensional square with the periodic boundary conditions. This can be done by generating pairs of coordinates (x i , y i ), i = 1, . . . , N, consisting of 2N independent random numbers uniformly distributed on the unit interval [0, 1] and connecting any two points i and j by an edge of the network if the distance between them is smaller than : For the periodic boundary conditions the coordinate differences are calculated as follows ∆x ij = min |x j − x i |, 1 − |x j − x i | and analogously for ∆y ij . The node degree distribution of the network obtained in this way follows the binomial law where a = π 2 is the area of a circle of radius . The mean degree distribution is k = (N − 1)a, and the variance σ 2 (k) = (N − 1)a(1 − a). When a is of the order of 1/N, then the distribution becomes Poissonian in the large N limit. The node degree distribution (1) is identical as for Erdős-Rényi random graphs [52]. The two classes of graphs are however completely different. In particular, the average clustering coefficient for the geometric random networks is C = 1 − 3 √ 3 4π ≈ 0.586503 [50], while for Erdős-Rényi random graphs it approaches zero like 1/N as N tends to infinity [52].

Agent-Based Implementation of SIR Dynamics
We use a discrete-time stochastic implementation of the SIR dynamics [3,12]. The network is populated with agents residing on its nodes. The population is divided into three classes of susceptible (S), infectious (I) and recovered (R) nodes, which describe the state of each agent at time t. The states change in the course of evolution according to epidemic rules which are implemented in the model in the form of a discrete time stochastic process. Time is counted in days from the outbreak of the epidemic. Initially, that is for time t = 0, one agent, or a few ones are infectious, while all others are susceptible. An infectious agent remains infective for τ days on average, and then it recovers. This is simulated in the model by assuming that the probability of remaining infective till the next day is q and of recovering 1 − q. The lifetime distribution of infectious state is given by a geometric law The mean lifetime of an infectious state is related to the probability q as follows which means that for the expected infectious period is τ days. We symbol . . . stands for expected value. Clearly, for τ 1 the probability distribution (2) can be approximated by P i (t) ≈ e −t/τ /τ. Once an infectious person recovers, he or she remains immune and healthy until the end of the SIR of evolution. Later we will modify the SIR dynamics by superimposing on it the death dynamics by modifying some of the rules described in this section. In particular, we shall assume that a recovered person may die with some probability and then reappear as a susceptible newborn. This means, in particular, that the R state may change to S with some probability. We shall discuss the death dynamics in the ensuing subsections. The resulting dynamics is similar to that used in MSIR models [3].
If an infectious node, a, is in contact with a susceptible one, b, the disease can be transmitted from a to b, if the contact is sufficient for disease transmission. Let p be a probability of transmission from a to b in one day. The probability p t of a transmission within t days is p t = 1 − (1 − p) t . The lifetime of an infectious state is a random variable (2) so the transmission probability for the whole infectious period is equal to the expected value p t = 1 − (1 − p) t . This yields for q given by (4). A node has on average k neighbours, so the number of infections generated by a single infected node, in a fully susceptible population, is This equation relates the basic reproduction number R 0 to the parameters p, τ and k of the model. The epidemic evolution is implemented in a synchronous way. This means that all states are updated simultaneously. States at time t + 1 are computed from states at time t. The following rules are used to update the states. If a node is recovered at time t, it remains recovered at time t + 1. If a node is infectious at time t it remains infectious at time t + 1 with a probability q. Otherwise, it changes to recovered. If a node is susceptible at t it changes to infectious with a probability p * . Otherwise, it remains susceptible. The probability p * , that a susceptible node becomes infectious is related to the transmission probability p, by the following relation p * = 1 − (1 − p) i * where i * is an effective number of infectious neighbours and i n is the number of infectious nearest neighbours of the node in the network, that is those which are connected to it by a direct edge. I is the total number of infected nodes in the network. The parameter α ∈ [0, 1] interpolates between the local and non-local (global) transmission modes. In the local transmission mode, that is for α = 0, i * is equal to i n , while in the non-local transmission mode, that is for α = 1, i * is proportional to all infectious nodes on the network k I/N. Later, we shall compare the results of local and non-local transmissions with the results for classic SIR models [3,12]. In the classic approach one usually uses the continuous time formalism. The epidemic evolution is described by a set of first order ordinary differential rate equations for the fractions of susceptible, infectious and recovered agents: is conserved during the evolution [3,12]. s(t) is a non-increasing function of time t and r(t) is a non-decreasing function. The infectious fraction, i(t) increases for t < t max and reaches a maximum for t = t max such that R 0 s(t max ) = 1. Indeed, as one can see from Equation (8), the derivative di/ds = −1 + 1 R 0 s changes sign when this condition is fulfilled. For t > t max the epidemic begins to die out and i(t) decreases from the maximum to zero: i(t) → 0 when t → ∞. The fraction of susceptible population for t → ∞ gives the level of herd immunity s(t) → s hi . The value s hi can be found from Equation (8). In particular, if i(0) is very close to zero and s(0) = 1 − i(0), then s hi is a solution to the equation ln s hi = R 0 (s hi − 1). This yields s hi ≈ 0.4172, 0.2032, 0.1074, 0.0595 for R 0 = 1.5, 2, 2.5, 3, respectively, to give some examples.
We use the following input parameters in the Monte-Carlo simulations of the epidemic on geometric random networks: the number of agents N, the mean node degree k , the basic reproduction number R 0 , the expected duration of the infectious period τ, the probability α of long-range transmissions. As an initial configuration, we choose I 0 randomly selected infectious nodes. The remaining nodes are susceptible. The probability to remain infectious till the next day is calculated from Equation (4). The probability of virus transmission from an infectious to a susceptible agent within one day is calculated from Equation (6) which gives An example of input values used in the simulations is N = 10 5 , k = 100, R 0 = 2.5, τ = 10, α = 0, I 0 = 5.

Modelling Background Conditions
In order to assess the impact of epidemics on death statistics, one also has to determine the death statistics and the background conditions in the absence of an epidemic. This is per se an interesting and very complex problem since it involves demographic factors, efficiency of healthcare systems, statistics of diseases, and many other factors. This is beyond the scope of this paper. We only model here basic factors to assess how death statistics change during a pandemic. The population is divided into classes according to health conditions. In the simplest version of the model we introduce two classes that correspond to healthy people and people with chronic diseases. We label the classes by H and C, respectively. The division is symbolic, but it allows the inclusion of statistical correlations between health conditions and mortality in simulations. This is modeled by choosing the mortality rate in the C class to be much larger than in the H class. The second important difference between the classes is that the death probability during epidemics increases faster in the C class than in the H class. The details are given in the next subsection where we discuss modelling of death statistics.
We assume that the size of the population is constant during the epidemic. The number of deaths is compensated by the number of newborns. This modifies the SIR dynamics that we described in a simplified version in the previous section. Denote the fraction of healthy people at time t by h(t) = H(t)/N, the fraction of chronically ill people by c(t) = C(t)/N and the fraction of deaths by We implement the population dynamics as a discrete time stochastic process (Markov chain) with the following evolution equation The matrix in this equation is a stochastic matrix. It describes the transition probabilities between the states H, C, D. The transition rates p DH and p DC add up to one p DH + p DC = 1, which means that the number of deaths is equal to the number of newborns. The parameter p DH is the probability that newborns are healthy at birth. For sake of simplicity, but without loss of generality, we additionally assume p HD = p DC = p CH = 0. The condition p HD = 0 means that the mortality rate of healthy people is zero or it is much smaller than the mortality rate of chronically ill people. The condition p DC = 0 means that a dead person is replaced by a healthy newborn. Thus the total size of the population is conserved. The condition p CH = 0 means that a chronically ill person does not become healthy again. Under these assumptions the last equation can be simplified to The transfer matrix has only two free parameters: β-the rate of becoming chronically ill and γ-the rate of dying. This stochastic process has a stationary state In our study we choose β and γ to reproduce the values d * = 2.3 · 10 −5 or d * = 2.74 · 10 −5 which correspond to the daily mortality rates in the USA and in Poland, as discussed in Section 1. We keep the ratio h * /c * = 3, so that the simulated population approximately consists of 75% people in the H class and 25% in the C class. For this choice, the paremeters of the transfer matrix (11) are and h * = 3 4 (1 − d * ) and c * = 1 4 (1 − d * ). We conclude this section with two remarks. Firstly, we have assumed that there is no direct transfer from H to D, from D to C and from C to H within one day, by setting p HD = 0, p DC = 0 and p CH = 0. One should note that the probabilities of transfers between these classes in two (or more) days are non-zero Secondly, the square or a higher power of the transfer matrix (11) is also a stochastic matrix. In principle, one can replace the original transfer matrix with any power of it, and interpret it as a daily transfer matrix. This will not change the stationary state. The stationary state is a left eigenvector of the transfer matrix associated with the eigenvalue 1 and it is identical for the transfer matrix (11) or any power of it. The transfer matrix (11) has three eigenvalues. The one which has the largest absolute value is λ 1 = 1 and the second largest is λ 2 ≈ 1 − β − γ. The eigenvalue λ 2 tells us about correlation of states at different times t, t . The correlation function decays exponentially as exp(−|t − t |/T). The correlation time T can be derived from λ 2 : T ≈ −1/ log(λ 2 ) ≈ 1/(β + γ). For the transfer matrix (11) T is of order 10 4 . By raising this matrix to the n-th power and interpreting the resultant matrix as a daily transfers matrix one can reduce the autocorrelation time from T to T/n.

Simulating Death Statistics during Epidemic
Let us begin this section by recalling the philosophy behind splitting the population into parts H and C. The mortality rate and the course of SARS for Covid-19 are known to be strongly correlated with age and co-existing diseases. Elderly people and people with chronic conditions die from Covid-19 SARS more frequently than young and healthy people. If one wanted to make the model very realistic one should divide the population into many age groups and, for each, collect good statistics on SARS frequency and mortality and implement these statistics into the model. This would make sense only if all other elements of the model were realistic. This is not the case in our study. The model we develop is minimalistic but it should of course implement all important factors, including the correlation between underlying diseases and SARS mortality. The split into two classes with distinct statistical properties is the simplest way of doing it. For example, we assume that the frequency, p H,sars , of SARS cases in H class is much smaller than the frequency, p C,sars , in the C class. For the sake of simplicity, we assume that the frequency p C,sars is an order of magnitude larger than p H,sars . The values p H,sars and p C,sars have to be consistent with the average SARS frequency which was previously assumed to be 1%: where h * ≈ 3/4 and c * = 1/4. In our simulations we use the following values p H,sars = 1/300 and p C,sars = 3/100, which give the correct average. For this choice, the frequency of SARS cases in the C class is almost ten times larger than in H. This is the first major difference between H and C classes. Another factor that plays an important role in the death statistics during epidemics is the fatality rate for SARS, which also should be significantly different for H and C. In the model, we distinguish four situations, labeled by C 0 , C 1 , H 0 , and H 1 : • We assume that the probabilities of dying from SARS are 1.0, 0.3, 0.9 and 0.1 for C 0 , C 1 , H 0 and H 1 , respectively. These values model a different course of SARS depending on co-existing diseases and access to a ventilator. They mean that respiratory ventilation increases the probability of staying alive from 0% to 70% for people with SARS in the C class, and from 10% to 90% for people with SARS in the H class.
In the simulations, as input paramaters, we use probabilities of dying within one day. They are related to the probabilities of dying in the whole period of infection by an equation identical to Equation (5) in which p t is interpreted as the probability of dying from SARS during the whole period and p is the probability of dying within one day. For τ = 10, the corresponding daily rates are 1.0, 0.041, 0.474, 0.011 for compartments C 0 , C 1 , H 0 , H 1 , respectively. During an epidemic, the number of people with SARS may easily exceed the number of ventilators available. In the simulations we set V = 27 (or V = 53) ventilators per 100, 000 people. These numbers are close to those for Poland (USA), as discussed in Section 1. A patient with SARS occupies a ventilator until he or she recovers or dies. In the model, this takes ten days (τ = 10) on average. So, if for some time there are more than 2.7 (5.3) new SARS cases a day per 100, 000 people in Poland (USA), the demand for ventilators will exceed the healthcare capacity.
The ventilator availability is simulated as follows. At any moment of time, the algorithm keeps track of the number of available ventilators. If this number is larger than zero, and there is a new SARS case, the number is decreased by one, and one SARS patient is moved between compartments C 0 to C 1 or H 0 to H 1 , respectively. The ventilator is occupied until the patient recovers or dies, in which case the number of available ventilators is increased by one. Initially, the number of ventilators is set to V per 100,000 people.
Another factor that has to be taken into account in assessing the epidemic total death toll is a lower efficiency of the healthcare system during epidemic [46][47][48]. This has an impact on the increase of deaths from non-Covid-19 SARS causes. The effect is significant in the group of people with oncological cardiovascular diseases [46][47][48], but also in the group of people who require continuous medical assistance. To estimate this effect, systematic statistical surveys should be carried out. Here we just assume that the number of deaths from other causes than those directly related to SARS increases by a factor 1 + x during an epidemic, where x is a few percent. In the model this is implemented by changing the value of the parameter d * from d * to d * (1 + x) and recalculating the parameters β and γ (13) of the Markov transfer matrix (11) for days when the number of infectious agents is I > 0.

Modes of Infection Transmission
In the model, the epidemic spreads on a geometric random network through local and non-local transmission modes. The non-local mode is selected with probability α, and the local one with 1 − α, as described before. For α = 1, the epidemic spreads by the classical SIR mean-field dynamics [3,12] which depends only on the node degree distribution, while α = 0 it follows a quasi-diffusive dynamics reflecting the geographic population distribution. In Figure 1 we show phase portraits for epidemics with different values of α on random geometric networks with N = 10 5 nodes. As one can see from Figure 1 the results of simulations for α = 1.0 are very well described the by the phase-portrait (8) of the classical SIR compartmental model [3,12]. The number of infectious agents it maximal at S max /N ≈ 1/R 0 ≈ 0.4 and the herd immunity is achieved for S hi /N ≈ 0.1 − 0.11, which is the place where the curve crosses the horizontal axis. This value is close to the mean-field prediction (8). The value of the basic reproduction number of the best fit to the theoretical curve given by the mean-field solution (8) is R 0 = 2.53. It differs by one percent from the value R 0 = 2.5 used in the Monte-Carlo simulations. The difference can be attributed to the fact that the classical mean-field dynamics is deterministic [3,12] and R 0 is a number, while in the simulations the dynamics is stochastic and R 0 is the mean value of a random variable. The variance of this random variable introduces some corrections to the effective value of R 0 .
The phase portrait starts to deviate from the mean-field solution when α decreases (see Figure 1). As shown in the right panel in Figure 1 the phase portraits for different simulations for α = 1 lie on top of each other and are consistent with the classical SIR solution. The curves for α = 0 have stochastic shapes and they differ from each other.
The herd immunity value S hi weakly depends on α (see Figure 1). The values of S hi /N ≈ 0.10 − 0.11 are almost identical for α = 1, and α = 0. What depends on α is the height of the curve which is a few times larger for α = 1 than for α = 0. This means that long-range social mixing significantly speeds up epidemic spreading. The effect is illustrated in Figure 2 where we compare dynamics of the epidemics for four different scenarios which differ by the basic reproduction number R 0 and the long-range social mixing parameter α. One can see that the spread of epidemic depends not only on the basic reproduction number R 0 but also on the long-range social mixing parameter α. Decreasing the parameter α models closing airports, schools, churches, sport arenas, etc., while decreasing the reproduction number R 0 models social distancing that is maintaining physical distance between people, reducing the frequency of personal contacts, wearing masks as well as disinfection, quarantine, isolation, etc. In the next section, we will evaluate the impact of these measures on mortality during the epidemic using Monte-Carlo simulations.
Let us make a couple of remarks to conclude this section. If one used, in the simulations, Erdős-Rényi random graphs with exactly the same node degree distribution (1), then one would observe the classical mean-field epidemic dynamics [3,12] independently of the value of α. There would be no distinction between the local and non-local transmission modes. Spatial distribution of nodes plays an important role in imitating geographic epidemic spreading. Epidemics spreading in classical random networks [53] are completely different than in geometric graphs, or more generally, in spatial networks, where it has a quasi-diffusive character [18].
The second remark regards the interpretation of results. The trajectories shown in Figure 2 represent single courses of epidemic for the given parameters. The model is stochastic and non-linear so trajectories for other time courses for the same parameters may look differently. For example, the epidemic may die out before it reaches say 1% of the whole population, because of a statistical fluctuation. We performed multiple runs to see how often it ends below the 1% threshold. The results for the four scenarios from Figure 2, are presented in the column A of Table 1. The column T shows the average duration time of epidemic, r hi is the immune fraction of the population at the end of the epidemic, and R hi /T is the average number of daily infections. Averages were calculated only from those cases that exceeded the 1% threshold. Table 1. Column A shows the percentage of simulated epidemics that expired before reaching 1%; T is the average duration of epidemics; r hi is the percentage of recovered people at the end of epidemic; and R hi /T is the average number of new infections per day for the four scenarios shown in Figure 2. We shall use the four scenarios in the next section to analyze excess deaths during epidemics. It is, therefore, useful to take a look at the last column of Table 1 which contains averages of quotients of R hi at T. The values correspond to the average numbers of new infections a day per 100,000 people for scenarios 1, 2, 3 and 4, respectively. Since the model assumes that the frequency of SARS cases is 1%, this means that one can expect c.a. 4.1, 2.1, 1.4 and 0.35 new SARS cases a day per 100, 000 people on average, and much more in the peak. These numbers should be compared with the healthcare system capacity, which is modeled by the number of available ventilators (or ICU beds) which are V = 27 for Poland and V = 53 for the USA per 100, 000 citizens, as discussed in Section 1. A ventilator is occupied on average for τ = 10 days, thus, the maximum capacity of the healthcare system to admit new SARS patients is V/τ = 2.7 or 5.3.

Assessing Mortality Rate for Different Scenarios
We are now going to compare excess death statistics for the simulated epidemics for six scenarios:

2.
R 0 = R 0 = 2.5 and α = 0.0. This simulates a suppression of virus transmission through reducing long-range social mixing.
R 0 = 1.5 < R 0 and α = 0.0. This simulates a quasi-lockdown. Both the local and non-local transmission modes are restricted.

5.
A quasi-lockdown for 300 days, as in item 4, and then do-nothing strategy, as in item 1.

6.
A quasi-lockdown for 600 days, as in item 4, and then do-nothing strategy, as in item 1.
As far as the mortality rates and the capacity of the health-care system are concerned the 'US' and 'PL' systems imitate the situation in the USA and in Poland, while 'US-V2' and 'PL-V2' simulate a hypothetical situation when the capacity of the healthcare systems would have been doubled in the two countries. The 'US-D' and 'PL-D' systems, in turn, simulate a situation when a pharmaceutic therapy would have reduced the number of SARS patients who require mechanical ventilation. The resulting effect of introducing an effective drug and doubling the number of ventilators is simulated by the configurations 'US-V2D' and 'PL-V2D'.
The six scenarios in those eight systems are studied for different values of the parameter x, which controls the increase of mortality from non-Covid-Sars causes [46][47][48]. We scan the range of x from 0% to 5%. The results are collected in Tables 2-9. Each entry corresponds to the average number of additional deaths after 1000 days per 100, 000 people, calculated from 100 independent simulations. The values in parenthesis represent statistical uncertainties. Only cases of epidemics that exceeded the 1% population threshold were included in the analysis. The resulting values should be referred to the expected number of deaths in 1000 days per 100,000 people in the absence of an epidemic, that is: 2740 in Poland and 2300 in the USA. When analyzing data in the tables, it is worth remembering that for scenario 4, the epidemic lasts longer than 1000 days (see Table 1). Table 2. Excess deaths 1000 days after the outbreak for 'US'.
x S1 S2 S3 S4 S5 S6 Let us, for illustration, present some results graphically. In Figure 3 we show an example of time evolution of the number of additional deaths during 1000 days after the outbreak in the system 'US-V2' for six different scenarios. The first four scenarios are shown in the left panel in Figure 3 and the remaining two in the right one. In the left panel, we additionally draw two reference curves representing the worse-case scenario when no SARS patients receive required medical assistance during an epidemic, and the scenario when there is no epidemic. For scenario 1, the number of daily new infections is large and the number of SARS cases exceeds the healthcare system capacity, so the number of daily SARS deaths is large. The effect manifest as a steep part of the mortality curve. The epidemic lasts a short time. For scenarios 2, 3 the epidemic lasts longer but the number of daily new infections is much lower. In effect, most of SARS patients obtain required medical attention, so the daily excess mortality rate is much lower than in scenario 4. In scenario 4, the epidemic spreads very slowly. The number of the new daily SARS cases is small, much below the healthcare system capacity. SARS patients are optimally treated, however, deaths from causes other than Covid are increasing due to protracted epidemiological restrictions. The graphs in the right figure show what happens when the lockdown lasts for 300 or 600 and then it is completely lifted. One can see that at the end of the studied period the total number of deaths is roughly the same as in the 'do-nothing' strategy, shown in the figure for reference.
As the next example, we compare in Figure 4, additional deaths 1000 days after the outbreak for all six scenarios in the systems 'PL-V2' and 'US-V2' which simulate hypothetical situations of the doubled capacity of the healthcare systems in Poland and the USA. The slope of the graphs increases with the duration of the epidemic as additional deaths from causes other than Covid are increasing with time. In particular one can see that the curves for scenarios 3 and 4 intersect for x close to 2%. In other words, these two strategies are comparable in this case. For economic reasons, strategy 3 is, however, much better than strategy 4, because it takes much less time. We also see in the figure that strategies 1 and 6 lead roughly to the same number of additional deaths.
In Figure 5 we compare the effect of doubling the healthcare system capacity, measured by the number of ventilators (or ICU beds) in scenario 1 ('do nothing strategy') and scenario 4 ('quasi-lockdown') in the USA. We use parameters as for the 'US', 'US-V2', 'US-D' and 'US-V2D' systems. In the left figure, we plot graphs for the 'do-nothing strategy'. We see that excessive mortality is approximately 23% and slowly varies with x. If the number of ventilators doubled, the excessive mortality would drop to around 15%. The introduction of a drug, that reduces the number of SARS cases requiring respiratory ventilation to 0.5%, would reduce the excessive mortality to approximately 8%, and if additionally, the number of ventilators doubled, to approximately 5%. The picture is completely different for scenario 4, as shown in the right figure. We see that the graphs for 'US' and 'US-V2' basically overlap, meaning that the doubling of the number of ventilators has no effect on mortality in this case. The same holds for 'US-D' and 'US-V2D'. Clearly, in scenario 4, the quantities of daily infections are so small that the healthcare system has a sufficient capacity. Additional ventilators are unnecessary in this case.

Discussion
We conducted a Monte-Carlo study of epidemic spreading on random geometric networks to assess the efficiency of non-pharmaceutic interventions in reducing the total number of surplus deaths during Covid-19-like epidemics. We discussed strategies based on social distancing and restricting long-range social mixing. They have different effects on epidemic spreading. Social distancing reduces the basic reproduction number R 0 to some effective reproduction number R 0 < R 0 . Restrictions on long-range social mixing reduce virus transmission between remote places. When long-range mixing is large, an epidemic spreads via mean-field dynamics. When it is small, it spreads via quasi-diffusive dynamics depending on the geographic population distribution. We studied the two modes of disease transmission here.
There are two sources of deaths which contribute to the total death toll during a Covid-19-like epidemic. One is related to Covid-SARS, and the other one to other diseases. The number of SARS deaths depends on the capacity of the healthcare system, which in the model is simulated by the number of available respiratory ventilators. If the daily number of new SARS cases exceeds V/τ, where V is the number of ventilators and τ is the number of days of using one ventilator for one SARS patient, some people with SARS will not be ventilated and will have lower survival probability. This effect was simulated in the model. If one assumes that there are 27 or 53 ventilators per 100,000 people, as in Poland or the USA, and τ is approximately 10 days, then V/τ = 2.7 or 5.3. As long as the number of SARS patients is below V/τ, then the number of deaths caused by SARS is maximally reduced. This effect can be achieved by slowing down the epidemic. On the other hand, the number of excess deaths from other causes may increase [46][47][48] with the epidemic duration so it is not beneficial to slow the rate of epidemic spreading too much. The optimal solution is to keep the number of SARS cases close to the capacity of the healthcare system, but not much below it.
We also showed that a strategy of maintaining the lockdown for some time and then releasing it by removing all restrictions has a similar effect on the number of deaths in the long term as if the do-nothing strategy was introduced right at the beginning. The deaths differ only by the time when they occur: in the do-nothing strategy the mortality is large at the beginning while in the other case it is large when the lockdown is released.
A strict lockdown makes sense only when one wants to gain time to increase the healthcare system capacity, for instance, buying new ventilators, increasing the number of ICU beds, training medical personnel or improving medical and epidemic procedures, or when an effective drug or vaccine is expected to be introduced in a short time. Otherwise, the optimal strategy is to keep the epidemic progress at the level that the number of SARS patients at any time is roughly equal to the capacity of the healthcare system. If the number of SARS cases is much larger than that, too many people will die of SARS. If it is much smaller than that and the epidemic will last too long, many additional people can die from cancer, cardiovascular diseases and chronic diseases, due to later diagnoses, later admissions for hospitalization and restricted access to health services [46][47][48].
Social distancing reduces the herd immunity level, see Table 1. This means that after lifting restrictions on social distance and restoring normal contacts between people, the percentage of immune people will be below the herd immunity level for the unrestricted system. The system will be unstable, in the sense that a new single infectious person may trigger a new outbreak. The situation is similar to a superheated liquid, where boiling may occur spontaneously at any time. For example, the total number of deaths in the simulated epidemic is comparable for scenarios 2 and 3 (see Figure 3), but the percentage of immune people at the end of the epidemic is 89% in scenario 2 and 59% in scenario 3 (see Table 1). The value 89% is close to the herd immunity level for R 0 = 2.5 while 59% is far below it. This means that the epidemic in scenario 3 can restart from the level 59% when the restrictions are lifted and a new infective person appears. This example shows that strategies reducing long-range social mixing bring better effects than introducing social distancing locally. They are, however, much more difficult to implement.
Let us compare the current Covid-19 mortality to typical mortality rates in Poland. Rates are quoted as daily deaths per 100,000 people. For example, in 2014 the daily mortality from all causes was approximately 2.71 including 0.69 from cancer [45,54]. The daily number of deaths registered as Covid-19 deaths between March the 5th and September the 27th, 2020, in Poland was approximately 0.031 [55]. According to the WHO data, the cumulative number of registered Covid infections in Poland in the quoted period was 227.2 per 100,000 people (approximately 2 per mille). The rate of spreading for Covid-19 is much slower than the simulated epidemics (see Table 1). Covid-19 has been spreading very slowly so far. If Covid-19 continued spreading at this rate the epidemic would need many years to end, unless an efficient vaccine is introduced. The number of registered cases is probably much smaller than the real number of Covid-19 cases because mainly suspected cases have been tested so far, so the statistics may be very biased. Assuming the number of cases is up to ten times underestimated, that would mean that at the end of September 2020 approximately 2% of people in Poland are immune to Covid-19, and 98% are still susceptible and face a Covid-19 infection. If the epidemic speeds up now too quickly, the effect can be as in scenario 5 or 6, discussed in this paper. According to the model, initial suppression of an epidemic does not reduce accumulated deaths over the long term, but extends the duration of the epidemic. Most of the European countries decided to suppress the Covid-19 in the first six months after the outbreak. Sweden took a different approach.
The comparison [56] shows that at the beginning there were relatively more deaths in Sweden than in other European countries, but this comparison does not take into account that the epidemic in Sweden is at a more advanced stage, which means that there are more people who are already immune to Covid-19. One has to wait with comparisons until the end of the epidemic.

Conclusions
Let us underline that the model developed in the paper does not attempt to simulate the Covid-19 pandemic but only to imitate some of its aspects. The basic assumptions are that immunity can be obtained only by infection and that reinfections are rare and can be neglected. Under these assumptions, the pandemic ends only after the herd immunity is achieved. The model has been constructed in a minimalistic way. The scale parameters of the model, like the number of ventilators and mortality rates simulated the real values. The conclusions drawn from the model can be treated qualitatively. Let us recall the main ones: • Strong suppression of an epidemic in the early stages does not significantly reduce the total number of deaths over the long term, but increases the duration of the epidemic; • In the absence of an efficient drug and a vaccine, the optimal strategy for reducing the total death toll for Covid-19-like epidemics, is to keep the number of new infections at a level where the number of SARS cases is as close as possible to the capacity of the healthcare system. • In the early stages of an epidemic, suppression should be only then implemented when one wants to gain time to increase the efficiency of the healthcare system or if the introduction of a drug or a vaccine is expected in a short time.
In contrast to the model, in the real world, it is very difficult to fine-tune the parameters that control the rate at which an epidemic spreads and to implement appropriate measures in society, without a vaccine.
Funding: This research received no external funding.