The potential impact of public health care denial on the transmission dynamics of COVID-19 in South Africa

Abstract

When the demand for public health care increases, governments often prioritize citizens over foreign nationals. In South Africa, limited resources and socio-economic inequalities pose unique challenges to epidemic control. The overcrowding and increasing demand for public healthcare have led to protests by some community groups, which have led to the denial of healthcare to migrants. Denying treatment to some infected individuals has the propensity to lead to an increase in the size of an epidemic. We introduce a novel epidemiological model that incorporates health care denial as a dynamic factor influencing the transmission of COVID-19. It incorporates healthcare denial as a key parameter influencing the progression and recovery rates of infections. The study presents a novel framework for understanding the intersection of healthcare access denial and the transmission dynamics of COVID-19. While much of the existing literature has focused on the direct effects of healthcare interventions on pandemic control, this research uniquely emphasizes the role that restricted access to healthcare services, whether due to policy decisions, resource shortages, or system inefficiencies, can exacerbate the spread of infectious diseases. The treatment class of the model is partitioned to account for individuals denied treatment at public healthcare facilities. Analytical results establish conditions for the existence and stability of both disease-free and endemic equilibria, with the basic reproduction number R₀ explicitly derived to quantify transmission potential under varying healthcare access scenarios. Sensitivity analysis reveals that increasing denial of care can significantly elevate R₀, resulting in higher infection peaks, prolonged epidemic duration and greater cumulative mortality. Numerical simulations further illustrate the non-linear relationship between treatment accessibility and outbreak severity. The findings highlight that equitable healthcare provision is not only a public health necessity but also a critical determinant for reducing the COVID-19 burden. Policy implications stress the integration of inclusive healthcare strategies to ensure epidemic resilience and minimize transmission risks, especially in vulnerable populations. Strategies that will accommodate every infected person who goes to the hospital for treatment should be adopted to reduce the disease burden.

Introduction

South Africa has been among the leading countries in the COVID-19 pandemic in Africa, often detecting new waves earlier than neighbouring states. ¹ During successive waves (e.g. driven by the Beta, Delta and Omicron variants), incidence rates, test positivity and mortality rose in step with more transmissible variants, even as testing capacity improved. ² Regionally in southern Africa, South Africa has functioned as an early warning signal: waves tend to emerge there first, and then spread to neighbouring countries on average several weeks later. ¹ At the global level, these trends mirror broader patterns: more transmissible variants causing sharper increases in cases and, over time, a decoupling of case numbers and death rates in places with rising immunity (through prior infection or vaccination).^2,3

These sharp increases led to overcrowding in public hospitals beyond the hospital's capacity. In 2024, it was reported that almost 84% of South Africans rely on the public health system for their medical needs. ⁴ The COVID-19 pandemic not only put further strain on the healthcare sector but also partially exposed the vulnerability of high-risk demographic groups to medical care. ⁵ A pandemic like COVID-19 can cause massive hospital and clinic overcrowding. In impoverished nations, healthcare services are sometimes insufficient to combat epidemic diseases. Due to the insufficient infrastructure to fight against pandemics in the health care system in South Africa, some patients are denied treatment in preference to others. South Africa, despite being the most developed country in Sub-Saharan Africa, has faced significant healthcare system challenges, especially in terms of resource availability, healthcare workers and capacity to handle a large number of COVID-19 cases. The country's public health system struggles with inequities in healthcare access, especially in rural or economically disadvantaged areas. Unlike in other sub-Saharan countries, where treatment denial is mainly due to overwhelmed healthcare facilities, which led to triaging of patients, with those in more critical conditions prioritized for treatment, In South Africa, access to medical treatment is prioritized for individuals in critical condition, with citizens receiving priority over non-citizens. The impact of treatment denial in each of these different countries is the same: a heightened risk of severe disease and death. Healthcare denial also likely increased community transmission, as individuals with untreated COVID-19 symptoms may have continued interacting with others, spreading the virus in communities, making the disease control very difficult. In South Africa, immigrants are frequently refused care due to their nationality and/or lack of valid identification^6–13 For instance, A two-year-old Zimbabwean boy died of poisoning after being turned away at Charlotte Maxeke Hospital in South Africa because his mother could not produce his birth certificate or pay Rs. 5,000. ¹⁴ Also, a 3-year-old girl from Uganda with cerebral palsy and epilepsy was denied hospital care in South Africa since administrators demanded payment before opening her file—even though children under six are entitled to free care by law citedoctors. Elfreda Spencer died in the UK in 2018 after being denied chemotherapy for a year because she could not afford the fees. ¹⁵ The denial of healthcare to immigrants in South Africa has been reported in several studies and media accounts^6–13; however, it is important to note that these accounts reflect observed or perceived incidents rather than an officially established or legally endorsed policy. To maintain scientific and ethical integrity, such claims are therefore presented as evidence of reported occurrences within the healthcare system, acknowledging that they may result from individual actions, administrative barriers or systemic challenges rather than a formal national directive. Refusing to treat a person or group of infected people for any reason merely makes the disease more likely to spread, especially when it's as contagious as COVID-19. Such acts could have disastrous results, particularly if a sizable number of people are excluded from treatment. However, over time, this has been disregarded in COVID-19 epidemic modeling.

In the last few years, researchers have used mathematical models to study the impact of various social, economic, and political factors on the transmission dynamics of COVID-19, investigating how certain COVID-19 control strategies affected the dynamics of the disease.^16–19 Eastman et al. ²⁰ created a mathematical model in 2021 to investigate how social behavior affected the pandemic's trajectory in the province of Ontario. Their strategy was predicated on a conventional SEIRD model with a variable transmission rate that is contingent upon population behavior. Their model demonstrated the necessity for attention as the economy reopens by indicating that the effective reproduction number, and consequently the pandemic's trajectory, is sensitive to the observance of public health regulations. Additionally, Rai et al. ²¹ assessed the effectiveness of social media advertisements in combating the coronavirus outbreak in India. They thought that increasing awareness among susceptible populations would alter public attitudes and behavior around COVID-19, reducing the likelihood of contracting the virus and, hence, slowing its spread. Their findings suggest that in order to reduce the prevalence of disease in India, government officials and health authorities should consistently raise awareness via the internet and social media platforms about the need to hospitalize symptomatic individuals and place asymptomatic individuals under quarantine. Wilk et al. ²² more recently modeled how government actions affected the COVID-19 epidemic. They addressed the issue of parameter non-estimability using data from a collection of European countries. They created a set of customized models, including thirteen distinct pandemic control methods, treating a subset of characteristics as universal and the others as country-specific. They then estimated the parameters of these models without making any assumptions beforehand.

In this study, we investigate the effects on disease dynamics of denying treatment to immigrants infected with COVID-19 in South Africa. We consider a scenario in which ill patients are turned away from healthcare services and continue to spread the infection within their communities and networks. Here, we look at a scenario where a group of individuals is denied health care based on nationality, non-documentation and insufficient resources. The study distinguishes two groups of patients seeking treatment: those who are accepted and those who are denied treatment. It uses a dynamical systems model as a tool to assess the possible effects of reintegrating the latter group of people into the community to ascertain how this influences the spread of the disease within the community.

The rest of the article is structured as follows: The article's introduction is provided in Section 1. Model formulation in Section 2, Section 3 presents the model analysis, Section 4 describes the model parameterization and sensitivity analysis, Section 5 presents numerical simulations, Section 6 describes the effects of treatment denial on the disease burden, Section 7 gives the policy implication of the results while Section 8 wraps up the work.

Method

Formulation of compartmental model

