A Full-Scale Agent-Based Model of Lombardy COVID-19 Dynamics to Explore Social Networks Connectivity and Vaccine Impact on Epidemic

COVID-19 outbreak is an awful event. However it gives to the scientists the possibility to test theories about epidemic. The aim of this contribution is to propose a individual-based model of Lombardy COVID-19 outbreak at full-scale, where full-scale means that will be simulated all the 10 millions inhabitant population of Lombardy person by person, in a commercial computer. All this to test the impact of our daily actions in epidemic, investigate social networks connectivity and in the end have an insight on the impact of an hypothetical vaccine.

COVID-19 · Epidemic model · Agent based · Individual based · vaccine · Social networks 1 Introduction The first alarm from COVID-19 has arrived from China [1], but one of the most dramatic outbreak was in Italy at the end of January [2]. This epidemic was the setting of a change in risk management: the use of mathematical modeling [3]. In this direction must be mentioned the paper of [ 4], where is explored the COVID-19 outbreak in agent-based model of about 750 thousands inhabitants of the city of Urmia, with the movement of agents approximated by places. On the same direction is the contribution [5], where is taken in account a model of transmission with a sub-sampled population of 9 thousands people living in Daegu. But there are also small-scale models as [6]. Must be also mentioned the model developed by University of Palermo [7] based on the paper [8]. The aim is to introduce a qualitative full-scale agent-based model able to reproduce the Lombardy COVID-19 dynamics, modeling its outbreak and descent, including as much real data as possible. Being Lombardy population about 10.06 millions makes this model very large scale in comparison with previous contributions. In second instance will be taken in account several alternative scenarios, in order to weight the events of those days. Then the model of social interaction used in epidemiological simulations, will be explored in terms of graph theory [9], in order to analyze it as a social network [10]. In conclusion will be drawn several conclusions about the impact of our habits in COVID-19 outbreak. 2 Methods 2.1 Model structure The key idea of the model is to create a three layers model, as can be seen in figure 1. The first layer is an agent based particle model of Lombardy. Every agent is an inhabitant of the region, making this model a full-scale model of the Italian region. Then we have 10.06 millions agents that move with a random walk [11]. This little detail is part of the novelty of the contribution, because (at the best of author’s knowledge) it is the first try to simulate a so large population individual by individual for this purpose. . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 14, 2020. ; https://doi.org/10.1101/2020.09.13.20193599doi: medRxiv preprint NOTE: This preprint reports new research that has not been certified by peer review and should not be used to guide clinical practice. A PREPRINT - SEPTEMBER 13, 2020 Environment + Agents Collision Detection Epidemic Model R S S I D Figure 1: Three layers structure of the model. On the agents model is built a collision detection algorithm, able to detect if two particles have distance less than a fixed value. However the large-scale reached by the model will require an ad hoc algorithm for this purpose, that will be a challenging problem. On the collision detection model is built the epidemic model, that is a SIRD (Susceptible-Infected-Recovered-Dead) model [12]. The model has been filled with many experimental findings on COVID-19 and some fitting parameters are fitted on experimental findings. In the whole analysis have been used as much open access contributions as possible, moreover the whole project is fully open source at the site [13]. 2.2 Agents model The model simulates the behavior of each inhabitant of Lombardy using the approximation of Random walk [11]. Then every particle moves with a random vector at every step. In particular the model runs at 6 Frames per days, that is a good frame rate considering that the scale time for Epidemic phenomena is usually months. The approach of random walk can appear unrealistic, but this hypothesis ha shown to be appropriate to model very large scale systems, like in gas thermodynamics [14]. In addition to the random walk, has been added a weak velocity field with a dependence of 1 r2 , where r is the distance between two particles, like in gravitational field, in order to aggregate particles. The drift speed of the particles is constant and selected with a Weibull distribution [15] with scale parameter 6 and shape parameter 1.5. The particles speeds has been adjusted through a multiplicative constant in order to make the average path length of a particle in a day about 43 Km, as suggested in the report [16]. 2.3 Environment model The setting of the simulation is Lombardy, then the environment model will be a closed 2D box with boundary shape following Lombardy border. In order to keep particles inside the region, has been introduced a bouncing condition at the border, so that a particle that tries to cross the border bounces backward. This condition is very popular in gas thermodynamics [11]. The initial conditions of particles in terms of position are distributed following the actual density of population of Lombardy, extracted from the image [17] with techniques of image analysis and then must be intended as an approximation of the real data. 2 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 14, 2020. ; https://doi.org/10.1101/2020.09.13.20193599doi: medRxiv preprint A PREPRINT - SEPTEMBER 13, 2020 lock down beginning lock down beginning lock down beginning Figure 2: Simulation of COVID-19 Outbreak. (Top-Left) Density of Population. (Top-Right) log10 of infected percentage per cell. (Bottom) From left to right: Infected number, Recovered number, Deaths number and Recover ratio ( Recovered Deaths ). Solid line is model simulation, dotted line is extracted data by [18] for Lombardy and vertical dotted blue line is the date of DPCM 9 March [19]. 2.4 Collision detection Starting from the assumption that the algorithm has been designed to run on a commercial computer in parallel (the one used has a AMD 3900X 12-cores and 64 GB of RAM) and in reasonable time (about 20 minutes of calculation for 14 days of simulation), the collision detection algorithm had a central role in the implementation of the algorithm. Because to find all the points with distance less than a constant in a set ofN points, in general has complexity order of N 2. In our model N ≈ 107, then the complexity is order of 1014, that’s an incredible number! Then the idea is to subdivide the Lombardy with a grid with thick about 20 Km and in every cell of the grid apply the collision detection, all this in parallel through cells. This multi-scale processing has lead to a speed up of the code, reducing at the same time the RAM used simultaneously in computation, making possible a simulation with 10 millions particles at the same time. However this approach neglects all the connections across the borders of the cells, the connections detected fraction (if the detection distance is about 1 Km) has been estimated computationally on a small scale experiment to be about 97.8% [13], that is a good percentage in a qualitative model. However this percentage can be increased using a multiple grids detection scheme, but it is beyond the aims of this contribution. 2.5 Epidemic model The Epidemic model is a SIRD (Susceptible-Infected-Recovered-Dead) model [3]. Most of the models available up to now are called population models [20]. A population model is a model where every node is modeled by a set of differential equations and it models a sub-population of a region. The number of people modeled by a single node can range from hundreds to millions. In our model, every node is a single person. So the model is not an ODE (Ordinary Differential Equations) model, but a stochastic agent-based model. Every node has 4 states: 1. Susceptible: it hasn’t already contracted the disease. It can be become Infected with a probability pI for each contact with an Infected; 2. Infected: it is infected and then it can infect Susceptible nodes. After E days will change its state to Recovered or to Dead, in particular it has probability pD to die and 1 − pD to be recovered; 3. Recovered: it has been healed from the disease and it can’t contract it either infect Susceptible nodes any more; 4. Dead: the node has died, so it can’t infect any more. 3 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 14, 2020. ; https://doi.org/10.1101/2020.09.13.20193599doi: medRxiv preprint A PREPRINT - SEPTEMBER 13, 2020 lock down beginning lock down beginning lock down beginning Figure 3: Simulation with social distancing. (Top-Left) Density of Population. (Top-Right) log10 of infected percentage per cell. (Bottom) From left to right: Infected number, Recovered number, Deaths number and Recover ratio (Recovered Deaths ). Solid line is model simulation, dotted line is extracted data by [18] for Lombardy and vertical dotted blue line is the date of DPCM 9 March [19]. 3 Results 3.1 Lombardy Outbreak First of all I want to point out that all the simulations are available in .avi format at the link [13], as well the MATLAB code to reproduce them. Then the first scenario is the Lombardy outbreak of March 2020 [2]. Our simulation starts 29 February 2020 and terminates 14 March 2020. The main realistic parameters are pI = 1 40500 (extracted from the report [21]) and pD = 0.3 (estimated from Worldometer [22] that has been cited also in [23]). The fitted parameters are the collision radius of 1 Km (that can appear enormous comparing to the 1 m distance suggested by WHO [ 24], but must be taken in account that there are 6 frames per day, then 1 Km is the radius of interaction of a person that stays in the same place for 4 hours) and the days of disease E = 7 ([25] suggests E < 10). The results of simulations can be seen in figure 2. The model is able to explain the experimental data until about 9 March 2020. What happened after this date belongs to history. Indeed 9 March was the day of the DPCM "Misure urgenti in materia di contenimento e gestione dell’emergenza epidemiologica da COVID-19" [19], a law in order to contain the COVID-19 outbreak. This law was the beginning of the Italian lock down. In this discrepancy between the data and the model, can be seen all the efforts of Italian people to block the Coronavirus. This simulation also warns us on what would have could happened. 3.2 Our habits make the difference This second scenario is inspired from the paper [26]. In this paper the authors show that keeping a 2 meters distance between people halves the risk of contracting COVID-19. Then the aim of this simulation is to simulate this kind of social distancing by halving pI in the model. The results, in figure 3, is that the COVID-19 (in this particular scenario) isn’t enough contagious to spread as in the experimental data. This simple simulation shows the striking role that we have in fighting this disease, indeed a simple action, like to keep a social distance, can make the difference between a virus under control and an epidemic. This simple fact has already been observed in experimental finding about Germany [27], where a synthetic method was used to estimate the contagion without the use of masks. On the same path has been performed a simulation of lock down, reducing the daily average kilometers done by a node from 43 to 5 and reducing the interaction distance from 1 Km to 100 m. The results of this simulation can be seen in figure 4. These simple assumptions were enough to defeat the virus. 4 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 14, 2020. ; https://doi.org/10.1101/2020.09.13.20193599doi: medRxiv preprint A PREPRINT - SEPTEMBER 13, 2020 lock down beginning lock down beginning lock down beginning Figure 4: Simulation of lock down. (Top-Left) Density of Population. (Top-Right) log10 of infected percentage per cell. (Bottom) From left to right: Infected number, Recovered number, Deaths number and Recover ratio ( Recovered Deaths ). Solid line is model simulation, dotted line is extracted data by [18] for Lombardy and vertical dotted blue line is the date of DPCM 9 March [19]. Figure 5: COVID-19 Outbreak simulation connectivity. (Top-Left) Density of Population. (Top-Right) log10 of group percentage per cell. (Bottom-left) Degree distribution of small group. (Bottom-Right) Daily degree distribution of small group. 5 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 14, 2020. ; https://doi.org/10.1101/2020.09.13.20193599doi: medRxiv preprint A PREPRINT - SEPTEMBER 13, 2020 Figure 6: Lock down simulation connectivity. (Top-Left) Density of Population. (Top-Right) log10 of group percentage per cell. (Bottom-left) Degree distribution of small group. (Bottom-Right) Daily degree distribution of small group. 3.3 Network Topology The impact of topology in an epidemic model is an hot topic [28] in the debate about social networks. For this reasons has been performed a test, one thousands particle has been chosen and then tracked in all the simulation to find the total number of connection with the whole population. In graph theory, the number of connection of a node is called degree [9]. With this test has been determined the degree distribution and the daily degree distribution (average degree per day) of this small group of people varying in time. However here will be presented just the final result, for the full simulation refer to [13]. The first scenario is the COVID-19 Outbreak scenario, figure 5. Can be seen that the distribution has an evident left tail (in contrast with the right tail of BA models [10]). Probably this fact emerges from the simulation time of 14 days (in contrast with the social networks, that usually takes years to be built). It is also interesting the scenario of lock down. In this scenario we have a fall of the connectivity from thousands average connections per day to order of hundred. This results shows the importance of the policy of lock down in COVID-19 containment. 3.4 The descent In this scenario will be taken in account the period between 31 May 2020 and 14 June 2020. In this period Italy is concluding its lock down and the active cases number is decreasing. For this simulation the kilometers per day has been arbitrary set to 15, because of the lack of additional information about mobility during this period. The probability of contagion has been halved to take in account the social distancing. The radius of interaction and the duration of the disease has been tuned to reproduce experimental data. In particular the value for radius of interaction is 300 m and as disease duration has been taken E = 35, so 5 weeks. This value (so higher in comparison to the one of the outbreak) could be motivated by a clinical protocol more accurate and by the queue created by the large number of infected, that could slow down the exams required to be declared recovered. The qualitative fitting can be seen in figure 7. 3.5 The vaccine In the previous scenario of descent has been tested the impact of a vaccination of the 70% of population, near to the 62% suggested in [29]. The agents-based models are perfect to test strategies on the individuals. The result of the simulation is a strong decreasing of infected trend, that is something unexpected in a simulation of 14 days. The results are shown in figure 8. 6 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 14, 2020. ; https://doi.org/10.1101/2020.09.13.20193599doi: medRxiv preprint A PREPRINT - SEPTEMBER 13, 2020 Figure 7: Simulation of descent. (Top-Left) Density of Population. (Top-Right) log10 of infected percentage per cell. (Bottom) From left to right: Infected number, Recovered number, Deaths number and Recover ratio ( Recovered Deaths ). Solid line is model simulation, dotted line is extracted data by [18] for Lombardy. number of vaccinations Figure 8: Simulation of vaccination. (Top-Left) Density of Population. (Top-Right) log10 of infected percentage per cell. (Bottom) From left to right: Infected number, Recovered number, Deaths number and Recover ratio ( Recovered Deaths ). Solid line is model simulation, dotted line is extracted data by [18] for Lombardy. 7 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 14, 2020. ; https://doi.org/10.1101/2020.09.13.20193599doi: medRxiv preprint A PREPRINT - SEPTEMBER 13, 2020 4 Conclusions This model shows the importance of our simple action in an epidemic setting. Indeed the virus behavior is just an emergent behavior of our habits [10]. However this contribution open a new way in social network analysis, where the graph theoretical approach is substituted by agents. It also open the way to more realistic epidemic models, where the hypothetical scenarios can be tested directly on the agents, without any ODE mediation. 5 Acknowledgement The author received no financial support for the research of this article and all the simulations have been performed on the personal computer of the author. For these reasons doesn’t exist any conflict of interest in publishing the source code of the project in [13]. However author thanks the department of Mathematics and Informatics of the University of Palermo and the Institute of Biophysics of National Research Council (CNR) for funding author’s PhD and for being a constant guide for his growth.

