Scientific Papers

A COVID-19 epidemic model with periodicity in transmission and environmental dynamics

1. Introduction

The coronavirus disease (COVID-19) pandemic is now a worldwide epidemic that has been rapidly growing from the onset and is caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) which has a significant morbidity and mortality estimate of 0.5–2% of confirmed cases [1]. The world experienced its first zoonotic human coronavirus, in 2002, the severe acute respiratory syndrome coronavirus (SARS-CoV) spread to 37 countries in 2012, and the second is the Middle East respiratory syndrome coronavirus (MERS-CoV), also spread to 27 countries [2]. COVID-19 is currently recorded in 230 countries as of November 2022, with more than 638 million cases and 6.8 million deaths. South Africa was leading the African continent with over 4 million cases and 102, 000 thousand deaths.

Symptoms of COVID-19 infection include breathing difficulty, fever, fatigue, a dry cough, and in severe cases, bilateral lung infiltration and others also developed non-respiratory symptoms, including vomiting, diarrhea, and nausea [24] similar to the symptoms of SARS-CoV and MERS-CoV.

According to WHO [5], the COVID-19 virus is primarily transmitted between people through respiratory droplets and contact routes. Droplet transmission occurs when an infected person sings, breathes, sneezed, coughed, and has nasal discharge [6]. Transmission can also occur through fomites in the immediate environment around the infected person [7]. Therefore, the COVID-19 virus can be transmitted by direct contact with infected people and indirect contact with surfaces in the immediate environment or with objects used on or by the infected person [7].

In the analysis done by Kampf et al. [8] on 22 types of coronaviruses, it revealed that human coronaviruses such as severe acute respiratory syndrome (SARS) coronavirus, Middle East respiratory syndrome (MERS) coronavirus or endemic human coronaviruses (HCoV) could persist on inanimate surfaces such as metal, glass, or plastic for up to 9 days, and this gives evidence of the survival of the pathogens in the environment to increase risk of infection again. In Gralinski and Menachery [3], it was reported that samples taken from the Huanan Seafood Market, a live animal and seafood wholesale market in Wuhan, were positive for COVID-19, which suggested that pathogens could be transmitted through the environmental reservoir [9, 10]. COVID-19 has also been found in the stool of some infected individuals, which may contaminate the aquatic environment [2, 11].

Evidence suggested that infected humans continued to shed virus into the environment as long as they were infectious [12]. The European Congress of Clinical Microbiology and Infectious Diseases [13] reported a COVID-19 patient who tested positive for 505 days until death. Furthermore, a report by Spanish researchers [14] described a 52-year-old man who shed the virus after 189 days of chemotherapy and a 64-year-old man also continued to shed the virus for 169 days after being infected [15]. Rahmani et al. [16] examined the length of time in an infected person with SARS-CoV-2 continued to shed the virus. It was discovered that the typical person continues to shed the virus for a month. However, some individuals’ bodies continued to discharge the virus for a long time.

López and Rodo [17], Liu et al. [18], and Wang et al. [19] reported that control measures were important factors in containing epidemics and reducing the transmission of the disease. During the early stages, different levels of control measures had different impacts on the transmission of COVID-19 [2022]. Strict controls were by closing public facilities, and some mild control measures included maintaining social distancing, taking temperature measurements, and wearing masks. All these control measures were implemented from the first wave to the fourth wave of COVID-19 in South Africa since there were no pharmaceutical inventions. However, during those times, most developed mathematical models did not consider vaccinations [9, 2325]. Recent mathematical models that have been developed include vaccination as a major control measure [26, 27].

Regardless of these control measures, seasonality is another important factor that influences epidemics. There has been much controversy regarding COVID-19, with questions over its transmission regarding seasonal patterns including other seasonal epidemics such as flu [9, 2830]. Some studies have attempted to establish the relationship between COVID-19 and seasonality by varying meteorological factors [9]. In a study by Matson et al. [31] and Liu et al. [18], it pointed out that the stability of the virus in the air or on surfaces is quick to respond to environmental conditions, including humidity, temperature, sunlight, and more. The virus is more stable at low-temperature and low-humidity conditions, whereas warmer temperatures and higher humidity give it a half-life. Moreover, experimental data and regional analysis by Yao et al. [32] showed that COVID-19 in a high-temperature environment had a lower survival rate and infections declined in summer. At 4°, 22°, and 37°, SARS-CoV-2 can persist for 14, 7, and 1 day, respectively, under laboratory conditions. However, at 56°, it persistence reduce drastically to 10 min [18, 33]. In addition, Chin et al. [33] found that at 4° the virus was stable and inactive when increased to 70°. Confirmed cases of COVID-19 were found to be concentrated at an absolute humidity of 3 to 10g/m3 at air temperatures of 5 to 15° in Huang et al. [34]. The following study and experiment also confirmed with similar findings, SARS-CoV-2 can remain viable and infectious in aerosols for hours at room temperature 21°—23° and a fixed relative humidity of 65% [35]. Simulated sunlight can rapidly inactivate SARS-CoV-2 suspended in either simulated saliva or culture media and dried on stainless steel coupons [18, 36].

