Evolutionary analysis of replicator dynamics about anti-cancer combination therapy

1.   Introduction

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) was first discovered in the city of Wuhan, Hubei province, China, in December 2019 ^[1] and has since spread all over the world with over $181.1$ million reported cases and over $3.9$ million deaths globally as of June 29, 2021 ^[2]. Following the outbreak of COVID-19, several non-pharmaceutical interventions (NPIs) such as physical distancing, isolation, hand washing, stay-at-home order, closing of schools and businesses, travel restrictions, among others, were implemented all over the world to limit the spread of the disease ^[3,4,5,6]. Despite implementing these NPIs, there is still a significant number of COVID-19 cases and deaths reported daily. On January 20, 2020, the World Health Organization (WHO) declared COVID-19 a public health emergency of international concern (PHEIC) ^[7] and a pandemic on March 11, 2020 ^[8]. The COVID-19 pandemic has proven to be more than just a medical emergency; But one that affects society. It has significantly affected our day-to-day activities such as socializing, entertainment, education, tourism, health, and business. The need to find an effective way of controlling the pandemic through integrated science cannot be overemphasized ^[9].

Since the outbreak of the SARS-CoV-2 virus, researchers have been studying its transmission dynamics. The virus is believed to be transmitted through two major pathways namely: the direct and indirect pathways ^{[10,11,12,13,14]}. The direct pathway, also known as human-to-human transmission, includes transmission through inhalation of virus-infected droplets, coughing, sneezing and having physical contact with infectious persons. Indirect transmission may occur when a susceptible individual comes in contact with a contaminated commonly shared surface or object ^[10,13]. These contaminations usually occur when an infectious individual touches, coughs or sneezes on the objects.

Several mathematical models have been developed to study the transmission dynamics of COVID-19 ^{[15,16,17,18,19,20,21,22,23,24,25,26,27,28]}. In ^[26], two differential equation models were developed to study the effect of the exposed or latency period on the dynamics of COVID-19 in China. Another model was developed in ^[28] to study the effect of super-spreaders on the dynamics of COVID-19 in Wuhan, China. These and many other mathematical models have provided a lot of insights into the dynamics of COVID-19. During the early stage of the COVID-19 pandemic, physical distancing, wearing of face masks, hand washing, and other non-pharmaceutical interventions (NPIs) were the primary interventions against the disease worldwide. Many mathematical models have been used to study the impact of these NPIs on the dynamics of COVID-19 ^[19,29]. This includes the susceptible-exposed-infected-quarantine-recovered (SEIQR) model developed in ^[25] to study the impact of physical distancing on the dynamics of COVID-19 in British Columbia (BC), Canada, and the model of ^[23] used to study the impact of different mitigation strategies on the dynamics of COVID-19 in Ontario, Canada. In ^[30], a {susceptible-exposed-infected-removed} (SEIR) model was designed on a temporal network, which evolves according to the activity-driven paradigm. Their model was used to study the spread of the disease as a function of the fraction of the population following public health measures. Their result shows that physical distancing and mask-wearing can effectively prevent COVID-19 outbreaks if adherence to both measures involves a substantial fraction of the susceptible population. The effect of timing of adherence to COVID-19 NPIs such as social distancing on COVID-19 in the United States was studied using an agent-based model in ^[31]. Their model was applied to Dane County, Wisconsin, the Milwaukee metropolitan area and New York City (NYC), and show that the timing of implementing and easing the social distancing measures has major effects on the number of COVID-19 cases. The results obtained in the works of ^[32,33] are also in the same line with the findings mentioned above. Many of these models focus on direct transmission of COVID-19 through in-person contacts despite the emphasis by the world health organization (WHO) on indirect transmission of the disease ^[13].

A few researchers have studied indirect transmission of COVID-19 ^{[20,21,22,34,35,36,37]}. In ^[35], a novel deterministic susceptible-exposed-infected-removed-virus-death (SEIRVD) model was developed to study the potential impact of both direct and indirect transmission of COVID-19 on the dynamics of the disease in Ontario, Canada. Their results showed significant increase in number of cases with indirect transmission. They also highlighted the importance of implementing additional preventive and control measures to minimize the spread of COVID-19 through both transmission routes. A five-compartment model for both direct and indirect transmission pathways was developed in ^[37], and used to study the impact of both transmission pathways on the dynamics of COVID-19. The mathematical model of ^[36] also used indirect transmission mechanisms to study the dynamics of COVID-19. They examined the conditions under which contaminated objects may lead to a significant spread of SARS-CoV-2 during and after lock-down using a SEIR model with the addition of a fomite term. They proposed that the addition of a fomite term will help to better understand the transmission dynamics of the virus and also in policy making. To the best of our knowledge, none of these studies have looked at the effect of adherence and non-adherence to COVID-19 NPIs on both direct and indirect transmissions of the disease.

In this study, we develop a susceptible-exposed-infectious-recovered-virus (SEIRV) compartmental model for studying the transmission dynamics of COVID-19. Our model considers both direct and indirect transmission of the disease following the approach of ^[35]. It has a similar structure to the model of ^[25] used to study the impact of physical distancing on the contact rates in British Columbia (BC), Canada, where BC population was divided into the groups practicing physical distancing and those not practicing physical distancing. Here, we stratify our population into two groups: those that adhere to all COVID-19 non-pharmaceutical interventions (NPIs) and those that do not adhere. We compute the control reproduction number and final epidemic size for our model, and study the effect of different parameters of the model on these quantities. This model is used to assess the impact of adherence to COVID-19 non-pharmaceutical interventions and indirect transmissions on the dynamics of the disease. Our study highlights the importance of adhering to the COVID-19 NPIs as instructed by the World Health Organization (WHO) in reducing infections and eventually controlling the spread of the disease.

