Estimates of the collective immunity to COVID-19 derived from a stochastic cellular automaton based framework.

In the context of the propagation of infectious diseases, when a sufficient degree of immunisation is achieved within a population, the spread of the disease is ended or significantly decreased, leading to collective immunity, meaning the indirect protection given by immune individuals to susceptible individuals. Here we describe the estimates of the collective immunity to COVID-19 from a stochastic cellular automaton based model designed to emulate the spread of SARS-CoV-2 in a population of static individuals interacting only via a Moore neighbourhood of radius one, with a view to analyze the impact of initially immune individuals on the dynamics of COVID-19. This impact was measured by comparing a progression of initial immunity ratio-the percentage of immunised individuals before patient zero starts infecting its neighbourhood-from 0 to 95% of the initial population, with the number of susceptible individuals not contaminated, the peak value of active cases, the total number of deaths and the emulated pandemic duration in days. The influence of this range of immunities over the model was tested with different parameterisations regarding the uncertainties involved in the model such as the durations of the cellular automaton states, the contamination contributions of each state and the state transition probabilities. A collective immunity threshold of 55 % ± 2.5 % on average was obtained from this procedure, under four distinct parameterisations, which is in tune with the estimates of the currently available medical literature, even increasing the uncertainty of the input parameters.

Introduction

The Coronavirus Disease 19 (COVID-19) is a contagious disease which can be transmitted via interactions between human beings, by means of particles emanated from the respiratory system of a contaminated person to the mouth or nose of another person. Contaminated individuals present symptoms like fever, body aches and shortage of breath, which are caused by the Severe Accute Respiratory Syndrome-Coronavirus (SARS-CoV-2) (Ezhilan et al. 2021). This virus was initially found in animals, and the origin of its transition to human beings is still unknown (Velavan and Meyer 2020).

Generally, the mass contamination caused by this virus began in December of 2019 in the province of Hubei, located in the People’s Republic of China, initally diagnosed as a type of pneumonia. In the subsequent months, the World Health Organisation (WHO) declared a pandemic state and authorities from each country began to adopt health and safety protocols trying to contain COVID-19’s transmission (Velavan and Meyer 2020). These authorities have been relying on different types of quarantine, from social distancing to lockdowns, so that the effects of these decisions during the last year and a half have led to worldwide decrease in economy, poverty raise and human lives (Ciotti et al. 2020). Consequently, many economic activities were affected and damaged (Ezhilan et al. 2021), prompting political and sanitary agents, and the population in general, to start asking when these restrictions could be removed without causing a new wave of infections.

Understanding when to start removing the restrictions involves essentially understanding when and how collective immunity can be achieved (Britton et al. 2020), that is, the individual protection from infection given by a sufficiently large proportion of immune individuals in a society, decreasing the infection rate, or the average number of individuals infected by one infected person, to values under 1.0 (Randolph and Barreiro 2020).

Studying the behaviour of diseases in general, specially in planning prevention measures for authorities can cost time and human lives. But it is possible to resort to theories and models related to epidemiological studies, to develop mathematical and computational approaches to help predict and test control strategies, as well as understand the impact of immune individuals in a disease. In this context, cellular automata based models constitute a paradigmatic example, as their use in disease spreading studies have been shown to be both viable and reliable (Mikler et al. 2005; Monteiro et al. 2020; White et al. 2007; Schimit 2021).

Cellular automaton (CAs) are discrete dynamical systems defined by a regular lattice of cells, where each one is defined by a finite number of discrete states. The dynamics of cellular automata is defined locally, with each cell having its state changed according to its current state and those of its neighbouring cells. And even simple individual cell behaviours, cellular automata may lead to arbitrarily complex global behaviours. On their part, stochastic cellular automata are the ones whose cells define their next state based on a probability distribution.

Here, a computational model based on stochastic cellular automata is proposed to aid in understanding the COVID-19 dynamics, with a focus on collective immunity. We describe the principles of a stochastic cellular automaton which allows to emulate the spreading of COVID-19 in a population, so as to obtain the proportion of the population that has to be vaccinated, or become naturally immunised, for collective immunity to be established. Even though one might argue that the regular topology of cellular automata might not be adequate to reliably support the dynamics of a disease spreading, this has been the case in various successful efforts in the literature (Yakowitz et al. 1990; Fuentes and Kuperman 1999; Schimit and Monteiro 2009; Holko et al. 2016; Schimit and Monteiro 2011; Ahmed and Agiza 1998; Sirakoulis et al. 2000). So, in tune with those, our main goal here is to probe whether such an abstract and minimal model could capture a fundamental aspect of the disease like the collective immunity.

