A mathematical model to study herbal and modern treatments against COVID-19

1 Introduction

In December 2019, a serious disease known as the novel Coronavirus 2 (SARS-Cov-2 or COVID-19) was reported in Wuhan Province in China. It was declared a world-wide pandemic by the World Health Organization (WHO) on March 17, 2020 [1], [2]. COVID-19 is transmitted when people breathe in air contaminated by small airborne particles containing the virus. The air is contaminated by infected people who may or may not exhibit symptoms of the disease. Infected individuals spread the virus when they cough, sneeze or talk in proximity with other people who are not protected. To mitigate the spread of COVID-19, barriers measures such as physical distancing, wearing of face masks, self-quarantine; massive screening tests; isolation of infected and the vaccination are prescribed. The effectiveness of some of these measures to overcome the disease may be questionable as COVID-19 remains persistent to date. As of June 30, 2023, there have been more than 767,726,861 of confirmed COVID-19 cases and around 6,948,764 of deaths globally [3].

Individuals infected with COVID-19 present symptoms which are similar to those of flu. These symptoms include fever, cough, tiredness, loss of taste or smell and sometimes, diarrhoea, sore throat, headache, aches and pains [4]. At a critical stage they may have chest pain, difficulty in breathing, loss of speech and mobility. However, many of the infected cases are asymptomatic or present benign symptoms. Some infected people in Africa who have not been tested relate the disease to contagious respiratory illness caused by influenza viruses and they resort to natural products/traditional medicine for remedy.

The world health organization defines traditional medicine as the sum of the knowledge, skills and practices based on the theories, beliefs and experiences in different cultures, whether explicable or not, used in the maintenance of health as well as in the prevention, diagnosis, improvement or treatment of physical and mental illnesses [5]. Africa has a long history of traditional medicine and practitioners that play a crucial role in providing health care to the population. Many centuries before the introduction of modern medicine, African people depended solely on traditional medicine to treat and prevent human and animal diseases. WHO recognizes traditional, complementary and alternative medicine of proven quality, safety and efficacy [6], [7], [8], [9], [10]. However, the organization insists on supporting scientifically-proven traditional medicine. This is the bone of contention between advocates for traditional medicine and those for modern medicine.

Though the world health organization discourages the use of herbal remedies for the treatment of COVID-19 disease because of safety concerns, many infected people in Africa continue to seek herbal/traditional treatment or combine traditional medicine with modern medicine as adjuvants. Moreover, some governments have adopted many local products to treat COVID-19 patients. This brings to mind the Madagascar’s COVID-Organics which sparked controversy over its use for lack of sufficient scientific evidence or approval by global health authorities like WHO. COVID-Organics is an extract made from a local plants. This treatment was endorsed by the President of Madagascar, Andry Rajoelina in April 2020. Other African countries adopted (though not officially) local remedies for cure and prevention of COVID-19. For instance in Cameroon, Archbishop Samuel Kleda came up with two products: Elixir Covid and Adsak Covid; Dr. J. E. Azombo discovered Fagaricine and M. P. M. Ndi Samba in partnership with Dr. Peyou developed local Africa herbs that serve as a cure for COVID-19. Dr. Marlys developed a product called Ngul Be Tara or the power of ancestors to treat this pathology [11].

Despite ongoing debate on COVID-19 treatment, many people in Africa are either ignorant or do not believe in the existence of COVID-19. While some mistake the virus for influenza, there are some who associate the disease to witchcraft or to supernatural causes. Patients from these categories will prefer either self-treatment using natural products or visit herbalists rather than going to hospitals. The treatment is cheap, affordable and readily available. A study by [5] revealed that a good proportion of university staff and students used traditional medicine in treating COVID-19 disease during the South African nationwide lockdown. The COVID-19 pandemic is far from over, but no specific treatment for COVID-19 has proven 100 % effective and hence any effective treatment be in modern (mainstream) or traditional or a blend of the two will be of great relief to humanity. Biomathematicians with their mathematical models have an important role to play in this endeavour [12]. This is the motivation behind our model which incorporates among others parameters for modern and traditional medicines.

