Multi-Strain Epidemic Models of Mutating Viruses with Airborne Transmission Based on Cellular Automata and Ordinary Differential Equations

Abstract

The emergence of SARS-CoV-2 virus mutations is a pressing concern in the context of the COVID-19 pandemic. In this paper, a Susceptible–Infected–Recovered (SIR) compartmental model has been formulated in terms of cellular automata and ordinary differential equations to simulate the transmission dynamics of a virus. This model accounts for the potential of new virus mutations to create additional infections of a multi-strain disease while also considering the impact of vaccination on disease control within a population. The basic reproduction number of the disease is derived, and the effect of virus mutations and vaccination rates is evaluated on a population over five years. The results of numerical simulations demonstrate the significant role of maintaining high vaccination rates in controlling the spread of the virus, even when assuming that all variants have similar illness characteristics and that a single shot of vaccine provides complete and lifelong protection against all strains. The findings underscore the necessity for countries to implement a consistent and high-level vaccination plan as soon as vaccines become available in order to mitigate the impact of a pandemic effectively. In conclusion, this study highlights the importance of taking into account the potential impact of virus mutations for controlling the COVID-19 pandemic. Furthermore, it emphasizes the critical role of vaccination in limiting the spread of the virus, and emphasizes the need to implement and maintain high vaccination rates as part of a comprehensive approach to managing the ongoing pandemic.

1. Introduction

The COVID-19 pandemic is an ongoing global public health crisis, which has been further complicated by the emergence of several mutations in the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). These variants have distinct genetic profiles and exhibit varying clinical and epidemiological features. In March 2020, the D614G mutation, characterized by increased infectivity, became the predominant strain worldwide, replacing the strain reported in Wuhan [1]. In September 2020, a new variant, known as Cluster-5 or $\mathsf{\Delta }$FVI-spike, emerged from mink farms in Denmark, resulting in the culling of approximately 17 million animals [2].

In December 2020, a variant of concern, SARS-CoV-2 VOC 202012/01, was identified in the United Kingdom. Preliminary reports suggested that this variant had increased transmissibility and mortality risk [3,4]. It has since been reported in over 62 countries [4]. In the same month, a new strain, 501Y.V2, was identified in South Africa, replacing the original strain in some provinces. Studies have demonstrated that this variant has higher viral loads, increased transmissibility, and reduced vaccine efficacy [3,5,6].

Another variant, named P.1, emerged in Brazil in November 2020, causing a second wave of infections in the country [7]. Recently, the Omicron variant was identified, and it has quickly become the predominant strain globally [8]. Notably, several new mutations have been detected in the spike protein of the virus, which is the primary target of current COVID-19 vaccines. The Beta, Gamma, and Delta variants have shown increased transmissibility and may be associated with higher rates of hospitalization and mortality [9].

The emergence and spread of these SARS-CoV-2 variants underscore the urgent need for continued surveillance and research to understand their clinical and epidemiological impact. Effective management of this pandemic requires a coordinated global response, including widespread vaccination and public health measures, such as physical distancing, wearing masks, and avoiding large gatherings. Furthermore, ongoing monitoring and genomic sequencing of the virus will be crucial in identifying new variants and developing strategies to combat their spread. Furthermore, without serious research on modelling such crises, control strategies aimed at combating mutating viruses or variations would inevitably fall short [10].

Hence, it is crucial for epidemiological modelling to account for novel variants of a disease in order to comprehend the dynamics of disease spread and furnish comprehensive solutions to public health officials. The necessity of a multidisciplinary approach for studying and managing virus variants usually begins with the use of mathematical modeling and differential equations, aiming to establish a framework for analyzing the interplay between viral mutations, transmission dynamics, and their impact on population-level infection patterns [11]. Within the ambit of this paper, we propose a disease transmission model that involves viral infection via human contact. This model accommodates the spontaneous emergence of new virus variants during the infection process, while also accounting for a spectrum of vaccination rates as a measure of disease control.

Previous works presented multi-strain disease models for a variety of objectives. In general, there is a need to assess the transmissibility of a new mutation, and whether this new mutation will be the majority of cases, invading a population and becoming endemic. Diseases transmitted both directly and environmentally [12], considering an age-structured model [13], formulated for mean-field and pair approximation approaches [14], and considering Malaria [15] and COVID-19 [16] are a few examples of multi-strain disease studies. Consider two strains of the same pathogen. If one strain in a host population may limit the growth or presence of another strain, we have competitive exclusion, which has also been taken into account due to important implications for disease dynamics and control strategies [15,17]. For example, if a highly infectious strain of a pathogen is introduced into a population that is already infected with a less infectious strain of the same pathogen, the highly infectious strain may outcompete the less infectious strain and become dominant. This could result in a more severe disease outbreak and potentially make control efforts more difficult.

In spite of the concerted efforts made by these models to generalize the outcomes across multiple strains prevalent in a population, the complexity of mathematical analysis directly corresponds to the number of strains involved. In the study conducted by Breban et al. [12], the stability of equilibria and the determination of the basic reproduction number were derived for multiple strains; however, the discussion was limited to the consideration of only two strains. Several other scholarly works have investigated scenarios involving two [17,18,19] and three [20] strains. Furthermore, certain studies have focused on evaluating the process by which the most successful strain would dominate the entire human population, thereby incorporating an evolutionary aspect into the investigation [13,15]. Consequently, there is a demand for a flexible model that accommodates the presence of any necessary number of strains, while simultaneously enabling a comprehensive mathematical analysis of the basic reproduction number.

