A self-organizing criticality mathematical model for contamination and epidemic spreading

We introduce a new model to predict the spread of an epidemic, focusing on the contamination process and simulating the disease propagation by the means of a unique function viewed as a measure of the local infective energy. The model is intended to illustrate a map of the epidemic spread and not to compute the densities of various populations related to an epidemic, as in the classical models. First, the model is constructed as a cellular automaton exhibiting a self-organizing-type criticality process with two thresholds. This induces the consideration of an associate continuous model described by a nonlinear equation with two singularities, for whose solution we prove existence, uniqueness and certain properties. We provide numerical simulations to put into evidence the effect of some model parameters in various scenarios of the epidemic spread.


1.
Introduction. Classical mathematical models for epidemic are based on systems of differential equations describing the evolution of various populations (infective, susceptible, immune) which may possible occur during the epidemic development (see e.g., [8], [11], [6], [14]). The purpose of this paper is to introduce a new model to predict the spread of an epidemic, focusing on the contamination process and simulating the disease propagation by the means of a unique function viewed as an indicator of the local infective energy. This model allows to illustrate a map in real time of the epidemic spread and it is not though to compute the densities of various populations related to an epidemic, as done by the classical models. The model is set up as a cellular automaton (CA) whose rules describe the local character of the contamination, the passing from a noninfective state into an infective state, the spread process, and the possible epidemic natural or induced (by a treatment) extinction. This CA simulates the local transmission as a selforganizing criticality process, meaning that there exists a critical phase at which the epidemic propagation starts and at which it ends, after a possible growth and decrease. We refer the reader e.g., to [26], [1], [3], [19] for detailed explanations of a self-organizing criticality process. After writing the discrete system of equations which rule the CA we introduce a continuous description expressed by a nonlinear diffusion equation with two singularities in the diffusion term.
Previous approaches of epidemics by CA models have been done e.g., in [4], [5], [13] and more recently in [25] and [20]. In the last article, the state of each cell stands for the fraction of the susceptible and infected individuals of the cell at a particular time step and this allows that one cell can be occupied by many individuals. A local transition linear function is introduced to connect the states of the neighbor cells and other parameters such as the virulence of the epidemic and the rate of recovered infected individuals. In these models equations referring to the densities of the various populations (infected, susceptible, recovered) involved are considered, while in our model the process is described by the means of only one function.
The second section is devoted to the description of the rules of the CA selforganizing criticality type model. Suggested by the mathematical relation describing the rules of the CA, an appropriate continuous model represented by a nonlinear partial differential equation with two singularities is introduced. In Section 3 we study the existence and uniqueness of the solution to the continuous model, by relying on an appropriate approximating problem without singularities. After proving existence and uniqueness in the approximating problem, by a passing to the limit technique we deduce the existence, uniqueness and other properties of the solution to the limit model. In Section 4, on the basis of simulations using the CA, starting from possible scenarios of contamination in different episodes of epidemics, we analyze the probability distribution of the size, area and duration of the events. Also, we present some graphics showing the dynamics of the process and compute the space extension of the disease and its stopping time. The last section gives a sketch of the numerical method for the integration of the continuous model and presents numerical simulations focusing on the influence exerted by some parameters of the model. The essential particularities of this novel model are summarized in Conclusions.
2.1. The CA model. Certain observed time duration and space spreading of some epidemics outbreak indicate that there exists a power-law dependence of frequenciessize of events, [21]. Such kind of behavior was thought to be specific to the dynamics of the complex system and was associated with the existence of the critical state, [26]. To illustrate the dynamics of very simple cell automata models, the basic sandpile model, forest fire or percolation were created. In spite of their simplicity they proved to be very useful in many applications, as crack propagation and earthquakes [9], [22], evolutionary biology [23] and epidemics, [10], [15], [24]. The basic sand-pile model uses three concepts that can be beneficial when one analyzes the epidemic spreading, namely the local equilibrium, the threshold activation and the diffusive character of the rules update.
We consider that during an epidemic individuals are separated in three classes: infected (I), which can transmit the disease, susceptible (S), which cannot transmit disease can be infected and roughly speaking inactive (R) which cannot neither transmit nor get the virus. In this last class we include immune, dead of the lack of individuals. We call these classes phases and their dynamics is followed in a 2D space domain Ω, represented as a lattice with square cells denoted (i, j). We call neighbor cells two cells that have a common side. One knows from the classical theory of SIR models (see e.g., [11], [16]) that there exists an endemic equilibrium where both the susceptibles and infected individuals exists and this equilibrium can be locally asymptotically stable or unstable. On this base we can associate to each cell a measure, denoted by p ij (t) describing the sensibility of cell to be disturbed from its endemic equilibrium and the infective action it may have upon the neighbor cells. These cells particularities are viewed between two thresholds which are intrinsic properties of each cell. Thus, let p ij R and p ij I be two positive numbers, associated to each cell, p ij I > p ij R , used to classify the phases as More exactly, we assume that each cell is occupied by many of individuals and that p ij represents the local fraction between the number of infected individuals and the total number of individuals in a cell. The cell (i, j) exhibits infection if at least an individual is infected and then its p ij ≥ p ij I . There is not yet infection, but it possibly can occur, if p ij R < p ij < p ij I . The state p ij = p ij R stands for cells through which the disease cannot be spread by various reasons. This lower threshold is assigned to void cells, to those occupied only by immune who has got immunity previously or expresses the fact that the type of contact allowed by this cell is not favorable to disease inoculation. The variation of p ij gives an information about the intensity of the infection (e.g., if p ij >> p ij I ) or the degree of susceptibility (e.g., if p ij ∈ (p ij R , p ij I )), meaning how much is the energy of a cell to infect another one and how much is the probability that a S cell to become infected.
The two local critical thresholds, p ij I and p ij R showing the range of extension of the infective energy of a cell, can be set on the basis of a statistical analysis of historical outbreaks. The latter threshold acts as an obstacle, the process of contamination being permitted only for the states having the value above p R .
We write down the basic assumptions for the model.
(i) We do not consider an age-structured population and we assume that there is no flux across the boundary. (ii) Virus can be transmitted from a cell (i, j) with the state I (having p ij ≥ p ij I ) only to a neighbor cell (k, l), upon contact, indifferently of the value p kl (greater or lesser than its p kl I ). This transmission will be evidenced by the increase of the value p kl . (iii) The susceptibles are not equally infected and the infected are not equally infectious, the transmission and reception of the disease depending on the level of the phase function p ij. Thus, the cell with p ij having a value lesser but close to p ij I can jump into the phase I within a shorter time that a cell having p ij closer to p ij R . (iv) There exists the possibility of viral state maintenance, that is the value of the state p ij ≥ p ij I may increase for a while, and eventually decrease, due to the fact that the viral load increases with the infection time and then decreases during a natural (or controlled) healing process (see e.g., [14]). This is taken into account in the model as a source function b ij (p ij (t)) which may have the following form and which implicitly contains the dependence on the infection time since b ij depends on p ij ≥ p ij I . (v) The disease dynamics allows the possibility of reinfection, because contacts that can result in disease transmission may be repeated events; this means that there exist t 1 and t 2 positive, such that if p ij (t) ≥ p ij I and if p ij (t + t 1 ) < p ij I , then one allows that p ij (t + t 1 + t 2 ) ≥ p ij I . (vi) The level p ij I is critical for the process, that is a cell is activated only if p ij (t) ≥ p ij I , moment at which the spreading process is triggered. Otherwise, the cell does not change its phase level if p ij (t) < p ij I . Next, we present the rules of the CA running in the lattice (i, j), in the time interval (t, t + 1) and considering Γ ij the class neighborhoods Γ ij of the cell (i, j). We specify again that a neighbor cell is that which has a common side with the cell (i, j).
A cell with p kl ≥ p kl I is called activated. (r1) The contact of a cell (k, l) with p kl (t) ≥ p kl I with a neighbor cell (i, j) having p ij > p ij R determines the increase of the value p ij with a rate a kl (t). (r2) A cell (i, j) with p ij (t) = p ij R remains always at its value p ij R even if it has a contact with a cell with p kl ≥ p ij I . (r3) After becoming active a cell can exhibit an increase of its p ij , possibly followed by a decrease, due to the assumption (iv). (r4) The process stops when the cells have p ij ≤ p ij I , or when a cell with p ij ≥ p ij I is surrounded by 4 cells with p kl = p kl R . Thus, these rules of the CA specify that at time t + 1 the value of the phase function is determined by where H is the Heaviside function H(r) = 1 for r > 0 and H(r) = 0 for r < 0.
2.2. The continuous model. By expanding the sum with respect to the neighbor cells in formula (1) we get the 2D discretized form of the Laplace operator −∆(a(t, x)H(p(t) − p I (x)), and this suggests to consider the following continuous form to which we attach an initial condition In the study of PDEs with jumps in coefficients one usually proceed to the filling in jumps, so that we shall consider in our model the graph Heaviside (instead of the discontinuous function) For mathematical compatibility we add homogeneous Dirichlet boundary condition and remark that it is satisfied if p ≤ p I a.e. on Σ.
We mention that previous existence results and solution properties for certain self-organizing criticality models have been obtained in [3], for the sand-pile model described by the equation y t − ∆(y(t) − y c ) ∋ 0, and in [19] for the nonautonomous sand-pile model represented by both with homogeneous Dirichlet boundary conditions.
In the next section we give a meaning to the continuous autonomous model (2) and study the existence and properties of its solution for time independent parameters.
3. Existence results for the continuous autonomous model. We consider Ω an open bounded subset of R 2 with a C 1 class boundary Γ := ∂Ω, T > 0, finite and let Q := (0, T ) × Ω, and Σ : Let a(t, x) = 1, p R , p I be constant and b not depending explicitly on x. Then, model (2)-(4) reads Let σ be positive, arbitrary small, σ << 1, and introduce the approximating problem where Denoting and taking into account that, as explained before, b(p) = 0 for p < p I and b p ≥ 0 for p ≥ p I , we can write (8) as We introduce j σ , the potential of j σ , by and so we have j σ (r) = ∂ j σ (r) for all r ∈ R, where ∂ j σ denotes the subdifferential of j σ . The function j σ is is continuous, monotonically increasing and concave. Therefore, (13) becomes Now, we set and see that its inverse is continuous, monotonically increasing and convex. Here, j σ −1 : R → R, With this notation, equations (15) lead to the system where Since j σ −1 is monotonically increasing and j σ −1 (R) = R, there exists r d such that j σ −1 (r d ) = d. Moreover, r d does not depend on σ, because since d > 0, by (18) it follows that r d = d. Then, We prove that (19)-(21) has a solution which will be used to define the solution to the equivalent system (15).

3.1.
Hypotheses. In the subsequent part we make the assumptions that there exist M > 0 and L b > 0 such that and |b(r) − b(r)| ≤ L b |r − r| , for any r, r ∈ R.
(28) A self-organized criticality mathematical model based on the sand-pile CA was studied in [3] using the theory of m-accretive operators in L 1 (Ω) approach and then a time difference scheme to get a more regular solution for data in L 2 (Ω). Here, we shall use m-accretive operators in the H −1 (Ω) approach.
3.2. Functional framework. We denote by V the Hilbert space H 1 0 (Ω) with the standard scalar product and by V ′ its dual H −1 (Ω), on which we consider the scalar product We observe that by (29) we have Let σ > 0 be fixed and consider the Cauchy problem dz dt The operator A σ is defined from V ′ to V ′ , with the domain and some η σ (x) ∈ β σ (z(x)) a.e. x ∈ Ω.
We note that a solution to (31) is a solution in the sense of distributions to (19)- (21).

STELIAN ION AND GABRIELA MARINOSCHI
Definition 3.2. Let p 0 ∈ L 2 (Ω). We call a solution to (5)-(7) a function such that which is obtained as the strong limit in L 2 (Q) by with z σ is the solution to (31).
First, we prove the existence of the solution z σ to (31).
To this end, we show that A σ,λ is quasi m-accretive on V ′ . Indeed, let λ ′ > 0 and calculate where η σ,λ (x) = β σ,λ (z(x)) a.e. x ∈ Ω and I is the identity operator. But, by (25) we have and so we get λ . In these calculations we used (29) and (30). Now, let ξ ∈ V ′ and show that the equation has a solution z ∈ D(A σ,λ ) = V. This equation can be still written as where Therefore, by denoting we get z = G(z).
We have to prove that G is a contraction on L 2 (Ω). First, we derive an estimate by multiplying by z − z scalarly in V ′ the difference of eqs. (41) written for the pairs (ξ, z) and (ξ, z). We get Hence, the operator In conclusion, (44) has a unique solution and this implies that (41) has a unique solution.
If p 0 ∈ L ∞ (Ω), then where with M the upper bound in (27).
Proof. Let p 0 ∈ L 2 (Ω). If p 0 ≥ p R a.e. in Ω, then z 0,σ = j σ (p 0 − p R ) = p 0 − p R ≥ 0 a.e. on Ω, and since the subset on which p 0 (x) < p R (where z 0,σ = p0−pR σ ) is of zero measure, we have Then, Problem (19)-(21) has a unique nonnegative solution given by Theorem 3.3, with the estimates (62), (63) turning out to be independent of σ. Then, on a subsequence {σ → 0} it follows that and we have the estimate obtained by (62) by passing to the limit and using the weakly lower semicontinuity of the norm. Moreover, since B(z σ ) ≤ M , by (19), we get that dzσ dt σ lies in a bounded subset of L 2 (0, T ; V ′ ) and on a subsequence It follows that z σ → z strongly in L 2 (0, T ; L 2 (Ω)) and z σ (t) → z(t) strongly in V ′ , uniformly in t ∈ [0, T ]. We also obtain that B(z σ ) → B(z) strongly in L 2 (Q).
The relation remains true for z(t) due to the strongly convergence. Thus, we get (64) as claimed.

4.1.
Numerical simulations for the CA model. We recall that the classical sand-pile model (see e.g., [26]) is represented by the equation The model we propose here differs essentially from this one due to two new rules: a) in each cell there exists another threshold p ij R smaller than the activation threshold p ij I which acts as a barrier against the disease transmission; b) the introduction of a production term b ij (p ij ) which maintains the infective state in a cell.
We also note that if b ij and p ij R are zero then the model reduces to the sand-pile approach.
The steady state dynamics of such a system can be characterized by the probability distributions related to the occurrence of relaxation clusters of a certain size, duration, and area. etc. In the critical steady state these probability distributions exhibit a power-law behavior. A first analyze we shall do aims at putting into evidence power-laws probability distributions. However, we mention that has a certain degree of uncertainty as explained in [18].
Renumbering the cells, and denoting A := {i|p i ≥ p I } the set of active indices the dynamics of the model reads where N (i) is the set of neighbors. We analyze three statistics of our automaton, namely, the size S (number of toppling events), the duration T (number of update sweeps needed until all sites are stable again) and the area A (number of activated cells involved in CA running) of an epidemic episode. We prepare a critical state P 0 of the automaton by initializing each cell of the automaton with a random number between (p I , 3p I ) and then stabilizing it following the classical rules of the sand-pile (73). Using P 0 we generate several cases of epidemic increasing by a unity a random selected cell of P 0 and then stabilizing it following the rules (74).
Scenario 1: An active cell contaminates its susceptible neighbors. If its initial level is higher than a threshold p D its level is set below p R (quarantine), if its initial level is smaller than p D then its level is increased by a quantity, due to the production mechanism, and then the level is decreased by a number equal with the number of the neighbors of the cell (the diffusion mechanism). The production function b(p i ) and supply function s(p i ) have the general form This statistics is given in Figure 1 ( for p R = 3) and Figure 2 ( for p R = 4) where the log of the distribution function (number of cases with S > s, T > t, A > a) is represented against log s, log t, log a, respectively, for the parameters: α = 1, p D = 8, p I = 5. The power-law function f is fitted using the appropriate points which may lead to relevant values of the probabilities P (S > s), P (T > t), P (A > a), respectively. Let us look at the left graphic in Figure 1. The relevant interval for log a in this figure can be considered between (2, 4.5) which produces log values for the probability distribution between (−2, 0), with corresponding probabilities in the interval (10 −2 , 1). In the calculus of the power-law function one can neglect the points with log a > 5 giving probabilities lower than 10 −2 . The occurring uncertainty in this analysis is due to the various possibilities of the extension of the fit region in each distribution, so that in principle it is not possible to obtain the exponents with a high accuracy (see also detailed comments in [18]).  The corresponding graphics are represented in Figure 3 for p D = 8, p I = 5, p R = 3 and in Figure 4 for p R = 4.  However, it is interesting to compare the exponents of the avalanche distributions ruled by the classical sand-pile with the exponents obtained by us using the rule (74) and to see if they are in the same universality class. The probability distributions of the three statistics S, A and T in the previous simulations are approximated as In Table 1 we compare our results with the exponents of the sand-pile (BTW) model ( [26]) calculated by Lübeck and Usadel, [18].
One observes that in the case of p R = 3 the results obtained in both scenarios are close to the BTW model, but they are quite different from the case p R = 4. We can comment that in the case p R = 4 there exists a very low probability to have large contamination area. In Figure 5 the snapshots of the area of contamination are presented. The blue spots in the left and middle pictures correspond to different episodes of epidemic and the pink spots correspond to a common area of many episodes. Each cluster represents the extension of the epidemic, that is the area occupied by the infected cells during an episode developed according the rules (74). The right picture represents the area of an avalanche in the classical sand-pile. One notes that the average area of the spots in the left picture, where p R = 3, is greater than that in the middle, where p R = 4. Instead of the equation (19) we consider the regularized equation given by with the boundary conditions and initial data z| Γ = 0; z| t=0 = z 0 (x, y) (the subscript σ is dropped). We build up a discrete model by using the finite volume method to approximate the space derivatives and the Heun's method for time integration [12]. Let N Figure 6. Epidemic evolution for m p = 1 : z(t, x, 0, 5) at various t (left); z(t, 0.5, 0.5) (middle); infected area a(t) (right) the function z(t, 0.5, 0.5) for t ∈ [0, 0.008]; at right: the infected area a developed with respect to time. These graphics indicate an initial extension of infected area followed by its reduction until extinction.
On the contrary, Figure 7, corresponding to m p = 100, indicates a continuous extension of infected area up to a steady state. Figure 7. Epidemic evolution for m p = 100 : z(t, x, 0, 5) at various t (left); z(t, 0.5, 0.5) (middle); infected area a(t) (right) 5. Conclusions. We introduce a new model intended to illustrate a map of an epidemic spread, focusing on the contamination process and simulating the disease propagation by the means of a phase function. The model starts from a cellular automaton exhibiting a self-organizing criticality process with two thresholds.
In the CA model (1) we focus especially on the contamination process, that is for the passing from S to I, in order to follow the epidemic propagation in time and space. The self-organizing character of the model resides in the terms involving the Heaviside function H(p kl (t) − p kl I ) which shows that only active cells with p kl (t) ≥ p kl I can initiate the transmission to neighbor cells. The other cells with p kl R < p kl (t) < p kl I are inactive and their state remains constant as long as in their neighborhood there is no active cell. Cells having p ij = p ij R preserve this value for all epidemic duration. In fact, this was the reason of the introduction of the class R and of the threshold p R , which acts as a barrier against the disease propagation.
The model is focused on the passing from S to I (disease transmission and spread) and less on the passing from I to S which eventually happens in a cell. Also, the model does not expressly focus on the immunity getting, that is the direct passing from I to R. In principle, one could allow the infected cells to reach the p ij R state since an instantaneous jump from a value greater than p ij I to the value p ij R might be possible, depending on the rates a ij (t) and b ij (t) (e.g., by diffusion when b ij = 0).