Many mathematical models are presented in literature to explain and predict evolution of the COVID-19 pandemic since its outbreak (see for instance [13], [14], [15], [16], [17] and references therein). Harjule et al. [13] presents a review of the various mathematical and statistical methods adopted so far to anticipate and analyse the outspread of the virus across the globe. In [17], a model is presented to evaluate the impact of control measures, such as quarantine and hospitalization on the spread of the corona virus. It is shown that when control measures are implemented, unreported symptomatic cases decrease faster than when control measures are neglected. The authors of the article [16] focus on the usefulness of contact tracing, quarantine and physical distancing in mitigating the intensity of COVID-19. Applying their model on data from Poland, they found out that contact tracing may be prevented in 50 % to 90 % of cases and that the effect of quarantine was limited by the number of undiagnosed cases. Furthermore, the disease could not be controlled without physical distancing measures. In the paper [15], the short-term forecasts of the COVID-19 pandemic in Cameroon was studied using a mathematical model. The model considers both direct transmission (human to human transmission) and indirect transmission due to the circulation of free Coronavirus in the environment. They discovered that, the former mode of transmission lead to more cases than the latter. Their results showed that, the number of COVID-19 cases will still increase and there is a necessity to improve the implementation of preventive measures. In [14], the influence of the presence of multi-strains of Coronavirus was investigated. The results confirmed that the alpha, beta and gamma variants are more transmissible than the original viral strain.

To the best of our knowledge and at the time of writing this paper, no author has investigated the role of traditional medicine in the control of the COVID-19 pandemic through mathematical modelling. In this work, we propose a two-group deterministic COVID-19 model which takes into account educational campaigns and the fact that people infected with COVID-19 may choose either modern (allopathic) medicine, traditional medicine or may combine the two modes of treatment. The model is analysed mathematically and numerically. In the case where infected people only recourse to modern medicine for their treatment, we prove that the disease-free equilibrium (DFE) is globally asymptotically stable (GAS) when the control reproduction number is less than one. When it is greater than one, we prove the existence of a unique endemic equilibrium, locally asymptotically stable (LAS). When people receive both modern and traditional treatment, we prove that the model has a backward bifurcation phenomenon. However, we derive a threshold T which ensures the elimination of the disease in some cases when its value is less than one. Numerical simulations are carried out to support mathematical analysis results. We fit the model to the daily cumulative cases of COVID-19 reported from January to February 2022 in South Africa and investigate the usefulness of educational campaigns, traditional and modern treatments. Sensitivity analysis of the model is performed to identify the most influential parameters on the dynamics of the model.

The rest of this article is partitioned as follows: In Section 2, the model and its quantitative analysis are presented in a comprehensive manner. Section 3 deals with the mathematical analysis and the illustration of some mathematical results. Model calibration and validation is presented in Section 4. The usefulness of educational campaigns and the place of each kind of treatment are presented in Section 5. Global sensitivity analysis in order to determine the most influential parameters is addressed in Section 6. The discussion is provided in Section 7.

2 Model formulation

2.2 Model description

Our model follows the flowchart presented in Figure 1. The susceptible population is recruited at rate π. Among these recruited individuals, a proportion τ is educated on COVID-19 while the remaining (1 − τ) is not educated. The two susceptible classes can be infected by people from classes I ₁, J, I ₂, C and T, at constant rates β, ν ₁ β, ν ₂ β, ν ₃ β and ν ₄ β, respectively. Thus, the force of infection following the standard incidence is given by:

λ ( t ) = β ( I 1 ( t ) + ν 1 J ( t ) + ν 2 I 2 ( t ) + ν 3 C ( t ) + ν 4 T ( t ) ) N ( t ) .

Susceptible individuals who are educated on COVID-19 follow WHO recommendations. Their behaviours will reduce their susceptibility by a factor ξ ≔ (1 − ɛ), where ɛ measures the efficiency of protective measures (or social distancing) taken by the individuals in compartment S ₁.

Figure 1:

Flow diagram for Model (1) and for the restricted Model (2).

Susceptible individuals once becoming infected arrive in compartments I ₁ and I ₂, depending on whether they fall in the category of educated on the virus or not, respectively. Individuals in compartment I ₁ may recover at rate γ ₁, decease at rate (μ + δ ₁) or seek treatment through modern medicine at rate α ₁. The parameter μ is the natural mortality rate while δ ₁ is the mortality due to SARS-CoV-2 virus.