In Sect. 2, the computational model based on the proposed cellular automaton is discussed, and is followed by a section on the methodology underlying the computational experiments carried out. In Sect. 4, the results obtained are presented and discussed, according to four distinct parameterisations setups, hereby called Experiments. Finally, in Sect. 5 the results obtained are discussed and compared to related works.

The reference model and its computational implementation

The reference model

During the ongoing COVID-19 pandemic, the following four conditions were observed in contaminated individuals, related to the duration and the intensity of the disease, according to COVID-19 related literature at the beginning of the pandemic (Ferguson et al. 2020; Lauer et al. 2020; Liu et al. 2020; Wölfel et al. 2020):

No symptoms Individuals who do not present any symptoms of the disease but can contaminate others through social interactions;

Moderate symptoms Describing the ones who present a few symptoms but quickly become immune, but still with the possibility of contaminating others through social interactions;

Severe symptoms Those who have to be hospitalised and present a moderate death risk;

Critical symptoms Individuals in need of hospitalisation in an Intensive Care Unit (ICU), maintained by forced ventilation, and presenting a high death risk.

Figure 1 presents the reference model we rely upon (due to Rahmé 2020), based on Ferguson et al. (2020), Lauer et al. (2020), Liu et al. (2020). It illustrates the four previous possible dynamics of the disease after contagion and divides COVID-19 in terms of the stages of the disease, as it progresses. On the left-hand side, the figure displays the proportions of the population following each dynamics, which can be viewed as a trail followed over time, from contagion to full recovery (or death), with their durations, and correlated to the severity of the symptoms. Also, as severity rises, the recovery is not a single destination but a branch with two possibilities, one of them being the death of an individual, which occur at the different death rates expressed on the rightmost part of the figure. After recovery from the symptoms, the figure also shows the periods where the remaining viruses in an individual can still contaminate others, while their antibodies combat the disease. Each trail shown in the figure has internal segments regarding incubation period, manifestation of symptoms, hospitalisation and recovery, with their respective durations; for example, the period associated with No Symptoms, which describes the incubation period of the disease, takes around five days to end.

Fig. 1

Timeline, progression and mortality of COVID-19, from Rahmé (2020), based on Ferguson et al. (2020), Lauer et al. (2020), Liu et al. (2020), Wölfel et al. (2020). This chart describes the general scheme supporting the present work. The values were taken by the disease references at the time the chart was created, in 2020

As patient zero starts interacting with other members of a population, the disease will start spreading and, after some time, every person will be contaminated and undergo COVID-19 effects. Some will become immune, others will die, and others will not even be contaminated. The dynamics proposed and described by this reference model are experienced individually.

The model can be couched in terms of a stochastic cellular automaton (White et al. 2007; Mikler et al. 2005), by representing the variety of possibilities elicited in Fig. 1 as the states of a cell of the CA, the probabilities therein as those governing the state transitions among the states and each day of the timeline expressed as a time step of the automaton. The details about how we go about it are described in the next section.

Representing COVID-19 dynamics with a stochastic cellular automaton

The states

Let us consider a closed population with individuals allocated over a two-dimensional lattice, each one interacting only within those in its direct neighbourhood of the eight surrounding individuals (i.e., in a Moore neighbourhood). Also, let us assume the population has fixed boundaries, assuming a no interactions state to missing neighbours. In other words, the interactions will only be among the cells in a Moore nighbourhood of unitary radius and, if a cell is at the boundary of the lattice, its neighbourhood will only include those inside the lattice.

Consider that one of the individuals is arbitrarily contaminated with SARS-CoV-2 and becomes patient zero for COVID-19, thus entering one of the trails mentioned above (no symptoms, moderate symptoms, severe symptoms or critical symptoms).