Many mathematical models have been developed since the appearance of COVID-19 disease. Most of these models are based on the significant role of human-to-human transmission as done by Chan et al. [37], Ojo et al. [23], and Babasola et al. [24]. However, several studies have given evidence of pathogen transmission through the environmental reservoir. Few mathematical models have been developed so far to consider the role of the human-to-environment transmission [9, 10, 38]. However, seasonal patterns expose the limitations of many recent COVID-19 models that do not incorporate seasonality. In this study, we consider South Africa, with two main climate seasons in a year. During these periods, there has always been a surge in infections. We present a mathematical model for COVID-19 that considers double periodic transmission pathways: human-to-human periodic transmission and periodic transmission through the environmental reservoir (human-to-environment). We divide our human population into five sub-populations: susceptible individuals, exposed individuals, asymptomatic individuals, symptomatic individuals, and recovered individuals. The transmission rate incorporates the influence of seasonality, control measures (vaccination), and the environment. We analyzed the basic reproduction number (using the next infection operator) and establish that it is a sharp threshold for the COVID-19 models with the periodic transmission. The analysis method for extinction and persistence results for periodic epidemic systems is inspired by the research done in Wang and Zhao [2]. We fit our mathematical model to data obtained from the National Institute for Communicable Diseases, South Africa [39], to estimate parameters and then simulate the efficacy driven by the vaccination parameter.

The study is organized as follows. In Section 2, we present the model and its properties. In Section 3, a detailed mathematical analysis is performed on the model developed. In Section 4, we explain numerical simulations, and conclude the study with some discussion and recommendations in Section 5.

2. Model formulation and properties

2.1. Model formulation

We define N(t) as the total human population at time t. The human population is divided into five classes: the susceptible class S(t), the exposed class E(t); those who have been exposed to the virus but have not yet been diagnosed as COVID-19 positive, the asymptomatic class Ia(t); people who have been diagnosed with COVID-19 but do not exhibit any clinical symptoms, the symptomatic individuals Is(t); people with clinical symptoms who have been diagnosed as COVID-19 positive and recovered class R(t); and people who have recovered from the virus. The human population has a natural mortality rate of μ. The total population is given as

N(t)=S(t)+E(t)+Ia(t)+Is(t)+R(t).    (1)

Ce(t) is the pathogen concentration in the environment at a given time (t). The removal rate of the virus is σ. The environmental reservoir is defined in this study as fomites in the immediate environment around the infected person and sheddings by infected humans either through droplets or feces. The relationship between humans and the environment is represented in Figure 1.

Figure 1. Schematic diagram for the presented model. The human population is divided into susceptible class (S), exposed class (E), asymptomatic class (Ia), symptomatic class (Is) and the recovered class (R). Thus, the population at any time t is N = S + E + Ia + Is + R. In the environmental population, Ce is the environmental pathogens. We have used the dashed lines with arrowheads to indicate contact and solid lines with arrowheads to indicate movement.

The dynamics of the susceptible population is given in Equation (2), Λ being the recruitment rate. The second term δ(t, Ia, Is, Ce) determines the rate of new infection, the parameter δ is the disease transmission rate and is time dependent.

dSdt=Λδ(t,Ia,Is,Ce)SμS,    (2)