The remaining part of this paper is organized as follows: The model is developed in Section 2. In Section 3, the control reproduction number of the disease is calculated using the next generation matrix approach. We also computed a final size relation for the epidemic in this section. In Section 4, we numerically study the effect of important model parameters on the control reproduction number calculated in Section 3. In addition, we present numerical simulations of the SEIRV model for different scenarios in this section. These simulations are used to comprehend the effect of direct and indirect transmission, with or without adherence on the dynamics of the disease. We conclude the paper in Section 5 with a brief discussion.

2.   Mathematical model

A compartmental SEIRV model of COVID-19 epidemic is developed and analyzed. The model considers both direct and indirect transmission routes and divides the population into two main groups: adherent and non-adherent populations. The adherent population include individuals in the population that adhere to all COVID-19 non-pharmaceutical interventions (NPIs), while the non-adherent population contains those that do not adhere to the NPIs. The compartments of the model for non-adherent population are defined as follows: susceptible $(S)$; exposed $(E)$; infectious $(I)$; and recovered/removed $(R)$. Corresponding compartments for the adherent population are given by $S_a$, $E_a$, $I_a$, and $R_a$, respectively. Note that individuals in the exposed compartments $E$ and $E_a$ are not infectious ^{[26,38,39,40,41]}. We assume that recovered individuals have permanent immunity from the disease (no reinfection). The compartment $V$ accounts for the virus shed on contaminated surfaces and objects by infectious individuals.

Figure 1 shows a schematic diagram of the model, where the solid black arrows indicate the direction of the flow of individuals between the compartments at the rates indicated beside the arrows. The dashed red arrows show virus shedding or contamination of surfaces by infectious individuals. The differential equations for non-adherent population are given by

where $\beta$ and $\beta_v$ are the direct and indirect transmission rates, respectively, for non-adherent population. The parameter $\sigma_a$, with $0 \leq \sigma_a \leq 1$ is used to model the reduction in virus shedding (or contamination of surfaces), susceptibility, and onward transmission of disease due to adherence to NPIs. Parameters $\gamma_a$ and $\gamma$ are the rates of movement from adherence to non-adherence and vice versa, respectively, while $\mu$ denotes the rate of transitioning from the exposed compartment $E$ to the infectious compartment $I$, and $r$ is the recovery rate. The corresponding equations for the adherent population are given by

where $S_a, E_a, I_a$ and $R_a$ denote susceptible, exposed, infectious, and recovered individuals in the adherent population, respectively. The parameters $\beta_a$ and $\beta_{av}$ are the direct and indirect disease transmission rates, respectively, for the adherent population. We assume that individuals in both adherent and non-adherent population transition from exposed to infectious compartments, and recover from the disease at the same rates, $\mu$ and $r$, respectively. In addition, we assume that individuals move from adherent to non-adherent and vice versa, at the same rates in the susceptible, exposed, infectious, and recovered compartments. The disease dynamics for both non-adherent (2.1a) and adherent (2.1b) populations are coupled to the dynamics of the virus in the environment or on contaminated surfaces given by

where $m$ and $m_a$ are the virus shedding (or surface contamination) rates for non-adherent and adherent populations, respectively. The parameter $\delta$ is the decay rate of SAR-CoV-2 in the environment and on surfaces, and $\tau$ is the environment cleaning/sanitization rate. We assume that our population is constant throughout the epidemic period so that $N = N_{\star}(t) + N_a(t)$ for all time $t \geq 0$, where $N_{\star}(t) = S(t) + E(t) + I(t) + R(t)$ and $N_a(t) = S_a(t) + E_a(t) + I_a(t) + R_a(t)$. In our SEIRV model (2.1), the compartments $S, E, I, R, S_a, E_a, I_a, R_a$ have dimension of individual, while the virus compartment $V$ has dimension of virus particle. The dimensions of the model parameters are given as follows:

where $[\lambda]$ represents dimension of $\lambda$. Note that $\sigma_a$ is a dimensionless parameter in our model.

It is important to mention that there are relationships between the parameters $\beta, \beta_v$, and $m$ for non-adherent population and $\beta_a, \beta_{av}$, and $m_a$ of adherent population, respectively. These relationships are used when computing the control reproduction number in Section 3.1 and final size relation in Section 3.2. We define the fraction of adherence in the population as $f = \gamma/(\gamma + \gamma_a)$, on which the dynamics of the ODE system (2.1) depends. Since adherence to COVID-19 NPIs reduces transmission rate for the adherent population, we equate the direct and indirect transmission rates for the adherent population to those of non-adherent population as follow

where $0 \leq \sigma_a \leq 1$. When $\sigma_a = 0$, we have $\beta_a = 0$, $\beta_{av} = 0$ and $\sigma_a I_a = 0$ in (2.1a) and (2.1b), which implies that adherence to NPIs is 100$\%$ effective in preventing infections in the adherent population, and that adherent individuals do not infect the non-adherent individuals. On the other hand, when $\sigma_a = 1$, we have $\beta_a = \beta$, $\beta_{av} = \beta_v$ and $\sigma_a I_a = I_a$, which implies that disease is transmitted at the same rate in both adherent and non-adherent populations. Table 1 shows the model variables and their description.