To understand the dynamics of COVID-19 transmission in South Africa, this study uses a mathematical model. COVID-19 outbreaks have persisted for over three years since 2020. Demographic effects will be incorporated into this study's model (the important parameters are π and µ, as described in Table 4). The whole population is divided into seven subpopulations. Individuals who are susceptible to the disease (S), those who have been exposed to the virus (E), those who are contagious and either asymptomatic or symptomatic (I_a) and (I_s) respectively, infected individuals who are given treatment (T_a), infected individuals who are denied treatment (T_d) and the recovered individuals (R). Individuals are recruited into the (S) compartment through birth or migration at a rate π. Susceptible individuals are exposed to the virus at a rate λ and transition to the (E) compartment. After a while, these exposed individuals become infectious at a rate σ. A proportion τ of the infectious individuals who are symptomatic move to the (I_s) compartment, while the remaining, (1−τ) who are asymptomatic move to the (I_a) compartment. It is assumed that asymptomatic patients are more likely to transmit COVID-19 than symptomatic ones because their infection often goes undetected. Unlike symptomatic individuals, who may isolate or seek treatment once they feel ill, asymptomatic carriers continue with normal activities, unknowingly spreading the virus. Since no control measures such as testing, isolation or treatment are applied to them, they play a significant role in sustaining and accelerating community transmission. A proportion ϕ of the individuals in the class (I_s) who seek health care are treated and hence moved to the class (T_a) at a rate α₁, while the remaining proportion (1−ϕ) are denied treatment and hence moved to the (T_d) compartment at a rate α₂. Individuals from the (T_a), (T_d) and (I_a) compartments who survive, recover and move into the (R) compartment at rates γ₁, γ₂ and γ₃ respectively. When individuals are denied treatment, their progression through the disease can be different compared to individuals who receive treatment. By partitioning the treatment class, the model can capture these differences in disease dynamics. For instance, those who are denied treatment may experience faster deterioration in health, leading to more severe symptoms or a higher probability of reaching critical stages of the disease. Denied treatment might also increase the likelihood of death or a longer time to recovery, which can influence overall mortality rates in the population. Conversely, those who receive treatment may have a slower progression, higher chances of recovery and reduced mortality, which can slow transmission. Denying treatment can also increase the effective infection rate. Partitioning the treatment class into treated and denied treatment generally allows for a more granular analysis of the system. The deceased compartment is excluded because it does not influence disease transmission dynamics. Once individuals die, they can no longer infect others or transition into other states of the model, making the compartment epidemiologically irrelevant to the spread of infection. This simplification reduces model complexity without compromising the accuracy of predictions related to transmission and control strategies. Tables 1 and 2 describe the related model parameters and state variables, whereas Figure 1 shows the model's flow diagram.

Figure 1.

The model diagram.

Table 1.

Model parameters and their meanings.

Parameter	Description	Unit
π	Recruitment rate	People/day
β	Effective human-to-human contact rate	day−1
µ	Rate of natural death	day−1
σ	The rate at which those exposed become contagious	day−1
γ 1	Recovery rate of treated individuals	day−1
γ 2	Recovery rate of those denied treatment	day−1
δ	Transition rate from asymptomatic to symptomatic class	day−1
α 2	Rate of treatment denial	day−1
ϕ	Proportion of individuals that are given treatment	dimensionless
α 1	Rate of treatment acceptance	day−1
η, η2, η3	Modification parameters	dimensionless

Table 2.

Model state variables and their meanings.

State variable	Description	Unit
S	Susceptible individuals	People
E	Exposed individuals	People
Is	Infected symptomatic individuals	People
Ia	Infected asymptomatic individuals	People
Ta	Individuals given treatment	People
Td	Individuals denied treatment	People
R	Recovered individuals	People

The mathematical model's equations are listed (and explained) in (1) below.

$$\begin{array}{rl}\dot{S}=& \pi -(\mu +\lambda )S,\\ \dot{E}=& \lambda S-(\sigma +\mu )E,\\ {\dot{I}}_a=& \sigma (1-\tau )E-({\gamma }_3+\delta +\mu )I_a,\\ {\dot{I}}_s=& \sigma \tau E+\delta I_a-({\alpha }_1\varphi +(1-\varphi ){\alpha }_2+\mu )I_s,\\ {\dot{T}}_a=& {\alpha }_1\varphi I_s-({\gamma }_1+\mu )T_a,\\ {\dot{T}}_d=& {\alpha }_2(1-\varphi )I_s-({\gamma }_2+\mu )T_d,\\ \dot{R}=& {\gamma }_3{I}_a+{\gamma }_2{T}_d+{\gamma }_1{T}_a-\mu R.\end{array}$$

(1)

with

$$\lambda =\beta (I_a+{\eta }_1{T}_a+{\eta }_2{I}_s+{\eta }_3{T}_d),$$

(2)

where $0<{\eta }_1<{\eta }_2<{\eta }_3<1$ .

Model analysis

In this section, we show that the model under study is mathematically and biologically meaningful, compute the basic reproduction number, and analyze the steady states. The R class is redundant; that is, individuals in this class do not contribute to disease transmission; hence, we omit it in the rest of the analysis.

Feasible region

We prove that the region Ω given below attracts all non-negative and bounded solutions of the model (1),

$$\mathrm{\Omega }=\left\{(S,E,I_a,I_s,T_a,T_d)\in {\mathcal{R}}_{+}^6:N(t)\le {\displaystyle \frac{\pi }{\mu }}\right\}.$$

To show that Ω is positively invariant, consider the first model equation,

$${\displaystyle \frac{dS}{dt}=\pi -(\lambda +\mu )S,}$$

which is solved using an appropriate integrating factor to obtain

$$S(t)=S_0{e}^{-\underset{0}{\overset{t}{\int }}(\mu +\lambda (s))ds}+\left[{\displaystyle \frac{\pi }{\mu }{e}^{-\underset{0}{\overset{t}{\int }}(\mu +\lambda (s))ds\times \underset{0}{\overset{t}{\int }}{e}^s(\mu +\lambda (\tau ))d\tau }}\right]>0.$$

For the solution of the second equation, suppose there is a time t^′ such that $E(t^{\mathrm{\prime }})=0,\dot{E}(t^{\mathrm{\prime }})<0$ and all the other variables are non-negative for 0 < t^′ < t. Then

$${{\displaystyle \frac{dE}{dt}}|}_{t=t^{\mathrm{\prime }}}=\lambda (t^{\mathrm{\prime }})S(t^{\mathrm{\prime }})>0,$$

which is a contradiction. Hence, E(t) > 0. Similarly, we can show that $I_s(t)>0,I_a(t)>0,$ T_d(t) > 0, and $R(t)>0$ .

To show that the solutions of Model (1) are bounded in $\mathrm{\Omega }$ , let $N=S+E+Ia+Is+Ta+Td,$ then,

$${\displaystyle \frac{dN}{dt}=\pi -\mu N-(1+\sigma )\tau E-{\gamma }_3{I}_a-{\gamma }_2{T}_d-{\gamma }_1{T}_a\le \pi -\mu N.}$$

Solving the differential inequality, we obtain

$$N(t)\le {\displaystyle \frac{\pi }{\mu }+\left(N_0-{\displaystyle \frac{\pi }{\mu }}\right)e^{-\mu t}.}$$

(3)

.

As $t\to \mathrm{\infty },$ the right-hand side of model (3) tends to $\frac{\pi }{\mu }$ . Thus, $N(t)\le \frac{\pi }{\mu }.$ The solutions of the equation (1) are thus bounded in $\mathrm{\Omega }.$