The present model is constructed utilizing Probabilistic Cellular Automata (PCA) and Ordinary Differential Equations (ODE) within the framework of the Susceptible–Infected–Recovered (SIR) compartmental model. The use of cellular automata has been documented in the literature as a viable tool for simulating COVID-19 [21,22,23], as well as variations of compartmental models, particularly the Susceptible–Exposed–Infected–Recovered (SEIR) model, which accounts for an exposed state before infection [24,25,26]. Indeed, compartmental models have been tailored to meet specific research needs, such as the investigation of the impact of the COVID-19 pandemic on healthcare systems utilizing states related to hospitalized individuals [23,26]; the assessment of the influence of migrant-exposed individuals in predicting COVID-19 trends [27]; the use of separate asymptomatic and symptomatic infected states to study the effects of social isolation [22,23,28]; the propagation of two SARS-CoV-2 strains, separating the respective exposed and infected states [18]; and even the study of an elementary cellular automaton evolution for a visible representation of important features of the protein of SARS-CoV-2 [29]. Furthermore, PCA and ODE have been employed in tandem to compare the pairwise approach for modelling the contact between individuals (PCA) to the mean-field concept of ODE, with ODE representing a mean-field approximation for PCA [14,30,31].

In summary, cellular automata possess specific constraints that researchers must consider. These models operate under the assumption of cellular homogeneity, neglecting subtle variations and individual-level heterogeneity, thereby potentially compromising their accuracy in capturing the intricate intricacies of real-world scenarios. Additionally, they may encounter difficulties when incorporating evolving population structures and behavioral adaptations during disease outbreaks. Nonetheless, despite these limitations, which are also prevalent in alternative mathematical and computational methodologies, cellular automata continue to hold significance in the study of disease propagation dynamics, serving as valuable complements to other epidemiological modelling approaches.

A broad range of mutation and vaccination rates are taken into account, yielding numerous scenarios. The compartmental model consists of three states: susceptible, infected, and recovered (for the cases with zero viral mutations). During the process of infection, the virus has a chance of $\mu$ to undergo a mutation, thus producing a new strain and a distinct type of infected state. With the PCA and ODE models described, the basic reproduction number is derived from ODE using the next-generation matrix process [32,33].

The model has been formulated based on specific assumptions, consisting of the homogeneity of disease features across all viral variants, the probability of reinfection among individuals who have already recuperated from the disease by contracting new strains, the administration of single-dose vaccines that guarantee durable immunity against all viral variants, immediate onset of immunity post-inoculation, and a consistent and abundant supply of vaccines. Therefore, this paper presents a numerical exercise that showcases the effectiveness of achieving elevated vaccination rates in mitigating disease burden and preventing the emergence of forthcoming variants that may pose a significant risk to public health.

In the next section, the model is fully described, with the results in Section 3. The discussion is presented in Section 4.

2. Methods

The probabilistic cellular automaton is structured on a square lattice with side n. The lattice consists of $n\times n=N_L$ cells that represent individuals in one of the disease states. In the context where only a single viral strain exists within the population, an individual in the susceptible state may undergo infection with a probability of $P_i$ per time step. Following infection, the individual has a probability of recovering and transitioning to a recovered state with a probability of $P_C$, as well as a probability of experiencing fatality due to the disease with a probability of $P_D$ per time step. Furthermore, a recovered individual may experience natural death with a probability of $P_N$. In the event of fatality, a susceptible individual replaces the deceased member to ensure a constant population size.

The probability of individuals becoming infected is a function of their neighborhood. Each individual is connected to C neighbors located within a Moore radius r per time step. The neighborhood alters at each time step, with C connections being randomly established using the following procedure. An individual has a likelihood of establishing a connection to a cell in layer i, with the probability $q_{i,r}=2(r+1-i)/r(r+1)$. By layer i, consider the cells that are at a Chebyshev distance of i. A random cell from this layer is selected, and C connections are established. Subsequently, the number of infected neighbors $v_I$ determines the probability of an individual being infected, given by $P_i\left(v_I\right)=1-e^{-k{v}_I}$, where k is a parameter associated with the disease’s infectivity. This approach has been employed in several previous works, such as [22,23,31]. When multiple disease strains are present in the population, the number of infected neighbors is counted separately, and the likelihood of becoming infected with each disease type is determined randomly to prevent any prioritization. It is assumed that an individual may only contract one strain of the disease at a time.

The introduction of new mutations into a model leads to the addition of new states. For example, in the case where two strains are circulating within a population, $S_X$ represents the state of individuals susceptible to all strains, while $S_{X-1}$ and $S_{X-2}$ denote the states of individuals susceptible to all strains except strain 1 and strain 2, respectively. These states represent individuals who have previously recovered from either strain 1 or strain 2. Infected individuals of strain 1 and strain 2 are denoted as $I_1$ and $I_2$, respectively. Recovered individuals are represented by the state $S_{X-1,2}$.

To illustrate the possible transitions between states, Figure 1 depicts a transition state diagram for a single strain, shown in blue with dark arrows, as well as the possibilities when three strains are present within the population, shown in light red with grey arrows. For instance, individuals who have recovered from strain 1 may become infected with either strain 2 or strain 3, and consequently, transition to either $S_{X-1,2}$ or $S_{X-1,3}$, respectively. In the event of individual mortality, replacement occurs with individuals who are susceptible to all strains and are represented by the state $S_X$. This mechanism maintains a constant population size.

