key: cord-0526158-o9qjgcz8 authors: Martucci, Roberta; Mascia, Corrado; Simeoni, Chiara; Tassi, Filippo title: Hospital management in the COVID-19 emergency: Abelian Sandpile paradigm and beyond date: 2021-02-23 journal: nan DOI: nan sha: 76433ac5207c99cf03a549973560a003ad0061a1 doc_id: 526158 cord_uid: o9qjgcz8 In this article, we propose a mathematical model -- based on a cellular automaton -- for the redistribution of patients within a network of hospitals with limited available resources, in order to reduce the risks of a local/global collapse of the healthcare system. We attempt at developing a conceptual tool to support making rational decisions relevant to the optimisation of the allocation of patients into accessible medical facilities. The strategy is based on a version of the Abelian Sandpile model for the Self-Organised Criticality, with the idea of testing the paradigm for the management of patients among the COVID-19 hospitals in Italian regions. In particular, we compare the novel proposal to the standard management of connections between hospitals, showing a number of advantages at a local and global level, by means of a reliable indicator function introduced for measuring the effectiveness of the allocation strategies. Since the beginning of 2020, people from all parts of the world are struggling with the COVID-19 coronavirus pandemic, which has led to a hundred million of infections and over 2 million deaths worldwide [20, 49] . The health emergency is not over yet and it will take months before being under control, possibly becoming endemic among most of the world population [32] . The current data suggest the paramount importance of acting promptly, and especially for trying to reorganize the medical facilities in order to avoid saturation of the critical care capacity in the various territories (refer to Figure 1 for the evolution of employment of the intensive care units in Italy). Indeed, several countries have been experiencing recurring epidemic waves, with a consequent worsening of the social and economic conditions, thus increasing dramatically the pressure on the healthcare systems. The emergence of infectious diseases has to be considered as an inherent property of human and animal populations, which cannot be generally avoided and/or foreseen: before the present COVID-19 emergency, SARS Most epidemics are characterised by outbreaks in localized geographical areas, and there may be no pharmaceutical solutions already available to safeguard the people's health when the disease begins to spread largely around. Moreover, the development of effective vaccines typically takes a rather long time, and therefore temporary alternative strategies must be implemented (quarantine, travel restrictions, activities closing, social distancing, face masks obligation, ...) However, the demand for specific healthcare facilities such as hospitals consistently increases, reaching and sometimes exceeding critical levels. For example, in the case of COVID-19, the saturation of intensive care units depends on the amount of patients requiring mechanical support for ventilation and proper equipment, so that any healthcare system can be overwhelmed if the number of severe cases becomes very high (refer to Figure 2 for the present-day situation of intensive care units in Italian regions). In that context, the rapid filling of the so-called COVID-19 hospitals in the areas most affected by the disease represents a serious danger while dealing with sudden emergencies at different scales (villages, cities, regions, countries). Hence, it is compelling to develop new arrangements of the healthcare system, for optimizing the overall structure especially in terms of an efficient distribution of patients among the most suitable facilities. In this article, we present innovative intervention options, by comparing the standard organisation of the Italian healthcare system with a novel proposal based on the Abelian Sandpile paradigm. We provide a systematic approach for an improved planning of the healthcare system by means of a mathematical model for the self-organisation of critical issues, which is capable to achieve some specific optimization inside the hospital management. In particular, we aim at reducing the risks of a local/global collapse of hospitals in times of crisis, while improving their functioning in normal times. More specifically, the basic model is grounded in the framework of cellular automata and postulates a large network of links between hospitals in a cooperative style of communication. Unfortunately, such connections are currently very limited, which drastically restricts the possibilities of reallocation for the supernumerary patients, and the strong urgency to enhance the hospital network is an important conclusion of the present analysis. The content of this article is organized as follows. In Section 2, we introduce the Abelian Sandpile paradigm, with a short description of its basic concepts and properties. Section 3 focuses on the adaptation of such paradigm to build up a novel proposal for the organisation of the healthcare system. In Section 4, we provide a gallery of examples aiming at illustrating how the proposed model works in the present context of COVID-19 epidemic. More specifically, we analyze two cases relevant to realistic applications: the central and peripheral outbreaks. The article is concluded by Section 5 which includes some generalisations of the model and (provisional) conclusions. ASM -Abelian Sandpile Model SOC -Self-Organized Criticality CAM -Cellular Automaton Model SRH -Sandpile with Redistribution to the Hub SID -Sandpiles with Internal Dissipation The notion of Self-Organized Criticality (SOC) was originally introduced by Bak, Tang and Wiesenfeld [4, 5] starting from a basic example proposed as a model for sandpiles (we refer to [3, 22, 37] for a general introduction to this broad subject). Since then, the concept has been expanded in many different directions, spanning from classical topics of physics (sandpile avalanches [27] , distribution of earthquakes [35, 47] , amplitude of solar flares [39] ) to less standard economic and socio-political contexts [2, 7, 9, 13, 24, 38, 42, 44, 52] , passing through computer networks and biological applications [1, 40, 53, 15] . At the same time, a huge effort to extend the mathematical tools to deal with theoretical questions has been made, thus contributing to drive SOC into an extraordinary crossroads of probabilistic approaches, graph theory, algebraic geometry, mathematical analysis and optimisation [6, 8, 10, 11, 12, 23, 28, 41] . 2.1. SOC and Sandpiles. The original application concerns with modeling of sandpiles, which are regarded as a manageable prototype of SOC, and fundamental contributions have been made by Dhar [17, 16, 18] notably in showing the crucial property of commutativity. On this account, the adjective Abelian has since been added to the technical terminology, leading to the actual meaning of Abelian Sandpile Model (ASM). Nevertheless, as usual in the most active fields of science, the terminology is not unique, and indeed similar topics are explored under different names (avalanches, chip-firing games, forest-fire models, parking functions, probabilistic abacus, Rotor-Router or Eulerian-walker models, ...) with more or less the same meaning. We refer to [14, 19, 29, 30, 36] for introductory presentations of the ASM and its various applications. We are particularly interested to the load balancing property, which denotes the method of distributing a set of tasks across a group of resources with the purpose of making their overall processing more efficient by balancing the workload of each operating unit. In the abstract formulation, units are represented by vertexes/nodes of a graph/network, with the corresponding connections represented by edges/links. The objective is to balance the loading process by allowing nodes to exchange particles with their neighbors through the incident edges. For the problem covered in this article, we build a specific ASM with critical height provided by the capacity of the medical facilities, and we propose to apply the Abelian Sandpile paradigm to achieve a methodical and efficient distribution of patients across the hospitals, in order to maintain the occupancy below the saturation level for dealing with eventual sudden and unexpected emergencies. A sandpile is also a type of Cellular Automaton Model (CAM), namely a discrete mathematical system fulfilling the following essential conditions: 1. the evolution takes place without external interventions; 2. the overall structural development depends only on local rules. To a large extent, the CAM is capable of simulating the dynamics of complex systems which are disposed to organize themselves through unstable critical states until reaching a stable configuration. In addition, the CAM implementation makes it possible to generate global coherent patterns starting from suitable local instructions, without any external supervisor in charge for understanding the process in its entirety [31, 43, 48] . In practical applications, each restricted portion of space contains a finite number of cells, which assume a finite set of (time-dependent) states. The initial configurationz 0 is the combination of a ground state configuration z 0 with some perturbation w, and for instancē After a given time interval ∆t > 0, the system typically comes to a new state z 1 , designated final configuration, which is determined by the changes of state of the single cell together with all the others, according to the preliminarily established series of (local) rules. Therefore, the evolution is described by means of a function Φ mapping some inputz 0 defined (from z 0 and w) over the graph/network Γ to an output z 1 , as specified by the simple formula z 1 = Φ(z 0 ). It is worth stressing that, in order to properly manage the patients between hospitals of the healthcare system, an optimal underlying ethicalcooperative paradigm has to be stipulated [34] . Sandpiles on a general network. We start by recollecting some basic definitions in graph theory [45, 25] . A graph (or network) is a pair Γ = (X, E) where X is a set whose elements are called vertexes (or nodes) and E is a set of paired vertexes called edges (or links). A rooted graph (or pointed graph) is a graph in which one vertex, denoted by x * , is distinguished as the root of the graph. In what follows, the root x * of the graph is also referred to as the hub, coherently with the subsequent meaning of the specific application to the healthcare system. Definition 2. An adjacent vertex y ∈ X of x ∈ X is a vertex which is connected to x by an edge, so that (x, y) ∈ E. The neighbourhood I(x) of a vertex x is the subgraph composed of the vertexes adjacent to x and all the edges connecting vertexes adjacent to x. The degree (or valency) of a vertex x, denoted by deg(x), is the number of edges which are incident to the vertex x. Let the graph Γ be connected, finite (with a finite number of vertexes denoted by x 1 , x 2 , . . . , x p ), simple (i.e. there are no loops connecting any vertex with itself) and undirected (i.e. the edges are bidirectional). We assume that a stock of identical particles is initially allocated at any vertex x i of the graph Γ. The configurations z 0 ,z 0 and z 1 , collecting the number of elements located at x i for any i, are natural-valued functions defined on the graph. The rules of the evolution are established to guarantee stability of the sandpile dynamics, starting from the definition of the height function given by z : X → N p (refer to Figure 3 for the graphical representation of a sandpile on a two-dimensional Cartesian grid). A stable configuration is a configuration of vertexes which are stable for any index i. A maximum stable configuration (or minimally stable state) is a configuration in which all the vertexes have a height z i of one unit lower than the corresponding threshold value deg(x i ). An almost stable configuration is a configuration of vertexes which are stable for any index i different from the root, and a maximum almost stable configuration is defined accordingly. In that framework, an unstable vertex is forced to topple over its neighbouring vertexes, thus causing a change of state in the entire configuration. A toppling (or firing) from the vertex x i is the mapping Ξ from the graph Γ to itself determined by the following (local) rules: and ∆ i→j is called the toppling matrix [17] . In principle, different stability criteria could be proposed involving values other than the degree deg(x i ) of the vertex x i , and modified redistribution criteria could also be considered, but we focus on the above definitions for the sake of simplicity in presentation. The addition of a particle to the sandpile is schematized by summing to the ground state configuration z 0 = (z 0 1 , z 0 2 , . . . , z 0 p ) the vector δ i with 1 at some fixed index i and zeros elsewhere. By iterating this procedure m times for possibly different choices of the index, the initial configurationz 0 is finally given by z 0 + w with w = δ i 1 + δ i 2 + · · · + δ im . It is a straightforward consequence that the operation of adding particles to a sandpile is associative. Moreover, the final configuration z 1 resulting from the evolution is independent of the order in which the topplings are performed, and therefore the operation of toppling is commutative [17] . As an illustrative case, we assume that deg(x i ) = deg(x j ) for two adjacent vertexes x i and x j . Then, we consider a sandpile in which both x i and x j are at their critical height, that is z i ≥ deg(x i ) and z j ≥ deg(x j ). According to the dynamics described above, a toppling from the vertex x i causes the vertex x j to become unstable, and viceversa, and the sandpile comes to a configuration in which the height z k for some k ∈ {1, 2, . . . , p} increases by a value i→k + j→k . Such procedure being symmetrical, by applying this argument repeatedly, we conclude that the same stable configuration is reached irrespective of whether x i or x j is toppled first, and in general regardless of the sequence in which unstable vertexes are toppled. 2.3. Two-dimensional Cartesian grids. We focus on the simple case of a two-dimensional Cartesian grid, with two standard schemes for adjacent vertexes given by the von Neumann neighbourhood, consisting of the four cells obtained moving one step towards North/East/South/West (refer to Figure 4 (left)), and the Moore neighbourhood, which includes also the four diagonal cells (refer to Figure 4 (right)). Of course, selecting a more general network provides a higher degree of flexibility and allows to closely represent the connections between various medical facilities of the healthcare system (refer to Figure 5 ). We consider a special grid of p = n 2 cells organized over a squared matrix of size n × n, where the particles are dropped randomly being allowed to stack on top of each other, so that a configuration of the sandpile is described by an element of the space of natural-valued matrices M n (N) := A ∈ R n×n : a ij ∈ N for any i, j = 1, 2, . . . , n . The position of each vertex x i (with its corresponding height z i ∈ N) is determined by the index i ∈ {1, 2, . . . , p = n 2 } according to the following For the von Neumann neighborhood, we have deg(x i ) = 4 for any index i associated with an element of the bulk of the matrix, and for instance Similarly, for the Moore neighborhood, we have deg(x i ) = 8 for any index i associated with an element of the bulk of the matrix, and for instance Example 1. We analyze the simple case n = 3 (p = n 2 = 9) with the von Neumann neighbourhood. By combining the perturbation w = 4 δ 5 with the ground state configuration z 0 = 0, we obtain the following initial and final configurations As expected, the threshold value deg(x 5 ) = 4 is reached at the central cell, thus inducing a destabilisation of the corresponding vertex, with the consequence that four particles are poured into the adjacent vertexes -the elements of the von Neumann neighbourhood-to compose a stable final configuration. As a matter of fact, the toppling from an unstable vertex changes the state of the adjacent vertexes, sometimes causing the appearance of instabilities at some other vertex. Then, subsequent topplings over the adjacent vertexes are generated, triggering a sequence of events which are evocative of sandpile avalanches and stopping only when all the cells return strictly below their threshold capacity. If the graph is infinite and connected, with a finite number of particles, all vertexes become stable after a finite number of topplings; moreover, it can be proven that the final (stable) configuration z 1 depends solely on the initial configurationz 0 , independently of the order in which the topplings are performed [17] . On the other hand, if the graph is finite, appropriate boundary conditions must be implemented: assuming that the particles are evacuated from the boundaries, an analogous result of stability holds [14] , but other types of boundary conditions may not guarantee the same property. During the operation of toppling, no particles are created as they are redistributed to neighbouring cells. For dissipative boundary conditions, the toppling can actually lead to the loss of particles if it occurs on a boundary cell, and we refer to this case as open boundary conditions. In this article, we pursue the idea of applying the Abelian Sandpile paradigm to the management of patients among the COVID-19 hospitals in Italy. The current organisation assigns responsibility for the healthcare system to the regions, which are restricted territorial bodies with their own statutes, powers and functions established by the Italian Constitution. We refer to the Lazio region when selecting (the order of magnitude of) the number of medical facilities relevant to the numerical simulations. More specifically, the healthcare system of the Lazio region is composed of about 100 hospitals, mostly located in the metropolitan area of Rome, which are structured within a network of reciprocal connections [33] . In that framework, the height function of the nodes (hospitals) indicates the number of hospitalised patients, and the links between different nodes represent the possibility of a direct exchange of patients (that is subject to constraints of geographical proximity, together with inherent organisational requirements of the connecting medical facilities). The node of the network exhibiting the largest number of links is the hub, which usually takes on a strategic function for the whole system, whilst the nodes with fewer links identify the hospitals with unfavorable geographical location and reduced capacity/functionality (each hospital having different resources, in the general case the threshold value changes from node to node). 3.1. Reassignment of outgoing particles. We have already discussed about the crucial role played by the boundary conditions associated with the ASM settled on finite networks, and indeed the existence of stable configurations is strongly influenced by the presence of a dissipation mechanism such as the open boundary conditions. Because the loss of particles through the edge of the network and the hypothesis of an infinite network are unrealistic assumptions for practical applications, we have to implement a suitable reassignment law for particles generating a critical height at some boundary cells. Hence, we are induced to redistributing the outgoing particles to the root of the network, namely the hub, which is typically a collecting node of the healthcare system. Therefore, we call the resulting process a Sandpile with Redistribution to the Hub (SRH), and we notice that actually the SRH is equivalent to the ASM when there are no topplings occurring at the edge of the network. In particular, for two-dimensional Cartesian grids with odd size, we incorporate the reassignment of outgoing particles from cells toppling at the border of the matrix by proposing that they are reallocated to the central cell. Furthermore, in case the threshold capacity is reached in several cells at the same time, we impose that the first cell to topple is always the hub, and then additional toppling is performed from the other nodes. In particular, the hub is allowed to topple only once, at the very beginning of the toppling procedure. We summarize the essential steps of the CAM workflow -which turns out to be also the implementation scheme of the numerical codes-as follows. 1. Initialisation. We choose a (stable) ground state z 0 = (z 0 1 , z 0 2 , . . . , z 0 p ) which satisfies the condition that z 0 i < deg(x i ) for any index i ; 2. Inflow. We combine the ground state z 0 with a perturbation w = (w 1 , w 2 , . . . , w p ) which indicates the number of new patients requiring hospitalization, and we obtain the initial configurationz 0 = z 0 + w ; 3. Hub toppling. If z 0 * ≥ deg(x * ), the hub performs a toppling towards the adjacent nodes of its neighbourhood I(x * ) ; 4. Additional toppling. If z 0 i ≥ deg(x i ) for some x i = x * , these nodes perform (a sequence of) topplings until reaching an almost stable configuration z 1 = Φ(z 0 ) ; 5. Iteration. We reinitialize the ground state with z 0 equal to z 1 , and we restart from 2. Comparison with the standard healthcare system management. In order to to compare the efficiency of the novel proposal with the standard organisation of the Italian healthcare system, we provide a reformulation of the current management of connections between hospitals in terms of the CAM paradigm. Perhaps surprisingly, the connection between medical facilities is presently mainly determined by their geographical proximity, corresponding to the basic adjacency of nodes. If a patient comes to a hospital where there are no places available, because the threshold capacity has been reached, then the single patient is reallocated to the nearest hospital with available places, and the reassignment is limited to the patients exceeding the threshold value. That being the case, the redistribution of patients is handled manually at each hospitalization, without foreseeing the possibility of repeated similar events which are instead highly probable during sudden and unexpected emergencies. We summarize the essential steps of the standard workflow as follows, by assuming the reassignement law from the edge of the network at the hub as in Section 3.1. 1-2. Initialisation/Inflow. We proceed as in Section 3.1 ; 3. Hub redistribution. If z 0 * ≥ deg(x * ), the hub reallocates only the exceeding patients to the adjacent medical facilities, starting from the less crowded ones (in case of equal crowding, a random choice is operated) ; 4. Additional redistribution. If z 0 i ≥ deg(x i ) for some x i = x * , the redistribution procedure is repeated for these nodes, until reaching an almost stable configuration z 1 = Ψ(z 0 ) ; 5. Iteration. We reinitialize the ground state with z 0 equal to z 1 , and we restart from 2. We illustrate the comparison by considering the following simple example. where Ψ denotes the evolution mapping of the grid according to the standard approach. Then, we reinitialize the ground state with Ψ(z 0 ) and the perturbation w is chosen as above, so that the configuration Ψ(z 0 ) + w is also unstable with respect to the central cell. By selecting randomly one of the neighbouring cells, another single particle is transferred, and after two additional iterations, the system reaches the final configuration On the other hand, the ASP paradigm suggests to transfer at the same time all the four patients initially allocated to the hub towards the von Neumann neighbourhood, to obtain Finally, under the assumption that new patients are always collected at the hub, the system can accept exactly 3 more particles from the additional iterations, before reaching a critical situation as before. Of course, it could be argued that, being the final configurations the same, there would be no obvious evidence to prefer the novel approach to the standard organization. However, the deciding factor is that we have unified within a single timestep the efforts for reallocating patients among the medical facilities, which is usually regarded as a source of stress for the overall healthcare system. Hence, we have left the structure time to reorganize without suffering a constant state of emergency, which is precisely the issue of the so-called predictive logistics. We collect numerical results from a selection of preliminary cases, which are nevertheless useful to understand the alternative methods of healthcare system management provided by the two approaches described in Section 3, and their inherent dynamics. We consider two-dimensional Cartesian grids, which are arranged into n × n (squared) matrices with the Moore neighbourhood, so that deg(x i ) = 8 for any index i ∈ {1, 2, . . . , p = n 2 }, and we assume the reassignment to the hub of outgoing particles from the edge of the network as in Section 3.1. We recall that each cell/node of the grid represents a hospital, and the height function reproduces the number of patients in the medical facilities, expressed as the percentage of capacity already achieved. We start by focusing on the possible outputs of a basic case consisting in a single hospital originally attaining its threshold capacity, which is located at the hub (central cell), with a (randomly chosen) stable ground state configuration. To highlight possibly critical situations, we boldmark all values above the threshold capacity, and we encircle the values 5, 6 and 7 dangerously close to the saturation level. Example 3. For n = 3, the hospitals are p = n 2 = 9, and we choose an initial configuration with all new patients coming to the hub x 5 , given by the perturbation w = 4 δ 5 , so that Because the hub has to manage a number of hospitalisations greater than its threshold capacity, some patients must be transferred to adjacent facilities. According to the standard strategy, only four patients need a different allocation to be found among the neighbourhing cells, ending at the (stable) configuration As it is clearly seen, an additional iteration with w = 4 δ 5 determines a new unstable/critical configuration at the central node (with again four patients in excess), demanding for a further reorganisation of the hub by means of another transfer operation. On the other hand, by employing the SRH approach, we obtain a different outcome given by such situation being more favorable in terms of load-balancing since an additional iteration with perturbation w = 4 δ * does not produce any unstable configuration. Moreover, there is still a certain amount of available places at the hub, and therefore the novel proposal solves the criticalities possibly generated by an emergency. Despite its simplicity, the above example suggests to attribute to each configuration the value of a given functional, aiming to provide an easily readable indicator of the effectiveness of the allocation strategy. This is necessarily a very delicate issue. As a first attempt, a possible choice is to count the total number of medical facilities attaining a given fraction of their threshold capacity, which are denoted by critical points. Such choice is actually quite questionable, since it does not take keep memory of the value of incoming patients at the beginning of the iteration. A different (rough) quantitative measure of the system management efficiency is which intends to evaluate the risk that a new patient comes to a given location, taking into account the previous history of the system (in this case, the effect of the perturbation w). For the case illustrated in the Example 3, the function F equals the state of the hub x 5 , so that F(w,z 0 ) =z 0 5 = 11, F w, Ψ(z 0 ) = Ψ(z 0 ) 5 = 7, F w, Φ(z 0 ) = Φ(z 0 ) * = 3 . The minimum value is achieved for the configuration Φ(z 0 ) corresponding to the SRH approach (or, equivalently, to ASM). Example 4. For n = 5, the hospitals are p = n 2 = 25, and we choose an initial configuration with all new patients coming to the hub x 13 , given by the perturbation w = δ 9 + δ 12 + 4δ 1 3 + 2δ 14 + δ 18 + δ 19 , so that , and the same observation previously done holds, that is such configuration likely determines an instability located at the hub after a subsequent process iteration, since the number of hospitalised patients is very close to the critical threshold. On the other hand, by employing the SRH approach, we obtain where the number of patients hospitalised at the hub is far away from the critical level. Indeed, by computing the indicator function F at the different outcomes, we deduce that since j w j = 10, and the smallest value is achieved by the configuration resulting from the SRH approach; thus, such strategy must be preferred to the standard one. Next, we explore a case where two subsequent topplings occur. Example 5. For n = 5, we choose the ground state and the initial configuration as follows, , in which the cell x 12 is unstable, and finally we arrive at the final configuration The matrices Ψ(z 0 ) and Φ(z 0 ) are represented in Figure 6 with different colors assigned to cells, corresponding to the degree of saturation achieved. Number Figure 6 . Model comparison from the Example 5: final configurations Ψ(z 0 ) (left) and Φ(z 0 ) (right). The colors attached to each cell/node represent its corresponding relative capacity (0=green, 1/2=yellow, 3/4=magenta, 4/5=red, 6/7=black) At the bottom of each subfigure, we report the counting of critical points, given by the number of cells with values 6 or 7. Although a quick glance can make us suspect that the left configuration is preferable, the situation is indeed different. In fact, first of all we distinguish the presence of two almost critical values 7 inside the configuration plotted on the left, against the presence of only one value 7 from the second approach. Furthermore, the hub -the facility most at risk taking into account recent history-has a value 7 in the first case, whilst a value 4 occurs in the second case, thus making the latter situation preferable. Finally, we compute the function F corresponding to the different configurations: since w/ j w j is everywhere zero except at the hub, where it is equal to 1, the values of F coincide with the values at the hub, which are given by suggesting again the advantages of configurations proposed by the SRH approach. 4.1. Simulation of multiple central outbreaks. We analyze the case with more than one hospital in the same restricted area (including the hub) which is concerned with new incoming patients. One might think of this event as the emergence of multiple outbreaks within some smaller district of a larger area. We consider a larger network, with n = 9, whose size is comparable to the capacity of the healthcare system in the Lazio region in Italy [33] . We consider the perturbation w = 2δ 32 +δ 33 +5δ 41 +2δ 42 and the ground state 1 2 1 4 6 2 3 6 2 3 2 5 2 1 3 2 4 3 3 3 1 5 6 2 3 1 3 6 1 5 3 2 5 2 3 4 1 3 2 1 6 3 4 5 5 1 4 1 2 2 5 1 2 3 4 1 3 4 5 6 2 5 3 4 2 3 5 2 2 6 3 1 1 6 5 2 4 4 2 so that j z 0 j = 250 and j w j = 10, corresponding to the initial configurationz 1 2 1 4 6 2 3 6 2 3 2 5 2 1 3 2 4 3 3 3 1 5 6 2 3 1 3 6 1 5 3 2 + 2 5+1 2 3 4 1 3 2 1 6+5 3+2 4 5 6 1 4 1 2 2 5 1 2 3 4 1 3 4 5 6 2 5 3 4 2 3 5 2 2 6 3 1 1 6 5 2 4 4 2 Then, according to the standard strategy, a final configuration is given by 1 2 1 4 6 2 3 6 2 3 2 5 2 1 3 2 4 3 3 3 1 5 6 2 3 1 Hence, the central subgraph around the hub is composed by the elements In the subsequent steps, any additional patients coming to the hub destabilises the configuration. On the other hand, by employing the SRH approach, we obtain The difference between Ψ(z 0 ) and Φ(z 0 ) is transparent regarding, specifically, the number of patients hospitalised in the hub. A further iteration with the same perturbation w would lead to a new critical situation for the hub in the first case, but it does not in the second. Figure 7 provides a representation of the two configurations Ψ(z 0 ) and Φ(z 0 . Counting the critical points in Figure 7 could give the (incorrect) impression that the standard approach has to be preferred with respect to the novel proposal. But considering the indicator function F applied toz 0 , Ψ(z 0 ) and as a consequence of the fact that memory of the dynamics -represented by wis now considered, and again the minimum value is achieved for the SRH strategy. Simulation of peripheral outbreaks. Next, we focus on the case of outbreaks occurring in the vicinity of a node located at the edge of the Cartesian grid. Specifically, we take n = 9 and the same ground state z 0 of the previous example, and an inflow matrix of patients w given by the submatrix 5 2 2 1 located at the top-right corner of the null matrix. Then, the initial configurationz 0 = z 0 + w is given bȳ 4 3 2 5 2 1 3 2 6 4 3 3 1 5 6 2 3 1 3 6 1 5 3 2 5 2 3 4 1 3 2 1 6 3 4 5 6 1 4 1 2 2 5 1 2 3 4 1 3 4 5 6 2 5 3 4 2 3 5 2 2 6 3 1 1 6 5 2 4 4 2 with the cell above the threshold capacity located at position (1, 8) . Therefore, four patients are in excess and require a different reallocation. Among others, a possible output of the standard strategy is the final configuration given by 1 2 1 4 6 2 5 7 4 3 2 5 2 1 3 4 6 4 3 3 1 5 6 2 3 1 3 6 1 5 3 2 5 2 3 4 1 3 2 1 6 3 4 5 On the other hand, by employing the SRH approach, we pass through an intermediate state  1 2 1 4 6 2 4 3 5 3 2 5 2 1 3 3 7 5 3 3 1 5 6 2 3 1 3 6 1 5 3 2 5 2 3 4 1 3 2 1 9 3 4 5 6 1 4 1 2 2 5 1 2 3 4 1 3 4 5 6 2 5 3 4 2 3 5 2 2 6 3 1 1 6 5 2 4 4 2 to reach the final configuration 1 2 1 4 6 2 4 3 5 3 2 5 2 1 3 3 7 5 3 3 1 5 6 2 3 1 3 6 1 5 4 3 6 2 3 4 1 3 2 2 1 4 4 5 6 1 4 1 3 3 6 1 2 3 4 1 3 4 5 6 2 5 3 4 2 3 5 2 2 6 3 1 1 6 5 2 4 4 2 . Figure 8 provides a representation of the two configurations Ψ(z 0 ) and Φ(z 0 , by indicating the level of saturation reached in each hospital. Finally, the indicator function F defined at the beginning of Section 4 takes the following values for the different configurations, F w,z 0 = 7.9, F w, Ψ(z 0 ) = 5.9, F w, Φ(z 0 ) = 4.4 , and this fact contributes to the advantages of the SHR approach. The mathematical model introduced in this article lays the foundation for an optimisation of the healthcare system management. A notable feature of the novel proposal is its scalability to various levels of description, and also the possibility of improving the experimental simulation by including In order to incorporate other relevant elements into the SHR model, such as recovery or (unfortunately) death of patients, we postulate an elimination mechanism inherent to the system, which is translated in mathematical formalism by considering the presence of some dissipation during the evolution. For instance, we describe the effect of removing particles from the network by a simple subtraction of a (randomly chosen, but possibly measurable) distribution ζ = (ζ 1 , ζ 2 , . . . , ζ p ) with the obvious constraint that 0 ≤ ζ i ≤ Φ(z 0 ) i for any index i = 1, 2, . . . , p . The modified model workflow is divided into the following steps. 1-2. Initialisation/Inflow. We proceed as in Section 3.1 ; 3-4. Hub/Additional toppling. We proceed as in Section 3.1 ; 5. Internal dissipation. We subtract the distribution ζ to the intermediate configuration z 1 , leading to the final configuration Φ(z 0 ) = z 1 − ζ ; 6. Iteration. We reinitialize the ground state with z 0 equal to Φ(z 0 ), and we restart from 2. We refer to this modified model as Sandpile with Internal Dissipation (SID). In the long run, after a large but finite number of iterations, a balance between the inflow and outflow steps has also to be incorporated in order to guarantee conservation of the total number of patients. In particular, it could be useful to add the hypothesis that where N is the total number of iterations, with w n and ζ n denoting the inflow and outflow contributions at the time-step n, respectively. Indeed, we notice that if the lefthand side is larger than the righthand side, the whole system risks to undergo a finite time collapse, by reaching its total capability -sum of the capacities of each single structure-in finite time. Other interesting features can be added to provide the model with a higher level of realism: for instance, the presence of some inertia to the transfer process [26] , giving preference to structures with a certain level of hospitalized patients [21] , or constrain additional bulk dissipation [46] . From the analysis developed in this article, we observe that the standard healthcare system management typically generates highly unstable and unbalanced configurations, where specific geographical areas with semi-empty hospitals alternate with others where all medical facilities are saturated, especially during sudden and unforeseen events like the spread of epidemics. Instead, following the SRH strategy for optimized management of connections between hospitals, seems to produce more sparse allocation of patients, which has to be considered as a preferable configuration in terms of load-balancing. On the other hand, it is crucial to improve the exchange of information and to provide decision-making tools to the local structures, in order to optimise the healthcare response in normal times and to avoid the collapse of individual hospitals in times of crisis. There are many relevant consequences in the socio-economic field: among others, we stress the riveting possibility of the automation of health protocols, meaning to build an application capable of learning something from the data autonomously, without receiving explicit instructions from the outside. Such conceptual experimentation could create a learning environment in which policy makers may gain a better understanding of how the system responds to their decisions, providing forecasts of potential different choices and strategies. We are aware that a paradigm shift is required and we hope to give a contribution to this respect. It is worth stressing that the agreement with realistic experimental data is presently very limited. However, the conceptual framework we have proposed in this article applies in principle to many different contexts, and these research directions are currently being explored. Self-organized criticality in living systems Self-organized criticality in a model of collective bank bankruptcies How nature works: The science of self-organized criticality Self-organized criticality: An explanation of the 1/f noise Self-organized criticality The steepest descent algorithm in Wasserstein metric for the sandpile model of self-organized criticality Self-organized criticality and stock market dynamics: An empirical study Self-organised criticality in stochastic sandpiles: Connection to directed percolation Modeling financial markets by self-organized criticality Controlling self-organized criticality in complex networks Self-organized criticality: Sandpiles, singularities, and scaling Metric features of self-organized criticality states in sandpile model Scaling laws and indications of self-organized criticality in urban systems Self-organized criticality model for brain plasticity Sandpiles and self-organized criticality Self-organized critical state of sandpile automaton models The Abelian sandpile and related models Theoretical studies of self-organized criticality On the effect of the drive on self-organized criticality Self-organised criticality -What it is and what it isn't Mean-field self-organized criticality Self-organized criticality of belief propagation in large heterogeneous teams Handbook of Graph Theory Crossover to self-organized criticality in an inertial sandpile model Optimization by Self-Organized Criticality Introduction to the sandpile model Theory of cellular automata: A survey Immunological characteristics govern the transition of COVID-19 to endemicity For whose benefit? The Biological and Cultural Evolution of Human Cooperation Self-organized criticality applied to natural hazards Sandpile models of self-organized criticality Power laws and self-organized criticality in theory and nature Workplace accidents and self-organized criticality Self-organized criticality in solar flares: A cellular automata approach Evidence for self-organized criticality in evolution Theoretical results for sandpile models of self-organized criticality with multiple topplings Self-organized criticality and the predictability of human behavior A brief history of cellular automata Mechanisms of self-organized criticality in social processes of knowledge creation Introduction to Graph Theory Proof of breaking of self-organized criticality in a nonconservative Abelian sandpile model Self-organized criticality A new kind of science Technological innovation, business cycles and selforganized criticality in market economies Self-organized criticality in a computer network model piazzale Aldo Moro 2 -00185 Roma (Italy) Email address: martucci.1701560@studenti.uniroma1.it (Corrado Mascia) Dipartimento di Matematica G. Castelnuovo, Sapienza Università di Roma, piazzale Aldo Moro 2 -00185 Roma (Italy) Email address: corrado.mascia@uniroma1.it (Chiara Simeoni) Laboratoire de Mathématiques J.A. Dieudonné CNRS UMR 7351, Université Côte D'Azur, Parc Valrose -06108 Nice Cedex 2 (France) Email address: chiara.simeoni@univ-cotedazur.fr (Filippo Tassi) Dipartimento di Matematica G. Castelnuovo, Sapienza Università di Roma The authors thank the Department of Mathematics G. Castelnuovo, Sapienza University of Rome, for hosting the electronic workshop COVID-19 calls for Mathematics -www1.mat.uniroma1.it/ricerca/convegni/2020/COVIDchiamaMATwhich has motivated this article. The authors thank prof. Ferdinando Romano (Department of Public Health and Infectious Diseases, Sapienza University of Rome) for the stimulating discussions about the logistic issues of the healthcare system organisation.The authors also thank prof. Stefano Finzi Vita (Department of Mathematics G. Castelnuovo, Sapienza University of Rome) for useful references on the mathematical theory of the Abelian Sandpile model. All authors read and approved the submitted manuscript. All authors have agreed to be personally accountable for the contributions, and to ensure that questions related to the accuracy or integrity of any part of this article are appropriately resolved and documented in the literature. All authors equally contributed to the conception and design of this article, the acquisition, analysis and interpretation of data, the drafting and revision of the manuscript.