Patients in compartment I ₂ seek traditional treatment at rate σ, recover at rate γ ₂, die at rate (μ + δ ₂). Some infected persons taking modern treatment combine it to traditional treatment as adjuvant (i.e. move to the class C) at rate η.

In compartment C, the death rate is (μ + δ _c ) and the recovery rate is γ _c . We assume that mass media disseminate information on the WHO guidelines at rate q. Meanwhile m _s , m _i , m _c and m _t measure the probability that a susceptible individual or an infected uneducated being influenced by the information received in the compartments S ₂, I ₂ and T. Thus due to sensitisation, m _s q, m _i q, m _c q and m _t q are the rates at which susceptible and infected individuals become aware and leave compartments S ₂, I ₂ and T to compartments S ₁, J and C, respectively. A complete list of parameters with definitions is presented in Table 1.

Table 1:

Parameters and their epidemiological interpretations.

Parameter	Epidemiological interpretation
β	Effective transmission rate of COVID-19 due to contact with individuals in compartment I 1.
ν 1, ν 2, ν 3, ν 4	Modification parameters for the infectiousness of individuals in compartments J, I 2, C and T, respectively.
π	Recruitment rate in the susceptible population.
τπ	Proportion of susceptible individuals recruited in the compartment S 1.
α 1	Exit rate from the compartment I 1 to the compartment J.
σ	Exit rate from compartment I 2 to the compartment T.
γ 1, γ 2, γ j , ω, γ c	Recovery rates of infected individuals in the compartments I 1, I 2, J and T, respectively.
δ 1, δ 2, δ j , ψ, δ c	Death rates due to COVID-19 in compartments I 1, I 2, J, T and C, respectively.
ɛ	Efficiency of protective measures taken by susceptible individuals in compartment S 1.
q	Spread rate of information in accordance with WHO guidelines.
m s	Probability that an unaware susceptible individual is influenced by information and moves to S 1.
m i (m c )	Probability that an uneducated infected individual is influenced by information and moves to J (C).
m t	Probability that infected who follow traditional medicine are influenced by information and move to compartment C.
η	Exit rate from compartment J to the compartment C.
μ	Natural mortality rate.

The model which follows is given by the ordinary differential system:

(1) S ̇ 1 ( t )   =   τ π + m s q S 2 − λ ( 1 − ε ) S 1 − μ S 1 , S ̇ 2 ( t )   =   ( 1 − τ ) π − λ S 2 − ( μ + m s q ) S 2 , I ̇ 1 ( t )   =   λ ( 1 − ε ) S 1 − ( μ + α 1 + γ 1 + δ 1 ) I 1 , I ̇ 2 ( t )   =   λ S 2 − ( μ + γ 2 + δ 2 + m c q + σ + m i q ) I 2 , J ̇ ( t )   =   α 1 I 1 + m i q I 2 − ( μ + γ j + δ j + η ) J , T ̇ ( t )   =   σ I 2 − ( μ + ω + ψ + m t q ) T , C ̇ ( t )   =   m c q I 2 + η J + m t q T − ( μ + γ c + δ c ) C , R ̇ ( t )   =   γ 1 I 1 + γ 2 I 2 + γ j J + ω T + γ c C − μ R .

Nonetheless, due to evolution and development in science, controversial views as regards some diseases within some societies are gradually changing for the better. In modern societies, the belief is that diseases are exclusively due to natural causes. Education in this context will be oriented towards means of prevention and not towards the choice of a mode of treatment, since people in such a society are naturally accustomed to going to hospital facilities when they are ill. Under this assumption, Model (1) becomes:

(2) S ̇ 1 ( t )   =   π − λ ξ S 1 − μ S 1 , I ̇ 1 ( t )   =   λ ξ S 1 − ( μ + α 1 + γ 1 + δ 1 ) I 1 , J ̇ ( t )   =   α 1 I 1 − ( μ + γ j + δ j ) J , R ̇ ( t )   =   γ 1 I 1 + γ j J − μ R .

For mathematical convenience, we define the following new parameters.