Steady states and basic reproduction number

The disease-free equilibrium point (DFE) E₀ is given by

$$E_0=\left({\displaystyle \frac{\pi }{\mu },0,0,0,0,0}\right).$$

The next-generation matrix approach was used to calculate the model reproduction number ²³ is given by

$$R_0={\displaystyle \frac{\pi \beta }{\mu }{\displaystyle \frac{[Q_3{Q}_4{Q}_5(1-\tau )+(\delta (1-\tau )+\tau Q_2)(Q_5(\varphi {\alpha }_1{\eta }_1+Q_4{\eta }_2)+(1-\varphi )Q_4{\alpha }_2{\eta }_3)]}{Q_1{Q}_2{Q}_3{Q}_4{Q}_5},}}$$

where $Q1=\sigma +\mu ,Q2=\delta +\gamma 3+\mu ,Q3=\alpha 1\varphi +(1\text{–}\varphi )\alpha 2+\mu ,Q4=\gamma 1+\mu ,$

$$Q5=\gamma 2+\mu .$$

By solving the system (4)

$$\begin{array}{rl}\pi -(\mu +\lambda )S=& 0,\\ \lambda S-(\sigma +\mu )E=& 0,\\ \sigma (1-\tau )E-({\gamma }_3+\delta +\mu )I_a=& 0,\\ \sigma \tau E+\delta I_a-({\alpha }_1\varphi +(1-\varphi ){\alpha }_2+\mu )I_s=& 0,\\ {\alpha }_1\varphi I_s-({\gamma }_1+\mu )T_a=& 0,\\ {\alpha }_2(1-\varphi )I_s-({\gamma }_2+\mu )T_d=& 0,\end{array}$$

(4)

we obtain the endemic equilibrium point E₁ given by

$${\mathcal{E}}_1=(S*,E*,I_a*,I_s*,T_a*,T_d*),$$

where,

$$\begin{array}{r}\begin{array}{cc}& S*={\displaystyle \frac{\pi }{\mu +{\psi }_5{I}_s*},E*={\psi }_3{\psi }_4{I}_s*}\\ & I_a*={\psi }_4{I}_s*,I_s*={\displaystyle \frac{\pi {\psi }_5-\mu {\psi }_6}{{\psi }_5{\psi }_6},}\\ & T_a*={\psi }_2{I}_s*,T_{d1}={\psi }_1{I}_s*.\end{array}\end{array}$$

$$\begin{array}{rl}{\psi }_1=& {\displaystyle \frac{{\alpha }_2(1-\varphi )}{{\gamma }_2+\mu },{\psi }_2={\displaystyle \frac{{\alpha }_1\varphi }{{\gamma }_1+\mu },{\psi }_3={\displaystyle \frac{{\gamma }_2+\delta +\mu }{\sigma (1-\tau )},{\psi }_4={\displaystyle \frac{{\alpha }_1\varphi +(1-\varphi ){\alpha }_2+\mu }{\sigma \tau {\psi }_3+\delta },}}}}\\ {\psi }_5=& \beta ({\psi }_4+{\eta }_1{\psi }_2+{\eta }_2+{\eta }_3{\psi }_1),{\psi }_6=(\sigma +\mu ){\psi }_3{\psi }_4\end{array}$$

By substituting the expressions for S^∗, E^∗ and $I_a*=$ , into the second equation of (4) and simplifying, we obtain the quadratic equation

$$(a_2{I}_s*+a_1)I_s*=0$$

(5)

where,

$$\begin{array}{rl}{a}_1=& \pi \beta ({\psi }_4+{\eta }_1{\psi }_2+{\eta }_2+{\eta }_3{\psi }_1-(\sigma +\mu ){\psi }_3{\psi }_4\\ a_2=& (\sigma +\mu ){\psi }_3{\psi }_4^2\end{array}$$

The solution $I_s*=0$ represents the DFE point while the solution of the equation

$$a_2{I}_{s1}+a_1=0,$$

(7)

represents the endemic equilibrium point. To get the maximum number of positive roots that the equation (7) can have, we apply Descarte's rule of signs as shown in Table 3. Clearly, the system under investigation has no endemic equilibrium point when a₁ and a₂ are both positive and negative and a unique endemic equilibrium point when a₁ and a₂ have opposite signs.

Table 3.

The maximum number positive roots of equation (5) that could exist.

a 2	a 1	Number of positive roots
+	+	0
+	−	1
−	+	1
−	−	0

Stability analysis of the equilibrium points

The DFE is locally stable when $R0<1$ following Theorem 2 in Reference. ²⁴

Theorem 1.

The endemic equilibrium is globally asymptotically stable in $\mathrm{\Omega }$ if $R0>1$ .

Proof.

If $R0>1$ , we have a unique endemic equilibrium. We consider a candidate Lyapunov function given by

$$\begin{array}{c}V=\left[S-S*-S*\ln \left({\displaystyle \frac{S}{S*}}\right)\right]+k_1\left[E-E*-E*\ln \left({\displaystyle \frac{E}{E*}}\right)\right]+k_2\left[I_a-I_a*-I_a*\ln \left({\displaystyle \frac{I_a}{I_a*}}\right)\right]\\ +k_3\left[I_s-I_s*-I_s*\ln \left({\displaystyle \frac{I_s}{I_s*}}\right)\right]+k_4\left[T_a-T_a*-T_a*\ln \left({\displaystyle \frac{T_a}{T_a*}}\right)\right]+k_5{(T_d-T_d*)}^2,\end{array}$$

where $k1,k2,k3,k4,k5$ are positive constants to be determined and the R class is redundant, hence omitted. At the endemic equilibrium, the partial derivatives with respect to each variable are

$$\begin{array}{r}\begin{array}{cc}& {\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }S}=\left(1-{\displaystyle \frac{S*}{S}}\right),{\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }I}=k_1\left(1-{\displaystyle \frac{E*}{E}}\right),{\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }{I}_a}=k_2\left(1-{\displaystyle \frac{I_a*}{I_a}}\right),}}}\\ & {\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }{I}_s}=k_3\left(1-{\displaystyle \frac{I_s*}{I_s}}\right),{\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }{T}_a}=k_4\left(1-{\displaystyle \frac{T_a*}{T_a}}\right),{\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }{T}_d}=2k_5(T_d-T_d*).}}}\end{array}\end{array}$$

The endemic equilibrium is clearly a critical point of V. The second derivatives are given by

$$\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }S}=\left(1-{\displaystyle \frac{S*}{S}}\right),{\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }I}=k_1\left(1-{\displaystyle \frac{E*}{E}}\right),{\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }{I}_a}=k_2\left(1-{\displaystyle \frac{I_a*}{I_a}}\right),}}}\\ {\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }{I}_s}=k_3\left(1-{\displaystyle \frac{I_s*}{I_s}}\right),{\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }{T}_a}=k_4\left(1-{\displaystyle \frac{T_a*}{T_a}}\right),{\displaystyle \frac{\mathrm{\partial }V}{\mathrm{\partial }{T}_d}=2k_5(T_d-T_d*).}}}\end{array}$$

The second derivatives of V are all positive at any point of $\mathrm{\Omega }$ , therefore, the Lyapunov function $V$ is concave up and the unique endemic equilibrium point is a minimum point of $V.$ We now show that $\frac{dV}{dt}\le 0$ . The time derivative of V is given by