δ(t,Ia,Is,Ce)=β(t)IaN+β(t)ηIsN+β(t)CeK50+Ce and β(t)={β^+β¯β^sin(2πtω),t0tt2(β^+β¯β^sin(2πtω))m(t),tt2.    (3)

The parameter β(t) is the effective contact rate of susceptible humans with the asymptomatic class Ia, symptomatic class Is, and the environment Ce, respectively. K50 gives the concentration of virus in the environment that yield 50% chance of infection with COVID-19. η measures the relative infectivity between Ia and Is, and it is measured between 0 ≤ η ≤ 1.

β^β̄ is the amplitude of the periodic oscillations when β̄=0 means there are no infections in the periodic function. Here, β^ is the baseline value or the time average, ω=3652=182.5 days since there are only two waves in a year and a time-varying vaccination parameter m(t), to formulate the influence of control measures. It gives the quantitative estimation of the control measure implemented in South Africa. We define our m(t)=emxtx, which falls between (0, 1], when m(t) = 1 then the epidemic spreads without any restrictions; however, m(t) can not be assumed to be one except when there are no any control measures, and when m(t) = 0, then there is no epidemic. mx controls the speed at which the vaccination reaches its maximum or minimum values, and we call it the vaccination efficacy in our study. Moreover, tx gives the start time of the vaccination.

For the exposed population, in Equation (4), the first term represents those who enter from the susceptible pool driven by the force of infection δ. The rate κ is the progression rate of exposed individuals to the asymptomatic and symptomatic classes.

dEdt=δ(t,Ia,Is,Ce)S(μ+κ)E.    (4)

dIadt=κpE(μ+α+γa)Ia,    (5)

The dynamics of the symptomatic individuals infected with COVID-19 is given in Equation (5). In the first term, p gives the proportion of humans that are moving from the exposed class to the asymptomatic class, α is the rate of transfer from the asymptomatic to the symptomatic class, and γa is the recovery rate for the asymptomatic class respectively.

The dynamics of the asymptomatic individuals infected with COVID-19 is given in Equation (6). In the third term, ψ is the disease-induced mortality rate and γs is the rate of recovery for symptomatic individuals.

dIsdt=κ(1p)E+αIa(μ+ψ+γs)Is,    (6)

The dynamics of the recovered individual are represented in Equation (7),

dRdt=γsIs+γaIaμR.    (7)

dCedt=ϵa(t)Ia+ϵs(t)IsσCe.    (8)

The environmental dynamics is given in Equation (8), where ϵa(t) is the shedding rate from the asymptomatic class, which is time-dependent, ϵs(t) is the shedding rate from the symptomatic class which is time-dependent, and the rate σ is is the removal of the virus from the environment, respectively. Here, ϵa(t) and ϵs(t) are given by

ϵa,s(t)=ϵ^a,s[1+ϵ¯a,ssin(2πtω)].    (9)

The amplitude of the periodic oscillations in ϵa, s(t) is given by ϵ^a,sϵ̄a,s and when ϵ̄a,s=0, no infections in the periodic function and baseline values or the time-average is ϵ^a,s. ω=3652=182.5 days since we have only two waves in a year.

We re-scale Equations (2)–(8), using the following substitutions:

s=SN,e=EN,ia=IaN,is=IsN,r=RN,x=CeKCe, and K˜=K50KCe.

The virus in the environment has a carrying capacity of Kce, Equations (2)–(8) changes into the following subsystems

dsdt=Λ˜δ˜(t,ia,is,x)sμsdedt=δ˜(t,ia,is,x)s(μ+κ)e,diadt=κpe(μ+α+γa)ia,disdt=κ(1p)e+αia(μ+ψ+γs)is,drdt=γaia+γsisμr,dxdt=ϵ˜a(t)ia+ϵ˜s(t)isσx,}    (10)


δ˜(t,ia,is,x)=β˜(t)ia+ηβ˜(t)is+β˜(t)xK˜+x and ϵ˜=ϵNKce.

We remove the fifth equation in Equation (10) since it is independent of the rest of the equations and rewrite the system as follows:

dsdt=Λ˜δ˜(t,ia,is,x)sμsdedt=δ˜(t,ia,is,x)s(μ+κ)e,diadt=κpe(μ+α+γa)ia,disdt=κ(1p)e+αia(μ+ψ+γs)is,dxdt=ϵ˜a(t)ia+ϵ˜s(t)isσx,}    (11)

2.2. Model properties

The human population and environmental reservoir are differentiable and periodic in time with a common period ω.

δ(t+ω,ia,is,x)=δ(t,ia,is,x) and ϵ˜a(t)ia+ϵ˜s(t)isσx=ζ(t,ia,is,x).

To make biological sense, we assume that the functions δ and ζ satisfy the following conditions for all t ≥0:

(P1) δ(t, 0, 0, 0) = ζ(t, 0, 0, 0) = 0, the disease-free equilibrium is unique and constant. Let E0 denote the disease-free equilibrium, which is given as

E0=(s,e,ia,is,x)=(Λ˜μ,0,0,0,0)    (12)

(P2) The force of infection, δ(t, ia, is, x) ≥ 0. This property ensures we have a non-negative for all time.