Let us represent the cell states of the cellular automaton by the notation , with t being the state duration in terms of the number of time steps of the automaton’s temporal evolution; id being just a state identificator; and standing for the type of state class it belongs to, out of the possibilities Susceptible, Infected and Removed. Notice that these are the same alternatives that define a SIR model of epidemic dynamics (Blavatska and Holovatch 2021).

Based on the previous notation, the reference model of Fig. 1 can now be expressed as follows. The durations of specific states with time dependent transition are obtained from the reference model, also being the mean values of probability distributions (further on this will be addressed in Sect. 3.3). These mean times are obtained from the time axis at the bottom of the chart expressed in days. For example, the duration 13 for state corresponds to duration of the magenta ICU/ventilation segment, just above the time axis in Fig. 1; the segment starts at day 8 and ends at day 21, thus resulting in the state duration of 13 days (to be translated to 13 times steps in the computational model).

S type state:

susceptible (): Individuals (cells) capable of being contaminated, during as long the neighbourhood contamination probability is enough to cause a change to the next state.

I type states:

no symptoms (): Absence of symptoms during the incubation period;

moderate symptoms (): Moderate symptoms manifestation;

severe symptoms (): Severe symptoms manifestation;

critical symptoms (): Critical symptoms manifestation;

hospitalised (): Hospitalisation period after severe symptoms;

ICU/ventilation (): Hospitalisation in intensive care unit, after the manifestation of critical symptoms;

contagious without symptoms (): Disease effects ceased, but it is still possible to contaminate others;

contagious with moderate symptoms (): Disease effects ceased, but it is still possible to contaminate others after moderate symptoms;

contagious with severe symptoms (): Disease effects ceased, but it is still possible to contaminate others after hospitalisation;

contagious with critical symptoms (): Disease effects ceased, but it is still possible to contaminate others after intensive care unit internation.

R type states:

deceased (): Individuals who did not resist the disease effects and were permanently removed from the population, thus no longer interacting;

immune (): Immunised individuals who cannot be contaminated again and do not contaminate others.

Contamination contributions

Notice that each state has its own contamination contribution, which indicates how much a cell can infect a neighbouring cell in the susceptible state, at each timestep. For example, an individual who did not have symptoms, after the incubation period (no symptoms state) would stay for some time contagious, even possibly unaware of it. For this state, the contamination contribution could be considered 100%. The following provides the whole picture we assumed:

susceptible, deceased and immune states were assumed to have 0% contamination contribution;

hospitalised and ICU/ventilation states were assigned 10% contamination contribution (a small rate, since the individuals in these states are under the more controlled conditions of a hospital);

The remaining states were considered with 100% contamination contribution.

State transitions

It is necessary to define and describe the state transitions between states. If any individual in a population gets infected, independently of the dynamics followed, he/she is capable of affecting the neighbourhood, with the contributions expressed in Sect. 2.2.2. In the same way, the neighbourhood is capable of afecting the individual, indicating that there is a correlation between the size of the neighbourhood and the probability of getting ill.

Accordingly, the transitions associated to the susceptible state are defined by the neighbourhood contamination probability , which is drawn from the average of the contamination contributions of each neighbour in a Moore neighbourhood of unitary radius (the sum of the contamination contributions defined by the current states of the eight neighbourhood cells divided by eight). For the other states, the state transitions are governed by the durations of the states (defined in days in the reference model and in time steps of the automaton).

So, we can define the state transitions as follows:The resulting state transition diagram of the model describing each probability and each state transition can be seen in Fig. 2.

Fig. 2

State transition diagram of the proposed cellular automaton. The durations of each state, superscripted in each state identification, are directly obtained from the reference model, but they really correspond in the computational experiments to the mean value of a probability distribution. At the time a state should transit to the next, this occurs with 100% probability, unless multiple different states might follow; in these cases, the transition probabilities are divided among the possible next states. States involved: susceptible (); no symptoms (); moderate symptoms (); severe symptoms (); critical symptoms (); hospitalised (); ICU/ventilation (); contagious without symptoms (); contagious with moderate symptoms (); contagious with severe symptoms (); contagious with critical symptoms (); deceased (); and immune ()

If a cell is in the susceptible state, the next state will be no symptoms with probability ;

The no symptoms state will remain for five time steps. The next possible states are contagious without symptoms with 30% probability, moderate symptoms with 55%, severe symptoms with 10% and critical symptoms with 5%;

