Cellular automata in the light of COVID-19.

Currently, the world has been facing the brunt of a pandemic due to a disease called COVID-19 for the last 2 years. To study the spread of such infectious diseases it is important to not only understand their temporal evolution but also the spatial evolution. In this work, the spread of this disease has been studied with a cellular automata (CA) model to find the temporal and the spatial behavior of it. Here, we have proposed a neighborhood criteria which will help us to measure the social confinement at the time of the disease spread. The two main parameters of our model are (i) disease transmission probability (q) which helps us to measure the infectivity of a disease and (ii) exponent (n) which helps us to measure the degree of the social confinement. Here, we have studied various spatial growths of the disease by simulating this CA model. Finally we have tried to fit our model with the COVID-19 data of India for various waves and have attempted to match our model predictions with regards to each wave to see how the different parameters vary with respect to infectivity and restrictions in social interaction.

Introduction

Epidemics and pandemics have a long story throughout human history. Recently human civilization has faced another pandemic named COVID-19. This pandemic has affected many countries through multiple waves. Total of 504,451,689 people have been infected worldwide and 6,222,430 people have died due to COVID-19 till 17 April 2022. In India 43,042,097 people have suffered and 521,781 have died as of 17/04/2022 due to this disease [1]. COVID-19 is caused by the virus which is named SARS COV-2. Multiple variants of this virus, like delta, omicron and many others makes it harder to control and predict its behavior. Recently another variant of COVID-19 named XE has been found [2]. Mathematical modeling helps us to understand the behavior of disease spread such that prevention and control strategies can be built. In addition, mathematical models can help us to find some inherent properties of the disease and nature of its spread.

There are many different types of models that have been used in the past to study various diseases. These models are mainly modified versions of the Kermack McKendrick SIR model which is based on a system of coupled ordinary differential equations [3]. Currently, the ODE-based models and statistical models are widely used in literature to model the temporal behavior of the spread of COVID-19 from different aspects. Most of these models have tried to analyze the spread of this disease and tried to predict its future behavior [4-7]. There are models which have proposed various intervention and vaccination strategies to prevent and control the spread of the disease [8-12]. Also, some authors have tried to predict different inherent properties of this pandemic like herd immunity and its chaotic nature [13-17]. These temporal models can give us much valuable information, however, most of these models assume that a population is homogeneously mixed and cannot describe any spatial behavior. To incorporate this spatial behavior, deterministic and probabilistic spatio-temporal models have been used in recent studies. Cellular automata (CA) is one such kind of spatio-temporal model.

Cellular automata (CA) has been used in many studies to model different aspects of epidemics. It has been widely used to model the disease spread of influenza and various vector-borne diseases, such as dengue [18-26]. A neighborhood condition is an important aspect in the CA. The most used neighborhood conditions are (i) Neumann’s neighborhood condition, (ii) Moore’s neighborhood condition, (iii) Extended neighborhood condition, and (iv) Random interactions. Coupled with these neighborhood conditions, various models such as SEIR, SEIRS, SEIRD, and SEIRQD have been studied with the help of CA to model the spatial growth of epidemics [19, 22, 27–32]. Currently, CA has gained a lot of momentum in the studies of COVID-19 [33-35]. Various advanced studies with Genetic algorithms and network models have been done for COVID-19 data [35-38].

In this paper, we have tried to model COVID-19 using cellular automata (CA) to find spatio-temporal behavior of it. We have also made some analysis to understand its behavior in different waves of the disease. A cellular automata (CA) model is represented in a square lattice and defined by some neighborhood and boundary conditions, the details of which are given in following sections.

This paper is arranged as follows: Sect. 2 consists of a detailed discussion of the model, neighborhood conditions, and probability of infection. The algorithm, results of the simulations have been shown in Sect. 3 and the data analysis is shown in Sect. 4. Finally, Sect. 5 consists of the conclusions of our model.

Mathematical model