Table 2 shows the parameters of the SEIRV model (2.1) with their descriptions and values. Some of these parameters were obtained from published articles, while others were either inferred, derived or varied. References are provided for parameters that were obtained from published articles. The direct and indirect transmission rates ($\beta_a$ and $\beta_{av}$), and shedding rate $m_a$ for adherent population were derived from those of non-adherent population, respectively, using the relations in (2.3).

3.   Model analysis

In this section, we calculate the control reproduction number and the final size relation for the epidemic using the SEIRV model (2.1). We also consider some limiting scenarios of the epidemic and derive the control reproduction number and final size relations for these scenarios from the general control reproduction number and final size relation, respectively.

3.1.   Control reproduction number

We compute the control reproduction number for our model (2.1) using the next generational matrix approach of ^{[42,43,44,45,46,47]} as used in ^[48]. The next generation matrix is a square matrix whose elements describes who is infecting whom in the population, and is defined as $\mathcal{M} = \mathcal{F}\mathcal{V}^{-1}$, where $\mathcal{F}$ and $\mathcal{V}$ are square matrices for the rates of new infections from the susceptible population and transfer between infected populations, respectively. The dominant eigenvalue of the next generation matrix is the basic (or control) reproduction number of the epidemic. See ^[42,43,47] for more details. To construct the $\mathcal{F}$ and $\mathcal{V}$ matrices, we consider the exposed and infectious compartments of the SEIRV model (2.1) given by,

Upon using the relations between the transmission rates for the adherent and non-adherent populations given in (2.3), (3.1) reduce to

Recall that the adherence fraction in the population is given by $f = \gamma/ (\gamma + \gamma_a)$. Therefore, at the beginning of the epidemic, $S(0) = (1-f)N = N_{\star}(0)$ and $S_a(0) = f N = N_a(0)$. Following the arrangement of the equations in (3.2) and the variables $\phi = (E, I, E_a, I_a, V)$, at the DFE = $(N_{\star}(0), 0, 0, 0, N_a(0), 0, 0, 0, 0)$, the matrix for {the rates of new infections from the susceptible population $\mathcal{F}$} is given by

Similarly, the Jacobian matrix for the rates of transfer between infected populations at the beginning of the epidemic is given by

Upon constructing the next generation matrix using (3.3) and (3.4), and computing its spectral radius using the symbolic computational software SageMath ^[74], we obtain that the control reproduction number satisfies

where $R_{D}$, which accounts for disease transmission through direct route among non-adherent and adherent individuals is given by

Here, $R_{D \star}$ and $R_{Da}$ are defined as

The second term in the control reproduction number (3.5a) accounts for indirect transmissions, and it is given by

where $R_{ID \star}$ and $R_{IDa}$ are the indirect transmissions contributed by the non-adherent and adherent populations, respectively, and are given by

The newly introduced variables $\chi$, $\chi_a$, $\Phi$, and $\Psi$ are defined as follows:

Suppose there is no adherent population, that is, {$f = 0$ and $\gamma = 0$ (consequently, $\sigma_a = 0$ and $\gamma_a = 0$), which gives $\Phi = \mu r$ and $\Psi = 0$}. Therefore, (3.5a) reduces to

where $\hat{R}_D = \beta/r$ accounts for infections cause through direct transmission, while $\hat{R}_{ID} = \beta_v m/ (r(\delta + \tau))$ accounts from those of indirect transmission. Observe that $\hat{R}_{D}$ will increase as both the transmission rate $\beta$ and {infectious period} $1/r$ increase. This shows that an infectious individual will spread the virus/disease more if the transmission rate of the {disease increases} and if it takes longer for them to recover from the disease. For $\hat{R}_{ID} = \beta_v m/ (r(\delta + \tau))$, infections increases as the shedding/contamination rate of infectious individuals increase and decrease as the virus decay rate $\delta$ and the environment cleaning or sanitization rate $\tau$ increases. Assuming a fraction $\tilde{f}$ of the population were already adhering to COVID-19 NPIs before the disease was introduced into the population (as a precaution before the first case arrives) and there is no movement from adherent to non-adherent and vice versa, that is, $\gamma = 0$ and $\gamma_a = 0$, but $\sigma_a \neq 0$, (3.5a) reduces to

where the parameter $\sigma_a$ with $0 \leq \sigma_a \leq 1$ account for the reduction in acquiring or transmitting the disease due to adherence to COVID-19 NPIs. Similar to (3.5a) and (3.6), the control reproduction number $\tilde{R}_c$ for this scenario also accounts for the contribution from both direct and indirect transmissions. Here, the contribution to $\tilde{R}_c$ due to direct transmissions is given by

where $\beta(1 - \tilde{f})/r$ accounts for infections through direct transmissions caused by those that do not adhere to the NPIs, while $\sigma_a^2 \tilde{f} \beta/r$ is the contribution from those that adhere to the NPIs. As one would expect, the transmissions caused by the adherent population will always be less than the infections caused by those that do not adhere if $\tilde{f} < 0.5$ since $0 \leq \sigma_a \leq 1$. For this scenario, the contribution from indirect transmission is given by

where $\beta_v (1-\tilde{f})m/[ r(\delta + \tau)]$ is the contribution from the non-adherent population and $\sigma_a \beta_v m_a \tilde{f}/[ r(\delta + \tau)]$ is from the adherent population. Note that $\tilde{f}$ is the fraction of the population that adhere to the COVID-19 NPIs at the beginning of the epidemic.

