**2. Model**

We use the cellular automata technique [21–24] to model disease propagation. The cellular automata technique is based on several assumptions, including:


The rule F defines the state *si* of the site *i* at time (*t* + 1) based on this site's state *si* at time *t* and the state of the sites in the *i*-th site's neighborhood N

$$s\_i(t+1) = \mathcal{F}(s\_i(t); \mathcal{N}(s\_i(t))).$$

We adopt the SEIR model for the simulation of disease spreading by probabilistic synchronous cellular automata on a square lattice with various neighborhoods N . Every agent may be in one of four available states: susceptible (S), exposed (E), infected (I) or recovered (R). The agents in E and I are characterized by different values of the probability of infection (*p*<sup>E</sup> and *p*<sup>I</sup> ) and different ranges of interactions (radii of neighborhoods *r*<sup>E</sup> and *r*<sup>I</sup> ). The considered neighborhoods (and their radii) are presented in Figure 1.

**Figure 1.** Sites in various neighborhoods N on a square lattice. (**a**) Von Neumann's neighborhood (*<sup>r</sup>* <sup>=</sup> 1, *<sup>z</sup>* <sup>=</sup> 4); (**b**) Moore's neighborhood (*<sup>r</sup>* <sup>=</sup> <sup>√</sup>2, *<sup>z</sup>* <sup>=</sup> 8); (**c**) neighborhood with sites up to the third coordination zone (*r* = 2, *z* = 12); (**d**) neighborhood with sites up to the fourth coordination zone (*<sup>r</sup>* <sup>=</sup> <sup>√</sup>5, *<sup>z</sup>* <sup>=</sup> 20); (**e**) neighborhood with sites up to the fifth coordination zone (*<sup>r</sup>* <sup>=</sup> <sup>2</sup> <sup>√</sup>2, *<sup>z</sup>* <sup>=</sup> 24).

Initially (at *t* = 0), every agent is in the S state, their age *a* is set randomly from a normal distribution with the mean value *a* = <sup>50</sup> × 365 days and dispersion *<sup>σ</sup>*<sup>2</sup> = <sup>25</sup> × 365 days. 'Patient Zero' in the E state is placed randomly at a single site of a square lattice with *L*<sup>2</sup> = 1002 nodes.

Every time step (which corresponds literally to a single day in the real-world) the lattice is scanned in typewriter order to check the possible agent state evolution:

• The susceptible agent may be infected (S→E) by each agent in state E or I present in his/her neighborhood N with radius *r*<sup>E</sup> or *r*<sup>I</sup> ≤ *r*<sup>E</sup> , respectively. We set *r*<sup>I</sup> ≤ *r*<sup>E</sup>

as we assume that infected (and aware of the disease) agents are more responsible than exposed (unaware of the disease) ones. The latter inequality comes from our assumption that the exposed agent is not careful enough in undertaking contacts with his/her neighbours while the infected agent is serious-minded and realizing the hazard of possible disease propagation and thus he/she avoids these contacts at least on the level assigned to exposed agents. The number and position of available neighbours who may infect the considered susceptible agent depend on the value of radius *r*<sup>E</sup> and/or *r*<sup>I</sup> as presented in Figure 1. The infection of the susceptible agents occurs with probability *p*<sup>E</sup> (after contacting with agent in state E) or *p*<sup>I</sup> (after contacting with agent in state I), respectively.