The transitions for the remaining states are defined just by the duration in time steps and can be better observed in Fig. 1. As stablished earlier, these durations are defined by probability distributions, as addressed in Sect. 3.3);

The transition to the immune or deceased states is also defined by the previous state durations but with a probability for each one. If the cell is in contagious without symptoms or contagious with moderate symptoms states, the probability of moving to immune state is 100%; if the cell is in hospitalised state, its probability to change to the immune state is 85% and to deceased is 15%; if the cell is in ICU/ventilation state the probability of turning to the immune state is 50% and to deceased is 50%.

Methodology of the computational experiments

Preliminaries

With the state transitions and probabilities defined, the proposed cellular automaton can be implemented and simulated, ending its iterations when a global convergence is achieved. This steady state condition is defined by a global configuration with all the cells in the states susceptible, deceased or immune. The initial configuration of the lattice is assumed to contain a randomly chosen single cell in the no symptoms state, surrounded by cells in the susceptible state or the immune state also randomly distributed (this will be better described in Sect. 3.2).

The characterisation of the dynamics of the temporal evolution of the cellular automaton can then be obtained by monitoring data such as:Consequently, the previous data provides what is needed to study the collective immunity effect related to COVID-19 in a population, specially the estimate of the necessary immunity ratio to achieve collective imunity.

The number of cells in susceptible state, representing the number of individuals not contaminated;

The number of cells in deceased state, representing the number of deceased individuals;

The number of cells in immune state, representing the number of recovered individuals;

The number of cells at the other states, representing the infected individuals.

Analysing collective immunity

With the computational model, the impact of immune individuals can then be studied by varying the initial global configuration of the lattice, which is done by supposing that a vaccination campaign has been performed, that leads to the increase of the initial number of immune cells; its effect is then monitored by four control variables, as mentioned below.

Initially, a lattice with cells in the susceptible state is created. Then, a single cell is randomly chosen to be in the no symptoms state, therefore becoming patient zero and a proportion of the remaining cells is randomly chosen to be in the immune state. In tune with the latter, we define immunity ratio of the lattice as the number of cells changed to the immune state initially, at time step zero, in comparison with the number of cells of the lattice. The rest of them stays in the susceptible state. Since the automaton is stochastic, a number of 10 simulations is performed, each one running until convergence; along the process, the variables below are monitored, and their values averaged over the 10 simulations:In the sequence, this procedure is repeated with increasing values of immunity ratio, i.e., the inital proportion of cells in the immune state. We rely on an immunity ratio range from 0% to 95%, with 5% steps. This procedure should be capable of allowing the analysis of the impact of the immunity ratio in the mentioned variables. Also, it enables the possibility of observing how each variable performs when immunity rises.

The proportion of cells in susceptible state at final convergence;

The maximum proportion of infected cells during a run;

The proportion of cells in deceased state at final convergence;

The number of time steps counted at final convergence.

Estimating collective immunity with different parameters

Aiming at strengthening the model’s flexibility, the contamination contributions, the durations of each state and the transition probabilities were defined as input parameters of the model, not constants. This approach was adopted to allow testing the impact of these parameters on the estimates of collective immunity. So, additionally to the procedure proposed, related to varying the initial immunity, four parameterisations are adopted to make a proper comparison of collective immunities obtained from each one.

The first one—that we refer herein as defining the computational Experiment I—is the same as the one we used previously in Lima and Balbi (2021) and considers the contamination contributions for each state as constants with the values proposed in Sect. 2.2.2. Also, the durations of each state follows a uniform distribution around the values shown in Fig. 2, here considered as means, .

The second setting—defining computational Experiment II—increases the degree of uncertainty of the automaton. The contamination contributions are then defined by a Gaussian distribution truncated by three standard deviations with the mean value equal to the contamination contribution (item 2.2.2). The deviation considered was 3.33%, so the probabilities of contamination for each state can variate by . The durations follow the same logic, being defined by a Gaussian distribution truncated by three standard deviations with the mean value equal to the duration given by Fig. 2. The standard deviation considered was 6.67%, so the probabilities of contamination for each state could vary by .