When vaccination is implemented in the population, $\gamma$ susceptible individuals of any type are vaccinated per time step. Following vaccination, individuals become immune to the disease and are henceforth impervious to infection. The only possible state transition is the occurrence of natural death, whereby the decreased individual is replaced by an individual of the susceptible state $S_X$.

The model was implemented using the C programming language and is available upon request to the author.

The ODE Model

Here, we propose an exercise to construct an ODE model from the PCA model in a systematic manner. Subsequently, we derive the basic reproduction number for a few cases, and extend it to the general scenario.

Let us begin by considering the case with a single variant. In this scenario, the ODE model can be expressed as follows:

$$\begin{array}{ccc}\hfill \frac{d{S}_{\left\{X\right\}}\left(t\right)}{dt}& =& -a_1{S}_{\left\{X\right\}}\left(t\right)I_1\left(t\right)+c{I}_1\left(t\right)+e{S}_{\{X-\{1\left\}\right\}}\left(t\right)+eV\left(t\right)-p{S}_{\left\{X\right\}}\left(t\right)\hfill \\ \hfill \frac{d{I}_1\left(t\right)}{dt}& =& a_1{S}_{\left\{X\right\}}\left(t\right)I_1\left(t\right)-b{I}_1\left(t\right)-c{I}_1\left(t\right)\hfill \\ \hfill \frac{d{S}_{\{X-\{1\left\}\right\}}\left(t\right)}{dt}& =& b{I}_1\left(t\right)-e{S}_{\{X-\{1\left\}\right\}}\left(t\right)-p{S}_{\{X-\{1\left\}\right\}}\left(t\right)\hfill \\ \hfill \frac{dV\left(t\right)}{dt}& =& p{S}_{\left\{X\right\}}\left(t\right)+p{S}_{\{X-\{1\left\}\right\}}\left(t\right)-eV\left(t\right)\hfill \end{array}$$

(1)

where $a_1$ is the infection rate constant; p is the vaccination rate constant; b is the recovering rate constant; c is the death rate constant due to the disease; and e is the death rate constant due to other causes.

As in the PCA model, here the total number of individuals is also constant, because $d{S}_{\left\{X\right\}}\left(t\right)/dt+d{I}_1\left(t\right)/dt+d{S}_{\{X-\{1\left\}\right\}}/dt+dV\left(t\right)/dt=0$, then $S_{\left\{X\right\}}\left(t\right)+I_1\left(t\right)+S_{\{X-\{1\left\}\right\}}\left(t\right)+V\left(t\right)=N$. Moreover, $a_1$, p, b, c, e can be estimated from simulations, since the ODE model is a mean-field approximation. Based on [31], the expressions that link these models are