3.2.   Final size relation

Next, we derive a relation for the final epidemic size of our SEIRV model (2.1). This final size relation gives the total number of cases of the disease in the population during the epidemic. We consider a scenario where there is no movement between adherent and non-adherent susceptible populations. In other words, we assume that susceptible individuals do not change their behaviour from adhering to the COVID-19 NPIs to non-adherence and vice versa. This way, the equation for the susceptible populations $S$ and $S_a$ in (2.1a) and (2.1b), respectively, reduce to

We assume that a fraction $f$, with $0 \leq f \leq 1$, of the susceptible population adhere to COVID-19 NPIs before a disease outbreak, so that $S(0) = (1-f)N = N_{\star}(0)$ and $S_a(0) = f N = N_a(0)$. In addition, we assume that at the beginning of the epidemic, {there is a small number $E_a(0) = E_{a0}$ and $I_a(0) = I_{a0}$ in the adherent population, with $E(0) = E_0$ and $I(0) = I_0$ in the non-adherent populations. All other compartments are initially empty}. Following the approach in ^{[48,65,66,67,68,69,70,71,72,73]} and using the ODE system (2.1) together with the assumptions for the susceptible population in (3.10), we derive the final size relation for non-adherent population as

where $\clubsuit_1 = R_{D \star} + \, R_{ID \star}$ and $\spadesuit_1 = R_{D a} + R_{ID a}$ with $R_{D \star}, R_{ID \star}, R_{D a}$ and $R_{ID a}$ as given in (3.5). The remaining variables $\diamondsuit_1$ and $\heartsuit_1$ are defined as follows

Similarly, the final size relation for adherent population is given by

where $\clubsuit_2 = \sigma_a\clubsuit_1, \spadesuit_2 = \sigma_a \spadesuit_1, \diamondsuit_2 = \sigma_a \diamondsuit_1$ and $\heartsuit_2 = \sigma_a \heartsuit_1$ with $\clubsuit_1, \spadesuit_1, \diamondsuit_1$ and $\heartsuit_1$ as given in (3.11). The final size relations in (3.11) and (3.12) for non-adherent and adherent population, imply that $S_\infty > 0$ and $S_{a \infty} > 0$, respectively. They give the relationship between the reproduction number $R_c$ and the final epidemic size. The total infected population over the course of the epidemic is given by $N_a(0) - S_{a \infty}$ and $N_{\star}(0) - S_\infty$ for the adherent and non-adherent population, respectively, which can be described in terms of the attack rates as $\left[1- \frac{S_{a \infty}}{N_a(0)} \right]$ and $\left[1- \frac{S_\infty}{N_{\star}(0)} \right]$ as in ^[67]. The final size relations in (3.11) and (3.12) can be written together in matrix form as

Also, from (3.11) and (3.12), we derive the following relation between the two final sizes

Since $0\leq \sigma_a \leq 1$, we have

When $\frac{S_{a\infty}}{S_{a0}} > \frac{S_\infty}{S_0}$, we conclude that the attack rate for the adherent population is smaller than that of the non-adherent population as a result of the reduction in the susceptibility and transmissibility of the adherent population. Observe that when $\sigma_a = 1$, which implies that there is no protection from adhering to NPIs, we have from (3.14) that $S_{a \infty} / S_{a0} = S_{\infty}/S_0$. This implies that the {attack rate} is the same in both adherent and non-adherent populations, and as a result, we expect the epidemic size to be the same in both populations. This result agrees with our numerical simulation shown in red in the top left panel of Figure 4. On the other hand, when $\sigma_a = 0$, which implies that adhering to NPIs is $100\%$ effective in preventing new infections, (3.14) reduces to $S_{a \infty} = S_{a0}$. This shows that there are no infections in the adherent population throughout the epidemic. Therefore, adherence to the NPIs during the pandemic is beneficial to reducing/containing the spread of the disease.

When the outbreak begins with no infectious individuals, that is, the outbreak begins through an indirect transmission, the final size relations in (3.11) and (3.12) reduce to

Similarly, when the outbreak begins with only direct transmission with an infective, (3.11) and (3.12) become

Now, consider an instance where there is no adherent population, so that $f = 0, \sigma_a = 0, \gamma = 0$, and $\gamma_a = 0$. For this scenario, we have $\clubsuit_1 = \hat{R}_c$, where $\hat{R}_c$ is as given in (3.6), $\spadesuit_1 = 0$, $I_a(0) = 0$, $N_{\star}(0) = N$, $\clubsuit_1 = \heartsuit_1$. Also, since $\sigma_a = 0$, we have $\clubsuit_2 = \spadesuit_2 = \diamondsuit_2 = \heartsuit_2 = 0$. This way, (3.11) and (3.12) reduces to

where $\hat{R}_c = \beta/r + \beta_v m/(r(\delta + \tau))$ is the control reproduction number, given in (3.6), for the scenario without adherent population. As expected, when there is no adherent population, the final size relation reduces to that of a single population, and depends on the control reproduction number. When a fraction $\tilde{f}$ of the population adhere to COVID-19 NPIs but there is no movement between the adherent and non-adherent populations throughout the epidemic, that is, $\gamma = \gamma_a = 0$ and $\sigma \neq 0$. The final size relation for non-adherent population given in (3.11) reduces to

where

Similarly, for the adherent population, the final size relation (3.12) reduces to