In this article, we have illustrated a Cellular automata (CA) model for epidemics and assumed the SEIR model as the base model. SEIR model stands for Susceptible–Exposed–Infectious–Removed. Here we have considered a square lattice, where each cell of the lattice is assumed to be a person. Each cell of the lattice can have the set of states, and these states are represented by the values . The updation of a cell’s state depends upon various conditions such as (i) the current state of the cell, (ii) the amount of time spent in the current state, and (iii) the current states of the neighbors. The main assumptions of our model are given below:

Every cell represents a person.

Only susceptible persons can interact with the other cells.

A removed person cannot be infected again.

For this CA model, we have assumed a periodic boundary condition. If a cell of ith row and jth column of a lattice is denoted by (i, j) then

One susceptible person can interact with a single person in each time step.

Neighborhood condition

Nearest neighborhood condition is a widely used concept in the literature. Here, it has been assumed that a particular cell can only interact with its nearest neighborhood cells. Such two famous neighborhood conditions are: (i) Neumann’s neighborhood condition and (ii) Moore’s neighborhood condition.

Figure 1 shows the two neighborhood conditions. Figure 1a shows Neumann’s neighborhood condition, where the nearest neighborhoods of any chosen (i, j) cell are the first neighborhood cells with respect to the chosen cell. Similarly, Fig. 1b shows Moore’s neighborhood condition. In this case, all first and second neighborhoods are treated as the nearest neighborhoods of the chosen (i, j) cell.

Fig. 1

Different neighborhood conditions. Blue cells represent the nearest neighborhoods of (i, j) cell

In this work, we have assumed that a cell can interact with any other cells depending on the probability of interaction () between them. Here we have assumed that the probability of interaction () of a cell (i, j) to any other cell varies inversely as a function of d (distance between two cells) in the form of a power law. Hencewhere n is the degree exponent and can have a value greater than zero. Here, we have assumed that the distance between two cells is not just the geometrical distance between two. It depends on the layer number (l). In Fig. 2, we have shown how layers are defined. In addition, it shows that a layer l contains the 8l number of cells. If we choose a lattice of size then the total number of layers in this lattice with respect to any cell is , when N is odd.

Fig. 2

Different layers of a lattice with respect to (i, j) cell

Hence, we can writeFor mathematical simplicity, we can assume . Hence, from Eq. 2 we can writeIf there are L number of layers then the above equation can be written aswhere . Hence, a person at the (i, j) cell can interact with any other cell of layer l with a probability . Thus, average interaction distance () can be defined asFigure 3a shows the variation of average interaction distance () with the degree exponent (n). From this figure, it can be found that the average interaction distance () is quickly decreases and saturates to unity as n increases. In addition, from Fig. 3b it can be seen that the average interaction distance () is approximately for exponents . Hence for the neighborhood condition is approximately similar to Moore’s neighborhood condition as discussed earlier and does not give any significantly different results.

Fig. 3

Plots of average interaction distance () with n and the total number of layers (L). a Plot of average interaction distance () with n by considering the total number of layers (L) = 50. b Plot of average interaction distance () with the total number of layers (L) for n=3,4, and 5

Probability of infection ()

Let, q denotes the disease transmission probability when a susceptible and an infectious person interact. The probability that a susceptible person will interact with any person at the layer l is . If the probability of finding an infectious person in that layer is then the probability that the susceptible person will be infected is . Hence, the probability of infection () of a susceptible person isAs, , from the above equation (Eq. 7) we can writeFrom the above equation (Eq. 8), we can say that the terms with small layer number (l) dominate the summation. Hence, the infection possibility of a susceptible person mainly depends on the infection situation around the person.

Thus in our model, instead of choosing a traditional neighborhood condition, where the degree of the interaction is fixed, we have assumed a model, where we can vary the degree of the social confinement by changing n (degree exponent). In addition, we have calculated the probability of infection () for this modified model.

