J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 Journal of the Nigerian Society of Physical Sciences Analysis of an HIV - HCV simultaneous infection model with time delay A. S. Hassan∗, N. Hussaini Department of Mathematical Sciences, Faculty of Physical Sciences, Bayero University, Kano, Nigeria. Abstract A novel mathematical delay model for simultaneous infection of HIV and hepatitis C virus is formulated and dynamically analyzed. Basic properties of the model are established and proved. Basic reproductive threshold is systematically calculated as the maximum of three subthreshold parameters. A disease free equilibrium is determined to be globally asymptotically stable for all values of the delay when the threshold is less than unity. However, when the threshold is greater than one, endemic equilibrium emerged which is shown to be locally asymptotically stable for any length of delay. Although the delay has no effect on stabilities of equilibria points, however, it is found to reduce the infectivity of the viruses as the length of the delay is increased. Epidemiological interpretations of the results and numerical simulations illustrating them are given. DOI:10.46481/jnsps.2021.109 Keywords: Delay, HIV/HCV simultaneous infection, reproduction number, equilibrium, stability. Article History : Received: 28 May 2020 Received in revised form: 27 June 2020 Accepted for publication: 18 July 2020 Published: 27 February 2021 c©2021 Journal of the Nigerian Society of Physical Sciences. All rights reserved. Communicated by: T. Latunde 1. Introduction Around 2.75 million people who have human immunodefi- ciency virus (HIV) are coinfected with hepatitis C virus (HCV) globally and and on average HIV-infected individuals are six time more likely to get HCV infection than HIV-uninfected [1]. Worldwide, HIV and HCV are public health challenges which affected several populations. Across the continents, there are about 37 million and 115 people infected with HIV and chronic HCV infections, respectively. Numerous mathematical models have been used to explore theoretical aspects of the transmission dynamics for superinfec- tion (strains or diseases never co-exist in a host) [2-4] and for co-infection (diseases can co-exist in the host) [5-15]. Most of ∗Corresponding author tel. no: +2348036289110 Email address: hashitu.mth@buk.edu.ng (A. S. Hassan ) the co-infection models, for simplicity, do not allow simulta- neous infection [5]. Owing to the growing empirical evidence of simultaneous infection (see, for example, [16-20]), which is now a major public health concern [16] and affects the epidemi- ology [21] and evolution [22] of infectious diseases, Zhang et al. [14, 15, 24] recently developed models which include si- multaneous infection. These models are only suitable to com- pletely curable diseases such as influenza. Alizon [5] studied the effect of co-infection where parasites compete. Of recent, some mathematical models for HIV/AIDS and HCV separately, are formulated and analyzed to explore impacts of the two dis- eases, see for example [25-28]. In Dhutta and Gupta [25], a mathematical model is pro- posed for the dynamics of HIV/AIDS with incorporation of weak CD4+ T cells. They computed the local stability of the infection-free and infection equilibria for the model when the 1 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 2 valuable reproduction number is less than and greater than one respectively. Furthermore, using Lyapunov’s second method and the geometric approach, they define novel conditions for the global stability of the equilibria. Further, Maimuna [27] formulated a mathematical model for the spread of HIV with an Anti-Retroviral Therapy (ART) intervention. The author es- tablished that dynamics of the model with no ART depends on the basic reproduction number. However, sensitivity analysis revealed that the basic reproduction number decrease with in- creasing number of infected humans who follow the ART treat- ment. Jia et al [26], constructed a dynamic model for HCV transmission and prevalence based on the reported data from China and determined the most influential parameters to eval- uate the effectiveness of control measures. Finally, in [28], the authors applied ideas from the ecology of infectious diseases to model the transmission of HCV in a population of injection drug users. The authors suggested that modelling HCV as an in- directly transmitted infection facilitates a more nuanced under- standing of disease dynamics. In addition, sensitivity analysis of parameters on the value of R0 was carried out and parame- ters related to an interaction with the environmental reservoir are found to be more influential. Time delay in epidemiology is mostly incorporated to repre- sent developmental stages, incubation period or wining of im- munity [29, 30]. It is an important component that cannot be neglected as it affects the dynamics of models causing insta- bilities of equilibria resulting in Hopf bifurcations [32, 31, 33]. In addition, an important threshold value (the basic reproduc- tion number), as can be seen later in this article, is expressed in terms of time delay [34, 30, 35]. Many delay models of HCV/hepatitis B virus (HBV) have been presented in the litera- ture to explore the effects of delay as incubation period. Baner- jee et al. [36] presented an intracellular time delay model for hepatitis C virus which has shown that the delay doesn’t affect the stability of steady-state, however, destabilized the infected equilibrium with resulting in Hopf bifurcation. Gourley et al. [37] considered the dynamical properties of HBV - delayed model with standard incidence formulation. They realized that the infected steady-state, expressed in terms of delay, is glob- ally stable regardless of the time delay length. In [33], Zhao and Xu presented a delay HCV model using Beddigton-DeAngels formulation and obtained global stabilities of equilibria when the threshold parameters satisfy certain conditions. In this study, we begin by developing a model using delay differential equation for HIV and HCV co-infection which in- cludes simultaneous transmission from person-to-person. Un- like in [5, 14, 15, 24], we introduce delay in order to capture the incubation periods of both HIV and HCV infections. In this research, the model is presented in Section 2 while the analysis is given in Section 3 . In Section 4, we present illustrative graphs for the dynamical properties of the model. Lastly, conclusion and summary of our findings are reported in Section 5. 2. Model Formulation The total human population at time t, N(t), is divided into four different classes of healthy (susceptible, (S (t))), HIV-only infected individuals (IH (t)), HCV-only infected individuals (IC (t)), individuals simultaneously infected with both HIV and HCV (IHC (t)), so that N(t) = S (t) + IH (t) + IC (t) + IHC (t). The susceptible population (S (t)) is increased by the recruit- ment of people into the community at a constant rate Π (all newly-recruited individuals are assumed to be susceptible to both infections) and by recovery from HCV (at a rate ψ). The population is decreased by infection with HIV only (at a rate λH ) or HCV only (at a rate λC ) or simultaneous infection with HIV and HCV (at a rate λHC ) (see, for instance, [19, 17] and ref- erences therein for the biology). However, the infections with HIV and HCV do not occur as soon as the infectives get into contact with the susceptible individuals. Rather, there are time lapses or delays, representing the latent periods, in which the infectives progress to fully infectious people. Thus we have IH (t − τ1)e−µτ1 representing the HIV only infectives that can only infect after the elapse of τ1, (the latent period for HIV) and e−µτ1 , the survival probability over the latent period. Simi- larly, IC (t −τ2)e−µτ2 represent the HCV only infectives that can infect after the elapse of τ2, (the latent period for HCV), e−µτ2 is the survival probability for natural death over the latent period. Lastly, IHC (t − τ3)e−µτ3 stand for the HIV and HCV infectives that can infect after the elapse of τ3, (the latent period for simul- taneous infection of HIV and HCV), while e−µτ3 is the survival probability over the latent period. Thus, λH = βH [ IH (t −τ1)e−µτ1 + ηIHC (t −τ3)e−µτ3 ] N(t) , λC = βC [ IC (t −τ2)e−µτ2 + ηIHC (t −τ3)e−µτ3 ] N(t) and λHC = βHC IHC (t −τ3)e−µτ3 N(t) . (1) In (1), βH , βC and βHC represent the effective contact rates for HIV, HCV and simultaneous HIV/HCV infections, respec- tively. The parameter η > 1 accounts for the assumed increase in infectiousness of dually-infected individuals due to high im- munosuppression caused by the simultaneous infection of HIV and HCV, compared to singly-infected with HIV or HCV. For mathematical simplification we assume that τ1 = τ2 = τ3 = τ. The susceptible population further decreased by natural death (at a rate µ; natural death occurs in all compartments at the same rate). dS dt = Π + ψIC − (λH + λC + λHC ) S (t) −µS (t). (2) The population of individuals infected with HIV-only (IH ) is generated by HIV infection following the effective contact (at the rate λH ) and by dually-infected individuals (IHC ) recov- ered from HCV (at a rate κψ; where the modification parameter 0 < κ < 1 models the reduced likelihood of dually-infected 2 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 3 person to recover from HCV, in comparison to those infected with HCV only). It is decreased by HCV infection at rate λC , HIV-related death at rate δC and natural death. dIH dt = λH S + κψIHC − θ1λC IC − (µ + δH )IH. (3) The population of individuals infected with HCV-only (IC ) is generated by HCV infection (at the rate λC ). It is decreased by HIV infection at rate λH ), recovery from HCV at rate ψ, HCV-related death at rate δH and natural death. dIC dt = λC S − θ2λH IC − (ψ + µ + δC )IC. (4) The class of individuals simultaneously infected (IHC ) is generated by HIV and HCV simultaneous infections at rate λHC , by HIV-only infected individuals (IH (t)) who are HCV infected at rate θ1λC and by HCV-only infected individuals (IC (t)) who are HIV infectious at rate θ2λH , where the parameter θi ≥ 1 i = 1, 2 captures the fact that singly-infected individuals have weak immune-system (due to illness) compared to wholly-susceptible individuals. The population suffers death due to the co-infection of both diseases (at a rate δHC ). dIHC dt = λHC S + θ1λC IH + θ2λH IC − (κψ + µ + δHC )IHC (5) The following system of equations describes the model based on the aforementioned assumptions and notations while the pa- rameters are presented in Table 1: dS dt = Π + ψIC (t) − (λH + λC + λHC ) S (t) −µS (t), dIH dt = λH S (t) + κψIHC (t) − θ1λC IH (t) − K1 IH (t), dIC dt = λC S (t) − θ2λH IC (t) − K2 IC (t), dIHC dt = λHC S (t) + θ1λC IH (t) + θ2λH IC (t) − K3 IHC (t), (6) with initial data S (t) = φ1(t) ≥ 0, IH (t) = φ2(t) ≥ 0, IC (t) = φ3(t) ≥ 0, IHC (t) = φ4(t) ≥ 0, for t ∈ [−τ, 0], (7) where φ1(t), · · · ,φ4(t) are continuous functions on the interval [−τ, 0] and K1 = µ + δH, K2 = ψ + µ + δC, K3 = κψ + µ + δHC . 2.1. Basic properties In this part, we present results for basic qualitative proper- ties of the model as follows. Theorem 1. The solution (S (t), IH (t), IC (t), IHC (t)) of the model (6), with initial data (7), exists for all time, t > 0 and is unique. Furthermore, the solution is nonnegative for all t > 0. Proof. For existence and uniqueness of solution, we start as follows: Let X(t) = (S (t), IH (t), IC (t), IHC (t)). The system (6), can be represented as dX dt = f (t, Xt), where Xt(θ) = X(t + θ), f is Lipschitz and continuous in X. It follows from Theorems 2.1, 2.3 in [38] that the model (6) has a unique solution (S (t), IH (t), IC (t), IHC (t)) satisfying the initial data (7). The positivity of solution is proved by the method of con- tradiction as shown below: Suppose the conclusion is not true. Under the given ini- tial function, there exists a time, t1 ∈ [0, +∞) where S (t) will changes sign, at least once, so that S (t) > 0 for all t ∈ (0, t1), S (t1) = 0 and S (t) < 0 for t > t1. Furthermore, IH (t) > 0, IC (t) > 0, IHC (t) > 0 for t ∈ (0, t1); or a time t2 such that IH (t) > 0 for all t ∈ (0, t2), IH (t2) = 0 when S (t) > 0, IC (t) > 0, IHC (t) > 0 for t ∈ (0, t2); or a time t3 such that IC (t) > 0 for all t ∈ (0, t3), IC (t3) = 0 when S (t) > 0, IH (t) > 0, IHC (t) > 0 for t ∈ (0, t3); or a time t4 such that IHC (t) > 0 for all t ∈ (0, t4), IHC (t4) = 0 when S (t) > 0, IH (t) > 0, IC (t) > 0 for t ∈ (0, t4). For the first case, consider the first equation in (6), thus dS dt (t1) = Π + ψIC (t1) − (λH + λC + λHC ) S (t1) −µS (t1), = Π + ψIC (t1), hence S (t1) = S (θ) ∫ t 0 exp (Π + ψIC (t1))du, > 0. This is a contradiction to the earlier assumption that S (t1) = 0. Therefore, t1 doesn’t exists, hence S (t) > 0 for all t > 0. Similarly for the second argument, consider the second equa- tion in (6): dIH (t2) dt = λH S (t2) + κψIHC (t2) − θ1λC IH (t2) − K1 IH (t2), = λH S (t2) + κψIHC (t2), hence IH (t1) = IH (θ) ∫ t 0 exp (λH S (t2) + κψIHC (t2))du, > 0. This also contradicts the earlier assumption that IH (t2) = 0. Therefore, there is no such time t2 exists. Hence IH (t) > 0 for all t > 0. Similar arguments can be applied to other two equations to show that IC > 0 and IHC > 0 for all t > 0. � Theorem 2. The biologically-feasible region Ω is positively- invariant for the model (6), where Ω = { (S, IH, IC, IHC ) ∈ R4+ : S + IH + IC + IHC ≤ Π µ } . 3 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 4 Table 1: Table 1: Description of model (6) parameters Parameter Interpretation Π Recruitment rate of susceptible humans by birth ψ Recovery rate of HCV λH HIV-only Infection rate λC HCV-only Infection rate λHC HIV and HCV simultaneous infection rate τ1, τ2 Latent period for HIV, HCV respectively βH Contact rate for HIV-only infection βC Contact rate for HCV-only infection βHC Contact rate for HIV and HCV simultaneous infection η > 1 Assumed increase in infectiousness of dually infected to single infection 0 < κ < 1 Modification parameter for reduced dually infected to single infection µ Natural death rate in all classes δH, δC, δHC HIV, HCV, HIV/HCV induced death rates θ1 ≥ 1,θ2 ≥ 1 Modification parameters for singly-infected to wholly susceptibles Proof. Adding the equations in (6) gives dN(t) dt = Π −µN(t) − [δH IH (t) + δC IC (t) + δHC IHC (t)], (8) so that, dN(t) dt ≤ Π −µN(t), (9) with positivity of IH, IC, IHC . It follows from (9), using Gron- wall lemma [39], that N(t) ≤ N(0)e−µ(t) + Π µ [ 1 − e−µ(t) ] . So that, if N(0) ≤ Π µ , then N(t) ≤ Π µ . Therefore, Ω is positively invariant. � 3. Existence and stability of equilibria At equilibrium, equations in (6) are equated to zero. In the absence of diseases, the DFE is obtained to be E0 = (S ∗, I∗H, I ∗ C, I ∗ HC ) = ( Π µ , 0, 0, 0 ) . (10) It should be noted that at equilibrium, IH (t−τ) = IH (t) = I∗H = 0, IC (t −τ) = IC (t) = I∗C = 0, and IHC (t) = IHC (t −τ) = I ∗ HC = 0. However, when there are diseases in the community, IH (t − τ) = IH (t) = I∗∗H , IC (t − τ) = IC (t) = I ∗∗ C and IHC (t) = IHC (t − τ) = I∗∗HC . Hence the following gives the implicit values for the endemic equilibrium: E∗∗ = (S ∗∗, I∗∗H , I ∗∗ C , I ∗∗ HC ), (11) where S ∗∗ = Π + ψI∗∗C λ∗∗H + λ ∗∗ C + λ ∗∗ HC + µ , I∗∗H = λ∗∗H S ∗∗ + κψI∗∗HC θ1λ ∗∗ C + K1 , I∗∗C = S ∗∗λ∗∗C θ2λ ∗∗ H + K2 , I∗∗HC = λ∗∗HC + θ1λ ∗∗ C I ∗∗ H + θ2λ ∗∗ H I ∗∗ C K3 , λ∗∗H = βH e−µτ [ I∗∗H + ηI ∗∗ HC ] N∗∗ , λ∗∗C = βC e−µτ [ I∗∗C + ηI ∗∗ HC ] N∗∗ , and λ∗∗HC = βHC e−µτI∗∗HC N∗∗ . (12) The local asymptotic stability of a given equilibrium of the model (6), is determined by linearizing the system about the equilibrium point [40, 41], and showing that all the roots of the transcendental polynomial have negative real parts. Thus linearizing the system (6) about the following variables: S (t), IH (t), IC (t), IHC (t), S (t −τ), IH (t −τ), IC (t −τ), IHC (t −τ) yields  dŜ (t) dt dÎH (t) dt dÎC (t) dt dÎHC (t) dt  = J  Ŝ (t) ÎH (t) ÎC (t) ÎHC (t) Ŝ (t −τ) ÎH (t −τ) ÎC (t −τ) ÎHC (t −τ)  , (13) where Ŝ (t) = S (t) − S ∗, ÎH (t) = IH (t) − I∗H , ÎC (t) = IC (t) − I ∗ C , ÎHC (t) = IHC (t) − I∗HC , J = [ J1 J2 ] , and J1 =  a1 a2 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15  , 4 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 5 J2 =  0 −βH e −µτS N − βC e−µτS N − e−µτ(βHη+βCη+βHC )S N 0 βH e −µτS N − θ1βC e−µτIH N e−µτ(βHηS−θ1βCηIH ) N 0 −θ2βH e −µτIC N βC e−µτS N e−µτ(βCηS−θ2βHηIC ) N 0 θ2βH e −µτIC N θ1βC e−µτIH N e−µτ(βHC S +θ1βCηIH +θ2βHηIC ) N  , with a1 = a2 − (λH + λC + λHC + µ), a2 = (λH + λC + λHC ) S N , a3 = −(βHη + βCη + βHC − λH − λC − λHC ) S N , a4 = −λH S N + λH + θ1λC IH N , a5 = −λH S N − θ1λC + θ1λC IH N − K1, a6 = −λH S N + θ1λC IH N , a7 = a6 + βHη S N − θ1βCη IH N + κψ, a8 = a9 + λC, a9 = −λC S N +θ2λH IC N , a10 = a9+θ2λH, a11 = a9+βCη S N−θ2βHη IC N , a12 = λHC S N +λHC−θ1λC IH N −θ2λH IC N , a13 = −λHC S N +θ1λC−θ1λC IH N − θ2λH IC N , a14 = λHC S N + θ2λH − θ1λC IH N − θ2λH IC N and a15 = −λHC S N −θ1λC IH N + θ1βCη IH N + θ2βHη IC N + βHC S N + θ2λH IC N − K3. Assuming X̂(t) = Dezt, where X̂(t) = (Ŝ , ÎH, ÎC, ÎHC ), is a solution with D = (d1, ...d4), be constants and z is a complex number. Equation (13) can be express as:[ J1 + e −zτJ2 ] Ĵ = 0, where Ĵ = [d1ezt d2ezt d3ezt d4ezt]T and 0 is 4 × 1 zero matrix. Therefore for nonzero solution, the transcendental equation at any equilibrium is given by: det(J1 + e −zτJ2) = 0. (14) However, before we determine the stabilities of equilibria, it is important to calculate the basic reproduction number of the model, denoted by R0. It can be recall that, susceptible hu- mans acquire HIV, HCV infections and HIV-HCV simultane- ously following effective contact with infected individuals with these diseases. It follows that, the number of HIV infectives generated by an infected human (near the DFE) is given by the product of the infection rate ( βHN∗ ), the probability that an infected human survives natural death (e−µτ) and the average duration of stay in the infected HIV class ( 1K1 ). Hence, the aver- age number of HIV infections generated by infected individual is given by RH = βH e−µτS ∗ K1 N∗ . (15) Similarly, the number of HCV infections generated by an infected individual with HCV (near the DFE), is given by the product of the infection rate ( βCN∗ ), the probability that an in- fected human survived natural death (e−µτ) and the average du- ration of stay in the infected HCV class ( 1K2 ). Therefore, the average number of HCV infected generated by infectious indi- vidual is given by RC = βC e−µτS ∗ K2 N∗ . (16) Furthermore, the number of HIV-HCV simultaneous infec- tions generated by an infected individual with the two diseases (near the DFE), is given by the product of the infection rate of infected HIV-HCV ( βHCN∗ ), survival probability (e −µτ) and the av- erage duration of stay in the class IHC ( 1 K3 ). Thus, the average number of HIV-HCV infections generated by infected individ- ual is given by RHC = βHC e−µτS ∗ K3 N∗ . (17) Finally, the maximum of the three sub thresholds in (15), (16) and (17) gives the value of R0, i.e. the average number of new infections generated by infected individuals (with HIV, HCV and HIV-HCV) is given by R0 = max { βH e−µτ K1 , βC e−µτ K2 , βHC e−µτ K3 } , = max{RH,RC,RHC}, (18) noting that at DFE, S ∗ = N∗ = Π µ . It is worth stating here that, the formulation of R0 in (18) is derived from the procedure in [12] used in a model that have more than one strains/diseases, to be the maximum of various sub thresholds for the infections. Moreover, similar approaches are also applied in [46, 42]. 3.1. Local asymptotic stability of DFE At DFE, we claim the stability properties as follows: Theorem 3. The DFE E0, is locally asymptotically stable (LAS) in Ω whenever R0 < 1 for all τ ≥ 0 and unstable if R0 > 1. Proof. At DFE, the transcendental equation (14) is reduced to G(z) = (z + µ)(z + K1 −βH e −τ(z+µ)) (19) ×(z + K2 −βC e −τ(z+µ)) (20) ×(z −βHC e −τ(z+µ) + K3) = 0. (21) Expanding equation (19), we have G(z) =z4 + P1z 3 + P2z 2 + P3z + P4 = e−zτ[Q1z 3 + Q2z 2 + Q3z + Q4], (22) where P1 = µ + K1 + K2 + K3, P2 = [µK1 + (µ + K1)K2 + µ + K1 + K2]K3(1 −RHC ), P3 = [(µ + K1)K2 + µK1 + µK1 K2]K3(1 −RHC ), P4 = µK1 K2 K3(1 −RHC ), Q1 = (βH + βC + βHC )e −µτ, Q2 = [µ(K2βH + βC ) + K1βC (1 −RH )e −zτ + βH + βC ](1 −RHC )K3e −µτ, Q3 = [µK1βC (1 −RH ) + (K2 + 1)βH + K2βH (1 −RC )e −zτ](1 −RHC )K3e −µτ, Q4 = µe −µτ[K1βC + K2βH (1 −RC )e −zτ]K3(1 −RHC ). First we determine the stability of E0 when τ = 0. From (19), when τ = 0, it is obvious that the eigenvalue −µ is negative, while βH − K1,βC − K2 and βHC − K3 are negative whenever RH < 1,RC < 1 and RHC < 1 or in general if R0 < 1. Hence, E0 at τ = 0 is stable when R0 < 1. 5 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 6 Further, we determine the stability of E0 when τ > 0. Here we investigate whether there is a root z = iy of (19), y ∈ R+, that may cross the imaginary axis and cause stability switches. Substituting z = iy in (19) or equivalently in (22), we find y that must satisfy G1(y) =| P1(iy) | 2 − |Q1(iy) | 2, (23) where P1(iy) = (y4 − P1iy3 − P2y2 + P3iy + P4), Q1(iy) = (−Q1iy3 − Q2y2 + Q3iy + Q4). However, (22) is implicit in e−zτ, therefore the process of separating real, imaginary parts and squaring (23) may not lead to explicit expression in y. Without loss of generality, one can see from (19) that G(0) = µK1 K2 K3[(1 −RH )(1 −RC )(1 −RHC )] > 0 if R0 < 1 and G(y) = +∞ as y → ∞. Therefore, when τ > 0, equation (19) has no positive root on (0, +∞), hence the DFE is LAS for all τ ≥ 0 if R0 < 1 and unstable if R0 > 1. The biological implication of this result is that the diseases can be completely eliminated whenever the initial population is close to the basin of attraction of the DFE. � 3.2. Global asymptotic stability of DFE To ensure complete eradication of the diseases in the popu- lation irrespective of the population size started with, we present and prove the following result. Theorem 4. The DFE E0, is globally asymptotically stable (GAS) in Ω whenever R0 < 1 for all delay τ ≥ 0. Proof. The global properties can be established using a Com- parison Theorem [23] and the approach in [43]. The equations of the three infected components in model (6) can be expressed using matrix-vector form as dIH (t) dt dIC (t) dt dIHC (t) dt  = (F −V1)  IH (t −τ) IC (t −τ) IHC (t −τ)  −V2  IH (t) IC (t) IHC (t)  − ( 1 − S N )  βH e−µτ 0 βHη 0 βC e−µτ βHη 0 0 βHC   IH (t −τ) IC (t −τ) IHC (t −τ)  , (24) where F =  βH e−µτ 0 0 0 βC e−µτ 0 0 0 βC  , V1 =  0 0 −κψ 0 0 0 0 −θ1λC −θ2λH  and V2 =  K1 + θ1λC 0 0 0 K2 + θ2λH 0 0 0 K3  are the matrices for new infection terms and transitions re- spectively. Since S ≤ N at every time t in Ω and by letting V = V1 + V2, it follows from (24) that dIH (t) dt dIC (t) dt dIHC (t) dt  ≤ (F −V)  IH (t) IC (t) IHC (t)  . (25) It can be shown that the spectral radius of FV−1 gives the value of R0 in (18). Furthermore, since the eigenvalues of the matrix (F − V) have negative real parts if R0 < 1 [43, 44], conse- quently, the linearized differential system (25) is stable when- ever R0 < 1. As a result, lim t→∞ (IH (t), IC (t), IHC (t)) → (0, 0, 0). (26) Therefore, from the first equation in (6) and substituting IH (t) = IC (t) = IHC (t) = 0 from (26), for any t, it follows that dS dt = Π −µS (t), (27) so that S (t) = Π µ . (28) Thus, lim t→∞ (S (t), IH (t), IC (t), IHC (t)) = ( Π µ , 0, 0, 0 ) . (29) By LaSalles invariance principle [45] E0 is the largest in- variance set in Ω, hence is GAS whenever R0 < 1. � 3.3. Stability analysis of endemic equilibrium Here, due to the complexity of transcendental equation, we give conditions for stability of any unique endemic equilibrium (EE), whenever it exists. From (12) the transcendental equation (14) is simplified to G2(z)= det  z − a1 −a2 + βH e −µτS N −a2 + βC e −µτS N −a3 −a4 z − a5 − βH e −µτS N −a6 θ1βC e −µτIH N a7 −a8 −a9 + θ2βH e −µτIC N z − a10 βC e −µτS N a11 a12 −a13 θ2βH e −µτIC N −a14 − θ1βC e −µτIH N z − a15  = 0, = P2(z) + Q2(z)e −zτ = 0, (30) where P2(z) and Q2(z) are polynomials of degree 4 and 3 respectively, with real coefficients and hence have no common imaginary roots whenever R0 > 1. Furthermore, substituting z = iy, as purely imaginary root, we define G2(y) =| P2(iy) | 2 − |Q2(iy) | 2, (31) which is a polynomial of degree 8 whose leading coefficient is positive. Suppose (31) has at least zero positive root, then the following stability result for E∗∗ can therefore be stated Theorem 5. If the polynomial G2(y) has (a) no positive root, then E∗∗ is locally asymptotically stable for all delays τ ≥ 0 whenever R0 > 1. 6 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 7 (b) at least one simple positive root, then as delay increases there will be n ∈ Z+ number of stability switches for fixed parameter values and the endemic equilibrium E∗∗ is locally asymptotically stable for 0 ≤ τ < τ∗ if R0 > 1, here, τ∗ = cot−1 ( − P2reQ2re +P2imQ2im −P2reQ2im +P2imQ2re ) y and the subscripts represents the real and imaginary parts of P2 and Q2. 3.3.1. Interior endemic equilibrium Here, we consider one of the interior equilibria with simul- taneous infection of HIV and HCV only. In the absence of HIV and HCV only infections, equation (12) will be reduced to E∗∗3 = (S ∗∗, 0, 0, S ∗∗(RHC − 1)), where S ∗∗ = ΠRHC e−µτβHC (RHC − 1) + µRHC , (32) provided RHC , e−µτβHC µ+e−µτβHC . Similarly, the transcendental equa- tion (30) is simplified to be G2(z) = P2(z) −Q2(z)e −zτ, where P2(z) = z 2 − ( p1 + p4)z + (p1 p4 − p2 p3); Q2(z) = −p5z + p1 p5 − p3 p5, (33) with p1 = − µK23 (RHC − 1) 2 e−µτβHC −µ; p2 = − µK3(RHC − 1) e−µτβHC ; p3 = µK23 (RHC − 1) 2 e−µτβHC ; p4 = − µK3(RHC − 1) e−µτβHC − K3; p5 = K3. When τ = 0, (33) yields G2(z) = z 2 − ( p1 + p4 + p5)z + ( p1 p4 − p2 p3) + ( p1 p5 + p3 p5), (34) so that after simplification, we have − ( p1 + p4 + p5) = K3(RHC − 1) βHC × [1 + K3(RHC − 1)] + µ, ( p1 p4 − p2 p3) + (p1 p5 + p3 p5) = δHC K 2 3 (RHC − 1) 2 + µK3(RHC − 1). (35) From (35), it can be seen that all the coefficients of G2(z) are positive whenever RHC > 1. Therefore, G2(z) is stable (all the roots have negative real parts). Furthermore, if τ > 0, we let z = iy be the root of G2(z) in (34) then from (31), separating real, imaginary parts and squar- ing gives F2(y) = y 4 + (−2A22 + A 2 11 − A 2 33)y 2 + A222 − A 2 44 = 0, (36) where A11 = p1 + p4, A22 = p1 p4 − p2 p3, A33 = p5 and A44 = p1 p5 + p3 p5. Here, − 2A22 + A 2 11 − A 2 33 = K23 (RHC − 1) 4 e−2µτβ2HC + 2µK3(RHC − 1)2 e−µτβHC + 2K33 (RHC − 1) 2 e−2µτβ2HC + 2µK33 (RHC − 1) e−µτβHC + K23 [I ∗∗2 HC − 1] S ∗∗2 + µ2 + K3, (37) and A222 − A 2 44 = (A22 − A44)(A22 + A44), = [ δHC (RHC − 1)2 e−µτβHC + K3(RHC − 1) + µ + µK3 ] [ K3(RHC − 1)2 e−µτβHC + µK3(RHC − 1) e−µτβHC + 2µK3 ] . (38) From (37) and (38), if RHC > 1, so also I∗∗2HC > 1, therefore F2(y) in (36) will have no positive root y that can cross the imaginary root as the delay is increased. Hence, the interior equilibrium E∗∗3 is absolutely stable for all delay, τ ≥ 0. Therefore, we claim the following result: Theorem 6. The endemic equilibrium E∗∗3 , (32), is locally asymp- totically stable for any delay, τ ≥ 0 when RHC > 1. 3.4. Analysis for impact of delay The impact of the time delay is analyzed using threshold approach to show whether increase (decrease) of the delay will have effect on the infectivity of the two diseases in the com- munity. We use the threshold quantity, R0, basic reproduction number, R0 = max{RH,RC,RHC}. (39) as a function of the delay τ. Since R0 = 1 is a threshold value for the infection in the population, we determine the critical delay value τcrit above (below) which the diseases have effect as follows: βHβCβHC e−3µτ K1 K2 K3 = 1, e−3µτ = K1 K2 K3 βHβCβHC , τ = τcrit = ln ( K1 K2 K3 βHβCβHC ) −3µ . (40) 7 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 8 0 500 1000 1500 0 0.5 1 1.5 2 2.5 x 10 4 Time (days) P op ul at io ns o f S , I H + I C + I H C S I H +I C +I HC (a) 0 500 1000 1500 0 5000 10000 15000 Time (days) P op ul at io ns o f S , I H , I C , I H C S I H I C I HC (b) 0 500 1000 1500 0 5000 10000 15000 Time (days) P op ul at io ns o f S , I H , I C , I H C S I H I C I HC (c) 0 500 1000 1500 0 0.5 1 1.5 2 x 10 4 Time (days) P op ul at io ns o f S , I H , I C , I H C S I H I C I HC (d) Figure 1: Numerical simulations of model (6), with parameter values from Table 2 with (a) the GAS of DFE, with βH = 0.002 = βC, βHC = 0.0025, τ = 3. In (b) and (c) βH = 0.002 = βC, βHC = 0.0025, τ = 1 and τ = 10 respectively (d) βH = 0.04 = βC, βHC = 0.045 and τ = 50. Table 2: Baseline values for the parameters of model Eq. 6 Parameter Baseline values References Π 300 day−1 Assumed ψ 6.8 × 10−5 day−1 [46] βH Variable Assumed βC Variable Assumed βHC Variable Assumed η 2 Assumed κ 0.013 Assumed µ 170 day −1 [46] δH 9.12 × 10−4 day−1 [47] δC 9.5 × 10−5 day−1 [46] δHC 9.32 × 10−4 day−1 Assumed θ1, θ2 2, 2 Assumed 8 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 9 0 200 400 600 800 1000 0 100 200 300 400 500 Time (days) In fe ct ed H IV p op ul at io n 20 40 60 200 300 400 (a) 0 500 1000 1500 0 0.5 1 1.5 2 x 10 4 Time (days) P op ul at io ns o f S , I H , I C , I H C S I H I C I HC (b) 0 500 1000 1500 0 2000 4000 6000 8000 10000 Time (days) C um ul at ve n ew c as es τ=1 τ=15 τ=25 τ=35 (c) Figure 2: Simulations of model (6), with parameter values as shown in Table 2, in (a) showing no oscillation with βH = 0.02 = βC, βHC = 0.025, τ = 30 (b) displaying the interior equilibrium, with βH = 0.0002 = βC, βHC = 0.025, RHC = 1.62, while in (c) the effect of time delay on the infectivity of the viruses with βH = 0.04 = βC, βHC = 0.045, using time delays τ = 1, 15, 25 and τ = 35. The rate of change of R0 with respect to τ is given by ∂R0 ∂τ = ∂ ∂τ max (RH,RC,RHC ) , = ∂ ∂τ ( βHβCβHC e−3µτ K1 K2 K3 ) , = −3µ βHβcβHC e−3µτ K1 K2 K3 . (41) Since the derivative in (41) is negative, it indicates that R0 decreases with increase of the time delay τ and vice varsa. This implies that, increasing the delay beyond the threshold value, with fix parameter values, will results in decreasing the number of infected individuals with the two diseases and those that are simultaneously infected. This result can be summarized as Proposition 1. The increase in time delay in the model (1) will reduces the infectivity of the two diseases whenever the delay is greater than the threshold value τcrit. 4. Numerical simulations Numerical experiments are conducted to illustrate the dy- namics as described in theoretical results and display the impact of delay in the model. Figure 1 (a) depicts the global asymp- totic stability for the DFE using parameter values from Table 2, except for βH = βC = 0.002, βHC = 0.0025, τ = 3 and different initial conditions so that R0 = 0.16. It can be seen that the sus- ceptible and total populations of infectives converges to Π µ and zero asymptotically, respectively. In Figures 1 (b)–(d), the local stability of endemic equilibrium is displayed using values from Table 2 except for in (b) τ = 1, βH = βC = 0.02, βHC = 0.025, (c) τ = 10 βH = βC = 0.002, βHC = 0.0025, and (d) τ = 30 , βH = βC = 0.04, βHC = 0.045, as proved in Theorem 5(a), so that R0 = 1.60, R0 = 1.42 and R0 = 1.45 respectively. In Fig- ure 2 (a), it is shown that the length of delay (τ = 30) doesn’t cause any oscillation in the stability of the endemic equilibria as observed in some delay models. Figure 2 (b), illustrates the stability if one of the interior equilibrium, E∗∗3 using parameter values in Table 2 while βH = βC = 0.0002, βHC = 0.025, so that RHC = 1.62. In this case, one or simultaneous infection will oc- cur. The effect of delay is illustrated in Figure 2 (c) to support Proposition 1 in which increasing the delay, τ from 1, 15, 25 to 35, caused remarkable decrease in the total number of infected individuals with HIV and HCV. 5. Concluding remarks In this work, a delay model of simultaneous infection of HIV and HCV is formulated and dynamically analyzed. The novelty (results) of the research is summarized and discussed below: (i) Basic properties (existence, boundedness and positivity of solution)of the model are stated and proved as Theo- rems 1 and 2. 9 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 10 (ii) The basic reproduction threshold is systematically ob- tained as the maximum of subthreshold values of the in- dividual viruses; a threshold parameter above which the viruses will persist in the population. (iii) A disease free equilibrium is found to be globally asymp- totically stable when the basic reproduction threshold is less than unity for any length of time delay. Under this case, the diseases will die out in the population whenever the basic reproduction number is brought and maintained below one irrespective of the initial populations started with. This is shown in Figure 1(a). (iv) However, if the basic reproduction threshold is greater than unity, endemic equilibrium exists which is shown to be locally asymptotically stable for all values of the incu- bation period (delay) as shown in Figures 1(b – d). This implies that whenever the initial population is within the basin of attraction of the endemic equilibrium and the ba- sic reproduction number is greater than one, the diseases will persist in the population no matter the length of de- lay. (v) Increasing the length of time delay is observed to be de- creasing the number of cumulative infectious people when the delay is above a critical value. It is worth remarking here that although the time delay has no effect on the sta- bility of equilibria however, it do affect the infectivity of the viruses as shown in the formulation of basic repro- duction number. The implication of this result is, if the incubation period can be extended by say intervention, the number of infected people will be reduced. (vi) Numerical simulations, using data from the literature are used to illustrate our results. Acknowledgments The authors are thankful to Department of Mathematics and Applied Mathematics, University of Pretoria, South Africa, where this research was initiated. We are equally grateful to the au- thorities of Bayero University, Kano for research leave, the anony- mous reviewers and the handling Editor, Tolulope Latunde (P.hD.) for their constructive efforts. References [1] World Health Organization: HIV and hepatitis coinfections. http:// www.who.int. Accessed: November, 2018. [2] C. Castillo-Chavez, J. X. Velasco-Hernandez, “On the relationship be- tween evolution of virulence and host demography”, Journal of Theoreti- cal Biology, 192 (1998) 437. [3] S. Levin, D. Pimentel, “Selection of intermediate rates of increase in parasite-host systems”, The American Naturalist, 117 (1981) 308. [4] M. A. Nowak, R. M. May, “Superinfection and the evolution of parasite virulence”, Royal Society London series B, 255 (1994) 81. [5] S. Alizon, “Parasite co-transmission and the evolutionary epidemiology of virulence”, Evolution, 67 (2013) 921. [6] R. M. Anderson, R. M. May, Infectious Diseases of Humans Dynamics and Control, Oxford University Press, Oxford, UK 1991. [7] R. A. Frederick, C. B. Robert, “The dynamics of simultaneous infections with altered susceptibilities”, Theoretical Population Biology. 40 (1991) 369. [8] R. M. May, M. A. Nowak, “Coinfection and the evolution of parasite virulence”, Royal Society London, seies B, 261 (1995) 209. [9] Z. Mukandavire, A. B. Gumel, W. Garira, J. M. Tchuenche, “Mathemati- cal analysis of a model for HIV-malaria co-infection”, Mathematical Bio- science and Engineering, 6 (2009) 333. [10] A. Y. C. Sanchez, M. Aerts, Z. Shkedy, P. Vickerman, F. Faggiano, “A mathematical model for HIV and Hepatitis C co-nfection and its assess- ment from a statistical perspective”, Epidemics, 5 (2013) 56. [11] E. W. Seabloom, P. R. Hosseini, A. G. Power, E. T. Borer, “Diversity and Composition of Viral Communities: Coinfection of Barley and Cereal Yellow Dwarf Viruses in California Grasslands”, The American Natural- ist, 173 (2009) 79. [12] M. van Baalen, M. W. Sabelis, “The dynamics of multiple infection and the evolution of virulence”, The American Naturalist, 146 (1995) 881. [13] A. Y. C. Sanchez, M. Aerts, Z. Shkedy, P. Vickerman, F. Faggiano, “A mathematical model for HIV and hepatitis C co-infection and its assess- ment from a statistical perspective”, Epidemics, 5 (2013) 56. [14] X-S. Zhang, K. F. Cao, “The impact of coinfections and their simulta- neous transmission on antigenic diversity and epidemic cycling of in- fectious diseases”, Biomed Research Internatioal, (2014), 375 doi : 10.1155/2014/375862. [15] X-S. Zhang , A. De Angelis, P. J. White, A. Charlett, R. G. Pebody, J. McCauley, “Co-circulation of influenza A virus strains and emergence of pandemic via reassortment: the role of cross-immunity”, Epidemics, 5 (2013) 20. [16] O. Balmer, M. Tanner, “Prevalence and implications of multiple-strain infections”, Lancet Infectious Diseases, 11 (2011) 868. [17] R. Ridzon, K. Gallagher, C. Ciesielski, E. E. Mast, M. B. Ginsberg, B. J. Robertson, C. C. Luo, A. DeMaria, “Simultaneous transmission of human immunodeficiency virus and hepatitis C virus from a needle-stick injury”, The New England Journal of Medicine, 336 (1997) 919. [18] S. Gorman, N. L. Harvey, D. Moro, M. L. Lloyd, V. Voigt, L. M. Smith, M. A. Lawson, G. R. Shellam, “Mixed infection with multiple strains of murine cytomegalovirus occurs following simultaneous or sequential infection of immunocompetent mice”, Journal of General Virology, 87 (2006) 1123. [19] G. Ippolito, V. Puro, N. Petrosillo, G. De Carli, G. Miche- loni, E. Magliano, “Simultaneous infection with HIV and hepati- tis C virus following occupational conjunctival blood exposure”, The Journal of the American Medical Association, 280 (1998), doi : 10.1001/ jama.280.1.28. [20] W. Liu, Z. D. Li, F. Tang, M. Wei, Y. Tong, L. Zhang, Z. Xin, “Mixed infections of pandemic H1N1 and seasonal H3N2 viruses in outbreak”, Clinical Infectious Diseases, 50 (2010) 1359. [21] S. Telfer, X. Lambin, R. Birtles, P. Beldomenico, S. Burthe, S. Paterson, M. Begon, “Species interactions in a parasite community drive infection risk in a wildlife population”, Science 330 (2010) 243. [22] P. Schmid-Hempel, Evolutionary Parasitology: The Integrated Study of Infections, Immunology, Ecology and Genetics, Oxford University Press, Oxford, UK, 2011. [23] H. L. Smith, P. Waltman, The Theory of the Chemostat, Cambridge Uni- versity Press, 1995. [24] X. Zhang, R. Pebody, D. De Angelis, P. J. White, A. Charlett, J. W. Mc- Cauley, “The possible impact of vaccination for seasonal influenza on emergence of pandemic influenza via reassortment”, PLoS ONE, 9 (2014) 12. doi : 10.1371/ journal.pone.0114637. [25] A. Dutta, P. K. Gupta, “A mathematical model for transmission dynam- ics of HIV/AIDS with effect of weak CD4+ T cells”, Chinese Journal of Physics, 56 (2018) 1045. [26] W. Jia, J. Weng, C. Fang, Y. Li, “A dynamic model and some strategies on how to prevent and control hepatitis c in mainland China”, BMC In- fectious Diseases, 19 (2019) 1. [27] D. A. Maimunah, “Mathematical model for HIV spreads control program with ART treatment”, Journal of Physics: Conference Series, 974 (2018) 1. [28] M. D. Miller-Dickson, V. A. Meszaros, S. Almagro-Moreno, C. B. Og- bunugafor, “Hepatitis C virus modelled as an indirectly transmitted in- fection highlights the centrality of injection drug equipment in disease 10 Hassan & Hussaini / J. Nig. Soc. Phys. Sci. 3 (2021) 1–11 11 dynamics”, Journal of Royal Society Interface, 16 (2019) 1. [29] G. Huang, Y. Takeuchi, “Global analysis on delay epidemiological dy- namic models with nonlinear incidence”, Journal of Mathematical Biol- ogy, 63 (2011) 125. [30] Y. N. Krychko, K. B. Blyuss, “Global properties of a delayed SIR model with temporary immunity and nonlinear incidence rate”, Nonlinear anal- ysis: Real world applications, 6 (2005) 495. [31] M. A. Safi, A. B. Gumel, “The effects of incidence functions on the dy- namics of a quarantine/isolation model with time delay”, Nonlinear Anal- ysis: Real World Applications, 12 (2011) 215. [32] Y. Wang, H. Wang, W. Jiang, “Stability switches and global Hopf bifurca- tion in a nutrient-plankton model”, Nonlinear Dynamics, 78 (2014) 981. [33] Y. Zhao, Z. Xu, “Global dynamics for a delayed Hepatitis C virus infec- tion model”, Electronic Journal of differential equations, 2014 (2014) 1. [34] G. Huang, Y. Takeuchi, W. Ma, D. Wei, “Global stability for delay SIR & SEIR epidemic models with nonlinear incidence rate”, Bulletin of Math- ematical Biology, 72 (2010). [35] S. Ruan, D. Xiao, J. C. Beier, “On the delayed Ross–Macdonald model for malaria transmission”, Bulletin of Mathematical Biology, 70 (2008) 4. [36] S. Barnerjee, R. Keval, S. Gakkhar, “Influence of intracellular delay on the dynamics of Hepatitis C virus”, International Journal of Applied and Computational Mathematics, 4 (2018) 1. [37] S. A. Gourley, Y. Kuang, J. D. Naggy, “Dynamics of a delay differen- tial equation model of Hepatitis B virus infection”, Journal of Biological Dynamics, 2 (2008) 140. [38] J. K. Hale, Introduction to Functional Differential Equations, Springer- Verlag, New York, 1993. [39] T. H. Gronwall, “Note on the derivatives with respect to a parameter of the solutions of a system of differential equations”, Annals of Mathematics, 20 (1919) 292. [40] E. Beretta, Y. Kuang, “Geometric stability switch criteria in delay differ- ential systems with delay dependent parameters”, SIAM Journal of Math- ematical Analysis, 33 (2002) 1144. [41] K. L. Cooke, P. van den Driessche, “On zeroes of some transcendental equations”, Funkcialaj Ekvacioj, 29 (1986) 77. [42] N. Hussaini, J. M.-S. Lubuma, K. Barley, A. B. Gumel, “Mathematical analysis of a model for AVL-HIV co-endemicity”, Mathematical Bio- science, 271 (2016) 80. [43] O. Y. Sharomi, Mathematical Analysis of Dynamics of Chlamydia tra- chomatis, PhD thesis, Department of Mathematics, University of Mani- toba, Canada, 2010. [44] P. van den Driessche, J. Watmough, “Reproduction numbers and sub- threshold endemic equilibria for compartmental models of disease trans- mission”, Mathematical Bioscience, 180 (2002) 29. [45] J. P. LaSalle, The Stability of Dynamical Systems, Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, 1976. [46] R. Shi, Y. Cui, “Global analysis of a mathematical model for Hepatitis C virus transmissions”, Virus Research. 217 (2016) 17. [47] Z. Mukandavire, W. G.Chiyaka, “Asymptotic properties of an HIV/AIDS model with a time delay”, Journal of Mathematical Analysis and Appli- cations 330 (2007) 933. 11