where $\tilde{\clubsuit}_2 = \sigma_a \tilde{\clubsuit}_a$ and $\tilde{\spadesuit}_2 = \sigma_a \tilde{\spadesuit}_1$ with $\tilde{\clubsuit}_1$ and $\tilde{\spadesuit}_1$ as defined in (3.20). Let $\tilde{S}_{0}$ and $\tilde{S}_{\infty}$ be the total susceptible population at the beginning and the end of the epidemic, respectively. Then at the beginning of the epidemic, we have

where $S_0$ and $S_{a0}$ are the initial susceptibles in non-adherent and adherent populations, respectively. Since there is no movement between the adherent and non-adherent population, the following conditions are satisfied at the end of the epidemic

where $S_{\infty}$ and $S_{a \infty}$ are the susceptibles left in non-adherent and adherent populations, respectively, at the end of the epidemic. Upon adding (3.19) and (3.21), and using (3.22) and (3.23), and $m_a = \sigma_a m$, we have

where $\hat{R}_c$ and $\tilde{R}_c$ are as given in (3.6) and (3.7), respectively.

In the next section, we present numerical computations of the control reproduction number (3.5) for different parameters of the model. We also present numerical simulations of the SEIR model (2.1), where we numerically computed final epidemic sizes for different scenarios.

4.   Numerical simulations

In this section, we study the effect of different parameters of the model on the control reproduction number $R_c$, given in (3.5). These results are presented as contour plots. In addition, we present numerical simulations of the SEIRV model (2.1) and investigate the effect of different parameters of the model on the disease dynamics and final epidemic size.

In Figure 2, we present contour plots of the control reproduction number (3.5) computed as a function of different parameters of the model (2.1). The top left panel of this figure shows $R_c$ as a function of the direct and indirect transmission rates. As one would expect, $R_c$ increases as both transmission rates increase. We also notice that the direct transmission rate $\beta$ has more effect on the control reproduction number than the indirect transmission rate $\beta_v$. This may be because indirect transmissions are mitigated by environment cleaning and virus decay rate. In the top right panel of the same figure, we have $R_c$ as a function of the adherence rate $\gamma$ and the direct transmission rate $\beta$. This result shows that the adherence rate has little effect on the control reproduction number when the transmission rate is low, but has more effect on it as the transmission rate increases. When the transmission rate is low, there is a lesser chance of getting infected, as a result, adherence to the NPIs will have little effect on the spread of infection. However, when transmission rate is high, adhering to the NPIs will significantly reduce disease transmission. Similar result is shown in the middle left panel of Figure 2 for the adherence rate $\gamma$ and the indirect transmission rate $\beta_v$. For this scenario, $R_c$ increases as the adherence to NPIs decreases even for a low indirect transmission rate, unlike the case of direct transmission (top right panel of Figure 2), where adherence rate does not affect $R_c$ when $\beta$ is small. The reason for this is that indirect transmission does not only depend on the transmission rate. It also depends on the environment cleaning/sanitization rate and the shedding and decay rate of the virus. If a lot of shedding is happening, or the decay rate of the virus is low, there may be a need to adhere more to the NPIs to reduce the spread of the disease.

For a fixed indirect transmission rate, $R_c$ increases as the virus shedding rate increases, although, if the transmission rate is low, $R_c$ remains constant irrespective of the shedding rate. Similarly, for a fixed shedding rate $m$ for non-adherent population, $R_c$ increases as $\beta_v$ increases, but remains constant when shedding rate is too low irrespective of the value of $\beta_v$ (see the middle right panel of Figure 2). This result shows that when infectious individuals do not shed much virus or contaminate surfaces, there will be no indirect disease transmission through these surfaces. Similar result is obtained when $R_c$ is computed as a function of the shedding rate $m_a$ for the adherent population and the indirect transmission rate $\beta_v$. In the bottom panel of Figure 2, we present $R_c$ as a function of the environment cleaning/sanitization rate and the indirect transmission rate $\beta_v$ (left panel), and the environment cleaning/sanitization rate and the shedding rate of non-adherent individuals $m$ (right panel). Both results show that $R_c$ decreases as the environment cleaning/sanitization rate increases, and increases with increase in the transmission and shedding $\beta_v$ and $m$, respectively.

Figure 3 shows the total infectious population $\big{(} I(t) + I_a(t) \big{)}$ over time for different rates of movement from adherent population to non-adherent and vice versa. These results are used to investigate the effect of adherence to COVID-19 NPIs on the dynamics of the disease. We consider four (4) different scenarios based on the rate of movement from adherent to non-adherent population $\gamma_a$ and the rate of moving from non-adherent to adherent population $\gamma$ for different values of the parameter $\sigma_a$, used to model the reduction in susceptibility, shedding of virus, and onward disease transmission for the adherent population. The first scenario shown in the top left panel of Figure 3 is for $\gamma = \gamma_a = 0$, where there is no movement of individuals between non-adherence and adherence populations. For this scenario, the epidemic final size is $60.19\%$ when $\sigma_a = 0.5$, $80.71\%$ when $\sigma_a = 0.75$, and $91.43\%$ when $\sigma_a = 1$. Recall that the epidemic final size is the total fraction of the population that becomes infected during the course of the epidemic. As expected, the largest epidemic occurs when there is no protection derived from adhering to the NPIs, that is, when $\sigma_a = 1$. Epidemic decreases as the protection from NPIs increases ($\sigma_a$ decreases). Similar results are observed for the remaining scenarios in this figure. When there is movement between the adherent and non-adherent population such that the rate of going from adherence to non-adherence is the same as that of going from non-adherence to adherence, specifically, when $\gamma = \gamma_a = 0.3$, the final epidemic size is $59.50\%$ for $\sigma_a = 0.5$, $81.20\%$ for $\sigma_a = 0.75$ and $91.43\%$ for $\sigma_a = 1$. The plots for this scenario are shown in the top right panel of Figure 3.