The agents' state modifications are applied synchronously to all sites. A single time step (from *t* to *t* + 1) of the system evolution described above is presented in Algorithm 1. The *L*2-long vector variables tmp\_pop[] and pop[] represent the current population (at time *t*) and the population in the next time step (*t* + 1), respectively. The presence of two such variables in a model implementation is caused by the synchronicity of cellular automaton. The *i*-th element of these vectors keeps information on the state (either S, E, I or R) of the *i*-th agent in the population. The age *a* (measured in days) for the *i*-th agent is kept in the *L*2-long vector age[] and its *i*-th element is either incremented (line 4 of Algorithm 1) or reset (lines 28 and 35 of Algorithm 1) in the case of removal. The random() function returns a real pseudo-random number uniformly distributed in [0, 1). The function *i*, *j* measures the Euclidean distance between agents *i* and *j*. The age-dependent daily death probability functions *fC*(·) and *fG*(·)—based on real-world data—are defined in the next paragraph in Equations (5) and (6). The numbers *n*S, *n*<sup>E</sup> , *n*<sup>I</sup> and *n*<sup>R</sup> of agents in various states must be initialized at *t* = 0 based on initial conditions. The cumulative number *n*<sup>D</sup> of deceased agents—earlier either in the exposed (E) or in the infected (I) state—is incremented in line 39 of Algorithm 1. The numbers *n*S, *n*<sup>E</sup> , *n*<sup>I</sup> , *n*<sup>R</sup> and *n*<sup>D</sup> are normalized to the total number *L*<sup>2</sup> of agents in the system before their return.

Based on American data on annual death probability [44] and assuming 365 days a year, we predict the daily death probability as:

$$f\_{\mathbb{G}}(a) = 184 \cdot 10^{-10} \exp(0.00023(a + 40259)),\tag{5}$$

where *a* is the agent's age expressed in days. The data follow Gompertz's exponential law of mortality [46,49].

Using the same trick, we estimate the probability of daily death for infected people of age *a* (expressed in days) as:

$$f\_{\mathbb{C}}(a) = \begin{cases} 5 \cdot 10^{-5} & \iff a \le 30 \text{ y}, \\ 2 \cdot 10^{-6} \exp(0.0003a) & \iff a > 30 \text{ y}. \end{cases} \tag{6}$$

These probabilities are calculated as the chance of death during SARS-CoV-2 infection (based on Polish statistics [50]) divided by (*τ*<sup>E</sup> + *τ*<sup>I</sup> ). We assume that infection lasts *τ*<sup>I</sup> ≈ 14 days and that an incubation process takes *τ*<sup>E</sup> ≈ 5 days [39]. Exponential fits (Equations (5) and (6)) to real-world data are presented in Figure 2.

```
Algorithm 1 Single time step in automaton
1: t ← t + 1
2: tmp_pop ← pop
3: for all i ∈ pop do
4: age[i] ← age[i]+1
5: if tmp_pop[i] = S then
6: for all j ∈ pop, j = i do
7: if tmp_pop[j] = I∧i, j ≤ rI∧ random()≤ pI then
8: pop[i]← E  S→E after contacting with I
9: nS ← nS − 1
10: nE ← nE + 1
11: break
12: if tmp_pop[j] = E∧i, j ≤ rE∧ random()≤ pE then
13: pop[i]← E  S→E after contacting with E
14: nS ← nS − 1
15: nE ← nE + 1
16: break
17: if tmp_pop[i] = E ∧ random() < 1/τE then
18: pop[i]← I  E→I
19: nE ← nE − 1
20: nI ← nI + 1
21: if tmp_pop[i] = I ∧ random() < 1/τI then
22: pop[i]← R  Recovered
23: nI ← nI − 1
24: nR ← nR + 1
25: if tmp_pop[i] = S ∨ tmp_pop[i] = R then
26: if random() < fG(age[i]) then
27: pop[i]← S  Removed. . .
28: age[i]← 0  . . . then Susceptible
29: if tmp_pop[i]=S then nS ← nS − 1
30: if tmp_pop[i]=R then nR ← nR − 1
31: nS ← nS + 1
32: else  tmp_pop[i] = E ∨ tmp_pop[i] = I
33: if random() < fC(age[i]) then
34: pop[i]← S  Removed. . .
35: age[i]← 0  . . . then Susceptible
36: if tmp_pop[i]=E then nE ← nE − 1
37: if tmp_pop[i]=I then nI ← nI − 1
38: nS ← nS + 1
39: nD ← nD + 1  cumulative number of deaths caused by infections
40: return(t, nS/L2, nE/L2, nI/L2, nR/L2, nD/L2)
```
**Figure 2.** Daily death probability *f*(*a*) for patients infected by the coronavirus (SARS-CoV-2, ×, *fC*(*a*), [50]) and natural death probability (+, *fG*(*a*), [44]).