Algorithm and simulations

Algorithm

Here, we have discussed the state updation algorithm of the SEIR model. As we have mentioned earlier, every cell’s state is denoted by a value (v) which is present in this set . The algorithm is given below:

Let, at time t there is a susceptible person at (i, j) cell. Therefore, the value of the (i, j) cell is and the probability of infection is . To find the infection possibility of the susceptible person, we will generate a uniform random number u between 0 and 1.

If, then the susceptible person is exposed and at time the state of the (i, j) cell will be changed from to .

else, at time the state of the (i, j) cell will be unchanged.

An exposed person () will remain exposed for number of days. After that, the person will be infectious and the state of the corresponding cell will be changed from to .

An infectious person will remain in this state for number of days. After that, the person will be removed (recovered or dead) and the state of the cell will be changed from to .

Simulation

In this part, we have done simulation of our model with . The values of the parameters and initial conditions that are used in the simulations are listed in the tables (Tables 1 and 2) below:

Table 1

Table for the parameter values that are used in the simulations

Description of the parameters	Parameters	Values of the parameters
Lattice size	N×N	101×101
Disease transmission probability	q	0.3
Latency period of the disease	τ_I	8 days
Removal period	τ_R	18 days

Table 2

Table for the initial conditions of the simulations

States	Initial values
S(0)	10200
I(0)	1
E(0)	0
R(0)	0

In Fig. 4, the first plot shows the temporal behavior of the exposed cases (E(t)) and the infectious cases (I(t)) of the epidemic. These temporal plots are averaged on 50 simulation samples. The rest plots of Fig. 4 are CA plots that represent the spatial evolution of disease spread. From these CA plots, we can hardly detect any clustering of the infected cases. This happens, because the average interaction distance, . Thus a susceptible person can be infected by an infectious person who is far away from the susceptible one.

Fig. 4

Plots of the temporal and spatial behavior of the disease spread for

In Fig. 5, the first plot again shows the temporal growth of the epidemic. It can be seen from the CA plots that for , clusters are formed. The reason behind this is the short average interaction distance. For , the average interaction distance, . In addition, from the temporal plot, we can see that the infection spread time is increased than in the case. This is also because of the short average interaction distance. For a short average interaction distance, only a few susceptible persons can interact with the infectious person. So, if most of those susceptible persons become infected then the infectious person cannot spread the disease further. Whereas, for case, an infectious person can interact with many susceptible persons. Hence, an infectious person can infect more people during the infectious period.

Fig. 5

Plots of the temporal and spatial behavior of the disease spread for

Figure 6 shows the evolution of the disease for . Here we can also clearly find the clusters. These clusters are more prominent than the case because of the less average interaction distance (). Value of the average interaction distance for is . Here we can see that the disease takes a longer time to fall for than for and . The reason behind this is the lower value of the average interaction distance which is discussed earlier.

Fig. 6

Plots of the temporal and spatial behavior of the disease spread for

Hence, from the above discussions, we can conclude that the clustering behavior of the disease spread depends on the average interaction distance () as well as on degree exponent n. In addition, the average interaction distance (), gives an average estimation of the number of susceptible persons who can interact with an infectious person which is represented by . Therefore, for a large (or small n) an infectious person can spread the disease to distanced region. Thus the infection period depends on the average interaction distance () and also on n.

Comparison with data

In this section, we have tried to fit our model with current COVID-19 data. Our model has four free parameters which are, (i) q : disease transmission probability, (ii) n : degree exponent, (iii) : mean latency period, (iv) : mean infectious period. We have optimized these free parameters for different waves of the COVID-19 pandemic in India. Here, we have considered each wave separately and normalized the active cases of each wave with the total number of infected cases in the respective wave. The data is taken from covid19india.org [39]. The date range of the different waves that we have considered here are given in Table 3.

Table 3

Date ranges for different waves