In the bottom panel of Figure 3, we show the results for $\gamma = 0.25$ and $\gamma_a = 0.3$ (left panel), and $\gamma = 0.3$ and $\gamma_a = 0.25$ (right panel). For the result in the left panel where the rate of movement from the adherent population to non-adherent ($\gamma_a$) is higher than that of going from non-adherence to adherence ($\gamma$), the epidemic final sizes are $64.86\%, 82.49\%$ and $91.43\%$ for $\sigma_a = 0.5, \sigma_a = 0.75$ and $\sigma_a = 1$, respectively. On the other hand, for $\gamma = 0.3$ and $\gamma_a = 0.25$ (bottom right panel of Figure 3), where the rate of movement from the non-adherence to adherence is higher than that of going from adherent population to non-adherent, the epidemic final size is $53.41\%$ for $\sigma_a = 0.5$, $79.83\%$ for $\sigma_a = 0.75$, and $91.43\%$ for $\sigma_a = 1$. We observe from the results in this figure that the epidemic size decrease as $\sigma_a$ decreases. This shows that fewer infections are occurring as people get more protection from adhering to the NPIs, since the benefit of adhering to NPIs increases as $\sigma_a$ decreases. In addition, the epidemic size when $\sigma_a = 1$ (no protection due to adherence to NPIs) is the same for any values of $\gamma$ and $\gamma_a$. For $\sigma_a = 0.5$ and $\sigma_a = 0.75$, the smallest epidemics were predicted for the scenario where the rate of movement from non-adherence to adherence $\gamma$ is higher than that of moving from adherence to non-adherence $\gamma_a$ (bottom right panel of Figure 3). Whereas the largest epidemics occur when the rate of adherence $\gamma = 0.25$ is lower than that of non-adherence $\gamma_a = 0.3$. These results show that there will be smaller epidemic if more people adhere to the NPIs.

The results in Figure 4 were obtained using the same parameters and initial conditions as those in Figure 3 but plotted in terms of the adherent and non-adherent populations. Presenting the results in this way enables us to study the dynamics of the two populations as the movement rates $\gamma$ and $\gamma_a$ are varied. We notice from the results in this figure that even though the final epidemic size is the same for all four scenarios when $\sigma_a = 1$ (as discussed earlier in Figure 3), the dynamics of the disease is different for each scenario. When $\gamma = \gamma_a$, the disease dynamics is the same for both adherent and non-adherent population when $\sigma_a = 1$ (see the red plots in the top panel of Figure 4). This is because the population is divided into two equal halves between the adherent and non-adherent population for these scenarios, and no protection is derived from adhering to the NPIs (since $\sigma_a = 1$). When $\sigma_a \neq 1$, for the scenario with $\gamma = \gamma_a = 0.3$, the dynamics of the disease in the two populations are the same irrespective of $\sigma_a$ (see top right panel of Figure 4). However, this is not the case for $\gamma = \gamma_a = 0$, even though $f = 0.5$ in both scenarios. For $\sigma_a = 0.5$ and $\sigma_a = 0.75$, there are more infections in non-adherent population than the adherent population. This is as a result of the protection derived from adhering to the COVID-19 NPIs. When $\gamma = 0.25$ and $\gamma_a = 0.3$, there are more infections in the non-adherent population than the adherent population for all values of $\sigma_a$. This is due to the fact that there are more people in the non-adherent population since the non-adherence rate is higher ($\gamma_a > \gamma$) and they do not get protection from the NPIs. The most interesting and counter-intuitive result is when $\gamma = 0.3$ and $\gamma_a = 0.25$ (bottom right panel of Figure 4). This result shows that there are more infections in the adherent population than non-adherent population. Recall that $\gamma$ is the adherence rate and $\gamma_a$ is the non-adherence rate. For $\gamma = 0.3$ and $\gamma_a = 0.25$, the initial fraction of the population that adhere to COVID-19 NPIs is $f \approx 0.5455$, which is larger than the population of those that do not adhere. As time goes on, more people adhere to COVID-19 NPIs than they stop adhering since $\gamma > \gamma_a$. This increasing the adherence population and consequently the adherent-susceptible and -infectious population. Therefore, causing more infections in the adherence population than in the non-adherent population. Note that even though there are more infections in the adherent population than the non-adherent population, the final epidemic size is still smaller for this scenario relative to that of the $\gamma = 0.25$ and $\gamma_a = 0.3$ scenario for $\sigma_a = 0.75$ and $\sigma_a = 0.5$, as discussed earlier and shown in the bottom left panel of Figure 4. The epidemic size is the same for both scenario when $\sigma_a = 1$.