$$\begin{array}{rl}{\displaystyle \frac{dV}{dt}=}& \left(1-{\displaystyle \frac{S*}{S}}\right){\displaystyle \frac{dS}{dt}+k_1\left(1-{\displaystyle \frac{E*}{E}}\right){\displaystyle \frac{dE}{dt}+k_2\left(1-{\displaystyle \frac{I_a*}{I_a}}\right){\displaystyle \frac{d{I}_a}{dt}}}}\\ & +k_3\left(1-{\displaystyle \frac{I_s*}{I_s}}\right){\displaystyle \frac{d{I}_s}{dt}+k_4\left(1-{\displaystyle \frac{T_a*}{T_a}}\right){\displaystyle \frac{d{T}_a}{dt}+2k_5(T_d-T_{d*}){\displaystyle \frac{d{T}_d}{dt},}}}\\ =& \left(1-{\displaystyle \frac{S*}{S}}\right)[\pi -\beta (I_a+{\eta }_1{T}_a+{\eta }_2{I}_s+{\eta }_3{T}_d)S-\mu S)\\ & +k_1\left(1-{\displaystyle \frac{E*}{E}}\right)(\beta (I_a+{\eta }_1{T}_a+{\eta }_2{I}_s+{\eta }_3{T}_d)S-Q_{1}E)\\ & +k_2\left(1-{\displaystyle \frac{I_a*}{I_a}}\right)(\sigma (1-\tau )E-Q_2{I}_a)+k_3\left(1-{\displaystyle \frac{I_s*}{I_s}}\right)(\sigma \tau E+\delta I_a-Q_3{I}_s)\\ & +k_4\left(1-{\displaystyle \frac{T_a*}{T_a}}\right)({\alpha }_1\varphi I_s-Q_4{T}_a)+2k_5(T_d-T_d*)({\alpha }_2(1-\varphi )I_s+Q_5{T}_d).\end{array}$$

(8)

At the endemic equilibrium, the system (1) yields the following:

$$\begin{array}{rl}\pi =& (\lambda *+\mu )S*,Q_1={\displaystyle \frac{\lambda *S*}{E*},Q_2={\displaystyle \frac{\sigma (1-\tau )E*}{I_a*}}}\\ Q_3=& {\displaystyle \frac{\sigma \tau E*+\delta I_a*}{I_s*},Q_4={\displaystyle \frac{{\alpha }_1\varphi I_s*}{T_a*},Q_5={\displaystyle \frac{{\alpha }_2(1-\varphi )I_s*}{T_d*}}}}\end{array}$$

Substituting the above expressions for the constants $\pi ,Q1,Q2,Q3,Q4$ and $Q5$ into (8), we have

$$\begin{array}{rl}& {\displaystyle \frac{dV}{dt}=\left(1-{\displaystyle \frac{S*}{S}}\right)[\beta I_a*S*\left(1-{\displaystyle \frac{S{I}_a}{S*I_a*}}\right)+\beta {\eta }_1{T}_a*S*\left(1-{\displaystyle \frac{T_{a}S}{T_a*S*}}\right)}\\ & +\beta {\eta }_2{I}_s*S*\left(1-{\displaystyle \frac{I_{s}S}{S*I_s*}}\right)+\beta {\eta }_3{T}_d*S*\left(1-{\displaystyle \frac{T_{d}S}{S*T_d*}}\right)+\mu (S*-S)]\\ & +k_1\left(1-{\displaystyle \frac{E}{E*}}\right)[\beta I_a*S*\left({\displaystyle \frac{S{I}_a}{S*I_a*}-{\displaystyle \frac{E}{E*}}}\right)+\beta {\eta }_1{T}_a*S*\left({\displaystyle \frac{T_{a}S}{T_a*S*}-{\displaystyle \frac{E}{E*}}}\right)+\beta {\eta }_2{I}_s*S*\left({\displaystyle \frac{I_{s}S}{I_s*S*}-{\displaystyle \frac{E}{E*}}}\right)\\ & +\beta {\eta }_3{T}_d*S*\left({\displaystyle \frac{T_{d}S}{T_d*S*}-{\displaystyle \frac{E}{E*}}}\right)]+k_2\left(1-{\displaystyle \frac{I_a}{I_a*}}\right)\left[\sigma (1-\tau )E*\left({\displaystyle \frac{E}{E*}-{\displaystyle \frac{I_a}{I_a*}}}\right)\right]\\ & +k_3\left(1-{\displaystyle \frac{I_s}{I_s*}}\right)\left[{\sigma }_{3}E*\left({\displaystyle \frac{E}{E*}-{\displaystyle \frac{I_s}{I_s*}}}\right)+\delta I_a*\left({\displaystyle \frac{I_a}{I_a*}-{\displaystyle \frac{I_s}{I_s*}}}\right)\right]+k_4\left(1-{\displaystyle \frac{T_a}{T_a*}}\right)\left[{\alpha }_1\varphi I_s*\left({\displaystyle \frac{I_s}{I_s*}-{\displaystyle \frac{T_a}{T_a*}}}\right)\right]\\ & -2k_5{T}_d*\left(1-{\displaystyle \frac{T_d}{T_d*}}\right)\left[{\alpha }_2{I}_s*\left({\displaystyle \frac{I_s}{I_s*}-{\displaystyle \frac{T_d}{T_d*}}}\right)\right].\end{array}$$

(9)

Let

$$\begin{array}{rl}x=& \frac{S}{S*},=\frac{E}{E*},z=\frac{I_a}{I_a*},\\ w=& \frac{I_s}{I_s*},u=\frac{T_a}{T_a*},v=\frac{T_d}{T_d*}.\end{array}$$

(10)

Substituting the expressions in (10) into, $\frac{dV}{dt}$ we have:

$${\displaystyle \frac{dV}{dt}=-\mu {\displaystyle \frac{{(S-S*)}^2}{S}+\beta I_s*f(x,y,z,w,u,v),}}$$

where,

$$\begin{array}{rl}f(x,y,z,w,u,v)=& {\psi }_{4}S*(1-xy)+{\eta }_1{\psi }_{2}S*(1-ux)+{\eta }_{2}S*(1-wx)+{\eta }_3{\psi }_{1}S*(1-vx)\\ & +k_1(1-y)[{\psi }_4(xz-y)+{\eta }_1{\psi }_{2}S*(ux-y)+{\eta }_{2}S*(wx-y)\\ & +{\eta }_3{\psi }_{1}S*(vx-y)]+k_2\sigma (1-\tau ){\psi }_3{\psi }_4(1-z)(y-z)\\ & +k_3(1-w)(\sigma {\psi }_3{\psi }_4(y-w)+\delta {\psi }_4(z-w))\\ & +k_4{\alpha }_1\varphi (1-u)(w-u)-2k_5{\alpha }_1(1-v)(w-v)\end{array}$$

We expand and group the coefficients of the same variable, set the terms containing variables with non-negative coefficients to zero to get rid of the positive, non-constant part of f, and solve for $k_1,k_2$ and $k_3$ . We obtain:

$$\begin{array}{rl}{k}_1=& 1,\\ k_2=& {\displaystyle \frac{\beta }{\sigma (1-\tau )}\left({\psi }_4{\psi }_3+{\displaystyle \frac{\sigma {\eta }_1}{Q_3{\psi }_2}}\right)S*,}\\ k_3=& \beta {\displaystyle \frac{\delta {\psi }_4{\eta }_2}{\sigma {\psi }_3}S*,}\\ K_4=& \beta {\displaystyle \frac{{\alpha }_1\varphi }{\sigma {\psi }_1+{\psi }_2},}\\ k_5=& \beta {\displaystyle \frac{{\alpha }_1{\eta }_1+{\eta }_2{\psi }_2}{{\eta }_2{\psi }_3{\psi }_4+\varphi }.}\end{array}$$