The third parameterisation—related to computational Experiment III—expands the range of values applied by the previous one, by doubling their variances. The contamination contributions were considered with a variation of and the durations were considered with a variation of .

The fourth parameterisation—related to computational Experiment IV—turns the variation of the remaining parameters involved in the reference model also as stochastic. So, it improves the previous experiment by defining the transition probabilities () from no symptoms state also as Gaussian distributions truncated by three standard deviations with the mean value equal do the probabilities shown in Fig. 2. The deviation considered as 3.33%, so the variance was .

Table 1 presents a briefing of each parameterisation, also with the range of possible values regarding the distribution and variance of each experiment.

Table 1

Parameterisations employed in the experiments

Parameter		Parameterisations
		Experiment I		Experiment II		Experiment III		Experiment IV
Name	Symbol	Range	PD	Range	PD	Range	PD	Range	PD
T	S0f(P_c)	∞	U	∞	G	∞	G	∞	G
	I15	[4,6]		[4,6]		[3,7]		[3,7]
	I25	[4,6]		[4,6]		[3,7]		[3,7]
	I33	[2,4]		[2,4]		[2,6]		[2,6]
	I43	[2,4]		[2,4]		[2,4]		[2,4]
	I710	[8,12]		[8,12]		[7,17]		[7,17]
	I512	[10,14]		[10,14]		[8,20]		[8,20]
	I613	[10,16]		[10,16]		[6,14]		[6,14]
	I811	[9,13]		[9,13]		[9,21]		[9,21]
	I94	[3,5]		[3,5]		[3,7]		[3,7]
	I105	[4,6]		[4,6]		[2,6]		[2,6]
	R11∞	∞		∞		∞		∞
	R12∞	∞		∞		∞		∞
P_c	S0f(P_c)	0%	n	0%	G	0%	G	0%	G
	I15	100%		[80%,100%]		[60%,100%]		[60%,100%]
	I25	100%		[80%,100%]		[60%,100%]		[60%,100%]
	I33	100%		[80%,100%]		[60%,100%]		[60%,100%]
	I43	100%		[80%,100%]		[60%,100%]		[60%,100%]
	I710	100%		[80%,100%]		[60%,100%]		[60%,100%]
	I512	10%		[0%,20%]		[0%,40%]		[0%,40%]
	I613	10%		[0%,20%]		[0%,40%]		[0%,40%]
	I811	100%		[80%,100%]		[60%,100%]		[60%,100%]
	I94	100%		[80%,100%]		[60%,100%]		[60%,100%]
	I105	100%		[80%,100%]		[60%,100%]		[60%,100%]
	R11∞	0%		0%		0%		0%
	R12∞	0%		0%		0%		0%
P_t	I15→I25	30%	n	30%	n	30%	n	[20%,40%]	G
	I15→I33	55%		55%		55%		[45%,65%]
	I15→I43	10%		10%		10%		[0%,20%]
	I15→I710	5%		5%		5%		[0%,15%]

Key to the table: T: Duration of the states (1 time step = 1 day); : Neighbourhodd contamination probability; : State transition probability; PD: Type of probability distribution; U: Uniform probability distribution; G: Gaussian probability distribution; Range: The interval from the values are chosen according to the distribution adopted, with notation [min, max]; n: Fixed value

Results

The proposed procedure was executed and monitored, comparing the average value of each variable between the ten simulations for each immunity ratio. Considering the range of immunity proposed and the simulations performed for each immunity, as well the four experiments, a total number of 800 simulations was obtained from the model. This implementation was written in Wolfram language, native of the software Mathematica.

Experiment I

Figure 3 displays the impact of the immunity ratio on the population main metrics (remaining susceptibles, infection peak, deceaseds and total time steps) considering Experiment I (Table 1). First, this impact can be seen in the number of susceptibles at the final convergence (in other words, the ones who did not get ill). In a first and brief look, it is predictable that more immune individuals means less susceptibe individuals, which describes a negative linear correlation. This effect becomes perceptible after around 60% immunity. Before 50%, the final number of susceptibles is near zero, and starts to increase as immunity rises.

Fig. 3

Remaining susceptibles, deaths and infection peak as a function of the immunity ratio considering Experiment I