π 1 = τ π , π 2 = ( 1 − τ ) π , ϕ j = μ + γ j + δ j + η , m 1 = m s q , m 2 = μ + m s q , ϕ 1 = μ + α 1 + γ 1 + δ 1 , m 3 = m i q , α 2 = m c q , ϕ c = μ + γ c + δ c , ξ = 1 − ε , ν = m t q , ϕ t = μ + ω + ψ + ν , a j = α 1 ξ ϕ 2 π 1 + m 3 ξ ϕ 1 π 2 , b j = α 1 ξ ϕ 2 ( m 1 π 2 + m 2 π 1 ) + m 3 ϕ 1 π 2 μ , a c = ξ ϕ t ϕ j ϕ 1 α 2 π 2 + η a j ϕ t + ξ ϕ j ϕ 1 ν σ π 2 , b c = μ ϕ t ϕ j ϕ 1 α 2 π 2 + η b j ϕ t + μ ϕ j ϕ 1 ν σ π 2 k n = ( ϕ j μ + α 1 μ + γ 1 ϕ j + γ j α 1 ) , ϕ 2 = μ + γ 2 + δ 2 + m c q + σ + m i q ,

Λ = μ + δ 1 + δ 2 + δ j + ψ + δ c , a r = γ 1 ξ π 1 ϕ c ϕ t ϕ j ϕ 2 + γ j a j ϕ c ϕ t + γ c a c + μ γ 2 π 2 ϕ c ϕ t ϕ j ϕ 1 ξ , b r = γ 1 ξ ( m 1 π 2 + m 2 π 1 ) ϕ c ϕ t ϕ j ϕ 2 + μ γ 2 π 2 ϕ c ϕ t ϕ j ϕ 1 + γ j b j ϕ c ϕ t + γ c b c , a n = ξ π 2 μ ϕ c ϕ t ϕ j ϕ 1 + ξ π 1 μ ϕ c ϕ t ϕ j ϕ 2 + a j μ ϕ c ϕ t + ξ σ π 2 μ ϕ c ϕ j ϕ 1 + a c μ + a r , b n = π 1 μ ϕ c ϕ t ϕ j ϕ 2 ϕ 1 + μ 2 π 2 ϕ c ϕ t ϕ j ϕ 1 + ξ ( m 1 π 2 + m 2 π 1 ) μ ϕ c ϕ t ϕ j ϕ 2 + ξ π 2 μ ϕ c ϕ t ϕ j ϕ 2 ϕ 1 + b j μ ϕ c ϕ t + μ 2 σ π 2 ϕ c ϕ j ϕ 1 + b c μ + b r , c n = μ ϕ c ϕ t ϕ j ϕ 1 ϕ 2 ( μ π 2 + ( m 1 π 2 + m 2 π 1 ) ) , a p = β ξ π 1 μ ϕ c ϕ t ϕ j ϕ 2 + β ν 2 ξ π 1 μ ϕ c ϕ t ϕ j ϕ 1 + β ν 1 a j μ ϕ c ϕ t + β ν 3 a c μ + β ν 4 ξ σ π 2 μ ϕ c ϕ j ϕ 1 , b p = β ξ ( m 1 π 2 + m 2 π 1 ) μ ϕ c ϕ t ϕ j ϕ 2 + β ν 2 μ 2 π 2 ϕ c ϕ t ϕ j ϕ 1 + β ν 1 b j μ ϕ c ϕ t + β ν 3 b c μ + β ν 4 μ 2 σ π 2 ϕ c ϕ j ϕ 1 .

Model (1) considers the dynamics of the human population split in different compartments. Thus for biological reasons it is necessary that these variables remain positive. This is guaranteed by the following propositions, which ensure the well-posedness of Model (1) and give the biologically feasible region.

Proposition 1

For ( S 1 ( 0 ) , S 2 ( 0 ) , I 1 ( 0 ) , I 2 ( 0 ) , J ( 0 ) , T ( 0 ) , R ( 0 ) ) ∈ R 8 \ { 0 } , System (1) has a unique local solution. Moreover, if

S 1 ( 0 ) > 0 , S 2 ( 0 ) > 0 , I 1 ( 0 ) ≥ 0 , I 2 ( 0 ) ≥ 0 , J ( 0 ) ≥ 0 , T ( 0 ) ≥ 0  and  R ( 0 ) ≥ 0 ,

then