Then,

$$\begin{array}{rl}f(x,y,z,w,u,v)=& S*\left(2-xv-{\displaystyle \frac{1}{x}}\right)+{\psi }_1{\eta }_1{\varphi }_{1}S*\left(2-{\displaystyle \frac{xvz}{y}-{\displaystyle \frac{1}{x}}}\right)\\ & +{\eta }_2{\psi }_{2}S*\left(2-{\displaystyle \frac{xwv}{y}-{\displaystyle \frac{1}{x}}}\right)+{\psi }_1\varphi S*\left(2-{\displaystyle \frac{xuv}{y}-{\displaystyle \frac{1}{x}}}\right)\\ & +k_2\sigma \left(1-{\displaystyle \frac{y}{z}}\right)+{\displaystyle \frac{k_3{\psi }_3}{{\beta }_1}\left(1-{\displaystyle \frac{y}{w}-w}\right)+{\displaystyle \frac{k_3{\alpha }_2{\psi }_2}{{\beta }_1}\left(1-{\displaystyle \frac{z}{w}-w}\right).}}\end{array}$$

We need to prove that $\frac{dV}{dt}\le 0$ . We already have the term $-\mu {\left(\frac{S-S*}{s}\right)}^2\le 0$ .

of $\frac{dV}{dt}$ . It is left for us to prove that $f(x,y,z,w,u,v)\le 0.$

Note that the expression $\left(2-xv-\frac{1}{x}\right)$ is less than or equal to zero by the arithmetic-geometric-mean inequality, with equality if and only if x = v = 1. Also, the expressions $\left(2-\frac{xvz}{y}-\frac{1}{x})\right),\left(2-\frac{xwv}{y}-\frac{1}{x}\right)$ and $\left(2-\frac{xuv}{y}-\frac{1}{x}\right)$ are less than or equal to zero by the arithmetic-mean geometric-mean inequality, with equality if and only if $x=y=u=w=v=r=z=1$ . After some tedious algebraic manipulations of replacing the constants $k_1,k_2,k_3,k_4$ and $k_5,$ similar conclusions can be drawn for the remaining expressions in f. Therefore, f is negative and will be equal to zero if $x=y=z=u=w=v=1$ . So V is positive definite at the endemic equilibrium and $\frac{dV}{dt}\le 0$ with equality if and only if $S=S*,E=E*,I_a=I_a*,I_s=I_s*,T_a=T_a*,T_d=T_d*.$ The only invariant set contained in Ω is the set containing only the endemic equilibrium point. This shows that each solution that intersects ${\mathbb{R}}_{+}^6$ limits to the endemic equilibrium. Therefore, by LaSalle's invariance principle, ²⁵ the endemic equilibrium is globally asymptotically stable on $\mathrm{\Omega }$ .

Parameterization

A thorough collection of parameter definitions, baseline values, and literature references can be found in Table 4. Table 4 contains some parameter values that were fitted or approximated using the COVID-19 data for South Africa from February 2021 to March 2022 (see Figure 2), while others were taken from the literature. With an average life expectancy of 62.8 years in South Africa, the demographic parameter, µ, is calculated as follows: $\frac{\pi }{\mu }=1/62.8=\mathrm{0.0145.}$

Figure 2.

Curve fitting of daily confirmed cases of the first wave in South Africa.

Table 4.

Model parameter values.

Parameter	Value/day	Source
π	2623.68	Estimated
β	0.52	Fitted
μ	0.0145	Estimated
σ	0.4	Estimated
γ1	0.54	Estimated
γ2	0.65	Estimated
γ3	0.25	Estimated
δ	0.125	Estimated
α2	0.24	Fitted
ϕ	0.65	Fitted
α1	0.48	Fitted
τ	0.68	Fitted
η1	0.005	Fitted
η2	0.09	Fitted
η3	0.15	Fitted

year. ²⁶ The other demographic parameter, π, is estimated as follows: We assume that $\frac{\pi }{\mu },$ the limiting total human population in the absence of the disease, is 60.14 million, meaning that $\pi =957,643.3$ per year, or 2623.68 per day, given that 60.14 million people were living in South Africa as of 2022,. ²⁷ According to the COVID-19 epidemiological literature,^28–32 the transition rates $\sigma$ and ${\gamma }_3$ are calculated prior to the infectious stages, the virus goes through a dormant stage in which it keeps replicating inside its new host until it reaches a stage when the host can infect others. Since the average latent period $\frac{1}{\sigma }$ is 2.5 days, ²⁸ it follows that $\sigma =0,4.$ Also, with an average infectious period of 5 days,^29,30 we have that $\frac{1}{{\gamma }_3}=5$ so that ${\gamma }_3=\mathrm{0.2.}$ Unlike the individuals in $I_a,$ those in $I_s,$ are symptomatic with symptoms ranging from mild or moderate coronavirus disease to severe cases that can last for a longer period of time. Also, the infectious symptomatic individuals who are denied treatment, in their pursuit of treatment, may take longer to recover, if ever they recover, than those who are given treatment. Therefore, we estimate the transition rates ${\gamma }_1$ and ${\gamma }_2$ such that ${\gamma }_3<{\gamma }_1<{\gamma }_2$ We thus assume that ${\gamma }_1=0.54$ and ${\gamma }_2=\mathrm{0.65.}$ It takes about 2–14 days for pre-symptomatic individuals to become symptomatic, ³³ which gives an average period of 8 days between the two stages. We thus set $\frac{1}{\delta }=8$ and obtain $\delta =0.125$ . The treatment denial rate ${\alpha }_2$ and the remaining parameters are estimated through the model fitting process. The model is fitted to COVID-19 data for South Africa (first wave), and the fitting findings verify that the model estimation and the real data fit each other well. This process involves determining the best values for the parameters by adjusting the parameter values of the model to make its predictions as close as possible to observed data. We use a MATLAB optimization algorithm with the least squares fitting method, which minimizes the difference between the observed data and the model's predicted values. The error function quantifies how much the model's predictions deviate from the observed data. The goal is to adjust the parameter values so that the error function (sum of squared error function) is minimized. The error function is used by the optimization process to determine the values of the parameters that minimize the discrepancy between the output of the model and the observed data. Once the optimization algorithm converges and provides parameter estimates, it is evaluated how well the model fits the data. This is observed in the size of the sum of squared error (SSE) and the fitting curve. The model parameters are estimated to be the parameter values that give the best fit and the smallest SSE.

Sensitivity analysis