[1] Chen Wang, Peter Horby, Frederick Hayden, and George Gao. A novel coronavirus outbreak of global health concern. The Lancet, 395, 01 2020. [2] Coronavirus, primi due casi in italia «sono due cinesi in vacanza a roma» sono ar- rivati a milano il 23 gennaio. https://www.corriere.it/cronache/20_gennaio_30/ coronavirus-italia-corona-9d6dc436-4343-11ea-bdc8-faf1f56f19b7.shtml . [3] Giulia Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, A. di Matteo, and M. Colaneri. Modelling the covid-19 epidemic and implementation of population-wide interventions in italy. Nature Medicine, pages 1 – 6, 2020. [4] Navid Mahdizadeh Gharakhanlou and Navid Hooshangi. Spatio-temporal simulation of the novel coronavirus (covid-19) outbreak using the agent-based modeling approach (case study: Urmia, iran). Informatics in Medicine Unlocked, 20:100403, 2020. [5] Son Woo-Sik et al. Individual-based simulation model for covid-19 transmission in daegu, korea. Epidemiol Health, 42(0):e2020042–0, 2020. [6] E. Cuevas. An agent-based model to evaluate the covid-19 transmission risks in facilities. Computers in Biology and Medicine, 121:103827 – 103827, 2020. [7] Covistat19. https://sites.google.com/community.unipa.it/covid-19/covistat19. [8] Vito M. R. Muggeo. Estimating regression models with unknown break-points.Statistics in Medicine, 22(19):3055– 3071, 2003. [9] G. Giacopelli, M. Migliore, and D. Tegolo. Graph-theoretical derivation of brain structural connectivity. Applied Mathematics and Computation, 377:125150, 2020. [10] Albert-László Barabási and Réka Albert. Emergence of scaling in random networks. Science, 286(5439):509–512, 1999. [11] K. Pearson. The problem of the random walk. Nature, 72:294–294, 1905. [12] Yilun Shang. Mixed si (r) epidemic dynamics in random graphs with general degree distributions. Applied Mathematics and Computation, 219(10):5042 – 5048, 2013. [13] Giuseppe Giacopelli. Cts-ext github. https://github.com/mrjacob241/CTS-Ext. [14] Günter Radons. The thermodynamic formalism of random walks: Relevance for chaotic diffusion and multifractal measures. Physics Reports, 290(1):67 – 79, 1997. [15] A Papoulis and Unnikrishna Pillai. Probability, random variables and stochastic processes. McGraw-Hill, 4th edition, November 2001. [16] Osservatorio unipolsai 2018. http://www.unipolsai.com/sites/corporate/files/pages_related_ documents/cs_osservatorio-unipolsai-2018.pdf. [17] Come siamo distribuiti nel territorio? https://www.tes.com/lessons/RP1R6l5ZB05CMA/ lombardia-settore-secondario-terziario-vie-di-comunicazione-e-popolazione . 8 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 14, 2020. ; https://doi.org/10.1101/2020.09.13.20193599doi: medRxiv preprint A PREPRINT - SEPTEMBER 13, 2020 [18] Covid-19 situazione italia. http://opendatadpc.maps.arcgis.com/apps/opsdashboard/index.html# /b0c68bce2cce478eaac82fe38d4138b1. [19] Dpcm 9 marzo 2020. https://www.gazzettaufficiale.it/eli/gu/2020/03/09/62/sg/pdf. [20] Aleksa Zlojutro, David Rey, and Lauren Gardner. A decision-support framework to optimize border control for global outbreak mitigation. 09 2018. [21] Rajiv Bhatia and Jeffrey Klausner. Estimation of individual probabilities of covid-19 infection, hospitalization, and death from a county-level contact of unknown infection status. medRxiv, 2020. [22] Worldometer covid-19 italy. https://www.worldometers.info/coronavirus/country/italy/. [23] Preeti Dhillon, Sampurna Kundu, Chander Shekhar, Usha Ram, Laxmi Dwivedi, Suryakant Dwivedi, Sayeed Yadav, and Sayeed Unisa. Case-fatality ratio and recovery rate of covid-19: Scenario of most affected countries and indian states, 04 2020. [24] Who advice for public. https://www.who.int/emergencies/diseases/novel-coronavirus-2019/ advice-for-public. [25] Duration of isolation and precautions for adults with covid-19. https://www.cdc.gov/coronavirus/ 2019-ncov/hcp/duration-isolation.html. [26] Derek Chu, Elie Akl, Stephanie Duda, Karla Solo, Sally Yaacoub, Holger Schünemann, Amena El-Harakeh, Antonio Bognanni, Tamara Lotfi, Mark Loeb, Anisa Hajizadeh, Anna Bak, Ariel Izcovich, Carlos Cuello-Garcia, Chen Chen, David Harris, Ewa Borowiack, Fatimah Chamseddine, Finn Schünemann, and Marge Reinap. Physical distancing, face masks, and eye protection to prevent person-to-person transmission of sars-cov-2 and covid-19: a systematic review and meta-analysis. The Lancet, 395, 06 2020. [27] Timo Mitze, Reinhold Kosfeld, Johannes Rode, and Klaus Wälde. Face masks considerably reduce covid-19 cases in germany. medRxiv, 2020. [28] Gianluca Manzo. Complex social networks are missing in the dominant covid-19 epidemic models. Sociologia, 14:31–49, 2020. [29] Hojeong Park and Songhee Kim. A study on herd immunity of covid-19 in south korea: Using a stochastic economic-epidemiological model. Environmental and Resource Economics, 07 2020. 9 . CC-BY 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.(which was not certified by peer review)preprint The copyright holder for thisthis version posted September 14, 2020. ; https://doi.org/10.1101/2020.09.13.20193599doi: medRxiv preprint