S 1 ( t ) > 0 , S 2 ( t ) > 0 , I 1 ( t ) ≥ 0 , I 2 ( t ) ≥ 0 , J ( t ) ≥ 0 , T ( t ) ≥ 0  and  R ( t ) ≥ 0 , ∀ t > 0 .

Proof

The right-hand side of System (1) is a vector valued function from R 8 into R 8 , which is continuously differentiable on R 8 \ { 0 } . Thus this function is locally Lipschitz. Consequently, by the Cauchy–Lipschitz theorem, the initial value problem for this system has a unique local solution on some interval [0, T].

Assume S ₂(0) > 0 and S ₁(0) > 0. By the second equation of System (1), one has,

S 2 ( t ) = S 2 ( 0 ) exp − ∫ 0 t ( λ ( s ) + m 2 ) d s + exp − ∫ 0 t ( λ ( s ) + m 2 ) d s × ∫ 0 t π 2 ⁡ exp ∫ 0 v ( λ ( u ) + m 2 ) d u d v > 0 .

By the first equation of Model (1), one gets,

S 1 ( t ) = S 1 ( 0 ) exp − ∫ 0 t ( λ ( s ) ξ + μ ) d s + exp − ∫ 0 t ( λ ( s ) ξ + μ ) d s × ∫ 0 t ( π 1 + m 1 S 2 ( s ) ) exp ∫ 0 s ( λ ( u ) ξ + μ ) d u d s > 0 .

To prove the positivity of the variables I ₁, I ₂, J, T, C and R, we use the Tangent theorem [19]. For this, we have to show that ⟨u(y)|g(y)⟩ ≤ 0, for y on each of the hyperplanes I ₁ = 0, I ₂ = 0, J = 0, T = 0, C = 0 and R = 0, u being the outer normal vector to the hyperplane and g(y) is the vector field defined by the remaining equations.

For the hyperplane I ₁ = 0, we have the inner product

〈 ( − 1,0,0,0,0,0 ) | λ ξ S 1 , λ S 2 − ϕ 2 I 2 , m 3 I 2 − ϕ j , σ I 2 − ϕ t T , α 2 I 2 + η J + ν T − ϕ c T , γ 2 I 2 + γ j J + ω T + γ c C − μ R 〉 = − λ ξ S 1 ≤ 0 .

Using the other hyperplanes, one can get similar results which prove that the variables I ₁, I ₂, T, C, J, R remain non-negative for non-negative initial conditions.□

Proposition 2

Assume

S 2 ( 0 ) ≤ π 2 m 2 , S 1 ( 0 ) ≤ π 1 m 2 + m 1 π 2 μ m 2  and  π Λ ≤ N ( 0 ) ≤ π μ .

Then

∀ t > 0 , S 2 ( t ) ≤ π 2 m 2 , S 1 ( t ) ≤ π 1 m 2 + m 1 π 2 μ m 2  and  π Λ ≤ N ( t ) ≤ π μ .

Proof

Let us assume that S 2 ( 0 ) ≤ π 2 m 2 . By the second equation of System (1), it is straightforward that

S ̇ 2 ( t ) ≤ π 2 − m 2 S 2 .

So according to the Gronwall lemma,

S 2 ( t ) ≤ π 2 m 2 + S 2 ( 0 ) − π 2 m 2 e − m 2 t .

That is

S 2 ( t ) ≤ π 2 m 2 ∀ t > 0 .

If one assumes in addition that S 1 ( 0 ) ≤ π 1 m 2 + m 1 π 2 μ m 2 , the application of the Gronwall lemma once more to the first equation of System (1) leads to

S 1 ( t ) ≤ π 1 m 2 + m 1 π 2 μ m 2 ∀ t > 0 .

Assume now that π Λ ≤ N ( 0 ) ≤ π μ . Adding all the equations of Model (1) gives,

N ̇ ( t ) = π − μ N − δ 1 I 1 − δ 2 I 2 − δ j J − ψ T − δ c C .

Therefore,

π − Λ N ≤ N ̇ ≤ π − μ N .

By the Gronwall lemma, one obtains

π Λ + N ( 0 ) − π Λ e − Λ t ≤ N ( t ) ≤ π μ + N ( 0 ) − π μ e − μ t .

That is