β˜(t)ia+ηβ˜(t)is+β˜(t)xK˜+x0, for all t0.


δ(t,ia,is,x)ia=β˜(t)0,δ(t,ia,is,x)is=β˜(t)0,δ(t,ia,is,x)x=β˜(t)(K˜+x)β˜(t)x(K˜+x)20.}    (13)

ζ(t,ia,is,x)ia=ϵ˜a(t)0,ζ(t,ia,is,x)is=ϵ˜s(t)0,}    (14)

ζ(t,ia,is,x)x=σ0.    (15)

The first set of equations in Property (13) ensures that the rate of new infection increases with both the infected human population and the pathogen concentration. The second set of equations in Property (14) ensures an increase in human infections and, therefore, a higher level of human contribution to the environmental virus, leading to a higher growth rate for the pathogens. The property (15) ensures that the rate of change of the pathogen concentration would be negatively related to itself [40].

(P4) δ(t, ia, is, x) and ζ(t, ia, is, x) are both concave for any t ≥ 0. This property ensures that the second partial derivative is concave.

D2δ=[2δia22δiais2δiax2δisia2δis22δisx2δxia2δxis2δx2]0 and D2ζ=[2ζia22ζiais2ζiax2ζisia2ζis22ζisx2ζxia2ζxis2ζx2]0,

and is negative semidefinite everywhere. This is a common assumption for non-linear incidence [41, 42]. In our model, this condition regulates δ(t, ia, is, x) as a biologically realistic incidence based on a consequence of saturation effects: when the number of the infective or the environmental pathogens concentration is high, the incidence rate will respond more slowly than linearly to the increase in ia, is, and x. Similar arguments hold for ζ(t, ia, is, x).

(P5) δ(t, 0, 0, x) > 0 if x > 0; ζ(t, ia, is, 0) > 0 if ia, is > 0.

This property implies that infection can begin solely through indirect transmission. Positive pathogen concentration can lead to a positive incidence even if ia, is = 0 at the beginning. Infected people will contribute to the increase of the virus in the environment by shedding even if x = 0.

3. Basic reproduction number and analysis

From Equation (11), the disease-free equilibrium (E0) is given as


By using Wang and Zhao’s [2] approach, we obtain

FE0(t)=[0β˜(t)ηβ˜(t)β˜(t)K˜000000000000] andVE0(t)=[w000κρx00yαz00ϵa(t)ϵs(t)σ]


w=κ(1p),x=(μ+κ),y=(μ+α+γa), and z=(μ+ψ+γs).

We denote Y(t, s) as the evolution operator of the linear ω−periodic system

dydt=V(E0)y.    (16)

For each s ∈ ℝ, the 4 × 4 matrix Y(t, s) satisfies Equation 17 I is a 4 × 4 identity matrix:

dY(t,s)dt=V(E0)Y(t,s),ts, and Y(s,s)=I    (17)

Thus, the monodromy matrix ΦV(E0) of Equation (16) yields Y(t, 0), t ≥ 0. Given a periodic environment, let the initial distribution of infectious individuals be ϕ(s) ∈ Cω, the rate of new infections caused by the infected humans who were introduced at time s be F(s)ϕ(s), and Y(t, s)F(s)ϕ(s) represent the distribution of those infected humans who were newly infected at the time s and remains in the infected compartments at time t for ts.

Let Cω be the ordered Banach space of all ω-periodic functions from ℝ to ℝ that possesses the maximal norm ||.|| and the positive cone


As a result, the distribution of cumulative new infections caused by all of the diseased people who were introduced at the previous time t is given by


The framework created by Zhang and Zhao [43] was expanded upon by Wang and Zhao [2] to incorporate epidemiological models in periodic environments. Let L : CωCω, the next infection operator is donated by L, and R0 := ρ(L) is the basic reproduction ratio, where L is the spectral radius. To solve R0 numerically, see Wang and Zhao [2] and Assan et al. [44, 45].

Lemma 1. [43] The following statements are valid for model (Section 2.1):

(i) R0 = 1 if and only if ρFV(ω)) = 1.

(ii) R0 > 1 if and only if ρFV(ω)) > 1.

(iii) R0 < 1 if and only if ρFV(ω)) < 1.

3.1. The stability of the disease-free equilibrium and existence of periodic solutions

3.1.1. The stability of disease-free equilibrium

Theorem 1. For system 11, if R0 < 1, then disease-free equilibrium E0 is globally asymptotically stable.

