Abstract
Mathematical models play an important role in epidemiology. The inclusion of a spatial component in epidemiological models is especially important to understand and address many relevant ecological and public health questions, e.g., when wanting to differentiate transmission patterns across geographical regions or when considering spatially heterogeneous intervention measures. However, the introduction of spatial effects can have significant consequences on the observed model dynamics and hence must be carefully analyzed and interpreted. Cellular automata epidemiological models typically rely on simplified computational grids but can provide valuable insight into the spatial dynamics of transmission within a population by suitably accounting for the connections between individuals in the considered community. In this paper, we describe a stochastic cellular automata disease model based on an extension of the traditional Susceptible-Infected-Recovered (SIR) compartmentalization of the population, namely, the Susceptible-Hospitalized-Asymptomatic-Recovered (SHAR) formulation, in which infected individuals either present a severe form of the disease, thus requiring hospitalization, or belong to the so-called mild/asymptomatic class. The critical transmission threshold is derived analytically in the nonspatial SHAR formulation, and this generalizes previously obtained theoretical results for the SIR model. We present simulation results discussing the effect of key model parameters and of spatial correlations on model outputs and propose an algorithm for tracking the evolution of infection clusters within the considered population. Focusing on the role of import and criticality on the overall dynamics, we conclude that the current spatial setting increases the critical transmission threshold in comparison to the nonspatial model.
1. Introduction
Infectious diseases have shaped the global population throughout history, and (whether through devastating epidemics or recurrent outbreaks of endemic diseases) they have been—and often still are—responsible for a large number of deaths worldwide. The ongoing COVID-19 pandemic has certainly brought significant attention to the important contribution that mathematical modeling can have in monitoring the spread of infectious pathogens across a population, in better understanding its underlying mechanisms, and in exploring the effects of possible control measures aimed at containing the transmission of the considered disease (and possibly eradicating it from the observed population).
While the foundations of modern epidemiology based on compartmental models were laid in the first few decades of the 20th century, thanks to the work of public health physicians such as Ross, Hamer, McKendrick, and Kermack [1], the last 100 years have seen incredible advances in the field that were made possible also by the advent of computers and modern-day computational resources. The rapid increase in computer power of the last few decades has in fact allowed us to tackle many problems that were previously considered intractable and to simulate numerically complex systems across the entire spectrum of scientific disciplines, testing various working hypotheses and making predictions in many different hypothetical scenarios. This is true also in epidemiological modeling, where modern computers and high-performance simulations allow the rapid analysis of the detailed models often required in order to provide quantitatively accurate predictions for the design of suitable public health policies.
Probably the most famous compartmental model in epidemiology is the SIR model [2]. In this formulation, individuals in the population belong to one of three mutually exclusive categories: susceptibles (individuals that have never been exposed to the pathogen causing the disease), infected individuals (the ones who have been exposed to the pathogen, becoming ill, and are capable of transmitting it), and recovered (the ones who have gone through the infection, recovering with a certain degree of immunity)—hence the SIR name for the model. The temporal evolution of the number of individuals in each disease-related category is described mathematically via a set of three ordinary differential equations involving three rates describing the probability of susceptible individuals becoming infected, of infected individuals to fight off the disease and recover, and of recovered individuals to become susceptible once again, that is, the so-called infection, recovery, and waning immunity rates, respectively.
This model and its many variations (SI, SIS, SEIR, etc.) give the basic framework for the vast majority of epidemiological models found in the literature, and as such, they have been analyzed in detail and extensively used to describe disease transmission dynamics in many contexts (e.g., [3–5] just to name a few). These classical models can be further refined in order to address disease-specific questions or according to the availability of data, for example, by including multiple strains of the same type of virus (see, e.g., [6–11]) and by subdividing the population into groups according to age or the risk of contracting the disease (e.g., [12–15] and references therein). Adding stochasticity is also an important aspect that might be key in correctly describing observed transmission patterns, especially when population size is small or in general when the number of infectious individuals is low [16, 17].
While explicitly including a spatial component in an epidemiological model could be an unnecessary complication to answer certain ecological and public health questions, some degree of spatial resolution is essential when wanting to differentiate transmission patterns across geographical regions [18] or when considering spatially heterogeneous intervention measures [19].
In this work, we consider the SHAR extension of the SIR compartmental model and develop a cellular automata stochastic formulation to explicitly account for spatial effects in the considered framework. We discuss the role of import and study the behavior of the system in various parametric settings, determining the critical transmission threshold analytically in the nonspatial formulation and exploring the effect of spatial correlations on the more general spatiotemporal model. Finally, we introduce the concept of infection clusters and investigate their evolution (in terms of the amount of clusters present in the system as well as their individual size) as a function of different parameter combinations.
2. The SHAR Model
The SHAR model is an extension of the well-known SIR framework in which the I class is divided into two groups labeled and , respectively: stands for individuals developing a severe form of the disease and likely being hospitalized, while refers to infected individuals who are asymptomatic or have a mild form of the disease. The model dynamics is described by the following reaction scheme:
In addition to the infection rate , the recovery rate , and the waning immunity rate (already present in the general SIR formulation), the SHAR approach includes two extra epidemiological parameters: the severity ratio that gives the fraction of infected individuals who develop severe symptoms and the scaling factor , which differentiates the infectivity of mild/asymptomatic infections with respect to the baseline infectivity of severe/hospitalized cases [20]. The value of can be tuned to describe different situations: a value of indicates that severe cases have larger infectivity than mild cases (which in the case of an infectious respiratory disease, for example, could be linked to enhanced coughing and sneezing), while would be used to describe the scenario in which asymptomatic individuals and mild cases contribute more to the spread of the infection (e.g., due to their higher mobility and possibility of interaction) than severe cases (that are more likely to be detected and isolated). Here, we also consider the import factor which refers to the possibility of susceptible individuals contracting the disease via an undetected infection chain which started outside the studied population [21, 22].
The stochastic SHAR epidemic with import is modeled as a time-continuous Markov process to capture population noise. The temporal dynamics for the probability of having at time an integer number of susceptible, hospitalized, and asymptomatic can be described by the following master equation [23]:
where the number of recovered individuals follows from the constant population size assumption .
From the definition of the mean values we can derive dynamic equations for the averages. For example, for the mean susceptible population, we derive
and then insert the master equation (2) on the right-hand side. By using the so-called mean-field approximation [24] (essentially assuming that the covariance between variables is negligible and hence approximating the second order moments as follows: and ), and by suitably incorporating boundary modifications of Equation (2) for the cases in which , we obtain the mean-field ODE system for the SHAR model
where we have not included the to simplify notation.
2.1. The SHAR Model with Import around the Critical Threshold
Let (where and could either be the variables of the deterministic SHAR model or the mean-field approximations of the corresponding stochastic system) and assuming abundance of susceptibles (i.e., ), then we can rewrite in vector form the equations for the dynamics of the infected classes of Equation (5) as follows:
where the matrix and the vector are defined, respectively, as
The steady-state of Equation (6) is nothing but , where the matrix can be explicitly computed via straightforward calculations as leading to
The eigenvalues of the matrix are
so that community transmission is controlled (and hence, the steady-state is an asymptotically stable equilibrium) for values of below the critical threshold
The nonlinear dependence of on and for a fixed value of is illustrated in Figure 1, highlighting that different parametric combinations could lead to the same value of , but the main underlying modeling assumptions in terms of the scaling parameter would be preserved (in fact, for , there are no combinations with giving the same value of as any combination of with ).
(a)
(b)
In these settings, letting , we observe that
thus generalizing the critical threshold result obtained for the SIR model with import [25]. Notice that, rather than considering the asymptotic number of infected , in order to generalize previous results on the SIR model, we had to suitably rescale by the parameter. However, can also be written in compact form as
with , and hence, its qualitative dependence on the import term remains unchanged.
Figure 2(a) shows 100 realizations (obtained by using the Gillespie algorithm [26]) of the stochastic SHAR model in the subcritical regime (). Figure 2(b)) shows the limiting case of community spreading at criticality , where we can see a linear growth of infections, while in Figure 2(c), an exponential increase of cases is observed for supercritical community spreading with .
(a)
(b)
(c)
This transition of the dynamics from subcritical to supercritical regimes is also illustrated in Figure 3, where the mean-field solution is plotted for different values of both in natural and log-log scales. Note that throughout the manuscript, we use days, labelled by , as time units for the proposed simulations. However, could be changed as to represent any relevant time unit of interest (which would be typically dictated by the time scale on which real-life data for model validation would be available).
(a)
(b)
3. The Spatial SHAR Model
In this section, we present the spatial version of the stochastic SHAR model. Let us consider a lattice consisting of sites, each of which is occupied by an individual and labelled by the indicators , according to its state. Each site, , can either have a susceptible (, ), hospitalized (, ), asymptomatic (, ), or recovered (). Notice that there is no need to explicitly introduce the indicator since in the present settings, . The state vector for the system is thus given by . Let be the lattice adjacency matrix with entries , if sites and are connected, and otherwise. From the matrix and any given state configuration, the following quantities of interest can be straightforwardly computed: the number of neighbors to site , , the number of infected neighbors to site , , and the force of infection to site , .
The adjacency matrix gives the “neighboring” structure of the lattice. Although different network structures can be considered [27], it is known that the particular choice of lattice structure does not make much difference on the “universality” property of the system [28]. In general, in two-dimensional lattices, interactions are considered with the nearest 4 (Von Neumann neighborhood) or 8 (Moore neighborhood) neighboring lattice sites, as illustrated in Figure 4 for a lattice. In our present study, we will consider Von Neumann neighborhoods and apply periodic boundary conditions so that for all .
(a)
(b)
(c)
(d)
The reaction scheme for the spatial SHAR model is
where the processes involving contagion between individuals within the considered population may only occur when sites and are connected. The master equation for the stochastic spatial SHAR model with import is given by
3.1. Critical Threshold
In Section 2, the critical threshold for the stochastic SHAR model was obtained analytically from the mean-field approximation, see Equation (11). In the spatial case, the problem is much more challenging and since we do not have an analytical expression for this threshold, we must rely on numerical methods, which are often computationally expensive.
To get a first insight into the spatial dynamics, let us consider the one-dimensional SHAR process in which each lattice site is connected to its nearest neighbors. Using periodic boundary conditions, we can essentially think about a ring lattice (a network in which each node is connected to its 2 nearest neighbors, i.e., for all ). Figure 5 shows three different realizations of the SHAR process on a lattice of size . Each pixel row of the given plots shows the spatial configuration of the considered one-dimensional system at a given time point, while time increases from top to bottom. In all cases, we start from an entirely susceptible population (green) with one hospitalized individual in the center of the lattice and run the simulation up to the final time . Hospitalized individuals are shown in red, and asymptomatic cases are yellow, while recovered are blue. The set of parameter values for these simulations is , , , , and , while in Figure 5(a), in Figure 5(b), and in Figure 5(c). We can see that for the smallest value of considered here, only a few infections occur (only asymptomatic for the particular realization shown in Figure 5(a)) before extinction, but for increasing values of , larger outbreaks take place and import also starts to play a role (see the small outbreak generated by import towards the end of the simulation shown in Figure 5(c)).
(a)
(b)
(c)
For the two-dimensional case, let us consider a population of individuals distributed in a square lattice, once again adopting a Von Neumann neighborhood with periodic boundary conditions. The population is assumed to be entirely susceptible at the beginning, except for one initial infected individual located in the center of the lattice. Parameter values are , , , , and , as in the simulations for the nonspatial model of Section 2. Figures 6 and 7 show four snapshots in time of a single realization of the spatial SHAR system for two different values of the infection rate, namely, and , respectively. In both cases, a first outbreak is generated by the infected individual initially considered in the lattice center. Then, the import term produces more isolated outbreaks of different sizes. While some of these outbreaks involve only a few individuals and quickly die out, with all individuals recovering, others may involve many more active infections (severe and/or mild and asymptomatic), spreading widely and eventually collapsing with neighboring clusters, leading to much larger outbreaks. An exponential growth of cases is to be expected when the community transmission is supercritical (see discussion below).
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
In order to decide whether there is critical percolation threshold, numerous stochastic realizations are run for different infection rates . Figure 8 shows the mean infected population of 100 stochastic realizations for 10 different values of in log-log scale. As for the nonspatial case, we observe that the mean number of infected individuals for low values of eventually levels off, approaching a constant infected population size, while for the largest values of , we see a very fast increase in the mean number of infected (an exponential growth is expected for values of above critical threshold). However, we notice that for (i.e., when coincides with the of the nonspatial model in the considered parametric setting), the mean number of infected no longer grows linearly in time but rather levels off, indicating a subcritical behavior for this particular parameter value. Hence, we can infer that spatial correlations for the considered model contribute to maintaining disease transmission under control for larger values of the infectivity . Nevertheless, giving a more precise estimate of the critical threshold for this spatial model would require considering much larger lattice sizes (in order for the “abundance of susceptible” assumption to remain valid for longer times) and longer simulations (with the corresponding heavier computational burden), thus allowing to disambiguate between values of corresponding to a super- or subcritical regime. Obtaining accurate estimates of the critical exponents for the proposed spatial model is out of the scope of the presented study but will surely be a topic of our future investigation.
3.2. The Effect of Import and Cluster Formation
In this subsection, we focus on two key concepts of epidemiological dynamics, namely, import and cluster formation. For the sake of simplicity, let us assume that the considered infectious disease is spread via human-to-human contact (i.e., ignoring vector-borne diseases and disease transmission across species). When no infections are present within a considered population, imports are responsible for triggering disease outbreaks. Imports typically take the form of external visitors who remain for a short period of time within the considered population (so that their presence does not change the population size) either introducing or reintroducing the pathogen, or of individuals belonging to the studied population who travel outside and return infected. A disease outbreak may become an epidemic but may also frequently die out due to the stochastic nature of the process. However, import of the pathogen from outside the population can prevent permanent extinction and generate new and sometimes unexpectedly large outbreaks.
In the case of infectious diseases, the development of useful testing and tracing strategies is fundamental. If we know that a particular individual has been infected and know her connections (in this context, her closest neighbors), then by looking at their state, we can effectively monitor the evolution of the outbreak and plan suitable containment measures accordingly. Fast and reliable testing strategies are clearly essential in obtaining a clear picture of the disease dynamics in the population of interest, and tracking individuals who enter from outside is particularly important, due to their previously mentioned key role in potentially triggering new outbreaks.
3.2.1. An Infection Cluster Computing Algorithm
Here, we describe an algorithm to monitor the evolution of an epidemic by tracking the number of infected individuals generated by each isolated outbreak. We assume that the population under study is often and regularly tested as to quickly identify new infections (either due to import or direct contact with an infected neighbor). For this purpose, we define the concept of an infection cluster, that is, a group of connected individuals who have been infected by the same index case and are either still infected or have recovered but have not lost their immunity yet.
Our proposed algorithm is coupled to the stochastic SHAR model solver described above and stores the cluster information in a cluster configuration matrix, . Rows of the matrix correspond to different lattice sites, and a new column is added to for every cluster configuration change, i.e., whenever the number of clusters present in the system or any of the cluster sizes changes. is initialized as a column vector of zeros having the size of the spatial lattice used for the simulation, and the label “1” is assigned to the row corresponding to the only infected individual present in the center of the lattice. A new isolated cluster is created when there are no infected neighbors to the newly infected site. Whenever this happens, our algorithm assigns to the new cluster a new label (an integer corresponding to the total number of clusters present in the system at that particular time point) and all future infections in the same cluster are then assigned the same label. When two clusters collide, we choose for the new infection the label of the larger cluster involved in the collision, but other choices (e.g., random selection) are possible. Moreover, in all the presented simulations, we assume so that the number of clusters and their individual size never decreases. Clearly, as soon as , this would be no longer the case and care should be taken in ensuring the proper labeling of new clusters. Notice that in order to monitor the number of clusters and their sizes, it is not necessary to store the state of the entire lattice for each time step. Instead, a new column of is only created when a new infection occurs. The pseudocode is given in Algorithm 1.
| ||||||||||||||||||||||||||||||||||||||
Figure 9 shows three different aspects of the monitored cluster information for two realizations of the spatial stochastic SHAR model with import obtained, respectively, with (left column) and (right column). Figures 9(a) and 9(d) give the cluster configurations at the end of the simulation, namely, at , for both realizations, and in both plots, a different colour is assigned to each cluster in order to aid their visualization. Figures 9(b) and 9(e) track the temporal evolution of the number of clusters in each simulation, while Figures 9(c) and 9(f) display changes in each individual cluster size over time.
(a)
(b)
(c)
(d)
(e)
(f)
By comparing these two individual realizations, we see that when is larger, many more clusters form in the time interval considered and their individual size is generally much bigger than for the other value of the considered parameter. In order to confirm this trend, we perform a simple statistical analysis over several realizations of the spatial SHAR model. In particular, we monitor infection cluster information (i.e., cluster size and number of clusters) of 500 realizations for both values of , plotting in Figures 10 and 11 the corresponding histograms. The obtained results confirm our expectations. Specifically, Figure 10 shows that clusters of much larger size form for the larger value of (which is to be expected due to the larger infectivity and hence faster transmission of the pathogen), while in Figure 11, we observe that, in average, a slightly larger number of clusters (i.e., new isolated outbreaks) forms for the larger value of due to its inflating effect on the net import rate of the system.
(a)
(b)
(a)
(b)
4. Conclusions
In this work, we have presented the stochastic SHAR model with import and its mean-field approximation. We have computed the critical threshold of the system under the assumption of abundance of susceptibles and the asymptotic number of infected for controlled community transmission as a function of the model parameters, thus generalizing previously obtained results for the SIR model with import. Then, we introduced a cellular automata extension of the SHAR framework accounting for spatial variation within the system and we presented the master equation for the proposed model. Modifications to the critical threshold behavior in the presence of spatial correlations have been discussed in Section 3.1, showing that, in the parametric setting considered, these correlations pushed the system critical threshold towards larger values, hence contributing to maintaining community transmission controlled for a larger range of values. Finally, in Section 3.2, we have focused on the role of import in triggering and maintaining an epidemic and proposed a way to track the number of infected individuals generated by each index case by looking at cluster formation and their evolution in time. We see the proposed algorithm as a way to mimic the contact tracing procedure that might be implemented by the health system of a given region in order to monitor disease transmission. To conclude the manuscript, we would like to mention some of the limitations of the presented work and comment on possible directions of future work. The proposed analysis is aimed at providing an overview of the potential offered by the SHAR framework in a spatially extended context, focusing in particular on tools such as the infection cluster computing algorithm to study relevant features of the resulting spatiotemporal dynamics. We explored the effect of varying the infectivity parameter on the model output but kept all other parameters fixed at some given (and typically arbitrary) value. A much more in-depth exploration of the parameter space remains to be performed. Clearly, according to the considerations made in Section 2.1, we expect variations in to affect the critical threshold of the system, and increasing values of to produce more and more isolated outbreaks, following our considerations of Section 3.2. However, studying the combined effect of multiple parametric changes (including different population sizes) will surely be a natural extension of this work. We remark also that throughout this manuscript, we always assumed and hence the case (which will likely produce much richer dynamical scenarios) should also be investigated. Moreover, this particular work did not focus on describing the dynamics of any specific disease and therefore specific parameter values (or ranges) were not validated against any real-life data. Nevertheless, the SHAR framework already captures many important features of infectious disease transmission, and as such, we expect it to be a very versatile tool for further disease-specific investigation (e.g., in the context of COVID-19, seasonal influenza, etc.).
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare no potential conflict of interests.
Authors’ Contributions
M.A. and N.S. devised the mathematical framework and supervised the development of this study. D.K. and N.C. designed the computational framework and conceived and performed the numerical simulations. All authors contributed to the analysis of the results and to the writing of the manuscript.
Acknowledgments
This research is supported by the Basque Government through the “Mathematical Modeling Applied to Health” Project, BERC 2018-2021 program and by the Spanish Ministry of Sciences, Innovation and Universities: BCAM Severo Ochoa Accreditation SEV-2017-0718. N.C. acknowledges also the ELKARTEK KK-2020/00008 funding scheme.