aa r X i v : . [ phy s i c s . s o c - ph ] J un Modelling Excess Mortality in Covid-19-like Epidemics
Zdzislaw Burda ∗ AGH University of Science and Technology, Faculty of Physics and Applied Computer Science,al. Mickiewicza 30, 30-059 Krakow, Poland; [email protected]
We discuss a stochastic model to assess cumulative excess deaths during Covid-19-like epidemicsfor various non-pharmaceutic interventions. The model simulates three interrelated stochastic pro-cesses: epidemic spreading, availability of respiratory ventilators and changes in death statistics.Epidemic may spread either locally or globally. The local mode simulates virus transmission throughcontacts in the vicinity of the place of residence while the global mode simulates virus transmissionthrough social mixing in public places, sport arenas, airports, etc, where many people meet, wholive in remote geographic locations. Epidemic is modelled as a discrete time stochastic process onrandom geometric networks. In the simulations we assume that the basic reproduction number is R = 2 . and the infectious period lasts ca. ten days. We also assume that the virus leads to severeacute respiratory syndrome in about one percent of cases, which in turn almost surely lead to res-piratory default and death, unless the patient receives an appropriate medical treatment supportedby respiratory ventilation. For other parameters, like mortality rate or the number of respiratoryventilators per million of inhabitants, we take values typical for developed countries. We simulatepopulations of − people. We compare different strategies: do-nothing, social distancing, re-duction of social mixing and lockdown, assuming that there is no vaccine and no efficient medicine.The results of the simulations show that strategies that slow down the spread of epidemic too muchare inefficient in reducing the cumulative excess of deaths. A hybrid strategy in which lockdown isin place for some time and is then completely released is inefficient as well. ∗ [email protected] I. INTRODUCTION
Mathematical and computer modelling have proved to be very useful tools for controlling existing infectious diseases[1–4] as well as for analysing and forecasting epidemics [5–7]. Modelling of infectious diseases and epidemics has along history [8–11]. The foundations of the contemporary theoretical epidemiology were laid by W.O. Kermack andA.G. McKendrick [12]. Today, theoretical epidemiology is a matured field of research [1–4].In the last decades the classical epidemic models have been reformulated in the framework of complex networksscience [13]. Complex networks [14–17] are very well suited to encode heterogeneity of spatial distribution [18]and mobility of population [19–21]. New techniques, which go beyond the classical mean-field approach, have beendeveloped and successfully applied to modelling of epidemic spreading in heterogeneous systems. Among these newtechniques it is worth mentioning degree-based mean-field theory [22, 23], spatial and mobility networks [19–21]and meta-population approach [24] were one can superimpose hierarchical transportation network on the populationdistribution in communities, cities, regions and countries, to differentiate between disease transmission modes in theregional and global scales. The models became very realistic. They can be fed with real-world data and used toforecast real world epidemics [6, 7, 20, 25–30].In this article, we are developing a model for Covid-19-like epidemics that lead in about one percent of cases toSevere Acute Respiratory Syndrome (SARS), causing respiratory failure and death [31]. In order to make the modelfully realistic for a specific epidemic, say Covid-19, it would require taking into account all medical, demographic,social, economic, geographic and any many other factors with thousands of parameters. These factors are still underactive investigation and need to be determined. That is why in this paper we talk about a class of similar epidemicsrather than a specific one. We focus on key factors that are important to assess the impact of various epidemiccontainment measures on reducing excess deaths. They are the following: • In the absence of a vaccine immunity can only be obtained through infection; • People who get infected become infectious for some time; • People who recover are immune to reinfection; • About one percent of all infections lead to SARS; • SARS occurrence is correlated with health conditions and the age of the infected person; • SARS leads almost surely to respiratory failure and death unless the person receives an appropriate medicalassistance supported by artificial ventilation; • Respiratory ventilation decreases the death probability; The death probability is correlated with general healthconditions of the patient; • The health care system has a limited capacity - in particular, the number of intensive care beds and mechanicalventilators is limited; • The mortality rate of other causes like cancer, cardiovascular diseases, diabetes and other chronic diseasesincreases during epidemic because of epidemic restrictions on medical procedures;The model simulates the spread of epidemics, the availability of respiratory ventilators for patients with SARS duringepidemic, as well as a stochastic process that shapes population reference statistics. These three interrelated stochasticprocesses are important for assessing the number of excess deaths during Covid-19-like epidemics and for selectingthe optimal strategy to reduce deaths.
II. RESULTSA. Modes of infection transmission
Epidemic spreading is simulated on geometric random networks [32, 33], see Section IV A for details. There aretwo different disease transmission modes: a local and a global one. Local transmission corresponds to geographicepidemic spreading through human-to-human contacts near the place of residence. This is modelled as transmissionsbetween neighbouring vertices of the network. The global mode corresponds to the disease transmission throughcontacts at places, like hospitals, cinemas, sport arenas, schools, universities, churches, airports, communicationsmeans, workplaces and many others, where people, who live in different geographic locations, meet. This effect is I N F E C T I O U S SUSCEPTIBLEMEAN FIELD1.00.10.010.0020.0 0 5000 10000 15000 20000 25000 30000 0 20000 40000 60000 80000 100000 I N F E C T I O U S SUSCEPTIBLEMEAN FIELD1.01.00.00.00.0
FIG. 1. (Left) Phase portraits of a simulated epidemic, with different values of the long-range social mixing parameter α , onrandom geometric random network. The simulations are carried out on networks with N = 10 nodes and the mean nodedegree h k i = 100 . The basic reproduction number used in the simulations is R = 2 . , the infectious period duration is τ = 10 and the mixing parameter is α = 1 . , . , . , . , . (from top to bottom). The results for α = 1 . are shown in symbolsand they are compared to a theoretical mean-field result (8) (solid line) going through the symbols. The value of the basicreproduction number in the mean-field result is R = 2 . . (Right) Two different simulations for α = 1 . (symbols) comparedto the mean-field result (solid line), and three different simulations for α = 0 . . implemented in the model as transmissions between randomly selected nodes, independently of their positions on thenetwork. We shall call this effect long-range social mixing. Long-range social mixing leads to outbreaks in remotelocations and it accelerates epidemic spreading.In the model there is a control parameter α which interpolates between these two modes of transmission. Theepidemic spreads by long range social mixing with a probability α or it spreads locally with a probability − α . For α = 1 the epidemic is described by the classical SIR mean-field dynamics [3, 12] which depends only on the nodedegree distribution, while for α = 0 - by a quasi-diffusive dynamics reflecting the geographic population distribution.In Fig. 1 we show phase portraits for epidemics with different values of α on random geometric networks with N = 10 nodes. See Section IV B for details. As one can see from Fig. 1 the results of simulations for α = 1 . are very welldescribed the by the phase-portrait (8) of the classical SIR compartmental model [3, 12]. The number of infectiveagents is maximal at S max /N ≈ /R ≈ . and the herd immunity is achieved for S hi /N ≈ . − . , which is seenin the chart as the place were the curve hits zero. This value is close to the mean-field prediction (8). The value ofthe basic reproduction number, R = 2 . , of the best fit to the theoretical curve given by the mean-field solution (8)differs by one percent from the value R = 2 . used in the Monte-Carlo simulations. The difference can be attributedto the fact that the classical mean-field dynamics is deterministic [3, 12] and R is a number, while in the simulationsthe dynamics is stochastic and R is the mean value of a random variable. The variance of this random variableintroduces some corrections to the effective value of R .The phase portrait starts to deviate from the mean-field solution when α decreases (see Fig. 1) that is when long-range social mixing gets smaller. As shown in the right panel in Fig. 1 the phase portraits for different simulationsfor α = 1 lie on top of each other and are consistent with the classical SIR solution. For α = 0 the curves differ fromeach other and change in a stochastic way.The herd immunity value S hi weakly depends on α (see Fig.1). The values of S hi /N ≈ . − . are almostidentical 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 isillustrated in Fig. 2 where we compare dynamics of the epidemics for four different scenarios which differ by the basicreproduction number R and the long-range social mixing parameter α . One can see that the spread of epidemicdepends on the basic reproduction number R as well as the long range social mixing parameter α . The parameter α can be reduced by closing airports, schooles, churches, sport arenas, etc. The reproduction number can be loweredto a value R ′ < R through social distancing that is maintaining physical distance between people, reducing thefrequency of personal contacts, wearing masks as well as disinfection, quarantine, isolation, etc.In the next section we shall assess effects of these actions on the total mortality during epidemic using Monte-Carlosimulations of Covid-19-like epidemics. I N F E C T I O U S TIME(DAY)SCENARIO 1SCENARIO 2SCENARIO 3SCENARIO 4 0 20000 40000 60000 80000 100000 0 200 400 600 800 1000 1200 1400 1600 1800 I MM UN E TIME(DAY)SCENARIO 1SCENARIO 2SCENARIO 3SCENARIO 4
FIG. 2. The charts show the dynamics of epidemic for four scenario: (1) R = 2 . , α = 1 . ; (2) R = 2 . , α = 0 . ; (3) R = 1 . , α = 1 . ; (4) R = 1 . , α = 0 . . We use an agent-based implementation of SIR dynamics on geometric random networksdescribed in Section IV B. The population size is N = 10 , the mean node degree is h k i = 100 , the expected duration of theinfectious period is τ = 10 in all four cases presented in the figures. (Left) The number of infectious agents I ( t ) as a function oftime t expressed in days from the beginning of the epidemic. (Right) The number of immune agents: I ( t ) + R ( t ) = N − S ( t ) ,where I ( t ) , R ( t ) and S ( t ) are the numbers of infectious, recovered and susceptible agents, respectively. For scenario 1, the herdimmunity level is reached in t = 199 days. In scenarios 2,3,4, the herd immunity levels: , , are reached in: t = 398 , , days, respectively. B. Assessing mortality rate for different scenarios
We make the following assumptions in the model to incorporate these factors: • We assume that in the the studied period there is no vaccine and no medicine available so one can only applynon-pharmaceutic interventions to mitigate epidemic spreading. To be specific let us assume that this periodis days. The non-pharmaceutic interventions are divided into two classes: social distancing and reducinglong-range social mixing. The first strategy is implemented in the model by reducing the reproduction number R to some effective value R ′ which is smaller than the basic reproduction number R ′ < R , while the secondone - by lowering the value of the long-range social mixing parameter α ; • The basic reproduction number is R = 2 . . The mean infectious period is τ = 10 days. • The correlations between the occurrence of SARS and health conditions are modelled by dividing the populationinto a class of healthy people and a class of people who have some health issues, chronic diseases, or otherproblems requiring treatment and regular contacts with a doctor. The two classes are labeled in the model by H and C , respectively. We assume that the percentage of people in the H class is p H = 75% and in the C class p C = 25% . We assume also that SARS occurs with a probability p H,sars = 1 / in the H class and with aprobability p C,sars = 3 / in the C class. This gives one percent on average for the whole population. • The probability of dying of SARS is about in the H compartment and in the C compartment. Theprobability drops to about in the H class and to in the C class, if the person receives intensive medicalcare supported by mechanical ventilation. • The number of available mechanical ventilators is per million inhabitants. This sets the upper limit on thenumber of people who can be simultaneously ventilated; • The mortality rate of other diseases than those caused by the virus is per million per day and it increases to per million per day during the epidemic. In other words, the background mortality increases by about during the epidemic.The model is implemented by superimposing SIR dynamics on a background stochastic dynamics, see Section IV B.The background dynamics models the basic health and mortality statistics in the population. In our model it isimplemented as a Markov chain which preserves the total size of the population by replacing deaths by newborns.The stationary state of the Markov chain reproduces the statistics of the H and C classes and the mortality rate ofcauses not related to the virus, see Section IV C for details. The health care system capacity is modelled by a single -200 0 200 400 600 800 1000 0 200 400 600 800 1000 EX C ESS D EA T H S TIME(DAY)SCENARIO 1SCENARIO 2SCENARIO 3SCENARIO 4WORST CASENO EPIDEMIC -200 0 200 400 600 800 1000 0 200 400 600 800 1000 EX C ESS D EA T H S TIME(DAY)SCENARIO 1SCENARIO 5SCENARIO 6
FIG. 3. The cumulative number of excess deaths for different scenarios in Monte-Carlo simulations of epidemic for populationof N = 10 agents on random geometric network. The number is calculated as a difference between the cumulative numberof deaths and . t , where . corresponds to daily mortality rate in the population of hundred thousand people if there is noepidemic. (Left) The top curve corresponds to the worst case scenario, that is none of SARS patients receives medical assistanceduring the epidemic. The four curves below correspond to the scenarios 1-4 given in the main text. They also correspondto the four scenarios shown in Fig. 2. The curve in the bottom is an example of the background mortality data generatedby Monte-Carlo in the absence of epidemic. In this case the number of excess deaths fluctuates around zero. (Right) Thecumulative number of excess deaths for scenarios 1,5,6. One can see that the strategy of maintaining lockdown for some timeand then replacing it by the do-nothing strategy, does not significantly reduce the total number of deaths as compared to thesituation when the do-nothing strategy is introduced at the beginning of the epidemic. number, which accounts for the number of available respiratory ventilators. The ventilators are used to ventilate SARSpatients. The problem occurs when the number of SARS cases exceeds the number of available ventilators, becausesome people will not receive a required treatment and thus will have a lesser chance to stay alive. A ventilator isoccupied by the SARS patient until he or she recovers or dies and only then it can be used for another SARS patient.It may take ten days on average, so the number of occupied ventilators saturates quickly if there are more than new SARS cases per million inhabitants per day.We are now going to compare the total death toll of the simulated epidemic for six scenarios:1. R ′ = R = 2 . and α = 1 . . This simulates do-nothing scenario. Epidemic spreads without any restrictions.2. R ′ = R = 2 . and α = 0 . . This simulates a suppression of virus transmission through reducing long-rangesocial mixing.3. R ′ = 1 . < R and α = 1 . . This simulates social distancing and lowering the transmission rate.4. R ′ = 1 . < R and α = 0 . . This simulates lockdown. Both, the local and global transmission modes arerestricted.5. Lockdown for days, as in item 4, and then do-nothing scenario, as in item 1.6. Lockdown for days, as in item 4, and then do-nothing scenario, as in item 1.Fig. 3 shows results of Monte-Carlo simulations of the accumulated excess death toll in the epidemic as a functionof time for a population of people for days, for the six different scenarios. The first four scenarios areshown in the left panel in Fig. 3 and the remaining two in the right panel. In the left panel we additionally drawtwo reference curves representing the worse case scenario when no SARS patients receive required medical assistanceduring epidemic, and the best case when there is no epidemic. The charts in the right figure show what happens whenlockdown is introduced right at the beginning of the epidemic, maintained for or days, and lifted afterwards.As a reference we also show the curve for scenario 1 which is identical as in the left figure.Comparing the scenarios we see that after days there are excess deaths in scenario 1, in scenario 2, in scenario 3 and in scenario 4, so the lockdown strategy might seem to be optimal. However when one checks thetotal number of surplus deaths after days, which is , , and in scenarios 1-4, respectively, one cansee that lockdown is much less advantageous than strategies 2,3. The total death toll for scenario 4 is comparable asfor scenarios 2,3 but in this case the epidemic is not yet dying out in contrast to the cases 2,3 where the epidemic isover (compare Fig. 2). For scenario 4 the accumulated excess deaths curve is still rising after days. One can SA R S C ASES / VE N T I L A T E D TIME(DAY)SCENARIO 1: SARS CASESSCENARIO 1: VENTILATEDSCENARIO 2: SARS CASESSCENARIO 2: VENTILATEDSCENARIO 4: SARS CASESSCENARIO 4: VENTILATED 0 5 10 15 20 25 30 35 40 0 200 400 600 800 1000 D A I L Y M O R T A L I T Y TIME(DAY)SCENARIO 1SCENARIO 2
FIG. 4. (Left) The number of SARS patients and the number of occupied intensive care beds equipped with ventilators. Inthe simulations presented in the figure, the number of ventilators is limited to . The difference between the two curvescorresponds to the number of patients who will not receive adequate medical assistance. The figure shows results of Monte-Carlo simulations for scenarios 1,2,4. (Right) The number of deaths per day in scenario 1 and 2. For scenario 1 one can seea clear peak during the fast spread of epidemic. For scenario 2 one can also see a peak but it is wider and lower, and it onlyslightly exceeds fluctuations of the background mortality. We do not plot curves for other scenarios because the figures wouldbe less transparent. also see that maintaining lockdown for some time and then completely lifting it as in scenario 5 or 6 does not helpin reducing the excess mortality as compared to the do-nothing scenario. The best results are achieved by slowingdown the epidemic to some extent as in scenarios 2,3. These two strategies save lives of . of the whole population,that is about five thousand people per million. In Fig. 4 we show the number of SARS patients and the number ofventilated SARS patients during epidemic for various scenarios. The difference between the two curves corresponds tothe number of SARS patients who do not receive required medical assistance. In the right panel we show the numberof deaths per day. The number is a sum of deaths of SARS and of other causes. The former is roughly proportional tothe difference between the number of all SARS patients and ventilated SARS patients while the latter is proportionalto the mortality rate of other causes which is by ten percent higher during the simulated epidemic, see Section IV Cfor details.So far we have presented results of simulations for N = 10 . We end this section by noting that one obtainsqualitatively the same results for larger populations. As an example we compare simulations for N = 10 and N = 10 on geometric random networks with the mean degree h k i = 100 . The initial condition for SIR dynamicsis I = 5 infectious agents placed on randomly selected vertices of the network of size N = 10 . All others agents(vertices) are initially susceptible. The corresponding initial condition is I = 50 for N = 10 . Also other parametersare correspondingly rescaled, for example the number of available respiratory ventilators is ten times larger, theexpected number of deaths per day is ten times larger etc, to keep these numbers proportional to the population size.The results of the comparison are shown in Fig. 5. In the left panel we compare SIR dynamics. The figure showsthe fraction of infected people as a function of time. We see that the curves obtained in Monte-Carlo simulations for α = 1 are almost identical for N = 10 and N = 10 , as expected. This is the mean-field regime. For α = 0 epidemicspreading depends on geometric details of the network. Epidemic spreading has a quasi-diffusive character in thiscase. For the initial conditions that we described, the average distance between locations where epidemic initiallyoutbreaks is the same in the two cases, so one can expect that also the duration of epidemic is comparable. Andindeed, this is what one can see in Fig. 5 on the left hand side. The right panel shows the number of all SARSpatients and of ventilated ones, for N = 10 . The curves display a similar pattern as the one for N = 10 shown inFig. 4, except that the statistical fluctuations are now smaller due to the law of large numbers. III. DISCUSSION
We have conducted a Monte-Carlo study of epidemic spreading on random geometric networks to assess efficiencyof non-pharmaceutic interventions in reducing the total number of surplus deaths during Covid-19-like epidemics. Wehave discussed strategies based on social distancing and restricting long-range social mixing. They have different effectson epidemic spreading. Social distancing reduces the basic reproduction number R to some effective reproductionnumber R ′ < R . Restrictions on long-range social mixing reduce virus transmission between remote places. When I N F E C T I O U S F R A C T I O N TIME(DAY)SCENARIO 1: 10E6SCENARIO 1: 10E5SCENARIO 2: 10E6SCENARIO 2: 10E5 0 100 200 300 400 500 600 700 0 200 400 600 800 1000 SA R S C ASES / VE N T I L A T E D TIME(DAY)SCENARIO 1: SARS CASESSCENARIO 1: VENTILATEDSCENARIO 2: SARS CASESSCENARIO 2: VENTILATEDSCENARIO 4: SARS CASESSCENARIO 4: VENTILATED
FIG. 5. (Left) Comparison of infectious fractions in the Monte-Carlo simulations of SIR dynamics for α = 1 and α = 0 ongeometric random networks with N = 10 and N = 10 nodes. Initially there are I = 5 infected nodes on random locationsof the network with N = 10 nodes and I = 50 infected nodes on random locations of the network with N = 10 nodes. Theresults for α = 1 are almost identical for N = 10 and N = 10 , so the two curves overlap. (Right) The number of SARS casesand the number of SARS ventilated cases for N = 10 . Compare with the left panel in Fig. 4. long range mixing is large, the epidemic spreads via mean-field dynamics, when it is small it spreads via quasi-diffusive dynamics depending on the geographic population distribution. We have studied the two modes of diseasetransmission.There are two sources of deaths which contribute to the total death toll during a Covid-19-like epidemic. One isrelated to SARS, caused by the virus, and the other one to other diseases. The number of SARS deaths depends onthe number of respiratory ventilators. If the daily number of new SARS cases exceeds V /τ where V is the number ofventilators and τ is the number of days of using one ventilator for one SARS patient, some people with SARS willnot be ventilated and will have less chance to survive. If one assumes that there are ventilators per million, and τ = 10 , then V /τ = 27 . As long as the number of SARS patients is below per million then the number of deathscaused by SARS is minimal. This effect can be achieved by slowing down epidemic. On the other hand the numberof excess deaths of other causes will increase proportionally to the epidemic duration so it is not beneficial to slowit down 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. In regions, where the number of available ventilators is very small and V /τ iscomparable to the excess mortality of other diseases the best solution is to do nothing. For example, in one assumesthat there are excess deaths of other causes during epidemic per day per million and there are less then V = 30 ventilators per million then the do-nothing approach is optimal. For the same reason it is not very beneficial to slowdown the epidemic to the extent that there are only a few SARS cases per day. In this case the epidemic will last along time and there will be many surplus deaths of other causes. We have also shown that a strategy of maintaininglockdown for some time and then releasing it by removing all restrictions has a similar effect on the number of deathsas if one introduced the do-nothing strategy right at the beginning. The deaths differ only by the time when theyoccur: in the do-nothing scenario the mortality is large at the beginning of the epidemic while in the other case it islarge when the lockdown is released. A strict lockdown makes sense only if one expects that an effective medicine ora vaccine are going to be approved in a short time, or if one needs some time to increase the capacity of the healthcare system by buying new respiratory ventilators and training people, etc. Otherwise, the optimal strategy is to keepthe epidemic progress at the level that the number of SARS patients at any time is roughly equal to the number ofrespiratory ventilators. If the number of SARS cases is much larger than that, too many people will die of SARS; ifit is much smaller than that, the epidemic will last too long and as a result too many people will die of other causes.For the simulated epidemic the strategies based on either social distancing or reducing long range social mixing(scenarios 2,3) give results close to the optimal one. Indeed, as one can see in Fig. 4 the number of SARS deathsin the peak of epidemic in scenario 2 does not significantly exceed the background mortality. In reality, however, itis much more difficult to tune parameters that control the spread of the epidemic and to find a way to implementappropriate actions in society.Social distancing reduces the effective reproduction number to a value R ′ < R below the basic reproductionnumber. The herd immunity level for R ′ is smaller than the herd immunity level fo R . This means that after liftingrestrictions on social distance and restoring normal contacts between people, the percentage of immune people will bebelow the herd immunity level. The system will be unstable, in the sense that a new single infective person may triggera new outbreak of epidemic. The situation is similar to superheated liquid, where boiling may occur spontaneouslyat any time. For example, the total number of deaths in the simulated epidemic is comparable for scenarios 2 and 3(see Fig. 3), but the percentage of immune people at the end of epidemic is in scenario 2 and in scenario 3,as one can see in Fig. 2. The value is close to the herd immunity level for R = 2 . while is far below it.This means that the epidemic in scenario 3 can restart from the level when the restrictions are lifted and a newinfective person appears. This example shows that strategies reducing long range social mixing bring better effectsthan social distancing. They are however much more difficult to implement.The model that we described in this paper can be developed in many directions. One can introduce a full MSEIRdynamics, one can expand the background dynamics by implementing age groups, each having own distribution ofhealth conditions and death statistics, one can introduce social mixing matrices for different public places where longrange social mixing takes place, one can use the meta-population approach to model in detail the long range socialmixing mechanisms and one can also introduce some temporal aspects. IV. METHODSA. Random geometric networks
Random geometric networks are constructed by the proximity rule [32, 33]. Two nodes are connected by anedge if they lie within the given distance from each other. The simplest example is a network constructed byconnecting randomly distributed points in a d -dimensional Euclidean space. One can apply this construction toimitate a geographical proximity network in two-dimensions to model a network of contacts between people. Forsake of simplicity we assume that the points are uniformly distributed on a two-dimensional square with the periodicboundary conditions. This can be done by generating pairs of coordinates ( x i , y i ) i = 1 , . . . , N consisting of N independent random numbers uniformly distributed on the unit interval [0 , and connecting any two points i and j by an edge of the network if the distance between them is smaller than ǫ : ∆ x ij + ∆ y ij ≤ ǫ . For the periodic boundaryconditions the coordinate differences are calculated as follows ∆ x ij = min ( | x j − x i | , − | x j − x i | ) and analogouslyfor ∆ y ij . The node degree distribution of the network obtained in this way follows the binomial law P ( k ) = (cid:18) N − k (cid:19) a k (1 − a ) N − − k , (1)where a = πǫ is the area of a circle of radius ǫ . The mean degree distribution is h k i = ( N − a , and the variance σ ( k ) = ( N − a (1 − a ) . When a is of order /N the distribution becomes Poissonian in the large N limit. Theaverage clustering coefficient of the network is large h C i = 1 − √ π ≈ . . [32]. This means that neighbours ofa node are directly connected with a very high probability. This feature of geometric random networks constitutesthe main difference to Erdőes-Rényi random graphs [34], for which the clustering coefficient is of order /N . This hasalso consequences for epidemic spreading on such networks. B. Agent-based implementation of SIR dynamics
We use a discrete time stochastic implementation of the standard SIR dynamics [3, 12]. The network is populatedwith agents residing on its nodes. The population is divided into three classes: susceptible (S), infectious (I) andrecovered (R) which describe the state of each agent at time t . The states change in the course of evolution accordingto epidemic rules which are implemented in the model in the form of a discrete time stochastic process. Time iscounted in days from the outbreak of the epidemic. Initially, that is for time t = 0 , one agent, or a few ones areinfectious, while all others are susceptible. An infectious agent remains infective for τ days on average, and then itrecovers. This is simulated in the model by assuming that the probability of remaining infective till next day is q andof recovering is − q . The life time distribution of infectious state is given by the following exponential law P i ( t ) = (1 − q ) q t − , t = 1 , , . . . (2)The mean life time of infectious state is related to the probability q as follows τ = h t i = ∞ X t =1 tP i ( t ) = 11 − q (3)which means that for q = τ − τ (4)the expected infectious period is τ days. Clearly, for τ ≫ the probability distribution (2) can be approximated by P i ( t ) ≈ e − t/τ /τ . Once an infective person recovers he or she remains recovered for the rest of the evolution. If asusceptible agent has a contact with an infective one and this contact is sufficient for disease transmission the agentbecomes infective. In the model there are two different transmission modes: a local one – between neighbouring nodeson the network, and a global one – between randomly selected nodes independently of their locations on the network.The probability of the global transmission is α and of the local one − α . For α = 0 the spread of epidemic is entirelyshaped by local properties of the network while for α = 1 only on average properties of the whole network. The globalmode is well captured by the classical mean-field dynamics [3, 12].Let p be a probability that the disease is transmitted from an infectious agent to a susceptible agent within oneday. The probability p t of a transmission within t days is p t = 1 − (1 − p ) t . The life time of an infectious state is arandom variable (2). The transmission probability for the whole life-time can be calculated as the expected value of h p t i = 1 − h (1 − p ) t i for the random variable t (2). This yields h p t i = τ p p ( τ − (5)for q given by (4). An agent has on average h k i neighbours, so the number of infections generated by a single infectivein a fully susceptible population is R = h k i τ p p ( τ − . (6)This equation relates the basic reproduction number R to the parameters of the model.The epidemic evolution is implemented in a synchronous way in the model which means that all states are updatedsimultaneously. States at time t + 1 are computed from states at time t . The following rules are used to update thestates. If a node is recovered at time t it remains recovered at time t + 1 . If a node is infectious at time t it remainsinfectious at time t + 1 with a probability q . Otherwise, it changes to recovered. If a node is susceptible at t it changesto infectious with a probability p ∗ . Otherwise it remains susceptible. The probability p ∗ of becoming infectious isrelated to the transmission probability p by the following relation p ∗ = 1 − (1 − p ) i ∗ where i ∗ is an effective numberof infectious neighbours i ∗ = (1 − α ) i n + α h k i IN . (7)In the local transmission mode, that is for α = 0 , the quantity i ∗ is proportional to the number of direct infectiousneighbours i n . In the global transmission mode, that is for α = 1 , the quantity i ∗ is proportional to all infectiousnodes on the network h k i I/N . When α changes from zero to one, i ∗ interpolates between the two limiting cases.In the classic approach one usually uses the continuous time formalism [3, 12]. The epidemic evolution is describedby a set of first order ordinary differential rate equations for the fractions of susceptible, infectious and recoveredagents: s ( t ) = S ( t ) /N , i ( t ) = I ( t ) /N , r ( t ) = R ( t ) /N . The epidemic outbreaks if s (0) R > . The quantity φ ( t ) = i ( t ) + s ( t ) − R ln s ( t ) = const (8)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 s ( t max ) = 1 .Indeed, as one can see from eq. (8), the derivative di/ds = − R 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 ) → when t → ∞ . Thefraction of susceptible population for t → ∞ gives the level of herd immunity s ( t ) → s hi . The value s hi can be foundfrom Eq. (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 ( s hi − . This yields s hi ≈ . , . , . , . for R = 1 . , , . , , respectively, to givesome examples.We use the following input parameters in the Monte-Carlo simulations of the epidemic on geometric randomnetworks: the number of agents N , the mean node degree h k i , the basic reproduction number R , the expectedduration of the infectious period τ , the probability α of long range transmission. As an initial configuration we choosea configuration which consists of I randomly selected infectious nodes. The remaining nodes are susceptible. Theprobability to remain infectious till the next day is calculated from Eq. (4). The probability of virus transmissionfrom an infectious agent to a susceptible one within one day is calculated from Eq. (6) which gives p = 1 τ (cid:16) h k i R − (cid:17) + 1 . (9)An example of input values used in the simulations is N = 10 , h k i = 100 , R = 2 . , τ = 10 , α = 0 , I = 5 .0 C. Modelling background conditions
In order to asses the impact of epidemics on death statistics one has also to determine the death statistics and thebackground conditions in the absence of epidemic. This is per se an interesting and very complex problem since itinvolves demographic factors, efficiency of health care systems, statistics of diseases, and many other factors. Thisis beyond the scope of this paper. Here we restrict ourselves to model only primary factors which may help tounderstand how the death statistics changes during epidemic. In the model this is implemented as follows. Thepopulation is divided into classes according to health conditions. In the simplest version of the model we introduceonly two classes which 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 including in the model statistical correlations between healthconditions and mortality. The mortality rate in the C class is significantly higher than in the H class. The secondimportant difference 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 section.We assume that the size of population is constant during the epidemic. The number of deaths is compensated bythe number of newborns. Denote the fraction of healthy people at time t by h ( t ) = H ( t ) /N , the fraction of chronicallyill people by c ( t ) = C ( t ) /N and the fraction of deaths by d ( t ) = D ( t ) /N . We have h ( t ) + c ( t ) + d ( t ) = 1 .We implement the population dynamics as a discrete time stochastic process (Markov chain) with the followingevolution equation ( h ( t + 1) , c ( t + 1) , d ( t + 1)) = ( h ( t ) , c ( t ) , d ( t )) p HH p HC p HD P CH p CC p CD p DH p DC p DD . (10)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 equalto the number of newborns. The parameter p DH is the probability that newborns will be healthy at the beginning.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 much smaller than the mortality rate of chronicallyill people. The condition p DC = 0 means that a death is replaced by a healthy newborn. The condition p CH = 0 means that a chronically ill person does not become healthy again. Under these assumptions the last equation can besimplified to ( h ( t + 1) , c ( t + 1) , d ( t + 1)) = ( h ( t ) , c ( t ) , d ( t )) − β β
00 1 − γ γ . (11)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 h ∗ = γβ + γ + βγ ,c ∗ = ββ + γ + βγ ,d ∗ = βγβ + γ + βγ . (12)We assume that d ∗ = 27 · − . The value d ∗ is adjusted to a typical mortality rate in European countries which is per million per day. In most of our simulations we have h ∗ ≈ . and c ∗ ≈ . , which means that the simulatedpopulation consists in of healthy people and in of people with chronic, cardiovascular, oncologic or otherserious diseases. The corresponding parameters of the transfer matrix (11), which reproduce the stationary state, are β ≈ d ∗ /h ∗ = 36 · − and γ ≈ d ∗ / (1 − h ∗ ) = 108 · − .We conclude this section with two remarks. 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 effectivetransfers rates between these sectors are non-zero already for two days − β β
00 1 − γ γ = (1 − β ) β (2 − β − γ ) βγγ (1 − γ ) γ (1 − γ )1 − β β . (13)1The second remark is that the square or a higher power of the transfer matrix (11) is also a stochastic matrix. Inprinciple one can replace the original transfer matrix by 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 associatedwith the eigenvalue and it is identical for the transfer matrix (11) or any power of it. The transfer matrix (11) hasthree eigenvalues. The one which has the largest absolute value is λ = 1 and the second largest is λ ≈ − β − γ . Theeigenvalue λ tells us about correlations of states at different times t , t ′ . The correlation function decays exponentiallyas exp( −| t − t ′ | /T ) . The correlation time T can be derived from λ : T ≈ − / log( λ ) ≈ / ( β + γ ) . For the transfermatrix (11) T is of order . By raising this matrix to the n -th power and interpreting the resultant matrix as adaily transfers matrix one can reduce the autocorrelation time from T to T /n . D. Simulating death statistics during epidemic
We assume that the infection leads to SARS with a probability p H,sars = 1 / in the H compartment and with aprobability p C,sars = 3 / in the C compartment. This yields on average p sars = h ∗ p H,sars + c ∗ p C,sars ≈ /
100 = 1% (14)where h ∗ ≈ / and c ∗ = 1 / , as described in the previous section. Because p C,sars is nine times larger than p H,sars this imitates a correlation between the occurence of SARS and the health conditions of the person. In the algorithmthis is implemented by moving new members of the H and C classes with the probabilities p H,sars and p C,sars to theSARS risk group. Members of this group will develop SARS after infection.Another factor which plays an important role in the death statistics during epidemics is the fatality rate for SARS.The rate depends on the age and the general health conditions of the patient, and on whether the patient is ventilatedor not when SARS develops. In the model we distinguish four situations, labeled by C , C , H , and H : • C : Patients with SARS from the C class who are not ventilated • C : Patients with SARS from the C class who are ventilated • H : Patients with SARS from the H class who are not ventilated • H : Patients with SARS from the H class who are ventilatedand we set probabilities of dying of SARS to . , . , . , . for the C , C , H and H groups, respectively. In thesimulations we use daily SARS mortality rates, which are related to the probabilities of dying in the whole period byan equation identical to Eq. (5) in which h p t i is interpreted as the probability of dying of SARS during the wholeperiod and p is the probability of dying within one day. Assuming that τ = 10 , the corresponding daily rates are . , . , . , . for the C , C , H , H compartments, respectively. This choice of parameters corresponds to theassumption that respiratory ventilation increases the probability of staying alive from to for people with SARSin the C class, and from to for people with SARS in the H class. During epidemic the number of people withSARS may easily exceed the number of ventilators available. We assume that there are V = 270 ventilators (intensivecare beds) per million inhabitants as in some European countries. A patient with SARS occupies a ventilator untilhe or she dies or recovers. This takes takes ten days on average. So one can expect that if there are more than new SARS cases per day, the demand for ventilators will exceed the health care capacity.The ventilator availability is simulated as follows. At any moment of time we keep track of the number of availableventilators. If this number is larger than zero, and there is a new SARS case, the number is decreased by one, andthe SARS patient is moved from the C to C or H to H compartment, respectively. The ventilator is occupieduntil the patient recovers or dies, in which case the number of available ventilators is increased by one. Initially, thenumber of ventilators is set to per million.Another factor which has to be taken into account in assessing the epidemic total death toll is a lower efficiencyof the healthcare system during epidemic. This has an impact on the increase of deaths of other causes which arenot directly related to SARS. The effect is more pronounced in the group of people with chronic diseases who requirecontinuous medical assistance, but it is also seen in the group of healthy people because regular medical proceduresare slowed down due to epidemic precautions in hospitals. One would have to conduct systematic statistical surveysto estimate this effect. Here we just assume that the number of deaths of other causes than those directly relatedto SARS increases from to per million per day during epidemic - roughly by . We implement this effectby changing values of the parameters β and γ in the transfer matrix (11) from β = 36 · − and γ = 108 · − to β = 40 · − and γ = 120 · − , during epidemic.2 ACKNOWLEDGMENTS
I thank Krzysztof Malarz for interesting discussions and assistance in preparing plots. [1] N. T. Bailey,
The Mathematical Theory of Infectious Diseases , 2nd ed. (Hafner: New York, 1975).[2] R. M. Anderson and R. M. May,
Infectious Diseases of Humans: Dynamics and Control (Oxford University Press: Oxford,UK, 1992).[3] H. W. Hethcote,
The mathematics of infectious diseases , SIAM Review , 599–653 (2000).[4] M. Y. Li, An Introduction to Mathematical Modeling of Infectious Diseases (Springer International Publishing AG: ChamSwitzerland, 2018).[5] N. M. Ferguson, et al.,
Strategies for containing an emerging influenza pandemic in southeast asia , Nature , 209–214(2005).[6] N. M. Ferguson, et al.,
Report 9: Impact of non-pharmaceutical interventions (npis) to reduce covid-19 mortality andhealthcare demand , Preprint at Spiral (2020).[7] S. Flaxman, et al.,
Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe , Nature10.1038/s41586-020-2405-7 (2020).[8] D. Bernoulli,
Essai d’une nouvelle analyse de la mortalité causée par la petite vérole et des avantages de l’inoculation pourla prévenir (Académie Royale des Sciences: Paris, 1760) pp. 1–45.[9] K. Dietz,
The first epidemic model: A historical note on p. d. en’ko , Austral. J. Statist. , 56–65 (1988).[10] W. Hamer,
Epidemic disease in england , Lancet , 733–739 (1906).[11] R. Ross, The Prevention of Malaria , 2nd ed. (Murray: London, 1911).[12] W. O. Kermack and A. G. McKendrick,
A Contribution to the Mathematical Theory of Epidemics , Proceedings of theRoyal Society of London Series A , 700–721 (1927).[13] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani,
Epidemic processes in complex networks , Rev.Mod. Phys. , 925–979 (2015).[14] A.-L. Barabási and R. Albert, Emergence of scaling in random networks , Science , 509–512 (1999).[15] R. Albert and A.-L. Barabási,
Statistical mechanics of complex networks , Reviews of Modern Physics , 47–97 (2002).[16] S. N. Dorogovtsev and J. F. F. Mendes, Evolution of networks , Advances in Physics , 1079–1187 (2002).[17] M. E. J. Newman, The structure and function of complex networks , SIAM Review , 167–256 (2002).[18] M. Barthélemy, Spatial networks , Physics Reports , 1–101 (2011).[19] G. Chowell, J. M. Hyman, S. Eubank, and C. Castillo-Chavez,
Scaling laws for the movement of people between locationsin a large city , Phys. Rev. E , 066102 (2003).[20] V. Colizza, A. Barrat, M. Barthélemy, and A. Vespignani, The role of the airline transportation network in the predictionand predictability of global epidemics , Proceedings of the National Academy of Sciences , 2015–2020 (2006).[21] D. Balcan, V. Colizza, B. Gonçalves, H. Hu, J. J. Ramasco, and A. Vespignani,
Multiscale mobility networks and the spatialspreading of infectious diseases , Proceedings of the National Academy of Sciences , 21484–21489 (2009).[22] R. Pastor-Satorras and A. Vespignani,
Epidemic spreading in scale-free networks , Physical Review Letters , 3200–3203(2001).[23] A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press:Cambridge, UK, 2008).[24] V. Colizza, R. Pastor-Satorras, and A. Vespignani,
Reaction-diffusion processes and metapopulation models in heterogeneousnetworks , Nature Physics , 276–282 (2007).[25] M. Bootsma and N. Ferguson, The effect of public health measures on the 1918 influenza pandemic in us cities , Proceedingsof the National Academy of Science , 7588–7593 (2007).[26] Cauchemez, S. et al,
Household transmission of 2009 pandemic influenza a (h1n1) virus in the united states , New EnglandJournal of Medicine , 2619–2627 (2009).[27] P. Bajardi, et al.,
Human mobility networks, travel restrictions, and the global spread of 2009 h1n1 pandemic , PLoS ONE , e16591 (2011).[28] E.H. Otete, et al.,
Parameters for the mathematical modelling of clostridium difficile acquisition and transmission: Asystematic review , PLoS ONE , 10.1371/journal.pone.0084224 (2013).[29] M. Chinazzi, et al., The effect of travel restrictions on the spread of the 2019 novel coronavirus (covid-19) outbreak , Science , 395–400 (2020).[30] A.J. Kucharski, et al.,
Early dynamics of transmission and control of covid-19: a mathematical modelling study , The LancetInfectious Diseases , 553 – 558 (2020).[31] J. T. Wu, et al., Estimating clinical severity of covid-19 from the transmission dynamics in wuhan, china , Nature Medicine , 506–510 (2020).[32] J. Dall and M. Christensen, Random geometric graphs , Phys. Rev. Lett. , 016121 (2002).[33] M. Penrose, Random Geometric Graphs (Oxford University Press: Oxford, UK, 2003).[34] P. Erdős and A. Rényi,
On random graphs , Publicationes Mathematicae6