The key objective was to create a 3-layer model (Figure 1). The first layer was an agent-based particle model for Lombardy. Every agent is an inhabitant of the region, making this model a full-scale model of Lombardy. Therefore, we have 10.06 million agents who move according to the random walk theory [18]. The random walk behavior must be intended as an approximation of the actual motion of people during the day; this approximation was introduced to reduce the amount of information required to run the model and is widely used in many fields of science (eg, ideal gas theory) [18]. The large number of agents simulated is part of the novelty of this study because (to the best of the author’s knowledge) it is the first to attempt to simulate such a large population individual by individual for this purpose.

A collision detection algorithm was built within the agent-based model to detect whether 2 particles have a distance less than a fixed value. However, the large scale reached by the model required an ad hoc algorithm for this purpose; this challenging problem was solved via a square cells algorithm that permits the code to run in parallel (thereby decreasing the computational complexity of the task).

An epidemic model was built within the collision detection model, that is, a susceptible-infected-recovered-deceased (SIRD) model [19]. The model was filled with many experimental findings on COVID-19, and some fitting parameters were tuned on experimental findings. The whole analysis used as much open-access research as possible; moreover, the entire project was fully open source and is available on GitHub [20].

The model comprised three different layers:

The agents model simulates the behavior of each inhabitant of Lombardy using the approximation of random walks [18]. The displacement of the particles follows the density of inhabitants in Lombardy (ie, publicly available data). Even if more accurate data on people displacement and movement could be used, privacy concerns may not permit the open-source and open-access distribution of this data. Per the random walk approach, every particle moves with a random vector at every step. The model runs at 6 frames per day, which is a good frame rate considering that the scale time for epidemic phenomena is usually months; however, this can be improved as discussed in the Conclusions section. The random walk approach can appear unrealistic, but this hypothesis has been shown to be appropriate to model very large-scale systems (eg, gas thermodynamics [21]). In addition to random walks, a weak velocity field with a dependence of 1/r² was added, where r is the distance between 2 particles, as in a gravitational field, in order to aggregate the particles. The drift speed of the particles is constant and selected with a Weibull distribution [22] with a scale parameter of 6 and a shape parameter of 1.5. The particles’ speeds were 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 a report by UnipolSai Assicurazioni [23].

The setting of the simulation is Lombardy, making the environment model a closed 2D box with a boundary shape following Lombardy’s borders. In order to keep the particles inside the region, a bouncing condition was introduced at the border, so that a particle that tries to cross the border bounces backward. This condition is very popular in gas thermodynamics [18]. The initial conditions of the particles in terms of position are distributed following the actual density of the population of Lombardy, extracted from UnipolSai Assicurazioni [23] via image analysis [24]; this is then intended as an approximation of the real data.

Starting with the assumption that the algorithm has been designed to run on a commercial computer in parallel (the one used in the study has AMD 3900X 12-cores and 64 GB of RAM) and within reasonable time (about 20 minutes of calculation for 14 days of simulation), the collision detection algorithm played a central role in the implementation of the algorithm. In order to find all points with a distance less than a constant in a set of N points, a complexity order of N² is generally needed. In our model N≈10⁷, the complexity order was 10¹⁴, which is a large number.

Next, Lombardy was subdivided into a grid, 20 km in dimension. Collision detection was applied to every cell of the grid, and every cell was assigned to a separate parallel job to run the computation in parallel through the cells. This multiscale processing allowed for the speeding up of the code, reducing the RAM used simultaneously in computation, which made possible a simulation with 10 million particles at the same time. This approach neglects all the connections across the borders of the cells, but this is beyond the aims of this study.

The creation of this algorithm was a challenging aspect of this study. The idea was to use matrix optimization in order to speed up the computation. The territory was subdivided into 20-km–long cells, and the cells in every frame were completely independent, with the supposition that, on average, every cell contains m people. In order to compute the distance between all nodes in the network, we had to compute the order of N² pairwise distances.

With this scheme, we had to only compute the order of m² distances for each block multiplied by the number of blocks (which is about N/m) that is an order of Nm. Considering m small in comparison with N, it can be said that the scheme has a complexity near the order of N (for large N and small m). However, determining in which cell a person is located was also challenging because of the large size of the population. For these reasons, a simple grid scheme was used to locate nodes inside the cells. We used the following idea—supposing a segment from 0 km to L_C=2 km with N_c=4 cells:

If, for example, the point p=0.6 km needed to be located, the formula used to calculate this would be id_p=ceil(N_cp/L_c). The result is 2, indicating the second cell. Applying this formula for the x-axis and y-axis allows the algorithm to locate people in the cells. Although this algorithm may appear to be simple, it requires few calculations to be computed, which can make a substantial difference when a large number of agents is concerned.

The epidemic model is an SIRD model [3]. Most of the models available up to now are called population models [25]. A population model is a model where every node is modeled by a set of differential equations; it models a subpopulation 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. The model is not an ODE model, but a stochastic agent-based model. Every node has four states:

The proposed model was compared with a classical SIRD model [26] fitted with a parameter exploration scheme on outbreak data (Figure 2). As seen in Figure 2, the results are comparable in terms of the rooted mean square error of the data: the SIRD model had an error of 150 for the infected, 71 for the recovered, and 18 for the deceased; and the proposed model exhibited an error of 535 for the infected, 58 for the recovered, and 34 for the deceased. This indicates that our model has comparable performance with the SIRD model (outperforming in the recovered), but it is not ODE mediated and is thus suitable for testing alternative scenarios. Moreover, in this paper, since most of the parameters are realistic, the model can be run for a general epidemic upon collecting the few parameters required (which in this case were all open access) and fitting the two parameters left. However, the model can be made more precise by adding additional realistic data, which most of the time are not fully open access; this, however, is out of the scope of this study.