Second, the infection peak is also affected by immunity, describing the peak value of active cases. There is a negative correlation between immunity ratio and the infection peak which comes near to zero when immunity achieves 65%. Last, the death rate behaves similarly to the infection peak, coming down close to zero at 65% of immunity.

Figure 4 gives an idea about the emulated pandemic duration according to immunity ratio. Apparently, ratios under 55% lead to a positive correlation, or a tendency of the pandemic to take more time to end as the immunity rises. In other words, the remaining number of susceptibles do not receive protection against contamination and the the disease takes longer to spread. The maximum average time steps taken to convergence are near 640 steps at 55%, starting to decrease for bigger immunities.

Fig. 4

Variation in the number of time steps necessary to final convergence as a function of the immunity ratio considering Experiment I

Experiment II

Figure 5 depicts the impact of immunity ratio on the population main metrics considering Experiment II (Table 1). The effect previously mentioned related to the decrease of susceptibles with the increase of immunes becomes perceptible after precisely 60% immunity. Before 50%, the final number of susceptibles is near zero, starting to increase as immunity rises.

Fig. 5

Remaining susceptibles, deaths and infection peak as a function of the immunity ratio considering Experiment II

The infection peak presents the same negative correlation between immunity ratio and the infection peak which comes near to zero when immunity achieves 60%. The death rate behaves similarly to the infection peak, approaching zero at 60% immunity rate.

Figure 6 presents the amount of time steps taken to convergence, with ratios under 55% presenting a positive correlation. The maximum average time steps taken to convergence are near 720 steps at 55%, starting to decrease for bigger immunities.

Fig. 6

Variation in the number of time steps necessary to final convergence, as a function of the immunity ratio, for Experiment II

Experiment III

Figure 7 presents the impact of immunity ratio on the population main metrics considering Experiment III (Table 1). The effect previously mentioned related to the decrease of susceptibles with the increase of immunes becomes perceptible after precisely 60% immunity rate. Before 45%, the final number of susceptibles is near zero, starting to increase as immunity rises, but rising faster than in the previous experiments.

Fig. 7

Remaining susceptibles, deaths and infection peak, as a function of the immunity ratio, for Experiment III

The infection peak shows the same negative correlation between immunity ratio and the infection peak, which comes near to zero when immunity achieves 60%. The death rate behaves similarly to the infection peak, reaching near zero at 60% immunity.

Figure 8 presents the amount of time steps taken to convergence, with ratios under 55% presenting a positive correlation. The maximum average time steps taken to convergence are near 820 steps at 55%, starting to decrease for bigger immunities.

Fig. 8

Variation in the number of time steps necessary to final convergence as a function of the immunity ratio, for Experiment III

Experiment IV

Figure 9 presents the impact of immunity ratio on the population main metrics considering Experiment IV (Table 1). The effect previously mentioned related to the decrease of susceptibles with the increase of immunes becomes perceptible after precisely 60% immunity rate. Before 50%, the final number of susceptibles is near zero, starting to increase as immunity rises.

Fig. 9

Remaining susceptibles, deaths and infection peak, as a function of the immunity ratio, for Experiment IV

The infection peak shows the same negative correlation between immunity ratio and the infection peak, which approaches zero when immunity achieves 60%. The death rate behaves similarly to the infection peak, getting close to zero at 60% immunity.

Figure 10 displays the number of time steps taken to convergence, with ratios under 50% showing positive correlation. The maximum average time steps taken to convergence are near 720 steps at 50%, starting to decrease for bigger immunities.

Fig. 10

Variation in the number of time steps necessary to final convergence as a function of the immunity ratio, for Experiment IV

Conclusions

A stochastic cellular automaton based framework was built to perform estimates of the collective immunity of COVID-19 in a given population of static individuals interacting socially only with their neighbours in a Moore neighbourhood of unitary radius. The main idea was to evaluate the impact of immune individuals, from zero to 95% of the population, in the dynamics of propagation of SARs-CoV-2, by observing the remaining susceptibles, the infection peak, the number of deaths and the time taken for the disease spread to end, thus estimating the collective immunity of COVID-19.

Additionally, the proposed model was implemented considering some input parameters regarding the duration of the infected states of COVID-19 and the contamination contributions for these states. Four different parameterisations were applyed to the model, which allowed to compare the collective immunities:So, considering the four experiments performed, the average collective immunity is 55% with a variance of . The estimation of the collective immunity seems to be suitable to this framework, since it is little affected by the input parameters.