Finding the important model parameters that have a significant effect on the fundamental reproduction number R₀ and some of the infectious classes is typically the goal of global sensitivity analysis. It is crucial because some actions pertaining to these critical characteristics have the potential to significantly impact the transmission patterns of the disease under study. Here, a sensitivity analysis is performed using the Partial Rank Correlation Coefficient (PRCC) and Latin Hypercube Sampling (LHS) approaches to determine the effect of related parameters on the fundamental reproduction number. PRCC is a widely used method for figuring out the nonlinear monotonic relationship between two variables. The study shows the PRCC values and related p values, which are used to evaluate the level of uncertainty in the model parameters. The qualitative relationship between the input and output variables is represented by the value and sign of the PRCC. The parameter that has the greatest impact on the output variable is typically thought to have the highest PRCC value and the lowest p value (p < 0.01). Relative to R₀, the four most significant parameters (i.e. those with the four greatest PRCC values) are β, α₁, α₂ and ϕ as shown in the Tornado plot in Figure 3. The contact rate β and the rate of treatment denial α₂ have a positive influence on R₀ (that is, as β and α₂ increases, the value of R₀ increases). R₀ is negatively impacted by the parameters α₁ and ϕ; that is, as these parameters rise, the value of R₀ falls. The most significant parameter with the highest PRCC value is α₂. Its strong positive correlation with R₀ means that the model predicts that increasing this parameter will increase the disease's transmissibility. This is crucial for understanding the potential severity of an outbreak. A positive correlation with R₀ means that as healthcare access decreases (or denial increases), the disease can spread more easily. This could happen due to delayed care, increased viral shedding, or longer infectious periods, leading to more secondary infections. It could also lead to more rapid disease spread and higher mortality, highlighting the importance of maintaining accessible healthcare during a pandemic. The PRCC value of α₂ shows that a small increase in the parameter value leads to a large increase in R₀, this means the disease is highly sensitive to changes in this parameter. In this case, interventions aimed at controlling this parameter could have a significant impact on controlling the spread of the disease. The result of the sensitivity analysis suggests that increasing healthcare access, improving care or preventing healthcare denial could be one of the most effective strategies for reducing transmission. We identify α₁ ϕ and α₂ as the critical parameters because they are not only influential but are also related to the method of intervention in this model, hence can greatly affect the pandemic.

Figure 3.

Tornado plot of model parameters against the model reproduction number showing the correlation coefficient of each parameter relative to R₀.

Numerical simulation

To evaluate the effect of treatment denial on COVID-19 transmission in South Africa, we run numerical simulations using the baseline parameter values listed in Table 4.

Figure 4(a) shows the number of susceptible individuals over time. Initially, the number of susceptible individuals remains relatively constant, as they represent the majority of the population before any substantial COVID-19 transmission occurs. As the disease spreads, a proportion of susceptible individuals move to the exposed compartment. In the baseline simulation, we observe a gradual decline in susceptible individuals as more people become exposed or infected. The rate at which susceptible individuals transition to the exposed state depends on the contact rate and the transmission probability. If healthcare denial exacerbates the disease spread, the number of susceptible individuals could decrease faster as infections increase. Figure 4(b) shows the changes in the population of exposed individuals. This compartment is typically a latent period, during which the virus is incubating. The graph for exposed individuals shows a sharp increase once the infection begins to spread. The exposed population rises as individuals transition from susceptible to exposed due to virus exposure. The rise could be gradual or sharp, depending on the intensity of community transmission and the rate of exposure. The exposure rate is largely driven by the transmission dynamics and interaction rates. If healthcare denial results in increased viral transmission due to delays in accessing care, the number of exposed individuals may rise more quickly. However, the overall duration of exposure will influence the length of time individuals remain in this compartment before moving to the symptomatic or asymptomatic states.

Figure 4.

Time series plots of the susceptible population in (a) and the exposed population in (b).

Figure 5(a) shows the number of asymptomatic individuals. These individuals can still transmit the virus to others, albeit at a lower rate than symptomatic individuals. Initially, the number of asymptomatic individuals might be low, with a steady increase as more individuals become infected. However, the growth of this population will depend on factors such as the infectivity of asymptomatic individuals and the rate of progression from exposed to infectious asymptomatic status. Healthcare denial can affect the number of asymptomatic individuals because a lack of testing or care may result in a delay in identifying and isolating asymptomatic cases, thus allowing more individuals to unknowingly spread the virus. In a baseline scenario with no mitigation, this group is likely to grow steadily over time as the disease continues to circulate.

Figure 5.

Time series plots of the infected asymptomatic population in (a) and the symptomatic population in (b).

On the other hand, Figure 5(b) shows the number of symptomatic individuals over time. These individuals are often more contagious than asymptomatic ones and tend to seek medical treatment. The number of symptomatic individuals will initially be small but will increase as more people are exposed to the virus and develop symptoms. In a baseline scenario with healthcare denial, a higher proportion of symptomatic individuals may experience delays in receiving treatment, which can result in longer infectious periods and higher mortality rates. Healthcare denial could lead to more individuals remaining untreated and possibly deteriorating into more severe forms of the disease. This could extend the duration of the symptomatic phase and increase the infectivity of individuals in this category, thus accelerating disease spread.

Figure 6(a) tracks the number of individuals who receive treatment for their COVID-19 infection, whether they are symptomatic or asymptomatic. These individuals are treated in public healthcare facilities. The number of individuals receiving treatment should steadily rise as the disease spreads, but treatment capacity limits the number of individuals who can access care. In the baseline scenario, overwhelming demand and healthcare system strain may increase the number of individuals who are adequately treated. The capacity of healthcare facilities, the availability of resources (e.g. ventilators, beds, healthcare workers), and treatment access policies all determine how quickly individuals move from the symptomatic compartment to the treatment compartment. Healthcare denial could slow the rate at which people receive treatment, leading to higher mortality and prolonged disease dynamics.

Figure 6.

Time series plots of population of individuals given treated (a), those denied treatment (b) and the recovered individuals in (c).

Figure 6(b) represents individuals who are symptomatic but are denied treatment either due to overcrowded hospitals, lack of healthcare access or policy restrictions. These individuals may continue to infect others and may face worse health outcomes. The number of individuals denied treatment increases over time as the healthcare system becomes more overwhelmed by COVID-19 cases. These individuals may experience worse health outcomes and can contribute to increased mortality rates. Treatment denial may also extend the duration of infection, allowing for greater viral spread. The number of individuals who have recovered from COVID-19, either through natural immunity or successful treatment, is shown in Figure 6(c). The number of recovered individuals will rise over time, with a steady increase as people recover from infection. However, if healthcare systems are overwhelmed and treatment is denied, the recovery rate may be slower, or individuals may experience longer recovery times.

To see the associated consequences, we vary the critical parameters, which are α₁, α₂ and ϕ. We first examine the variations of COVID-19 progression with changes in the treatment denial rate α₂ and the treatment acceptance rate α₁, as presented in Figure 7. The other parameters remain the same as the baseline values in Table 4. The graph in Figure 7(b) shows how the number of infected individuals evolve over time in response to varying rates of treatment denial. In this simulation, the treatment denial rate refers to the proportion of symptomatic individuals who do not receive medical care due to overcrowded healthcare systems, policy restrictions or insufficient healthcare resources. As the treatment denial rate increases, we see a significant rise in the number of infected individuals over time. This is because individuals who are denied treatment tend to remain infectious for longer periods, both asymptomatically and symptomatically, thereby increasing the opportunity for further transmission within the community. Additionally, untreated symptomatic individuals may experience more severe disease progression, increasing morbidity and potentially leading to higher mortality, which further exacerbates transmission dynamics. Also, as the denial of treatment increases, infected individuals who do not receive care may have longer infectious periods, thus increasing the chances of spreading the virus. This will likely cause the number of infections to continue rising even in the face of existing public health efforts to control the spread. The rate of treatment denial directly impacts how quickly infected individuals progress to either recovery or severe disease. When treatment is denied, the likelihood of severe complications increases, which may not only impact individual health outcomes but also slow down the overall recovery rate in the population. Hence, delayed treatment results in higher overall infections due to prolonged transmission.

Figure 7.

Graphs of α₁ variations in (a) and α₂ variations in (b).

Key Takeaway: Increasing the rate of treatment denial contributes to a sharp rise in the number of infected individuals, particularly in the symptomatic category, as more people remain untreated and continue to spread the virus for longer periods.