$$\begin{array}{cc}\hfill a_1& \simeq \frac{\mathsf{\Delta }{I}_1{\left(t\right)}_{S_{\left\{X\right\}}\to I_1}}{S_{\left\{X\right\}}\left(t\right)I_1\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill p& \simeq \frac{\mathsf{\Delta }V{\left(t\right)}_{S_{\left\{X\right\}}\to V}+\mathsf{\Delta }V{\left(t\right)}_{S_{\{X-\{1\left\}\right\}}\to V}}{(S_{\left\{X\right\}}\left(t\right)+S_{\{X-\{1\left\}\right\}}\left(t\right))\mathsf{\Delta }t}\hfill \\ \hfill b& \simeq \frac{\mathsf{\Delta }{S}_{\{X-\{1\left\}\right\}}{\left(t\right)}_{I_1\to S_{\{X-\{1\left\}\right\}}}}{I_1\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill c& \simeq \frac{\mathsf{\Delta }{S}_{\left\{X\right\}}{\left(t\right)}_{I_1\to S_{\left\{X\right\}}}}{I_1\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill e& \simeq \frac{\mathsf{\Delta }S{\left(t\right)}_{S_{\{X-\{1\left\}\right\}}\to S_{\left\{X\right\}}}+\mathsf{\Delta }S{\left(t\right)}_{V\to S_{\left\{X\right\}}}}{(S_{\{X-\{1\left\}\right\}}\left(t\right)+V\left(t\right))\mathsf{\Delta }t}\hfill \end{array}$$

(2)

where

• $\mathsf{\Delta }{I}_1{\left(t\right)}_{S_{\left\{X\right\}}}\to I_1$ is the increase in $I_1\left(t\right)$ due to new infection cases from $S_{\left\{X\right\}}$ per time step;
• $\mathsf{\Delta }V{\left(t\right)}_{S_{\left\{X\right\}}\to V}$ and $\mathsf{\Delta }V{\left(t\right)}_{S_{\{X-\{1\left\}\right\}}\to V}$ are the increase in V due to vaccination of $S_{\left\{X\right\}}$ and $\mathsf{\Delta }{S}_{\{X-\{1\left\}\right\}}$ per time step;
• $\mathsf{\Delta }{S}_{\{X-\{1\left\}\right\}}{\left(t\right)}_{I_1\to S_{\{X-\{1\left\}\right\}}}$ is the increase in $S_{\{X-\{1\left\}\right\}}$ due to recovering of $I_1$ per time step;
• $\mathsf{\Delta }{S}_{\left\{X\right\}}{\left(t\right)}_{I_1\to S_{\left\{X\right\}}}$ is the increase in $S_{\left\{X\right\}}$ related to death by the disease of $I_1$;
• $\mathsf{\Delta }S{\left(t\right)}_{S_{\{X-\{1\left\}\right\}}\to S_{\left\{X\right\}}}$ and $\mathsf{\Delta }S{\left(t\right)}_{V\to S_{\left\{X\right\}}}$ are the increase in $S_{\left\{X\right\}}$ due to deaths from another cause of $S_{\{X-\{1\left\}\right\}}$ and V per time step.

In this system, the disease-free equilibrium (DFE) corresponds to a stationary solution satisfying $S^{\prime }=I_1^{\prime }=I_2^{\prime }=S_{X-\left\{1\right\}}^{\prime }=S_{X-\left\{2\right\}}^{\prime }=d{S}_{X-\{1,2\}}^{\prime }=V^{\prime }=0$ represented by $(S_X^{*},I_1^{*},S_{X-\left\{1\right\}}^{*},V^{*})=(e/(p+e),0,0,p/(p+e))$. This point in the state space is of a particular interest since it is used to calculate the basic reproduction number using the next-generation matrix [33] method. From Equation (1), the infectious subsystem is given by

$$\begin{array}{cc}\hfill \{\frac{d{I}_1\left(t\right)}{dt}& =a_1{S}_{\left\{X\right\}}\left(t\right)I_1\left(t\right)-b{I}_1\left(t\right)-c{I}_1\left(t\right)\hfill \end{array}$$

(3)

For this subsystem, the Jacobian matrix in the disease-free equilibrium $(S_X^{*},I_1^{*},S_{X-\left\{1\right\}}^{*},V^{*})$$=(e/(p+e),0,0,p/(p+e\left)\right)$ is

$$\begin{array}{c}\hfill J_{(I_1,I_2)}=\left[\begin{array}{c}{a}_{1}e/(p+e)-(b+c)\end{array}\right]\end{array}$$

(4)

Then, the Jacobian matrix is decomposed as $T+\mathsf{\Sigma }$, where T is the transmission part, containing the production of new infected, and $\mathsf{\Sigma }$ is the transition part, containing changes in the states [32], only for the infected states. With these matrices, we calculate the dominant eigenvalue (spectral radius) $\rho$ of the matrix $-T{\mathsf{\Sigma }}^{-1}$. Therefore,

$$\begin{array}{c}\hfill T=\left[\begin{array}{c}{a}_{1}e/(p+e)\end{array}\right]\end{array}$$

(5)

and

$$\begin{array}{c}\hfill \mathsf{\Sigma }=\left[\begin{array}{c}-(b+c)\end{array}\right]\end{array}$$

(6)

The eingenvalue of $-T{\mathsf{\Sigma }}^{-1}$ is $a_1.e/(b+c).(p+e)$. Then, $R_0$ is given by

$$\begin{array}{ccc}\hfill R_0& =& \frac{a_{1}e}{(b+c)(p+e)}\hfill \end{array}$$

(7)

as analyzed in [31] for the case $p=0$. If $R_0>1$, the disease is endemic, and thus present in the permanent regime of the system; if $R_0<1$, the disease is simply extinguished from the population.

If we consider two variants in the population, the ODE system presented in Equation (1) can be expanded to

$$\begin{array}{ccc}\hfill \frac{d{S}_{\left\{X\right\}}\left(t\right)}{dt}& =& -a_1{S}_{\left\{X\right\}}\left(t\right)I_1\left(t\right)-a_2{S}_{\left\{X\right\}}\left(t\right)I_2\left(t\right)+\hfill \\ & & c(I_1\left(t\right)+I_2\left(t\right))+\hfill \\ & & e(S_{\{X-\{1\left\}\right\}}\left(t\right)+S_{\{X-\{2\left\}\right\}}\left(t\right)+S_{\{X-\{1,2\left\}\right\}}\left(t\right)+V\left(t\right))-p{S}_{\left\{X\right\}}\left(t\right)\hfill \\ \hfill \frac{d{I}_1\left(t\right)}{dt}& =& a_1{S}_{\left\{X\right\}}\left(t\right)I_1\left(t\right)+a_{12}{S}_{\{X-\{2\left\}\right\}}\left(t\right)I_1\left(t\right)-(b+c)I_1\left(t\right)\hfill \\ \hfill \frac{d{I}_2\left(t\right)}{dt}& =& a_2{S}_{\left\{X\right\}}\left(t\right)I_2\left(t\right)+a_{21}{S}_{\{X-\{1\left\}\right\}}\left(t\right)I_2\left(t\right)-(b+c)I_2\left(t\right)\hfill \\ \hfill \frac{d{S}_{\{X-\{1\left\}\right\}}\left(t\right)}{dt}& =& -a_{21}{S}_{\{X-\{1\left\}\right\}}\left(t\right)I_2\left(t\right)+b_{11}{I}_1\left(t\right)-e{S}_{\{X-\{1\left\}\right\}}\left(t\right)-p{S}_{\{X-\{1\left\}\right\}}\left(t\right)\hfill \\ \hfill \frac{d{S}_{\{X-\{2\left\}\right\}}\left(t\right)}{dt}& =& -a_{12}{S}_{\{X-\{2\left\}\right\}}\left(t\right)I_1\left(t\right)+b_{21}{I}_2\left(t\right)-e{S}_{\{X-\{2\left\}\right\}}\left(t\right)-p{S}_{\{X-\{2\left\}\right\}}\left(t\right)\hfill \\ \hfill \frac{d{S}_{\{X-\{1,2\left\}\right\}}\left(t\right)}{dt}& =& b_{12}{I}_1\left(t\right)+b_{22}{I}_2\left(t\right)-e{S}_{\{X-\{1,2\left\}\right\}}\left(t\right)-p{S}_{\{X-\{1,2\left\}\right\}}\left(t\right)\hfill \\ \hfill \frac{dV\left(t\right)}{dt}& =& p(S_{\left\{X\right\}}\left(t\right)+S_{\{X-\{1\left\}\right\}}\left(t\right)+S_{\{X-\{2\left\}\right\}}\left(t\right))-eV\left(t\right)\hfill \end{array}$$

(8)

where the new constants are $a_{12}$ as the infection rate of $S_{\{X-\{1\left\}\right\}}$, generating an infected individual $I_1$, and $a_{12}$ as the infection rate of $S_{\{X-\{1\left\}\right\}}$, generating an infected individual $I_1$. The total number of individuals is also constant, because $S_{\left\{X\right\}}^{\prime }+I_1^{\prime }+I_2^{\prime }+S_{\{X-\{1\left\}\right\}}^{\prime }+S_{\{X-\{2\left\}\right\}}^{\prime }+V^{\prime }=0$, then $S_{\left\{X\right\}}\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+S_{\{X-\{1\left\}\right\}}\left(t\right)+S_{\{X-\{2\left\}\right\}}\left(t\right)+V\left(t\right)=N$. Note that $b_{11}+b_{12}=b$ and $b_{21}+b_{22}=b$. From the $I_1$ individuals that are cured from the disease, part is recovered only for strain one, and part is recovered for both strains; the same is true for $I_2$ individuals. The new estimations of the constant rates from PCA are

$$\begin{array}{cc}\hfill a_2& \simeq \frac{\mathsf{\Delta }{I}_2{\left(t\right)}_{S_{\left\{X\right\}}\to I_2}}{S_{\left\{X\right\}}\left(t\right)I_1\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill a_{12}& \simeq \frac{\mathsf{\Delta }{I}_1{\left(t\right)}_{S_{\{X-\{2\left\}\right\}}\to I_1}}{S_{\{X-\{2\left\}\right\}}\left(t\right)I_1\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill a_{21}& \simeq \frac{\mathsf{\Delta }{I}_2{\left(t\right)}_{S_{\{X-\{1\left\}\right\}}\to I_2}}{S_{\{X-\{1\left\}\right\}}\left(t\right)I_2\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill p& \simeq \frac{\mathsf{\Delta }V{\left(t\right)}_{S_{\left\{X\right\}}\to V}+\mathsf{\Delta }V{\left(t\right)}_{S_{\{X-\{1\left\}\right\}}\to V}+\mathsf{\Delta }V{\left(t\right)}_{S_{\{X-\{2\left\}\right\}}\to V}}{(S_{\left\{X\right\}}\left(t\right)+S_{\{X-\{1\left\}\right\}}\left(t\right)+S_{\{X-\{2\left\}\right\}}\left(t\right))\mathsf{\Delta }t}\hfill \\ \hfill b_{11}& \simeq \frac{\mathsf{\Delta }{S}_{\{X-\{1\left\}\right\}}{\left(t\right)}_{I_1\to S_{\{X-\{1\left\}\right\}}}}{I_1\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill b_{12}& \simeq \frac{\mathsf{\Delta }{S}_{\{X-\{1,2\left\}\right\}}{\left(t\right)}_{I_1\to S_{\{X-\{1,2\left\}\right\}}}}{I_1\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill b_{21}& \simeq \frac{\mathsf{\Delta }{S}_{\{X-\{2\left\}\right\}}{\left(t\right)}_{I_2\to S_{\{X-\{2\left\}\right\}}}}{I_2\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill b_{22}& \simeq \frac{\mathsf{\Delta }{S}_{\{X-\{1,2\left\}\right\}}{\left(t\right)}_{I_2\to S_{\{X-\{1,2\left\}\right\}}}}{I_2\left(t\right)\mathsf{\Delta }t}\hfill \\ \hfill c& \simeq \frac{\mathsf{\Delta }{S}_{\left\{X\right\}}{\left(t\right)}_{I_1\to S_{\left\{X\right\}}}+S_{\left\{X\right\}}{\left(t\right)}_{I_2\to S_{\left\{X\right\}}}}{(I_1\left(t\right)+I_2\left(t\right))\mathsf{\Delta }t}\hfill \\ \hfill e& \simeq \frac{\mathsf{\Delta }S{\left(t\right)}_{S_{\{X-\{1\left\}\right\}}\to S_{\left\{X\right\}}}+\mathsf{\Delta }S{\left(t\right)}_{S_{\{X-\{2\left\}\right\}}\to S_{\left\{X\right\}}}+\mathsf{\Delta }S{\left(t\right)}_{V\to S_{\left\{X\right\}}}}{(S_{\{X-\{1\left\}\right\}}\left(t\right)+S_{\{X-\{2\left\}\right\}}\left(t\right)+V\left(t\right))\mathsf{\Delta }t}\hfill \end{array}$$

(9)

with the new expressions being $\mathsf{\Delta }{I}_2{\left(t\right)}_{S_{\left\{X\right\}}}\to I_2$ as the increase in $I_2\left(t\right)$ due to new infection cases from $S_{\left\{X\right\}}$ per time step; $\mathsf{\Delta }{I}_1{\left(t\right)}_{S_{\{X-\{2\left\}\right\}}}\to I_1$ is the increase in $I_1\left(t\right)$ due to new infection cases from $S_{\{X-\{2\left\}\right\}}$ per time step; $\mathsf{\Delta }{I}_2{\left(t\right)}_{S_{\{X-\{1\left\}\right\}}}\to I_2$ is the increase in $I_2\left(t\right)$ due to new infection cases from $S_{\{X-\{1\left\}\right\}}$ per time step; $\mathsf{\Delta }V{\left(t\right)}_{S_{\left\{X\right\}\to V}}$ is the increase in V due to vaccination of $\mathsf{\Delta }{S}_{\left\{X\right\}}$ per time step; $\mathsf{\Delta }V{\left(t\right)}_{S_{\{X-\{1\left\}\right\}}\to V}$ is the increase in V due to vaccination of $\mathsf{\Delta }{S}_{\{X-\{1\left\}\right\}}$ per time step; $\mathsf{\Delta }V{\left(t\right)}_{S_{\{X-\{2\left\}\right\}}\to V}$ is the increase in V due to vaccination of $\mathsf{\Delta }{S}_{\{X-\{2\left\}\right\}}$ per time step; $\mathsf{\Delta }{S}_{\{X-\{1,2\left\}\right\}}{\left(t\right)}_{I_1\to S_{\{X-\{1,2\left\}\right\}}}$ is the increase in $S_{\{X-\{1\left\}\right\}}$ due to recovering of $I_1$ previously cured from strain two per time step; $\mathsf{\Delta }{S}_{\{X-\{2\left\}\right\}}{\left(t\right)}_{I_2\to S_{\{X-\{2\left\}\right\}}}$ is the increase in $S_{\{X-\{2\left\}\right\}}$ due to recovering of $I_2$ per time step; $\mathsf{\Delta }{S}_{\{X-\{1,2\left\}\right\}}{\left(t\right)}_{I_2\to S_{\{X-\{1,2\left\}\right\}}}$ is the increase in $S_{\{X-\{1,2\left\}\right\}}$ due to recovering of $I_2$ previously cured from strain one per time step; $\mathsf{\Delta }{S}_{\left\{X\right\}}{\left(t\right)}_{I_2\to S_{\left\{X\right\}}}$ is the increase in $S_{\left\{X\right\}}$ related to death by the disease of $I_2$; and $\mathsf{\Delta }S{\left(t\right)}_{S_{\{X-\{2\left\}\right\}}\to S_{\left\{X\right\}}}$ is the increase in $S_{\left\{X\right\}}$ due to deaths from another cause of $S_{\{X-\{2\left\}\right\}}$ per time step.

In this system with two variants, the disease-free equilibrium (DFE) corresponds to a stationary solution satisfying $S_X^{\prime }=I_1^{\prime }=I_2^{\prime }=S_{X-\left\{1\right\}}^{\prime }=S_{X-\left\{2\right\}}^{\prime }=S_{X-\{1,2\}}^{\prime }=V^{\prime }=0$ represented by $(S_X^{*},I_1^{*},I_2^{*},S_{X-\left\{1\right\}}^{*},S_{X-\left\{2\right\}}^{*},S_{X-\{1,2\}}^{*},V^{*})=(e/(p+e),0,0,0,0,0,p/(p+e))$. From Equation (8), the infectious subsystem is given by

$$\{\begin{array}{cc}\frac{d{I}_1\left(t\right)}{dt}& =a_1{S}_{\left\{X\right\}}\left(t\right)I_1\left(t\right)+a_{12}{S}_{\{X-\{2\left\}\right\}}\left(t\right)I_1\left(t\right)-(b+c)I_1\left(t\right)\\ \frac{d{I}_2\left(t\right)}{dt}& =a_2{S}_{\left\{X\right\}}\left(t\right)I_2\left(t\right)+a_{21}{S}_{\{X-\{1\left\}\right\}}\left(t\right)I_2\left(t\right)-(b+c)I_2\left(t\right)\end{array}$$

(10)

For this subsystem, the Jacobian matrix in the disease-free equilibrium $(S_X^{*},I_1^{*},I_2^{*},S_{X-\left\{1\right\}}^{*},$$S_{X-\left\{2\right\}}^{*},S_{X-\{1,2\}}^{*},V^{*})=(1,0,0,0,0,0,0)$, is

$$\begin{array}{c}\hfill J_{(I_1,I_2)}=\left[\begin{array}{cc}{a}_{1}e/(p+e)-(b+c)& 0\\ 0& a_{2}e/(p+e)-(b+c)\end{array}\right]\end{array}$$

Decomposing the Jacobian matrix as $T+\mathsf{\Sigma }$,

$$\begin{array}{c}\hfill T=\left[\begin{array}{cc}{a}_{1}e/(p+e)& 0\\ 0& a_{2}e/(p+e)\end{array}\right]\end{array}$$

and

$$\begin{array}{c}\hfill \mathsf{\Sigma }=\left[\begin{array}{cc}-(b+c)& 0\\ 0& -(b+c)\end{array}\right]\end{array}$$

The eingenvalues of $-T{\mathsf{\Sigma }}^{-1}$ are: $a_{1}e/(b+c)(p+e)$ and $a_{2}e/(b+c)(p+e)$. Then, $R_0$ is given by

$$\begin{array}{ccc}\hfill R_0& =& \underset{1\le i\le 2}{max}\left({\displaystyle \frac{a_{i}e}{(b+c)(p+e)}}\right).\hfill \end{array}$$

(11)

Expanding the system to more than two strains can become laborious. However, due to the $R_0$ calculation by next-generation matrix, the basic reproduction number will always have the form presented before. Therefore, considering $\sigma$ strains in the population, we have

$$\begin{array}{ccc}\hfill R_0& =& \underset{1\le i\le \sigma }{max}\left({\displaystyle \frac{a_{i}e}{(b+c)(p+e)}}\right).\hfill \end{array}$$

(12)

In the next section, the results are presented.

3. Results

The cellular automata is configured with a square lattice of size $n=1000$, and $N_L=$ 1,000,000 cells. At the beginning of the simulation, the states are homogeneously distributed over the lattice. To compare PCA and ODE models with one variant, one simulation in PCA is run with $S_{\left\{X\right\}}\left(0\right)=0.9998$, $I_1\left(0\right)=0.0002$, $k=0.3$, $P_C=0.2$, $P_D=0.2$, $P_N=0.1$, and $\gamma =\mathrm{10,000}$ vaccinations per time step (from time step $t=52$). After the system reaches the permanent regime, Equation (2) is used to estimate ODE constant rates. Figure 2 contains the time evolution of the PCA and ODE approaches.

For two variants, a PCA simulation is run with $S_{\left\{X\right\}}\left(0\right)=0.9996$, $I_1\left(0\right)=0.0002$, $I_2\left(0\right)=0.0002$, $k=0.15$, $P_C=0.2$, $P_D=0.2$, $P_N=0.1$, and $\gamma =5000$ vaccinations per time step (from time step $t=52$). From PCA, the ODE approach with two variants is simulated using the constant rates from Equation (9). Figure 3 contains the time evolution of the PCA and ODE approaches.

Note that there is a good agreement between the approaches. Other parameters have been tested, always returning such similar results for PCA and ODE.

If we assume that a single time step of the simulation corresponds to one week of real time, and if the average duration during which an infected individual can spread the virus is roughly three weeks, we can estimate the probability of cure as $P_C=1/3$. Additionally, we consider the probability of an infected individual dying from COVID-19 as 1%, denoted as $P_D=0.01$ [34,35,36]. Given the life expectancy of Brazil, which is approximately 76.3 years [37], we assume $P_n=1/(76.3\times 52)$. To achieve an $R_0$ value within the range of $[1,6]$, we adjust the values of $C=8$, $r=16$, and $k=0.15$ [36,38,39].

A summary of the variables and parameters of the PCA and ODE models is presented in

Table 1. Description of variables and parameters used for PCA and ODE models. * Note that C, r, and k were chosen to have an $R_0$ value within the range of $[1,6]$; we adjust the values to $C=8$, $r=16$, and $k=0.15$.

Variable/Parameter	Description	Value	Reference
n	lattice side	1000	-
$N_L$	total number of cells	1,000,000	-
$P_i$	probability of infection	-	-
k	infectivity parameter	0.15	* [36,38,39]
$P_C$	probability of cure	1/3	[34,35,36]
$P_D$	probability of dying due to disease	0.01	[34,35,36]
$P_N$	probability of dying due to other causes	$1/(76.3 \times 52)$	[37]
C	number of contacts per individual	8	* [36,38,39]
r	maximum radius for contact in the lattice	16	* [36,38,39]
$S_{\left\{X\right\}}\left(t\right)$	the number of individuals susceptible to	-	-
	all variants at time t
$S_{\{X-1\}}\left(t\right)$	the number of individuals susceptible to	-	-
	all variants except one at time t
$I_1\left(t\right)$	the number of infectious individuals at	-	-
	time t who were infected by strain 1
$I_2\left(t\right)$	the number of infectious individuals at	-	-
	time t who were infected by strain 2
$V\left(t\right)$	the number of vaccinated individuals at time t	-	-

.

Prior to analyzing the impact of mutations, Figure 4 presents the temporal evolution of the normalized cumulative deaths resulting from the disease, as a function of two key variables, k and vaccine efficacy $v_e$ (Figure 4a and Figure 4b, respectively). The vaccine efficacy utilized in all subsequent simulations is fixed at 100%; nevertheless, this figure enables consideration of two additional scenarios corresponding to lower vaccine efficacies. It is noteworthy that as k increases, so does the total number of deaths. The increase in deaths from $k=0.0095$ to $k=0.0100$ highlights how even a minor increase in k may significantly impact disease transmission as a result of increased mutational frequency. The simulations were executed under the conditions of $\gamma =5000$ vaccinations per time step, $v_e=1$, $\mu =1.5\times {10}^{-5}$, and $k=0.01$, except where indicated otherwise.

Once the parameters have been defined, it becomes feasible to assess the mutation and vaccination rates and generate simulation scenarios for the COVID-19 pandemic over a period of five years. The model posits that the mutation leading to the emergence of a new disease variant occurs during the infection process. In particular, when an individual is infected, there exists a probability $\mu$ that a new strain and disease variations may arise. Given that the vaccines were not immediately available, we assume that $\gamma$ non-infected individuals are vaccinated per time step after 52 time steps of the simulation have elapsed. Consequently, it becomes possible to evaluate the impact of the values of $\mu$ and $\gamma$ on the course of the pandemic’s evolution.

In order to accomplish this objective, we have conducted ten simulations for a duration of 260 time steps, encompassing all possible combinations of $\mu$ and $\gamma$ within the range of $\mu =1\times {10}^{-6}, 2\times {10}^{-6}, \cdots , 100\times {10}^{-6}$ and $\gamma =1000, 2000, \cdots ,$ 30,000, resulting in a total of 12,000 simulations. Our analysis has primarily centered on five outputs, namely:

• The number of successful strains by the conclusion of the simulation. Specifically, those strains that have caused more than 1000 cases will be considered in the final results.
• The average number of new infection cases per time step at the end of the five-year simulation period will be calculated across all strains.
• The total number of deaths caused by all strains during the simulation.
• The duration of the pandemic, measured in weeks, specifically the time taken for zero cases per time step to occur.
• The basic reproduction number at the conclusion of the five-year simulation period.

These results, in function of $\mu$ and $\gamma$, are graphically represented in Figure 5. Notably, our findings underscore that a low vaccination rate coupled with a high mutation frequency results in a significant increase in multiple epidemiological metrics, including the number of successful mutations, the average count of new cases five years following the commencement of the pandemic, the total mortality rate, and the time required for the disease to be eradicated. However, it is pertinent to mention that the duration of the simulation was limited to 260 time steps, and hence the plateau depicted in Figure 5d signifies that the disease was not eradicated by the end of the simulation period.

The final result to be presented pertains to the basic reproduction number at the culmination of the five-year period. As depicted in Figure 6, $R_0$ is illustrated as a function of the vaccination rate and mutation, utilizing the average value for the last ten temporal steps. It is noteworthy that for a vaccination rate exceeding 5000 vaccinations per week, the contagion will no longer disseminate within the population beyond a five-year span. This is an important result, substantiating that persistent and uninterrupted vaccination measures can culminate in the eradication of an epidemic.

In conclusion, our study underscores the importance of vaccination rates in controlling the spread of disease, particularly in the context of high mutation frequency. The simulation results highlight the significant impact of low vaccination rates and high mutation frequency on various epidemiological metrics, and the importance of maintaining sustained vaccination efforts in controlling the spread of disease.

4. Discussion

In the present study, a novel model has been proposed that utilizes probabilistic cellular automata and ordinary differential equations to simulate the transmission dynamics of a virus that may generate several active strains within the host population due to mutations. The basic reproduction number of the disease has been derived by employing the next-generation matrix process. Moreover, a vaccination strategy has been evaluated to mitigate the epidemic and confront the virus mutations over a period of five years. Although the model entails some strong assumptions, the findings underscore the critical significance of a consistent and comprehensive vaccination plan in curtailing the emergence of new mutations and ultimately eradicating the pandemic.

In this study, the COVID-19 pandemic served as a reference point for investigating a novel multi-strain disease model. Notably, the PCA flexibility, which is endorsed by the ODE approach, obviates the problem of high-dimensional analytical analysis commonly encountered in other publications in the field, where typically two [17,19] or three [20] strains are considered. The model presented herein is transferable to pathogenic diseases exhibiting drug resistance [17,40]. Although the model assumes equal epidemiological characteristics for all strains, the calculated basic reproduction number can be applied to strains possessing diverse transmissibility, illness severity, and mortality.

The world has undergone a considerable experience in combating the COVID-19 pandemic. Since December 2019, numerous mutations of the SARS-CoV-2 virus have dominated populations both locally and globally, exhibiting rapid dispersion, even during periods of lower international travel. The results presented herein exhibit two positive outcomes: vaccination reduces the number of successful SARS-CoV-2 mutations, and the incidence of COVID-19 cases and deaths decreased within the considered period, irrespective of the mutation rate. Additionally, Figure 5d illustrates a remarkable finding: regardless of the virus’s mutation rate, vaccinating 30,000 individuals (in our lattice) per time step (week) can end the pandemic within a year of vaccine availability. Since vaccination commenced after one year in the simulation (52 weeks), the pandemic could be over by March 2022, for all mutation rates. However, the pandemic’s end could be postponed by lower vaccine efficacy, individuals refusing to vaccinate in some countries, and mutations that evade vaccines. The significance of an early vaccination program has been substantiated in prior studies [20]. Furthermore, multi-strain models that incorporate vaccination strategies, assuming a two-dose regimen for achieving complete immunity, have demonstrated that one strain tends to prevail, ultimately leading to disease eradication when a substantial proportion of the population has received both vaccine shots, as outlined in the chosen parameters [16].

According to a recent study conducted by Oliveira et al. [41], the pandemic in Brazil could be effectively terminated within a year if a vaccination rate of 14 million individuals per week was achieved from the vaccine availability. Notably, the study’s findings were predicated upon the administration of a vaccine with a 50% efficacy rate, specifically Sinovac Biotech’s Coronavac, which was initially employed in Brazil’s vaccination campaign. In the current study, we have considered the use of a vaccine with a 100% efficacy rate, and our calculations have been based on a weekly dose administration of 30,000, which translates to 6 million doses for a population of approximately 200 million people in Brazil. It is imperative to acknowledge that this difference in the efficacy of the vaccine and the corresponding dose administration may affect the estimated timeline for terminating the pandemic.