The first scenario was the Lombardy outbreak of March 2020 [2]. Our simulation began on February 29, 2020, and terminated on March 14, 2020. The main realistic parameters were p_I=1/40,500 (extracted from Bhatia and Klausner [27]) and p_D=0.3 (estimated from Worldometer [28], which has also been cited by Dhillon et al [29]).

The fitted parameters have a collision radius of 1 km. This can appear very large compared to the 1-m distance suggested by the World Health Organization [30]; however, when taken into account that there are 6 frames per day, then 1 km is the radius of the interaction of a person who remains in the same place for 4 hours and the duration of the disease (in days) E=7 (the Centers for Disease Control and Prevention suggests E<10 [31]).

The results of the simulations can be seen in Figure 3. The model was able to explain the experimental data until approximately March 9, 2020. On this day, the Decree of the President of the Council of Ministers (DPCM) implemented measures to contain the COVID-19 outbreak [32], which included the beginning of the lockdown in Italy. This discrepancy between the data and the model was the result of the collective effort of the Italian populace to protect itself against SARS-CoV-2. Therefore, the simulation serves to provide a warning about what could have happened.

The second scenario was inspired by Chu et al [34], who showed that maintaining a 2-m distance between people halved the risk of contracting COVID-19. Thus, we aimed to simulate this kind of social distancing by halving p_I in the model. The results (Figure 4) showed that COVID-19 (in this scenario) was not contagious enough to spread as in the experimental data. This simulation demonstrated the striking role of a simple action like social distancing in fighting COVID-19 and highlighted the difference between a virus under control and a disease of epidemic proportions. This simple fact has already been observed in experimental findings in Germany [35], where a synthetic method was used to estimate the spread of the contagion without the use of masks.

We also performed a lockdown simulation, reducing the daily average kilometers traveled by a node from 43 km to 5 km and reducing the interaction distance from 1 km to 100 m. The results of this simulation can be seen in Figure 5. According to the model, these simple actions were enough to control the virus.

The impact of topology in an epidemic model is a popular topic [30,36] in the debate on social networks. Thus, we performed a test: 1000 particles were chosen and then tracked across all simulations to find the total number of connections (ie, contact between particles) made within the whole population. In graph theory, the number of connections of a node is called a degree [16]. This test allowed us to determine the degree distribution and the daily degree distribution (average degree per day) of this small group of people across time. However, only the final result is presented (the full simulation is available on GitHub [20]). The first scenario was the COVID-19 outbreak scenario (Figure 6).

It can be seen that the distribution has an evident left tail (in contrast with the right tail of the Barabási-Albert models [17]). This was probably due to the simulation time of 14 days (in contrast with human social networks, which usually take years to be built). The lockdown scenario was also interesting. In this scenario, we observed a decline in connectivity from thousands of average connections per day to hundreds (Figure 7). This shows the importance of lockdowns in COVID-19 containment.

This scenario took into account the period between May 31, 2020, and June 14, 2020. During this period, Italy concluded its lockdown, and the number of active cases was decreasing. For this simulation, the kilometers per day was set arbitrarily to 15 km because of the lack of additional information on mobility during this period. The probability of contracting the contagion was halved to account for social distancing. The radius of interaction and the duration of the disease were tuned to reproduce the experimental data. The value for the radius of interaction was 300 m and disease duration was 5 weeks (E=35). This value (which is higher in comparison to that of the outbreak) could be influenced by a clinical protocol more accurately and by the queue created by the large number of infected people, which could slow down the tests required to declare recovery. The qualitative fitting can be seen in Figure 8.

Using the previous scenario of a decline in cases, we tested the impact of vaccinating 70% of the population, similar to the 62% suggested by Park and Kim [37]. The agent-based models are suitable for testing strategies like vaccination at the individual level. The result of the simulation was a strong decrease in infections, which was unexpected in a simulation of 14 days. The results are shown in Figure 9.

This model demonstrated the importance of people’s actions in an epidemic setting. Indeed, the behavior of the virus was indicative of our own habits [17]. The agent-based model proposed here has shown great flexibility in simulating alternative scenarios; in contrast, although ODE models [6] are faster than the proposed model, they are not suitable for this task.

The model proposed is more computationally expensive than ODE models, which require the calculation of few differential equations to simulate large populations. In general, such algorithms are also faster than agent-based models. The proposed model, however, allows for the interpretation of complex parameters.

This study has explained the behavior of the COVID-19 outbreak in Lombardy and has validated the herd immunity threshold obtained with different techniques [37], even if the 62% proposed by Park and Kim [37] is less than the 70% proposed in this study. This contribution also provides a new methodology in social network analysis, where the graph theoretical approach is substituted by agents. It also paves the way to more realistic epidemic models, where hypothetical scenarios can be tested directly on the agents, without any ODE mediation.

This work provides a novel, efficient, and low-demanding (in terms of computational resources) population model. Many features remain to be introduced in the model, like an age-dependent virus model, the ability to introduce an age parameter in the model or a more precise spatial simulation based on big data, and the ability to simulate the habits of the population. In conclusion, future work could be done to increase the number of frames per day, thereby improving the performance of the agents.