Figure 7(a) illustrates the relationship between the rate of treatment acceptance and the number of infected individuals over time. As the treatment acceptance rate increases, we observe a sharp reduction in the number of infected individuals, particularly in the symptomatic category. With more individuals receiving prompt treatment, viral load in the body can be reduced more quickly, leading to shorter infectious periods and lower transmission rates. Additionally, more people will recover more rapidly, thus decreasing the overall pool of infectious individuals. When treatment is readily available, infected individuals who receive appropriate care will be more likely to recover faster, reducing the duration of their infectious period. Consequently, this leads to fewer opportunities for the virus to spread, which results in a decrease in overall infections. The availability of healthcare resources, such as ICU beds, ventilators, medications and healthcare staff, is critical to ensuring that individuals who are symptomatic can be treated effectively. The efficiency of the healthcare system in treating patients will influence how quickly the infected individuals recover and stop being a source of transmission. A high treatment acceptance rate also supports the control of severe outcomes, thus lowering both hospitalization rates and COVID19-related mortality.

Key Takeaway: As the rate of treatment acceptance increases, the number of infected individuals will decrease significantly. This is due to faster recovery times, reduced viral transmission and an overall improvement in public health outcomes.

As observed on the graph, when the denial rate reaches 80%, the peak number of infected individuals will be higher than when the denial rate is at the baseline value. The higher the rate of treatment denial, the more the creation of new infections. The reverse is the case when more infected individuals are given treatment. When the treatment acceptance rate reaches 85%, the peak number of infectives drops significantly compared to when the acceptance rate is at its baseline value. This may be due to the fact that infected individuals under treatment are kept in confined and controlled environments with very limited contact with susceptible individuals. Hence, the creation of new infections due to physical contact is limited.

Figure 8 illustrates an investigation of the impact of treatment denial on COVID-19 infection. While α₂ fluctuates, the other parameters stay the same as the baseline values in Table 4. It shows that the number of daily cases and epidemic magnitude both rise in tandem with the treatment denial rate α₂. The number of new infections significantly declines when the value of α₂ drops as shown by the red and blue curves in Figure 8. With the treatment denial rate reaching 5% and below, the curve quickly reaches its peak and flattens quickly as well. The peak number of cases is also lower compared to when the treatment denial rate is high. Decreasing the treatment denial rate, therefore, can largely decrease the creation of new infections and general disease burden. Generally speaking, Denying certain individuals treatment against COVID-19 can greatly affect the spread of the disease. The shaded region in Figure 8 represents the impact of treatment denial.

Figure 8.

The impact of treatment denial.

Figure 9(a) shows that as the treatment acceptance rate increases (i.e. a greater proportion of symptomatic individuals receive treatment), the basic reproduction number R₀ decreases. This is because treating symptomatic individuals leads to shorter infectious periods, reduced viral load, and consequently less opportunity for transmission to others. As more individuals receive care, the virus has fewer opportunities to spread, thus reducing R₀. When a higher proportion of symptomatic individuals are accepted for treatment, the healthcare system is more efficient at managing the disease, reducing the infectiousness of individuals over time. The basic reproduction number decreases as treatment is provided more widely and effectively, leading to fewer new cases generated by infected individuals. When either the treatment acceptance rate or the proportion accepted for treatment is low, symptomatic individuals either remain untreated or are not treated effectively. This leads to longer infectious periods and higher transmission potential, which in turn causes the basic reproduction number to rise. In such cases, the healthcare system is less effective in curbing the virus's spread, resulting in exponential growth of cases. Key Takeaways: Increasing treatment acceptance and proportion of individuals accepted for treatment consistently leads to lower R₀, which is indicative of a slower transmission rate.

Figure 9.

Contour plots of α₁ in (a) and α₂ in (b).

On the other hand, Figure 9(b) illustrates how changes in both the rate of treatment denial and the effective contact rate impact the basic reproduction number R₀ of COVID-19. The basic reproduction number is a critical epidemiological metric that measures the potential of a disease to spread in a susceptible population, with values above 1 indicating the potential for sustained outbreaks. Understanding the interaction between treatment denial and the effective contact rate is vital to evaluating the overall impact of healthcare denial on disease transmission dynamics. By visualizing these relationships, we can better understand how barriers to treatment and social behaviors interact to influence the potential for disease spread. When the treatment denial rate is high, a larger proportion of symptomatic individuals go untreated. This results in longer infectious periods, allowing individuals to spread the virus to more people. As treatment is withheld, the population has fewer opportunities to recover or reduce their infectiousness, leading to an increase in the basic reproduction number R₀. Essentially, high treatment denial enables the virus to circulate freely in the population for longer periods. When the effective contact rate is high, it indicates that individuals come into contact with others more frequently, thereby increasing the likelihood of transmission. This, in turn, raises the basic reproduction number R₀. A high contact rate, combined with high treatment denial, creates a situation where infection spreads more rapidly because there are both more opportunities for transmission and fewer opportunities for infected individuals to recover or reduce their infectiousness through treatment. When the rate of treatment denial and the effective contact rate are low, the spread of the virus can be controlled more effectively. Low treatment denial means that more symptomatic individuals receive treatment, leading to shorter infectious periods and a lower R₀. Low contact rates can also reduce the number of new infections generated, lowering the basic reproduction number R₀. Key Takeaways: Increased treatment denial and increased effective contact rates both contribute to a higher basic reproduction number R₀, suggesting that denial of treatment, coupled with frequent contact between individuals, will lead to faster disease transmission and larger outbreaks. Conversely, reducing treatment denial and lowering the effective contact rate can significantly reduce the R₀ and slow the spread of the disease.

Effects of treatment denial on disease burden

Outbreaks of infectious diseases often overwhelm hospitals, as seen in South Africa during COVID-19, where overcrowding reduced hospitals’ effectiveness.^34,35 With limited capacity, hospitals must either continue admitting patients, risking overcrowding and reduced care quality, or deny treatment to some,^36,37 which raises infection rates in the community. Both choices involve trade-offs: overcrowding increases in-hospital deaths and infections, while denial fuels community spread. A time-series analysis compared two scenarios: (ϕ = 1) all patients treated, see Figure 10(b), and (ϕ = 0.65) some denied treatment, see Figure 10(a). Results showed higher daily cases and fewer recoveries when treatment was denied. Despite compromised care in overcrowded settings, treating all patients remains the better strategy. The basic reproduction number was lower when all were treated (R₀ = 1.5) than when some were denied care (R₀ = 2.6), indicating higher disease burden and treatment costs under denial.

Figure 10.

Graphs of daily new cases for the two scenarios in (a) graphs of recovered population for the two scenarios in (b).

The results show that, although the health care quality will be compromised due to overcrowding, accepting to treat every infected individual is still a much better strategy than denying some them treatment. The basic reproduction number when ϕ = 1 and when ϕ = 0.065 are approximated at 1.5 and 2.6, respectively. This is indicative of a greater disease burden and, consequently, a higher treatment cost when some of the infected individuals are denied treatment compared to when all infected individuals are given treatment.

Implications of the results to public health policy in South Africa

Health care denial, which refers to situations where individuals or populations are unable to access needed health services, whether due to financial, logistical, systemic or policy-related barriers, has significant and far-reaching implications for public health policies. The denial of health care affects not only the individuals directly impacted but also the overall health of a population, and, consequently, the functioning of the healthcare system and economy of a country. Our findings indicate that public health care denial leads to

• Increased Disease Transmission and Burden on Public Health Systems. When individuals are denied or delayed access to care, they may continue to be infectious for longer periods, increasing the likelihood of disease transmission.

Policy implication: Countries need to implement policies ensuring universal access to healthcare, especially during outbreaks, to prevent rapid disease transmission. Failure to do so can result in large-scale public health crises.

• Worse health outcomes, including increased mortality. Denied or delayed treatment for serious diseases can lead to preventable deaths.