Experiment I led collective immunity to lie between 50% and 65% immunity rate. The peak of remaining susceptibles was found at 60% but the major shift in their increase began at 50%. The infection peak and the death rate approached zero at 65%. The big picture, in this case, shows collective immunity of 57.5% on average;

Experiment II entailed collective immunity in the slightly smaller range between 50% and 60% immunity rate. The peak of remaining susceptibles was found at 60%, with the major shift in their increase also beginning at 50%. The infection peak and the death rate came down close to zero at 60%. Overall, this case also points collective immunity of 55% on average;

Experiment III induced collective immunity to stay between 45% and 60% immunity rate. The peak of remaining susceptibles was found at 60% but the major shift in their increase started earlier than in the previous experiments, at 45%. The infection peak and the death rate went close to zero at 60%. Once again, the present case also makes it evident the onset of collective immunity of 52.5% on average.

Experiment IV induced collective immunity to stay between 50% and 60% immunity rate. The peak of remaining susceptibles was found at 60%, with the major shift in their increase also beginning at 50%. The infection peak and the death rate went close to zero at 60%. So, the present case also makes it evident the onset of collective immunity at 55% on average.

A difference between the parameterisations adopted was related to the maximum number of time steps taken to final convergence. In terms of collective immunity, all four led to a peak duration at 55% with a big decrease at 60%, but as randomness increased in the parameterisations, the emulated pandemic took more time to end, except for the last experiment, which presented a minor duration than the previous one. This could be explained by the combinations between lower contamination contributions and lower state durations, which would make it difficult for the virus to spread, while maintaining the main characteristics of the pandemic (infection peaks, death rates and maximum remaining susceptibles), but for the last experiment an additional explanation could be that the variations in the transition probabilities from the no symptoms state somehow compensate the variations in durations and contamination contributions. So, the estimation of the pandemic duration is more affected by the input parameters, thus rendering its estimation an unsuitable application to this framework.

The calculated collective immunity at 55% rate is close to the estimated value of 60% or little more found in the literature (Aguas et al. 2020; Jones 2020), also very close to the results in Lima and Balbi (2021), whose value was 60%. But it is important to remember that the collective immunity of COVID-19 is becoming more uncertain as new variants of SARS-CoV-2 appear (García-García et al. 2022).

So, the stochastic cellular automaton architecture defined showed a fair degree of robustness regarding the parameter values we tested. In order to further probe our approach, it is appealing to perform additional experiments considering larger variations for the contamination contributions and durations, as well some degree of asymmetry in the Gaussian distributions to improve the model sensibility in terms of time steps taken to final convergence. On the other hand, the model could be improved to cope with new variants of SARS-CoV-2 infection capability by modifying the cells neighbourhood size and distribution.

Another aspect that can be improved in the model is changing the immune () state from a sink state to a time-dependent state, considering that COVID-19 immunity has been proven not to be permanent, so that re-infections are possible, approximating the approach to a SIRS model (Fuentes and Kuperman 1999). This could be addressed by creating new state transitions back to the susceptible () state.

So, considering all above, it is fairer and more realistic to consider the results we presented as more related to the original, pre-delta variants of the SARS-CoV-2. However, as we argued, the model itself can be easily modified so as to account also for the more recent and important variants, such as omicron (Callaway 2022; Willyard 2022).

The chosen number of simulations per immunity ratio, 10, also seemed few to average the sthocastic simulations, but this choice is directly related to the computational effort to perform all the procedure (10 simulations for each immunity ratio), which takes almost four hours to be completed. More simulations to average the outputs will certainly require more computational capability or optmization of the procedure, such like a larger immunity interval, for example 10% instead of 5%.

Finally, the cellular atomaton has regularities which decreases its reliability in emulating the spreading of COVID-19, such like the size and distribution of the neighbourhood and the immobility of the cells. This could be adapted by performing random connections between cells, instead of the rigid Moore neighbourhood with unitary radius, like small-world type networks (Schimit and Monteiro 2009); and representing the mobility of the individuals in the population, thus increasing the parameters of the model (Sirakoulis et al. 2000). These are the directions we are taking.