Proof. Let (s, e, ia, is, x) be a non-negative solution to system (11). To complete the proof, it is sufficient to show this non-negative solution tends to E0 as t → +∞. The first equation of system (11) with 0 ≤ ns gives

s˙Λ˜μs    (18)

Hence, for any ε > 0, when tε > 0; when t > tε, we have ss0 + ε, and gives the following inequalities

dedt=δ(t,ia,is,x)(s0+ε)(μ+κ)e,diadt=κpe(μ+α+γa)ia,disdt=κ(1p)e+αia(μ+ψ+γs)is,dxdt=ϵa(t)ia+ϵs(t)isσx,}    (19)

Mε gives the coefficient matrix of the auxiliary system of system (19)

Mε=(0εβ˜(t)εηβ˜(t)εβ˜(t)K˜000000000000).    (20)

If R0 < 1, it is known from Lemma 1 which was originally stated in Wang and Zhao [2] that ρFV(ω)) < 1. We choose ε > 0 small enough giving ρFV+M(ω)) < 1. It can be concluded from Lemma 2.1 in Zhang and Zhao [43] that there exists a positive, ω−periodic function f¯(t)=(e¯(t),i¯a(t),i¯s(t),x¯(t)) such that f^(t)=eΘtf̄(t) is a solution to Equation (19), where Θ=1ωln(ρ(ΦFV+Mε(ω))). Here, ρFV+Mε(ω)) < 1 ⇒ Θ < 0, which implies f^(t)0 as t → +∞. Thus, the zero solution of Equation (19) is globally asymptotically stable. For any non-negative initial value, there is a sufficiently large M. Using the comparison theorem [46], we get f(t)Mf^(t),t > 0. Thus, we obtain e(t) → 0, ia(t) → 0, is(t) → 0, x(t) → 0 as t → +∞. By the theory of asymptotic autonomous system [47], we get


as t → +∞. Hence, if R0 < 1, the disease-free equilibrium E0 is therefore globally asymptotically stable.

3.2. Disease persistence


  X:​​​=Ω, X0:​​​={(s,e,ia,is,x)X:e(t) > 0,ia(t) > 0,is(t) > 0,x(t) > 0},X0:​​​=X\X0.

Consider P : XX, to be the Poincaré map associated with Equation (11); that is, P(z0) = u(ω, z0), ∀z0X, where ω is the period. u(t, z0) is the unique solution of Equation (11) with u(0, z0) = z0 = (s(0), e(0), ia(0), is(0), x(0)). We see that Pn(z0)=u(nω,z0),n0.

Lemma 2. For Equation (11), if R0 > 1, then there exist a ν > 0 such that, for any z0 = (s(0), e(0), ia(0), is(0), x(0)) ∈ X0 with z0E0ν, we have limn+supd[Pn(z0),E0]ν.

Proof. Since R0 > 1, Lemma 1 implies that E0 is unstable; then, ρFV(ω)) > 1. Take ε1 > 0 small enough such that ρFVMε1(ω)) > 1, where

Mε1=(0ε1β˜(t)ε1ηβ˜(t)ε1β˜(t)K˜000000000000).    (21)

Contrary, now suppose that the limn+supd[Pn(z0),E0]<ν for some z0X0. Without the loss of generality, we assume that d[Pn(z0),E0]<ν,n ≥ 0. By the continuity of the solution with respect to initial value, we have


From the periodicity of the system, for ε1 > 0, there exists tε1 such that, for all t > tε1, there holds


dedt=δ(t,ia,is,x)(s0+ε1)(μ+κ)e,diadt=κpe(μ+α+γa)ia,disdt=κ(1p)e+αia(μ+ψ+γs)is,dxdt=ϵa(t)ia+ϵs(t)isσx.}    (23)

Let now consider the Equation (23); we conclude that there exists a positive, ωperiodic function f¯(t)=(e¯(t),i¯a(t),i¯s(t),x¯(t)) such that f^(t)=eΘ1tf̄(t) is a solution of the Equation (23), where Θ1 = (1/ω)ln(ρFVMε1(ω))). Here, ρFVMε1(ω)) > 1 ⇒ Θ1 > 0, which implies that, for non-negative integer n, f^(nω)+ as n → +∞. For any non-negative initial value, there is a sufficiently small m > 0. Using the comparison theorem [46], we have f(t)mf^(t), > 0. Thus we obtain as follows:

e(t)+,ia(t)+,is(t)+,x(t)+ as t+.    (24)

This contradicts, thus proof end.