Policy implication: Policies promoting early intervention, accessibility to critical care, and affordable healthcare can significantly reduce preventable deaths and improve the overall health of the population.

• When healthcare access is denied or limited, people who would otherwise have received preventive care, early treatment, or routine health services may delay seeking help until their conditions worsen. This leads to increased demand for urgent and emergency care, placing a burden on healthcare systems.

Policy Implication: Policymakers should ensure that health systems are adequately funded and well-resourced to handle both routine and emergency healthcare needs. Additionally, policies aimed at preventive care, health education, and early intervention can reduce the strain on hospitals by decreasing the need for urgent care.

• In the event of a pandemic or epidemic, health care denial has direct consequences for controlling the spread of infectious diseases. People who are denied care are less likely to seek timely testing, treatment, or quarantine, which leads to increased disease transmission. Policy Implication: During pandemics, policies should prioritize equitable access to healthcare, including testing, treatment, and preventive measures like vaccines. Governments must ensure that public health interventions are available to everyone, especially the most vulnerable populations.

The denial of healthcare access is a critical issue that has far-reaching consequences for both individuals and society at large. Its impact on public health policies is profound and multifaceted, touching on aspects such as disease transmission, health equity, healthcare system capacity, economic outcomes, and social stability. For effective public health policy, countries must prioritize the expansion of universal health coverage (UHC), ensure accessibility to quality healthcare for all, and implement strategies to reduce health disparities. This is especially important in times of crisis, such as during pandemics, when timely and equitable access to healthcare is essential for controlling the spread of infectious diseases and minimizing their societal and economic impacts.

Conclusion

This analysis provided key evidence on the effects of treatment denial on the COVID-19 disease burden in South Africa. Our model parameters and specifications were estimated based on the population structure of South Africa, which has a median age of 27.6 years (compared with 38 years in the USA) and thus are likely to reflect epidemic scenarios in low- and medium-income nations with age distributions that are equally youthful. We choose the two scenarios studied in this work (where every infected person is given treatment and the scenario where some are denied treatment) based on their relevance to settings with limitations in formal healthcare infrastructure.

To assess the impact of treatment denial, especially in the regions where treatment denial is very rampant, we simulated the model with variations of the treatment denial rate and treatment acceptance rate. An overall increase in the peak value is observed for the infected population as the treatment denial rate crosses 80%, and the reverse result is observed when the treatment acceptance rate crosses 80%. Also, an increase in the size of the outbreak is observed in Figure 8 as the treatment denial rate increased from 5% to 40%. The epidemic curve quickly flattens as the treatment denial rate falls below 5%. Although the healthcare quality will depreciate due to overcrowding, where every infected individual is given treatment, the results from this study point to the fact that treatment must be made accessible to everyone if the epidemic is to be controlled. The two scenarios (ϕ = 1 and ϕ = 0.65) were further analyzed by projecting the number of daily new cases and the number of recovered individuals in each scenario. An increase in the basic reproduction number from 1.5 when ϕ = 1 to 2.6 when ϕ = 0.65 indicated a greater disease burden and a consequent increase in the size of the outbreak when some individuals are denied treatment compared to when every infected individual is given treatment. A fall in the peak of the epidemic curve and an increase in the population of recovered individuals is also observed when ϕ = 1 compared to ϕ = 0.65. Therefore, treatment denial is not good for the control and management of the COVID-19 epidemic and must be abolished if disease containment is attained.

The impact of treatment denial on COVID-19 dynamics can vary significantly between South Africa and other low-income and resource-poor countries, given differences in healthcare infrastructure, government policies, disease burden, and the availability of resources.

Model limitations and potential extensions for future work

While the model in this work provides valuable insights, it is important to recognize its limitations. These limitations stem from the assumptions made, the data available, and the complexity of real-world dynamics that might not be fully captured by the model. Here are some potential limitations of the model:

Simplifying Assumptions about Health Care Denial Homogeneity of Denial Impact: The model assumes that the impact of healthcare denial on individuals is uniform across the population, which may not be the case. For example, vulnerable groups such as the elderly, individuals with comorbidities, or marginalized populations (e.g. low-income groups, ethnic minorities) may experience more severe outcomes from treatment denial compared to others. These heterogeneous effects might not be fully captured.

The model was intentionally designed with certain limitations in mind, acknowledging that it does not capture all potential co-morbidities and conditions that could affect disease outcomes (e.g. HIV, diabetes, hypertension). The decision not to incorporate HIV coinfection into the model, despite its high prevalence in South Africa, is likely due to a combination of factors, including model focus and complexity, data availability, computational burden and the intention to provide a clearer understanding of healthcare denial's impact on COVID-19 transmission dynamics. While HIV co-infection is a significant factor in understanding the severity of COVID-19, the authors have chosen to prioritize simplicity and clarity in the current model, with plans to address this limitation in future research.

Granularity of Denial: The model treats healthcare denial as a binary variable (either treatment is provided or not), but in reality, denial could occur in varying degrees or be context-dependent. For example, partial treatment, delayed care or restrictions in access due to geographical or financial factors could all impact disease transmission differently.

Data Availability and Calibration Lack of Comprehensive Data: Accurate model calibration relies on comprehensive data on healthcare denial, transmission dynamics and individual health outcomes. During a pandemic like COVID-19, data on healthcare denial is often sparse, especially in low-resource settings like South Africa. This can lead to inaccurate model predictions if the data used to calibrate the model is incomplete or of poor quality.

Limited Historical Data: COVID-19 is a relatively new disease, and there may be limited historical data on how healthcare denial directly impacts the transmission dynamics of such a disease. As a result, the model may rely on approximations or assumptions based on previous pandemics, which may not perfectly reflect COVID-19's specific characteristics.

Over-Simplified Disease Transmission Dynamics Assumptions about Transmission Rates: Disease transmission is influenced by numerous factors, such as population density, behavioral factors (e.g. social distancing, mask-wearing) and public health interventions (e.g. lockdowns, vaccination). The model did not fully account for these additional factors, potentially leading to unrealistic estimates of transmission dynamics when healthcare denial is introduced.

Since treatment denial only takes place in public health facilities in South Africa, it just takes into account these facilities. Given that the high prevalence of HIV and non-communicable diseases in South Africa did impair COVID-19 outcomes when compared to resource-rich countries, we cautiously hypothesized that HIV had no moderating influence on the severity of the pandemic in that country.

Further work that can be done on the model is the evaluation of the cost implications of denying immigrants treatment given that those who would have been denied treatment have the potential of spreading the disease to the non-immigrant population. The model could be extended to include dynamic healthcare system capacities that change in response to the epidemic's progression. This would involve modeling hospital beds, intensive care units (ICUs), medical personnel and treatment availability as limited resources that fluctuate based on infection rates, healthcare interventions and resource allocation policies.

Future work could focus on optimal allocation of healthcare resources during a pandemic. For example, simulating triage protocols or the impact of prioritizing treatment for high-risk individuals could reveal important policy insights.

Extending the model to consider how medical supply shortages (e.g. ventilators, medications) influence treatment denial and disease outcomes could improve the realism of the model. A potential extension is to model how socioeconomic factors (e.g. income, education, access to insurance, living conditions) influence treatment denial and disease outcomes. This could include: Income and Insurance Status: Modeling how people with different income levels or insurance status experience varying degrees of healthcare access.

The model could be expanded to simulate how individuals change their behavior in response to healthcare denial. For example: Self-Treatment or Alternative Care: People denied care might resort to self-medication, alternative health practices or delayed care, which could influence disease progression and transmission rates.