In Figure 5, we investigate the effect of the indirect transmission rate $\beta_v$ and the shedding rate $m$ for non-adherent population on the dynamics of the disease. When the shedding rate is $m = 0.02$ (left panel of Figure 5), the final epidemic size is $46.66\%$ for $\beta_v = 0.15$, $48.50\%$ for $\beta_v = 0.35$, and $50.25\%$ for $\beta_v = 0.55$. As the shedding rate increases to $m = 0.125$ (right panel of Figure 5), the final epidemic size for $\beta_v = 0.15$ is $53.44\%$, for $\beta_v = 0.35$ is $62.11\%$, and $68.87\%$ for $\beta_v = 0.55$. This shows that the final epidemic size increases, as the transmission rate $\beta_v$ increases. When the shedding rate is low, there is not much difference in the epidemic size as the transmission rate $\beta_v$ increases. However, for a higher virus shedding rate or surface contamination rate, there is a significant difference in the final epidemic size as $\beta_v$ increases. This is because when infectious individuals do not shed much virus or contaminate shared surfaces, there will be no disease transmission even if the transmission rate (implicitly $R_c$) is high. On the other hand, when a lot of shedding or surface contamination is happening, the disease will spread more easily leading to bigger epidemics as $\beta_v$ increases. This results agree with the $R_c$ computation shown in the middle left panel of Figure 2. A similar result was obtained for the shedding $m_a$ for the adherent population (result not included).

Next, we study the effect of different initial conditions on the dynamics of the diseases. We consider two main scenarios: when the epidemic starts mainly through direct transmissions and when it starts mainly through indirect transmissions. For the scenario where the epidemic starts through direct transmission, we assume that there are only twenty infected individuals in the population at the beginning of the epidemic; five in the exposed and another five in the infectious compartments of the adherent population. Similar infected population is specified for the non-adherent population, that is, $E(0) = 5$, $E_a(0) = 5$, $I_a(0) = 5$ and $I(0) = 5$, and $V(0) = 1$, while the remaining compartments are empty. On the other hand, when the epidemic starts mainly through indirect transmissions, the initial conditions used are $E(0) = 1$, $E_a(0) = 1$, $I_a(0) = 1$, $I(0) = 1$, and $V(0) = 10$, with the remaining compartments empty. The results for both scenarios are shown in Figure 6 for different values of the movement parameters $\gamma$ and $\gamma_a$, where the solid curves are for the case where epidemic starts mainly through direct transmissions, and the dashed curves are for when the epidemic starts mainly through indirect transmissions. We observe from these results that for all the cases considered, the total infectious population increases more rapidly when direct transmissions dominate at the beginning of the epidemic, compared to when it starts mainly through indirect transmissions. However, the final epidemic size is the same irrespective of the initial conditions. This result agrees with the $R_c$ computation in the top left panel of Figure 2 that shows that the control reproduction number is more sensitive to the direct transmission rate $\beta$ compared to the indirect transmission rate $\beta_v$.

The results in Figure 7 show the effect of the environment cleaning or sanitization rate on the disease dynamics. When $\beta_v = 0.15$, the final epidemic size for $\tau = 0.1$ is $53.44\%$, $51.26\%$ for $\tau = 0.3$, and $50.00\%$ when $\tau = 0.5$. As the indirect transmission rate increases to $\beta_v = 0.5$, the final epidemic size also increases with $67.33\%$ for $\tau = 0.1$, $62.38\%$ for $\tau = 0.3$, and $59.23\%$ for $\tau = 0.5$. Overall, the final epidemic size decreases as the rate of cleaning/sanitizing the environment increases. When the indirect transmission rate $\beta_v$ is small, there is not much difference in the final epidemic size as $\tau$ increases, but for a higher transmission rate $\beta_v$, an increase in the environment sanitizing rate significantly decreases the epidemic final size. These results are consistent with the $R_c$ computation in the bottom left panel of Figure 2, where there is no significant difference in the control reproduction number as $\tau$ increase for a small value of $\beta_v$, and a significant decrease in $R_c$ as $\tau$ increases for a higher indirect transmission rate.

5.   Discussion

We have developed and analyzed a SEIRV model for studying the dynamics of COVID-19. Our model divides the population into two groups: adherent and non-adherent, and considers both direct and indirect transmission of the disease. Individuals in the adherent population are assumed to adhere to all COVID-19 non-pharmaceutical interventions (NPIs) such as hand washing, physical distancing, wearing face masks, and avoiding large gatherings. In contrast, those in the non-adherent population do not adhere to these NPIs. By adhering to the NPIs, those in the adherent population have reduced their susceptibility, transmissibility and infectiousness. {By direct transmission, we mean disease transmission that occurs when a susceptible individual comes directly in contact with the saliva and other respiratory droplets expelled by an infectious individual when they talk, sing, cough or sneeze}, while indirect transmission refers to the infections that occur through contaminated shared surfaces and objects such as door handles, utensils, etc. Individuals can move from adherent to non-adherent and vice versa at some rate in each human compartments.

We calculated the control reproduction number and final size relation for our model and studied the effect of the model parameters on the calculated control reproduction number. These results are presented as contour plots in Figure 2 and show that the direct transmission rate has more effect on the control reproduction number than the indirect transmission rate. We believe this is because indirect transmission depends on other factors such as shedding rate and environment cleaning rate. Similarly, when we studied $R_c$ with respect to the adherence rate $\gamma$ and the direct transmission rate $\beta$, we notice that adherence to NPIs has little effect on $R_c$ when the transmission rate is small. As the direct transmission rate increases, our results show that adhering to the NPIs would significantly reduce $R_c$. On other hand, adherence to NPIs is beneficial irrespective of the values of the indirect transmission rate. We also showed that the virus shedding rate has little effect on the control reproduction number if the indirect transmission rate is too small. However, as the transmission rate increases, an increase in shedding rate increases the control reproduction number significantly.