π Λ ≤ N ( t ) ≤ π μ ∀ t > 0  when  π Λ ≤ N ( 0 ) ≤ π μ .

□

From Propositions 1 and 2, one deduces the following result about the well-posedness of Model (1).

Proposition 3

The Model (1) is a dynamical system on the biologically feasible domain

Ω : = ( S 1 ( t ) , S 2 ( t ) , I 1 ( t ) , I 2 ( t ) , J ( t ) , T ( t ) , C ( t ) , R ( t ) ) ∈ R + 8 / S 2 ( t ) ≤ π 2 m 2 ,   S 1 ( t ) ≤ π 1 m 2 + m 1 π 2 μ m 2 , π Λ ≤ N ( t ) ≤ π μ .

3 Mathematical analysis of Model (1)

In order to highlight the influence of misconceptions in the dynamics of Model (1), we first analyse Model (2). The dynamics of this restricted model is summarized in the following theorem.

The proof of this theorem is presented in Appendix B. We illustrate in Figure 2 the local stability of the endemic equilibrium for Model (2) when R 0 W T > 1 . The stability of the DFE for R 0 W T < 1 will be illustrated later in numerical simulations using the estimated parameters.

Figure 2:

Stability of the endemic equilibrium E* for R c W T > 1 .

Figure 2 considers the value γ _j  = 0.0037; δ _j  = 0.0000952. The other parameters are as in Table 3. With these parameters R c W T = 2.9471 > 1 .

3.1 Determination of equilibria for Model (1)

It is straightforward to show that Model (1) has a unique equilibrium E ₀ = (S ₁₀, S ₂₀, 0, 0, 0, 0, 0, 0), with

S 10 = m 2 π 1 + m 1 π 2 μ m 2  and  S 20 = π 2 m 2 .

Let E * = S 1 * , S 2 * , I 1 * , I 2 * , J * , T * , C * , R * be a non-frontier equilibrium for Model (1). From the set of parameters defined above, one can easily get

(3) S 2 * = π 2 λ * + m 2 , S 1 * = π λ * + m 1 π 2 + m 2 π 1 λ * + m 2 ( λ * ξ + μ ) , I 1 * = λ * ξ S 1 * ϕ 1 = ξ π 1 λ * 2 + ξ ( m 1 π 2 + m 2 π 1 ) λ * ϕ 1 λ * + m 2 ( λ * ξ + μ ) , I 2 * = λ * π 2 ϕ 2 λ * + m 2 , J * = α 1 I 1 * + m 3 I 2 * ϕ j = a j λ * 2 + b j λ * ϕ j ϕ 1 ϕ 2 λ * + m 2 ( λ * ξ + μ ) , C * = α 2 I 2 * + η J * + ν T * ϕ c = a c λ * 2 + b c λ * ϕ j ϕ c ϕ t ϕ 1 ϕ 2 λ * + m 2 ( λ * ξ + μ ) , T * = σ π 2 λ * ϕ t ϕ 2 λ * + m 2 , R * = a r λ * 2 + b r λ * μ ϕ 1 ϕ 2 ϕ j ϕ t ϕ c λ * + m 2 ( λ * ξ + μ ) , N * = S 1 * + S 2 * + I 1 * + I 2 * + J * + C * + T * + R * = a n λ * 2 + b n λ * + c n μ + ϕ c ϕ t ϕ j ϕ 1 ϕ 2 λ * + m 2 ( λ * ξ + μ ) .

Set

P * = β I 1 * + β ν 2 I 2 * + β ν 1 J * + β ν 3 C * + β ν 4 T * .

Similar computations show that

P * = a p λ * 2 + b p λ * .

Thus

λ * = P * N * = a p λ * 2 + b p λ * a n λ * 2 + b n λ * + c n .

Therefore

λ * = 0  or  a n λ * 2 + ( b n − a p ) λ * + c n − b p = 0 .

The value λ* = 0 leads to the frontier equilibrium. The other value of λ* verifies the equation

(4) a n c n λ * 2 + ( b n − a p ) c n λ * + 1 − b p c n = 0 .

Multiplying both sides of Equation (4) by c _n /a _n , one gets

(5) λ * 2 + ( b n − a p ) a n λ * + c n a n 1 − b p c n = 0 .

we note that