Lemma 3. The following equations are established:

M:​​={z0X0:Pn(z0)X0,n0}    (25)

={(s0,0,0,0,0)X:s00}.    (26)

Proof. We know that


We prove by contradiction as follows:


Suppose that,

z0=(s(0),e(0),ia(0),is(0),x(0){z0X0:Pn(z0)X0,n0}\{(s0,0,0,0,0)}.    (27)

We assume that e(nω) > 0 without loss of generality. From the general solution of Equation (11), we know that

e(t) > 0,ia(t) > 0,is(t) > 0,x(t) > 0.

This means that, (s(t), e(t), is(t), is(t), x(t)) ∉ ∂X0. This contradicts with (s(0), e(0), ia(0), is(0), x(0)) ∈ ∂X0. Therefore, the equation is established.

Theorem 2. If R0 > 1, then there exist ε* > 0 such that any solution (s(t), e(t), ia(t), is(t), x(t)) of Equation (11) with initial value z0 = (s(0), e(0), ia(0), is(0), x(0)) ∈ X0 satisfies

limt+infe(t)ε*,limt+infia(t)ε*,limt+infis(t)ε*,limt+infx(t)ε*,    (28)

and at least one positive periodic solution is admissible in Equation (11).

Proof. We have proved that {Pn}n0 is uniformly persistent with respect to (X0, ∂X0). For any z0X0 from the first equation of Equation (11), it follows that

s(t)=e0t(δ(s˜,ia(s˜),is(s˜),x(s˜))+μ)ds˜[s(0)+Λ˜(0te0s˜1(δ(s˜,ia(s˜),is(s˜),x(s˜)+μ)ds˜ds˜1)].    (29)

Then, s(t) > 0, ∀t > 0. As generalized to non-autonomous equations [48], the irreducibility of matrix (30) implies that e(t) > 0, ia(t) > 0, is(t) > 0, x(t) > 0, ∀t > 0 where

M˜(t)=((μ+k)β˜(t)ηβ˜(t)β˜(t)K˜κp(μ+α+γa)00κ(1p)α(μ+ψ+γs)00ϵa(t)ϵs(t)σ).    (30)

X and X0 are positively invariant. We see that ∂X0 is relatively closed in X. E0 of Equation (11) is globally asymptotically stable. Lemma (3) means that E0 is a unique fixed point of P in M. In addition, E0 is an isolated invariant set in X, and Ws(E0)X0=ϕ. Every orbit in M approaches E0 and E0 is acyclic in M. By Zhao [49], it follows that {Pn}n0 uniformly persists with respect to (X0, ∂X0) and the solutions of Equation (11) uniformly persists with respect to (X0, ∂X0); that is, if R0 > 1, there exist ε* > 0 such that any solution (s(t), e(t), ia(t), is(t), x(t)) of Equation (11) with initial values z0 = (s(0), e(0), ia(0), is(0), x(0)) ∈ X0 satisfies

limt+infe(t)ε*,limt+infia(t)ε*,limt+infis(t)ε*,limt+infx(t)ε*    (31)

Furthermore, (s*(0),e*(0),ia*(0),is*(0),x*(0))X0 is fixed point for P. Furthermore, there exists some t̄[0,ω] such that s*(t̄) > 0. If it is not the case, s* ≡ 0. s*(t̄)0, for all t ≥ 0 by the periodicity of s*(t̄). From the first equation of Equation (11), we get 0=Λ~ > 0, which is a contradiction. Thus, t[t̄,t̄+ω], we obtain

s*(t)=et¯t(δ(s˜,ia*(s˜),is*(s˜),x*(s˜))+μ)ds˜[s*(t)¯+Λ˜(t¯tet¯s˜1(δ(s˜,ia*(s˜),is*(s˜),x*(s˜))+μ)ds˜ds˜1)] > 0.    (32)

When s*(t) > 0, ∀t ≥ 0 then s*(t) is periodic. Similarly, e*(t) > 0,ia*(t) > 0,is*(t) > 0,x*(t) > 0. Therefore, (e*(t),ia*(t),is*(t),x*(t)) is a positive ω−periodic solution of Equation (11).

4. Numerical simulation

We applied our model to study the COVID-19 epidemic in South Africa. We used the new case data published daily by the National Institute for Communicable Diseases [39]. These data are made up of the daily reported new cases, the recovered, disease-induced deaths, and cumulative cases for all South African provinces.

We fitted the constructed model to new cases from South Africa to illustrate our mathematical results. The total population of South Africa in the year 2020 was 59.31 million according to Worldometer [50]. The mortality rate μ was obtained by using the life expectancy of South Africa, which was found to be 64.3 years. The proportion of exposed individuals p is in the interval (0, 1).

We conduct numerical simulations for an epidemic period starting from 11th April 2021 to 11th August 2022 in Table 1. We set time to be in days and period (ω = 182.5) days due to the number of waves in a year.

Table 1. Time and days for the waves used in the estimation of initial conditions.

To get the initial condition for asymptomatic individuals that was not given in the collected data, we used the study done by Kleynhans et al. [51]. Their studies reveal that the exact number of COVID-19 infections in Africa may be 97% higher than the number of confirmed reported cases. The rest of the initial conditions of state variables are taken from the data and are given below in Equation (33).

Third wave:      E0=16972, Ia0=16282, Is0=438, R0=890429, Ce=1000,Fourth wave:    E0=16000, Ia0=15296, Is0=598, R0=1331581, Ce=1500,Fifth wave:      E0=30000, Ia0=28744, Is0=889, R0=890429, Ce=1800,}    (33)

Initial conditions for susceptible and the asymptotic infectives are given by S0 = N − (E0 + Ia0 + Is + R0) and Ia0=973×Is0, respectively.

We considered the periodic transmission rate

δ˜(t,ia,is,x)=β˜(t)ia+ηβ˜(t)is+β˜(t)xK˜+x,    (34)

β~(t) is given as a piece-wise function below,

β˜(t)={β^[1+β¯sin(2πt182.5)],t0tt2β^[1+β¯sin(2πt182.5)]m(t),tt2.    (35)


t0=11th of April 2021, t1=29th of October 2021,    t2=19th of April 2022,t3=19th of August 2022.

The parameter β is positive and denotes the maximum value of the transmission rate. Time-varying vaccination parameter m(t). We define our m(t)=emxtx, which falls between (0, 1], mx is the speed of the vaccination (called the vaccination efficacy) and tx is the start time of the implementation of the vaccine. For our numerical simulations, we fitted our model to the data starting from the third wave to the fifth wave, see Assan and Nyabadza [45] for first and second wave fitting. We added the damping to the fifth wave because that is when mandatory vaccination was implemented in South Africa. We used the function curve-fit from the python module scipy.optimize to fit our data which uses non-linear least squares to fit data to a functional form, see Figure 2.

Figure 2. Fitting our mathematical model to the third, fourth, and fifth wave of data in piece-wise form. New cases in South Africa for the third, fourth, and fifth waves. The third wave from 11th April to 29th October 2021, the fourth wave from 30th October 2021 to 19th April, 2022 and the fifth wave from 20th April, 2022 to 19th August, 2022. The blue dots denote the new cases from the data, and the red lines denote our mathematical model prediction. We can see that in each wave our model constructed fitted well for the given parameters in Table 2.

In Table 2, we discuss in detail the values of the parameters used in Equation (11). When fitting Model (11) to the South African data obtained, we divided days according to the number of waves we have, starting with the third wave which occurred on the 400th day. When comparing the periodic transmission rate β~(t), the fifth wave has a low transmission rate; this is due to the mandatory vaccination that was introduced by the government, which reduces the disease-induced death rate drastically. Parameter values used for plotting are within the range from the third to the fifth wave, see Figures 36.

Table 2. Parameter estimation of the third, fourth, and fifth wave in South Africa in days.

Figure 3. Proportion of infected humans for different vaccination efficacy; mx = 0.9, 0.5, and 0.2 and with other parameters made constant as in Table 2. The infection persists and a periodic solution with ω = 182.25 days forms after a long transient.

We hypothetically assume carry capacity K = 106, as a result of difficulties in estimating such a parameter, the removal rate of the virus is σ=130 per day, and Lambda = 14, 000 and based on the model properties in Section 2.2, 0<B^,ϵ^a,s<1. We also chose B̄ and ϵ̄a,s to be 10. Figure 3 gives the efficacy of the vaccination mx in the disease transmission. When mx is carried out to 0.9 (the green wave), the disease transmission reduces drastically, but the disease persists after a long transient with a periodic solution ω = 182 days.

Figure 4 gives the periodic threshold of R0 by varying the efficacy of the vaccination to see its effect on infection transmission. Figure 5 gives the maximum or minimum values of the vaccination speed (vaccination efficacy), which is (0, 1.5] and remains constant (disease dies out) from 1.5 and the upper bound does not make sense after mx = 1.5.

Figure 4. Gives the periodic threshold of R0 for various mx in model (2). When mx = 1, then R0 < 1, and if mx = 0, then R0 > 1.

Figure 5. Maximum and minimum values for mx in model (2). mx = (0, 1.5] and does not make sense after 1.5 for the disease dies out completely.

Figure 6 gives two different initial conditions for the virus in the environment. The long-term behavior showed the same patterns, with infection moving toward a positive solution and a long oscillating transient, respectively, showing that a little shedding of the virus in the environment can sustain infection in the environment for the long term (see the Supplementary material for sample algorithm for some of the figures).

Figure 6. Pathogens density for different initial conditions (500 and 1, 000) over a period of 5 years of disease persistence and with the rest of the parameters being constant as defined in Table 2. The virus persists and a periodic solution with ω = 182.25 days forms after a long transient.

5. Conclusion and discussion

We presented a general non-autonomous COVID-19 model in a periodic environment. Seasonally variational factors have been incorporated into the incidence function δ and the pathogen function ζ. We proposed multiple periodic transmission pathways for COVID-19, human-to-human transmission, and human-to-environment transmission in this study. A detailed analysis of our model was conducted, and publicly reported data from the National Institute for Communicable Diseases was used for the study in South Africa.

Using the next infection operator, R0 = ρ(L), it was demonstrated that the disease-free equilibrium is globally asymptotically stable if R0 < 1. However, when R0 > 1, the disease persists uniformly, indicating that there is at least one positive periodic solution.

Our numerical simulation results demonstrate the application of our model to the COVID-19 outbreak in South Africa. Our mathematical model fits the reported data well, and through data fitting, we obtained some of our parameter values for other simulations. The South Africa data for COVID-19 show evidence of a seasonal pattern in the new cases, where peaks appear twice a year, that is in December to January and July to August due to various festive activities during these seasons. During the first four waves of COVID-19 in South Africa, non-pharmaceutical interventions were used as control measures. In this study, we proposed a control measure (vaccination), seasonality, and the role of the environment explaining the differences between the epidemic curves, starting from the third wave. A time-varying vaccination parameter m(t) was defined with the assumption that m(t) rises or falls exponentially. We implemented this in the fifth wave when the implementation of vaccination was mandatory in the South African population. We set m(t) to be between (0 − 1]; when m(t) = 1, then the epidemic spreads without any restrictions. However, m(t) can not be assumed to be one except when there are no control measures, and when m(t) = 0, then there is no epidemic. mx is the speed of the vaccination (vaccination efficacy).

We varied the efficacy of the vaccination on the infected population on a scale of 0.2, 0.5, and 0.9. This was done separately to see their effects on the disease transmission rate. It was observed that when vaccination has 0.5 efficacy, it slightly reduces the infection rate since half of the population will be aware of the government and policymakers’ prevention strategies. When vaccination is carried out to 0.9 or a mandatory vaccination is carried out effectively, the transmission rate decreases substantially, reducing the infected population. In presenting the results of R0, we varied mx while keeping other parameters fixed, and our result indicates that when mx = 1, then R0 < 1, and when mx = 0, then R0 > 1. From our result, we found out that the parameter mx = (0, 1.5] remains constant after 1.5, that is, the disease dies out completely.

Our simulation results clearly show that individuals continue to shed the virus back into the environment as long as individuals are infected. Thus, an increased transmission rate causes more susceptible individuals to become infected. This study attempts to model COVID-19 by taking seasonality, control measures (vaccination), and the influence of the role of the environment. The findings have implications for developing initiatives and regulations aimed at reducing COVID-19 transmission in the presence of seasonality.

The model presented in the study is not without limitations. After vaccination, the periodic nature of the data somewhat vanished. The modeling of such a decline requires functions that reflect damped oscillations, which will be a valuable addition to the model. The model does not include aspects of hospitalization and intensive care, which were central to the pandemic. The inclusion of these would surely make the model more robust. The model does not capture heterogeneity observed in the disease, especially aspects of age, COVID-19 variants, long COVID-19, co-morbidities, and obesity. These were important determinants of disease progression and outcomes. Despite these limitations, the model presents some useful findings in the modeling of the different waves of COVID-19.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

All authors contributed to the formulation of the model, its review and analysis, and drafting and editing of the manuscript.


The Organization for Women in Science for the Developing World (OWSD) and Swedish International Development Cooperation Agency (SIDA) funded the study.


The authors are grateful to the Department of Mathematics and Applied Mathematics, University of Johannesburg for the production of this paper.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at:


Source link