Furthermore, we studied the effect of some model parameters on the disease dynamics. We showed that the disease dynamics is different when there is no movement between adherent and non-adherent population, compared to when the movement between the two groups happens at the same rate. When there is no movement between the groups, there are less infections in the adherent population compared to the non-adherent population as a result of protection from adhering to NPIs. On the contrary, when movement happens between the two groups at the same rate, the level of infection in the two populations is the same, simply because the entire population becomes homogeneous as a result of the movement between adherent and non-adherent population happening at the same rate. When the non-adherence rate is higher than the adherence rate, there are more infections in the non-adherent population relative to the adherent population. In a situation where the adherence rate is higher than the non-adherence rate, the infection in the adherent population is higher compared to the non-adherent population, although, one would have expected the infections in the adherent population to be low because of the protection from adhering to the NPIs. This counterintuitive observation is due to the movement parameters $\gamma$ and $\gamma_a$, and the parameter $\sigma_a$ used to model protection in the adherent population due to NPIs. Observe from the results in Figure 4 that we have more infections in the adherent population only when the rate $\gamma$ of adhering to COVID-19 NPIs is larger than the non-adherence rate $\gamma_a$. For this scenario, the initial fraction of the total population adhering to the COVID-19 NPIs is larger than the non-adherence population, and as the epidemic continues, we have more people moving to the adherence population that they are becoming non-adherent. Therefore, increasing the adherent susceptible and infected populations, and consequently leading to a larger epidemic in this population compared to the non-adherent population. As the protection derived from adhering to the NPIs increases, fewer adherent individuals are getting infected as shown in the bottom left panel of Figure 4. Overall, even though there are more infections in the adherent population when the rate of movement from non-adherent to adherent is higher than that of moving vice versa, the final epidemic size is smaller for this scenario than the other scenarios. This result shows the benefits of adhering to the COVID-19 NPIs in reducing the epidemic size.

Also, we examined the effect of different initial conditions on the dynamics of the diseases. We considered when the epidemic starts mainly through direct transmissions and when it starts mainly through indirect transmissions. For the scenario where the epidemic starts through direct transmissions, we assume that only twenty infected individuals are in the population at the beginning of the epidemic with a unit of virus on contaminated surfaces. On the other hand, when the epidemic starts mainly through indirect transmissions, we assume that there are only four infected individuals in the population and ten units of viruses on contaminated surfaces at the beginning of the epidemic. These initial conditions are carefully selected to model the dominance of one transmission route over the other at the beginning of the epidemic. Our results show that the epidemic increases more rapidly when it starts mainly through direct transmissions, compared to when it starts through indirect transmissions, although the final epidemic size is the same irrespective of the initial conditions. This result agrees with the result obtained when we studied the effect of the direct and indirect transmission rates on the control reproduction number. Finally, we study the effect of environment sanitization on disease dynamics. Our results show that environmental sanitization reduces the prevalence of COVID-19 and consequently reduces the final epidemic size.

An immediate and interesting future direction in this project is to fit the model to reported cases of COVID-19 in a specific location and use region-specific parameters. By fitting the model to the cases data, important parameters that are not immediately available can be estimated from the data. The model can then be used to make forecasts/predictions and test the effect of implementing different intervention strategies. Another important future direction is to include vaccination compartments into the model. The COVID-19 vaccines became available and have been administered since early 2021. As of July 6, 2021, over 2.9 billion doses of the vaccine have been administered ^[2]. To properly capture the dynamics of COVID-19 in any population starting from January 2021, vaccination needs to be included in the model. It would be worthwhile to stratify the population by age so that the model can be used to answer age related questions such as those of school reopening and age-specific vaccination plans. There are currently four SARS-CoV-2 variants of concern (VOC) that are believed to be in circulation worldwide. These VOC are believed to be more contagious, leading to more severe sickness than the initial SARS-CoV-2 virus ^[75]. An interesting extension of our model would be to incorporate these variants of concerns.

Our model is prone to some limitations. This includes our assumptions that the population is well-mixed. We know contact rates vary from one person to another depending on their age group and activity level. Also, mixing patterns are different for individuals of different age groups. Using homogeneous mixing, we assume that everyone in the population has the same contact rate and mixing pattern. Another limitation of our model is that we have clustered all COVID-19 non-pharmaceutical interventions (NPIs) into one group, assuming that any individual that adheres to one of them will adhere to all. This may not be the case in reality. Some individuals may adhere to only a few of the NPIs. It would be nice to distinguish between the NPIs and study the effect of different NPIs on the disease dynamics. Similarly, we have assumed that the rate of movement from adherence to non-adherence and vice versa are constant over time. In a real-world scenario, these rates may change from time to time and may be affected by government policies and their implementation of the NPIs. {Further limitation to our study includes lack of stochasticity. An interesting future work may include adding uncertainties and parameter estimation into this modeling framework.} Despite these limitations, our model can capture the overall dynamics of the COVID-19 epidemic considering infections transmitted through direct and indirect routes, and with and without adherence to the COVID-19 NPIs.

Acknowledgements

This research is funded by Canada's International Development Research Centre (IDRC) and the Swedish International Development Cooperation Agency (SIDA) (Grant No. 109559-001). J.D. acknowledges that this research was conducted during her maternity leave.

Conflict of interest

The Authors declare that there is no conflict of interest.