( b n − a p ) a n = c n ( b n − a p ) + a n b p a n c n − b p c n .

Thus Equation (5) is equivalent to the equation

(6) λ * 2 + ( K − R c ) λ * + C 1 − R c = 0 ,

with

R c = b p c n , C = c n a n and  K = c n ( b n − a p ) + a n b p a n c n .

The discriminant Δ ( R c ) of Equation (6) is equal to

Δ ( R c ) = R c 2 + R c ( − 2 K + 4 C ) + K 2 − 4 C .

Set

R 1 = K − 2 C − 2 C ( 1 − K ) + C 2  and  R 2 = K − 2 C + 2 C ( 1 − K ) + C 2 .

Whenever they exist, R ₁ and R ₂ are the roots of the equation Δ ( R c ) = 0 . According to [20], we have the following theorem:

Theorem 2

The following statements hold.

1. Model (1) always has a unique disease-free equilibrium E ₀ = (S ₁₀, S ₂₀, 0, 0, 0, 0, 0, 0).

2. If R c > 1 , Model (1) has a unique positive equilibrium E*.

3. If K > 1, Model (1) has no positive equilibrium when R c ≤ 1 . That is Equation (6) exhibits a forward bifurcation.

4. If K ≤ − 4 C then R ₂ < 0 and Model (1) has two positive equilibria when R c < 1 , that is Equation (6) exhibits a full backward bifurcation.

5. If − 4 C < K < 1 then 0 < R ₂ < 1 and Equation (6) exhibits backward bifurcation. That is Model (1) has:

   1. No positive equilibrium if R 0 < R 2 .

   2. One positive equilibrium if R 0 = R 2 .

   3. Two positive equilibria if R 0 ∈ ( R 2 , 1 ) .

R c is the number of secondary infections produced by an index case in a completely susceptible population, during its entire infectious period in spite of the treatments and the educational campaigns implemented.

The existence of two positive equilibria in some cases whenever the control reproduction number is less than one shows that our model may undergo a backward or a full backward bifurcation phenomenon. The existence of backward bifurcation makes the control of the disease more difficult. The requirement that R c should be less than one is no longer sufficient for the elimination of the disease. More efforts need to be made to reduce and maintain R c to a value less than R c c ≔ R 2 in order to ensure the elimination of the disease. Figure 3 shows the bifurcation diagrams produced by MATLAB. Firstly, Figure 3(a) indicates that a stable disease–free equilibrium co-exists with a stable endemic equilibrium, when R c c < R c < 1 . This supports the presence of the backward bifurcation proven. On the other hand, Figure 3(b) illustrates the case where the disease-free equilibrium loses its stability while the endemic equilibrium is stable for the values of R c greater than one. This highlights the forward bifurcation phenomenon. The presence of a full backward bifurcation shows that, in some cases, two positive equilibria may always co-exist with the DFE. In this case, the global elimination of the disease is not possible since one endemic equilibrium is always stable when R c < 1 . The bifurcation diagram exhibiting the existence of a full backward bifurcation phenomenon is given in Figure 3(c).

Figure 3:

Bifurcation diagrams. (a) Is plotted for: μ = 10.000814; β = 0.0501; ɛ = 0.7396; τ = 100.03799; m _s  = 0.1123; m _i  = 0.1021; q = 0.3548; α ₁ = 100000000.0277; m _c  = 0.0238; σ = 0.068; η = 0.0327; m _t  = 0.02208; γ ₁ = 0.07571; γ ₂ = 0.00000612; γ _j  = 0.150; ω = 0.0197; γ _c  = 0.0488; δ ₁ = 0.1765; δ ₂ = 0.0116; δ _j  = 0.03608; ψ = 0.0152; δ _c  = 0.102; π = 10.003; ν ₁ = 10.000609; ν ₂ = 1.4014; ν ₃ = 10.7224; ν ₄ = 10.9503. (b) Is plotted for: μ = 100.000814; α ₁ = 0.0277; ω = 0.0000197; ν ₃ = 10.7224; ν ₄ = 1.9503 (the other parameters are same as in (a)). (c) Is plotted for: ɛ = 0.07396; m _s  = 1.1123; q = 5.3548; α ₁ = 10000.0277 (the other parameters same as in (a)).