Waves	Start date	End date
First wave	30-Jan-2020	16-Feb-2021
Second wave	17-Feb-2021	31-Oct-2021

To fit the model with the data we have optimized the sum of squared errors (SSE): : fraction of the active cases from the data. : fraction of the infectious cases from the model.

The results of the best fit parameter values are given below.

We have optimized COVID-19 data with a lattice space. Figure 7a, b shows the fitted model along with data for the first and the second wave.

Fig. 7

Model fitting of the COVID-19 data of India a Model fitting of the first wave. b Model fitting of the second wave

Here we can see that the model fits reasonably well with the data. Also from Table 4, we see that the disease transmission probability (q) increases and the degree exponent (n) decreases in the second wave as compared to the first wave. Increment of q represents that in the second wave of the disease it was possibly more infectious and spread faster than in the first wave. Also decrement of n indicates that the interaction between infectious and the susceptible population spread out to a larger distances in the second wave as compared to the first wave. This is possibly because of the COVID protocols which is relaxed much more in the initial phase of the second wave as compared to the first wave.

Table 4

Fitting parameter values for different waves

Waves	q	n	τ_Idays	τ_Rdays
First wave	0.1950	1.9310	5	11
Second wave	0.2406	1.3449	8	10

Conclusions

In this section, we have summarized the main features and results of our model. The cellular automata (CA) is a very common tool to model a disease spread and has been used extensively in literature for studying different systems. In this paper, we have modeled the CA by proposing a new neighborhood criteria. Usually in earlier studies, neighborhood condition is such that the neighborhood of a lattice cell is always fixed. Whereas in our model, rather than choosing a specific neighborhood condition, we assume that a lattice cell can interact with any other cell at distance d with a certain probability which is called interaction probability () (Eq. 5). We have assumed that the interaction probability () is a function of the distance (d) and has a form of inverse power law with degree exponent n. Here, exponent n is a very important parameter as it enables us to tune the social confinement of our model. With this newly defined neighborhood criteria, we have calculated various relations such as average interaction distance () and the probability of infection () to understand and represent our model properly.

From Fig. 3a, we can see that the average interaction distance () deceases and saturation is reached to as n increases. Therefore, for a higher n, a person can mainly interact with nearest neighbors. However, for a smaller value of n, a person can interact with the distant neighbors. Hence higher values of n represents higher social confinement and viceversa. Also we want to mention that for exponents, , the average interaction distance, . Thus the values do not give us any new results.

In the simulation section, we have studied the temporal and spatial behavior of our model for different degree exponents (n). As n increases, the disease spread becomes slower and it is more clustered. This happens because of the decrease in the average interaction distance () or in the other words increase in the social confinement with increasing values of n. Thus the disease transmission probability (q) and the degree exponent n regulates the speed of the disease spread.

Also, we have compared our model with COVID-19 data of India for different waves. We have first normalized the active cases of a wave with the total number of infected cases in that wave. Then we have optimized the sum of squared errors of the fraction of infectious cases (Eq. 9) to get the best fit result with the data. Here all simulations are done on a lattice. We have found that the disease transmission probability (q) increases in the second wave than the first wave. This means that the disease is more infectious in the second wave than the first wave. In addition, the degree exponent (n) decreases in the second wave. This implies that the decrement in the COVID-19 restrictions (or decrement in the degree of social confinement) at the initial time of the second wave played a significant role to the faster spread of the disease. Our model fits the peak of the waves well, however, fall in the data at the end of both waves and plateauing at the end of the second wave have some mismatch with our fitted model. This possibly indicates that our model needs to be modified to incorporate these aspects which we will look at in our future works.

Modeling spread of a disease is a very complicated process, since several factors have to be considered. Non-uniform distribution of the population and economic situation of the regions are two major factors which affects the disease spread. In future we want to look at these complex aspects by refining this model to find the behavior of the disease spreading more accurately. We would also like to study these possibilities in the context of COVID-19.