Review

Novel Approaches to Pediatric Cancer: Immunotherapy

  • From the early 20th century, immunotherapy has been studied as a treatment modality for cancers, including in children. Since then, developments in monoclonal antibodies and vaccine therapies have helped to usher in a new era of cancer immunotherapeutics. However, efficacy of these types of therapies has been limited, mostly in part due to low tumor immunogenicity, cancer escape pathways, and toxicities. As researchers investigate the cellular and molecular components of immunotherapies, mechanisms to improve tumor specificity and overcome immune escape have been identified. The goal of immunotherapy now has been to modulate tumor escape pathways while amplifying the immune response by combining innate and adaptive arms of the immune system. Although several limiting factors have been identified, these recent advances in immunotherapy remain at the forefront of pediatric oncologic therapeutic trials. Immunotherapy is now coming to the forefront of precision treatment for a variety of cancers, with evidence that agents targeting immunosuppressive mechanisms for cancer progression can be effective therapy [1-3]. In this review, we review various types of immunotherapy, including the cellular biology, limitations, recent novel therapeutics, and the application of immunotherapy to pediatric oncology.

    Citation: Payal A. Shah, John Goldberg. Novel Approaches to Pediatric Cancer: Immunotherapy[J]. AIMS Medical Science, 2015, 2(2): 104-117. doi: 10.3934/medsci.2015.2.104

    Related Papers:

    [1] Churni Gupta, Necibe Tuncer, Maia Martcheva . Immuno-epidemiological co-affection model of HIV infection and opioid addiction. Mathematical Biosciences and Engineering, 2022, 19(4): 3636-3672. doi: 10.3934/mbe.2022168
    [2] Shengqiang Liu, Lin Wang . Global stability of an HIV-1 model with distributed intracellular delays and a combination therapy. Mathematical Biosciences and Engineering, 2010, 7(3): 675-685. doi: 10.3934/mbe.2010.7.675
    [3] Georgi Kapitanov, Christina Alvey, Katia Vogt-Geisse, Zhilan Feng . An age-structured model for the coupled dynamics of HIV and HSV-2. Mathematical Biosciences and Engineering, 2015, 12(4): 803-840. doi: 10.3934/mbe.2015.12.803
    [4] Nawei Chen, Shenglong Chen, Xiaoyu Li, Zhiming Li . Modelling and analysis of the HIV/AIDS epidemic with fast and slow asymptomatic infections in China from 2008 to 2021. Mathematical Biosciences and Engineering, 2023, 20(12): 20770-20794. doi: 10.3934/mbe.2023919
    [5] Nicolas Bacaër, Xamxinur Abdurahman, Jianli Ye, Pierre Auger . On the basic reproduction number $R_0$ in sexual activity models for HIV/AIDS epidemics: Example from Yunnan, China. Mathematical Biosciences and Engineering, 2007, 4(4): 595-607. doi: 10.3934/mbe.2007.4.595
    [6] Kazunori Sato . Basic reproduction number of SEIRS model on regular lattice. Mathematical Biosciences and Engineering, 2019, 16(6): 6708-6727. doi: 10.3934/mbe.2019335
    [7] Georgi Kapitanov . A double age-structured model of the co-infection of tuberculosis and HIV. Mathematical Biosciences and Engineering, 2015, 12(1): 23-40. doi: 10.3934/mbe.2015.12.23
    [8] Xue-Zhi Li, Ji-Xuan Liu, Maia Martcheva . An age-structured two-strain epidemic model with super-infection. Mathematical Biosciences and Engineering, 2010, 7(1): 123-147. doi: 10.3934/mbe.2010.7.123
    [9] Andrew Omame, Sarafa A. Iyaniwura, Qing Han, Adeniyi Ebenezer, Nicola L. Bragazzi, Xiaoying Wang, Woldegebriel A. Woldegerima, Jude D. Kong . Dynamics of Mpox in an HIV endemic community: A mathematical modelling approach. Mathematical Biosciences and Engineering, 2025, 22(2): 225-259. doi: 10.3934/mbe.2025010
    [10] Wenshuang Li, Shaojian Cai, Xuanpei Zhai, Jianming Ou, Kuicheng Zheng, Fengying Wei, Xuerong Mao . Transmission dynamics of symptom-dependent HIV/AIDS models. Mathematical Biosciences and Engineering, 2024, 21(2): 1819-1843. doi: 10.3934/mbe.2024079
  • From the early 20th century, immunotherapy has been studied as a treatment modality for cancers, including in children. Since then, developments in monoclonal antibodies and vaccine therapies have helped to usher in a new era of cancer immunotherapeutics. However, efficacy of these types of therapies has been limited, mostly in part due to low tumor immunogenicity, cancer escape pathways, and toxicities. As researchers investigate the cellular and molecular components of immunotherapies, mechanisms to improve tumor specificity and overcome immune escape have been identified. The goal of immunotherapy now has been to modulate tumor escape pathways while amplifying the immune response by combining innate and adaptive arms of the immune system. Although several limiting factors have been identified, these recent advances in immunotherapy remain at the forefront of pediatric oncologic therapeutic trials. Immunotherapy is now coming to the forefront of precision treatment for a variety of cancers, with evidence that agents targeting immunosuppressive mechanisms for cancer progression can be effective therapy [1-3]. In this review, we review various types of immunotherapy, including the cellular biology, limitations, recent novel therapeutics, and the application of immunotherapy to pediatric oncology.


    In the last 15 years, US deaths due to the opioid epidemic have quadrupled from nearly 12,000 in 2002 to more than 47,000 in 2017 [1]. In October 2017, the US Department of Health and Human Services declared the opioid crisis a national public health emergency [2]. The increase in injection drug use and reduction of behavioral inhibition have also contributed to the spread of infectious diseases, particularly HIV [3]. The two epidemics – opioid and HIV – are intertwined and modeling them in tandem will lead to understanding their interdependence [4].

    The HIV epidemic, which begun in 1981, has been modeled extensively both on immunological level and on epidemiological level. On within-host level typically models involve healthy and infected CD4 cells and the virus [5,6,7,8,9]. On the between-host scale, the very basic models include susceptible, infected and AIDS classes [10,11]. Multi-scale models of HIV, mostly of the nested type [12], have also received significant attention in the recent years [13,14,15,16].

    Modeling of the opioid epidemic is more recent and picked up with the increasing importance of the substance abuse disorder. Early models [17] assume opioid use can be modeled similarly to infectious disease spread and that assumption has been used in current models of heroin use [18,19,20,21]. The transition to opioid use typically starts from prescription drug misuse and modeling of that has also been drawing attention lately [22].

    Despite the importance of the HIV and the opioid epidemics and the clear interdependence between the two, relatively little modeling has been done at the interface of the two epidemics. The within-host modeling of the interplay between opioid and HIV has first drawn attention. Based on experiments in monkeys, Vaidya et al. [23,24,25,26] model the within-host dynamics of HIV and an opioid. The between host dynamics has been addressed even less. Duan et al. [27] models the two epidemics on population level and finds that the best control strategy should reduce the probability of opioid affected individuals getting HIV and should target the drug abuse epidemic. In this paper, we investigate the interplay of the two epidemics with a multi-scale model that incorporates both a within-host component and a between host component.

    The sexual contacts leading to HIV are not homogeneous across the population. Typically, few individuals in the population partake in a lot of contacts, while most of the members of the population partake in few contacts. This type of heterogeneity of contacts in HIV transmission is modeled by scale-free networks. Modeling infectious disease dynamics on networks has been also drawing attention in the recent years, resulting in a book devoted to this topic [28]. We use here a network modeling approach introduced in [29] which gives a closed form model equations [30]. In this modeling framework of networks, ODE and age structured models have been investigated but the only multi-scale model on networks in closed form seems to be [16]. Here we expand the modeling framework developed in [16] to include the opioid epidemic thus considering for the first time a multi-scale network model of two diseases. This model is in closed form and we are able to perform analysis on it.

    In Section 2 we introduce the multi-scale network model of HIV and opioid which consists of within-host component, between-host component and linking functions used to connect the two scales. In Section 3 we discuss the existence and stability of the disease-free equilibrium and compute an explicit form of the reproduction numbers. In Section 4 we focus on the existence and stability of the semi-trivial equilibria. In Section 4 we also compute the invasion numbers of the two epidemics. In Section 5 we perform simulations and consider different scenarios with differing parameter values. In Section 6 we summarize our results.

    We modify a well-known within-host model of HIV by explicitly including the opioid drug concentration $ C(\tau) $ and its impact on the average susceptibility of target cells. The model takes the following form.

    $ {dTdτ=sdTTk(C)ViT, dTidτ=k(C)ViTδiTi,dVidτ=NvδiTicVi,dCdτ=ΛdcC,
    $
    (2.1)

    Here $ T $ are the target cells, $ T_i $ are the infected target cells and $ V_i $ is the virus (HIV). Target cells are produced at rate $ s $ and cleared at rate $ d_T $. Infected cells die at a rate $ \delta_i $, releasing $ N_v $ viral particles at death. The clearance rate of the virus is denoted by $ c $. Opioid is taken at doses $ \Lambda $ is degraded at rate $ d_c $. Infection rate of target cells by HIV in the presence of opioid is given by

    $ k(C) = k_0+\frac{k_1 C(\tau)}{C_0 +C(\tau)}, $

    where $ C_0 $ is the half saturation constant, $ k_0 $ is the transmission coefficient in the absence of opioid and $ k_1 $ is a maximal increase in infection rate due to opioids.

    A shortcoming of this model is that it does not consider multiple infection routes and drug resistant viral strains. This problem could have been remedied if a model such as the one formulated in [31] was used instead, as the base model. But because of the network structure, the system is already complicated, and addition of strains and infection route would further complicate it. We believe addition of strains and infection route would not give any different insights with the network structures.

    To introduce the model, we define a complex network with size (number of nodes) $ N $ where each node is either occupied by an individual or vacant. The states for the epidemic transmission process on the network are divided into vacant state $ E $, susceptible state $ S $, opioid state $ U $, infected state $ V $, co-affected state $ i(\tau, t) $ and AIDS state $ A $. Nodes change states at rates to be introduced, and HIV transmission is governed by the network connections. A vacant state becomes susceptible state at the recruitment rate. Susceptible, opioid, infected, co-affected and AIDS states can change their state into a vacant state at natural death rate $ \mu $ or at disease-induced death rates $ d_u $, $ d_v $, $ d_i(\tau) $, $ d_a $, respectively. A susceptible state can be infected with HIV and change into an infected state, or can become opioid-dependent and change into opioid state. HIV and opioid states can get co-affected by adding the other disease. An opioid state or co-affected state can move to a susceptible or HIV-infected states respectively due to treatment denoted as $ \delta $. HIV-infected or co-affected states can move to the AIDS state at rates $ \gamma_v $ and $ \gamma_i(\tau) $.

    For an epidemic network, degree of a node is the number of contacts the node has with other nodes. We assume that the network contacts are HIV type contacts. That is an edge between any two nodes represents a potential for transmission of HIV, either through sexual contact or intravenous drug usage. For a network with maximal degree $ n $, the average network degree is given by

    $ < k > = \sum\limits_{k = 1}^n kp(k)\, , $

    where $ p(k) $ is the probability that a randomly chosen node has degree $ k $. That is, $ < k > $ is a constant, when $ n $ is pre-specified. This is a standard notation for average degree. For further background, [32] is a good source for network science. Empirical studies suggest that many real-life HIV networks have scale-free degree distribution $ p(k)\sim k^{-\eta} $, where $ 2 < \eta < 3 $ [33,34]. The conditional probability $ p(j | k) $ that a node with degree $ k $ is connected to a node with degree $ j $, is given by

    $ p( j | k) = \frac{j p(j)}{ < k > }\, . $

    Basically, $ p(j | k) $ gives probability that an individual with $ k $ contacts is connected to an individual with $ j $ contacts. We assume contacts that lead to opioid addiction are homogeneous (transmission of opioid addiction can be between two nodes with the same probability). With the structure of a complex network and infection age, let $ S_k(t) $, $ U_k(t) $, $ V_k(t) $, $ A_k(t) $, be the number of susceptible, opioid-dependent, HIV infected and AIDS nodes respectively with degree $ k $ at time $ t $ for $ k\in \{1, 2, \dots, n\} $. Let $ i_k(\tau, t) $ be the density of co-affected nodes of degree $ k $ at time $ t $ and with infection age $ \tau $. Then we can formulate the following network multi-scale model of HIV and opioid epidemics.

    $ {dSk(t)dt=Λkλu(t)Sk(t)kλv(t)Sk(t)μSk(t)+δUk(t),dUk(t)dt=λu(t)Sk(t)kqvλv(t)Uk(t)(μ+du+δ)Uk(t),dVk(t)dt=kλv(t)Sk(t)quλu(t)Vk(t)(μ+dv+γv)Vk(t)+δ0σ(τ)ik(t,τ)dτ,ik(t,τ)dt+ksik(t,τ)dτ=(μ+di(τ)+γi(τ)+δσ(τ))ik(t,τ),ksik(t,0)=kqvλv(t)Uk(t)+quλu(t)Vk(t)dAk(t)dt=γvVk(t)+0γi(τ)ik(t,τ)dτ(μ+da)Ak(t).
    $
    (2.2)

    The total population size of degree $ k $ is $ N_k(t) = S_k(t)+U_k(t)+V_k(t)+I_k(t) $ where $ I_k(t) = \int_{0}^{\infty}i_{k}(t, \tau)d\tau $. The force of HIV infection, $ \lambda_v(t) $ takes into account the heterogeneous mixing in the network:

    $ λv(t)=1<k>nk=1kp(k)λkv(t) where λkv(t)=βkv1Vk(t)+0βkv2(τ)ik(t,τ)dτNk(t).
    $
    (2.3)

    Thus $ \lambda_v^k(t) $ denotes the force of infection from a node with degree $ k $, where $ \beta^k_{v_1} $ is the infection rate from $ V_k(t) $ (HIV-infected node with degree $ k $) per effective contact, and similarly, $ \beta^k_{v_2}(\tau) $ is the time-since-infection dependent infection rate from $ i_k(\tau, t) $ (co-affected node with degree $ k $) per effective contact. Then the force of infection $ \lambda_v(t) $ for the heterogeneous network model is obtained by summing over all the degrees of the network $ \lambda_v^k(t) $ times the probability that the node with degree $ k $ is linked to the node with degree $ j $. This is the average force of infection from each contact, and is therefore multiplied by $ k $ in the equations (2.2) for HIV transmission to nodes of degree $ k $. The force of opioid addiction, $ \lambda_u(t) $ is then given by,

    $ λu(t)=nk=1λku(t) where λku(t)=βuUk(t)+0ik(t,τ)dτNk(t).
    $
    (2.4)

    Since the opioid contacts are homogeneous within the network, the force of addiction is obtained by summing over all the force of addictions of degree $ k $, $ \lambda_u^k(t) $. The within-host model remains the same as in (2.1) and the linking functions are given below.

    To link the within-host and between host models, we use data [35] to determine the form of the linking of the transmission coefficient $ \beta^{k}_{v_2}(\tau) $ [36] to the viral load. Fitting to the data [36], we obtain the following function for $ \beta^{k}_{v_2}(\tau) $:

    $ \beta^{k}_{v_2}(\tau) = \frac{\beta^k_0 V^r_i(\tau)}{B+V^r_i(\tau)}, $

    where $ r\approx 1 $. Further, we use the suggested functions in [37] to link the remaining $ \tau $-dependent rates:

    $ d_i(\tau) = d_0\left(T(0)-T\right) +d_1, \qquad \gamma_i(\tau) = \gamma_0\left(T(0) -T\right), \qquad \sigma(\tau) = \sigma_0, $

    where $ \beta_0^k, B, d_0, d_1, \gamma_0, \sigma_0, $ are constants. Disease-induced death rate $ d_i $ and transition to AIDS rate $ \gamma_i $ do not depend explicitly on the viral load because the viral load is high during the acute HIV phase but these rates are low during this same stage.

    Equilibria are time-independent solutions of the system and often determine its long-term behavior. We find the equilibria of this system by setting the derivatives with respect to $ t $ equal to zero.

    $ {Λkλu(t)Sk(t)kλv(t)Sk(t)μSk(t)+δUk(t)=0,λu(t)Sk(t)kqvλv(t)Uk(t)(μ+du+δ)Uk(t)=0,kλv(t)Sk(t)quλu(t)Vk(t)(μ+dv+γv)Vk(t)+δ0σ(τ)ik(t,τ)dτ=0,ksik(t,τ)dτ=(μ+di(τ)+γi(τ)+δσ(τ))ik(t,τ),ksik(0)=kqvλv(t)Uk(t)+quλu(t)Vk(t).
    $
    (3.1)

    At the disease-free equilibrium $ U_{k}(t) $, $ V_{k}(t) $ and $ i_{k}(t, \tau) $ are zero for all $ k $. So, the disease-free equilibrium is given by

    $ \varepsilon^{0} = \left(\frac{\Lambda_{1}}{\mu}, 0, 0, 0, \cdots, \frac{\Lambda_{n}}{\mu}, 0, 0, 0\right). $

    To determine the stability of the disease free equilibrium and to find the basic reproduction numbers of HIV infection and opioid addiction, denoted by $ \mathcal{R}_{u} $ and $ \mathcal{R}_{v} $ respectively, we linearize the system around the disease-free equilibrium. We take $ S_{k}(t) = S_{k}^{0} + x_{k}(t) $, $ U_{k}(t) = u_{k}(t) $, $ V_{k}(t) = v_{k}(t) $, $ N_{k}(t) = N_{k}^{0}+n_{k}(t) $ and $ i_{k}(t, \tau) = y_{k}(t, \tau) $. Then linearizing the system (2.2) takes the following form

    $ {dxk(t)dt=λu(t)S0kkλv(t)S0kμxk(t)+δuk(t),duk(t)dt=λu(t)S0k(μ+du+δ)uk(t),dvk(t)dt=kλv(t)S0k(μ+dv+γv)vk(t)+δ0σ(τ)yk(t,τ)dτ,yk(t,τ)dt+ksyk(t,τ)dτ=(μ+di(τ)+γi(τ)+δσ(τ))yk(t,τ),ksyk(t,0)=0
    $
    (3.2)

    where,

    $ \lambda_u(t) = \sum\limits_{k = 1}^n \beta_u \frac{u_k(t)+ \int_0^\infty y_k(t, \tau)d\tau}{N_k^{0}}, $
    $ \lambda_v(t) = \frac 1{ < k > }\sum\limits_{k = 1}^n kp(k) \frac{\beta^k_{v_1}v_k(t)+ \int_0^\infty \beta^k_{v_2}(\tau)y_k(t, \tau)d\tau}{N_k^{0}}, \quad \text{and} \quad N_k^0 = S_k^0 = \frac{\Lambda_k}{\mu}. $

    We look for solutions of the form $ x_{k}(t) = x_{k0}e^{\lambda t} $, $ u_{k}(t) = u_{k0}e^{\lambda t} $, $ v_{k}(t) = v_{k0}e^{\lambda t} $, $ y_{k}(t, \tau) = y_{k}(\tau)e^{\lambda t} $ and obtain the following eigenvalue problem,

    $ {(λ+μ)xk0+λu(t)S0k+kλv(t)S0kδuk0=0,(λ+μ+du+δ)uk0λu(t)S0k=0,(λ+μ+dv+γv)vk0kλv(t)S0kδ0σ(τ)yk(τ)dτ=0,ksyk(t,τ)dτ+λyk(τ)=(μ+di(τ)+γi(τ)+δσ(τ))yk(τ),ksyk(0)=0
    $
    (3.3)

    Solving the fourth equation of (3.3) we get

    $ yk(τ)=yk(0)π(τ)eλτks=0
    $
    (3.4)

    where

    $ \pi(\tau) = e^{- \frac 1{k_s}\int_0^\tau \mu+d_i(\xi)+\gamma_i(\xi)+\delta\sigma(\xi)d\xi}. $

    Substituting (3.4) in the second equation of (3.3) we have,

    $ uk0=S0k(λ+μ+du+δ)βunj=1uj0N0j.
    $
    (3.5)

    We notice that all $ u_{k0} $ have the same sign. Summing both sides from $ 1 $ to $ n $ after dividing with $ N_{k}^{0} $, we get

    $ nk=1uk0N0k=βu(λ+μ+du+δ)nk=1S0kN0knj=1uj0N0j
    $
    (3.6)

    If $ \sum_{k = 1}^{n}\frac{u_{k0}}{N_{k}^{0}} = 0 $ then $ u_{k0} = 0 $ for all $ k $ which is not the case. So cancelling $ \sum_{k = 1}^{n}\frac{u_{k0}}{N_{k}^{0}} $ from both sides of the equation we obtain,

    $ 1=βu(λ+μ+du+δ)nk=1S0kN0k=nβu(λ+μ+du+δ)
    $
    (3.7)

    since $ N_k^0 = S_k^0 = \frac{\Lambda_k}{\mu} $. We define

    $ Ru=nβuμ+du+δ.
    $
    (3.8)

    Substituting (3.4) in the third equation of (3.3) we have,

    $ vk0=kS0k(λ+μ+dv+γv)1<k>nj=1jp(j)βjv1vj0N0j.
    $
    (3.9)

    We note that all $ v_{k0} $ have the same sign. Multiplying both sides by $ \frac 1{ < k > }\sum_{k = 1}^n kp(k) \frac{\beta^k_{v_1}}{N_k^{0}} $ and summing from $ 1 $ to $ n $ we obtain,

    $ 1<k>nk=1kp(k)βkv1vk0N0k=1<k>nk=1k2p(k)βkv1S0kN0k(λ+μ+dv+γv)1<k>nj=1jp(j)βjv1vj0N0j.
    $
    (3.10)

    If $ \frac{1}{ < k > }\sum_{k = 1}^n kp(k) \frac{\beta^k_{v_1}v_{k0}}{N_k^{0}} = 0 $, then $ v_{k0} = 0 $ for all $ k $ which is not the case. So cancelling $ \frac 1{ < k > }\sum_{k = 1}^n kp(k) \frac{\beta^k_{v_1}v_{k0}}{N_k^{0}} $ from both sides we get,

    $ 1=1<k>nk=1k2p(k)βkv1S0kN0k(λ+μ+dv+γv)=1<k>nk=1k2p(k)βkv1(λ+μ+dv+γv).
    $
    (3.11)

    So we define

    $ Rv=1<k>nk=1k2p(k)βkv1(μ+dv+γv).
    $
    (3.12)

    Now we can prove the following theorem,

    Theorem 1. If $ \max\{\mathcal{R}_{u}, \mathcal{R}_{v}\} < 1 $ then the disease free equilibrium is locally asymptotically stable. If $ \max\{\mathcal{R}_{u}, \mathcal{R}_{v}\} > 1 $, then the disease-free equilibrium is unstable.

    Proof. Suppose

    $ G1(λ)=nβuλ+μ+du+δ,G2(λ)=1<k>nk=1k2p(k)βkv1(λ+μ+dv+γv).
    $
    (3.13)

    Then we notice that $ \mathcal{G}_{1}(0) = \mathcal{R}_{u} $, $ \mathcal{G}_{2}(0) = \mathcal{R}_{v} $, $ \lim_{\lambda \to \infty}\mathcal{G}_{1}(\lambda) = 0 $, $ \lim_{\lambda \to \infty}\mathcal{G}_{2}(\lambda) = 0 $.

    We claim that if $ \max\{\mathcal{R}_{u}, \mathcal{R}_{v}\} < 1 $ then the disease free equilibrium is locally asymptotically stable, that is all the solutions of $ \mathcal{G}_{1}(\lambda) = 1 $ and $ \mathcal{G}_{2}(\lambda) = 1 $, have negative real parts. To show this, we proceed by way of contradiction. Suppose one of the equations $ \mathcal{G}_{1}(\lambda) = 1 $ and $ \mathcal{G}_{2}(\lambda) = 1 $ has a solution $ \lambda_{0} $ with $ \Re{(\lambda_{0})} \geq 0 $. Then,

    $ 1 = |\mathcal{G}_{1}(\lambda_{0})|\leq |\mathcal{G}_{1}(\Re{\lambda_{0}})|\leq |\mathcal{G}_{1}(0)| = \mathcal{R}_{u}, $

    or

    $ 1 = |\mathcal{G}_{2}(\lambda_{0})|\leq |\mathcal{G}_{2}(\Re{\lambda_{0}})|\leq |\mathcal{G}_{2}(0)| = \mathcal{R}_{v}. $

    This is a contradiction. Hence $ \varepsilon^{0} $ is locally asymptotically stable when $ \max\{\mathcal{R}_{u}, \mathcal{R}_{v}\} < 1 $. If $ \max\{\mathcal{R}_{u}, \mathcal{R}_{v}\} > 1 $, let us assume, without loss of generality, $ \mathcal{R}_{u} > 1 $. Then $ \mathcal{G}_{1}(\lambda) = 1 $ has a positive solution, $ \lambda^{*} $. Thus, the disease free equilbrium is unstable.

    In this section, we prove the existence and stability of the two boundary equilibria $ E_1^* $ and $ E_2^* $ corresponding to opioid addiction and HIV transmission in a single population respectively. To obtain $ E_1^* $ we let $ S_k(t) = S_{k1}^* $, $ U_{k}(t) = U_{k1}^* $ and $ V_{k}(t) = 0 $ and $ i_{k}(t, \tau) = 0 $ for all $ k $, i.e., $ E_1^* = (S_{11}^*, U_{11}^{*}, 0, 0, \cdots, S_{n1}^{*}, U_{n1}^{*}, 0, 0) $. We get the following equations

    $  ΛkμSk1Sk1λu(U1)+δUk1=0 Sk1λu(U1)(μ+du+δ)Uk1=0 λu(U1)=βunk=1Uk1Nk1
    $
    (4.1)

    From the second equation of (4.1) we get,

    $ S_{k1}^* = \frac{\mu+d_u+\delta}{\beta_{u}}\frac{U_{k1}^{*}}{\sum\limits_{j = 1}^{n}\frac{U_{j1}^{*}}{N_{j1}^{*}}}. $

    Summing both sides of the equation from $ k = 1 $ to $ k = n $, and dividing by $ N_{k1}^{*} $ we get,

    $ nk=1Sk1Nk1=μ+du+δβu=nRu
    $
    (4.2)

    We know $ \frac{S_{k1}^{*}+U_{k1}^{*}}{N_{k1}^*} = 1 $, and so

    $ nk=1Sk1+Uk1Nk1=n
    $
    (4.3)

    i.e.,

    $ \sum\limits_{k = 1}^{n}\frac{U_{k1}^{*}}{N_{k1}^{*}} = n\left(1-\frac{1}{\mathcal{R}_{u}}\right) $

    Plugging in the values in the first equation of (4.1) and solving for $ U_{k1}^{*} $ we get

    $ Uk1=ΛkμRu1+μ+du=Λk(μ+du)(1μ(1Ru)(μ+du))
    $
    (4.4)

    Now when $ \mathcal{R}_{u} < 1 $, $ \sum_{k = 1}^{n}\frac{U_{k1}^{*}}{N_{k1}^{*}} < 0 $, which implies $ U_{k1}^{*} < 0 $ for at least one $ k $. When $ \mathcal{R}_{u} > 1 $, $ U_{k1}^{*} > 0 $. Thus, $ E_{1}^{*} $ exists if and only if $ \mathcal{R}_{u} > 1 $.

    To obtain $ E_2^* $ we let $ S_k(t) = S_{k2}^* $, $ V_{k}(t) = V_{k2}^* $ and $ U_{k}(t) = 0 $ and $ i_{k}(t, \tau) = 0 $ for all $ k $, i.e., $ E_2^* = (S_{12}^*, 0, V_{12}^{*}, 0, \cdots, S_{n2}^{*}, 0, V_{n2}^{*}, 0) $. We get the following equations

    $  ΛkμSk2kSk2λv(V2)=0, kSk2λv(V2)(μ+dv+γv)Vk2=0, λv(V2)=1<k>nj=1jp(j)βjv1Vj2Nj2.
    $
    (4.5)

    From the second equation of (4.5) we get

    $ kS_{k2}^{*} = \frac{\mu+d_{v}+\gamma_{v}}{\lambda_{v}(V_{2}^{*})}V_{k2}^{*}. $

    Dividing both sides by $ N_{k2}^{*} $ and multiplying with $ \frac 1{ < k > }\sum_{k = 1}^n kp(k)\beta^k_{v_1} $ and adding from $ 1 $ to $ n $ we get,

    $ 1<k>nk=1k2p(k)βkv1Sk2Nk2=(μ+dv+γv).
    $
    (4.6)

    We know $ \frac{S_{k2}^{*}+V_{k2}^{*}}{N_{k2}^*} = 1 $, and so

    $ 1<k>nk=1k2p(k)βkv1(Vk2+Sk2)Nk2=(μ+dv+γv)Rv,
    $
    (4.7)

    which gives us,

    $ \frac 1{ < k > }\sum\limits_{k = 1}^n k^2p(k)\frac{\beta^k_{v_1}V_{k2}^{*}}{N_{k2}^{*}} = (\mu+d_{v}+\gamma_{v})(\mathcal{R}_{v}-1). $

    So when $ \mathcal{R}_{v} < 1 $, $ E_{2}^{*} $ does not exist. Now

    $ N_{k2}^* = V_{k2}^*+S_{k2}^{*} = V_{k2}^*+V_{k2}^*(\frac{\mu+d_{v}+\gamma_{v}}{k\lambda_{v}(V_{2}^{*})}), $

    and so we obtain,

    $ Vk2Nk2=11+μ+dv+γvkλv(V2),
    $
    (4.8)
    $ λv(V2)=1<k>nk=1kp(k)βkv111+μ+dv+γvkλv(V2)=1<k>nk=1kp(k)βkv1kλv(V2)kλv(V2)+(μ+dv+γv).
    $
    (4.9)

    $ \lambda_{v}(V_{2}^{*}) = 0 $ is a solution to (4.5) which gives the disease free equilibrium. Then for a boundary equilibrium, $ \lambda_{v}(V_{2}^{*}) > 0 $ is a root of $ f(\lambda_{v}) $, where

    $ f(\lambda_{v}) = \frac 1{ < k > }\sum\limits_{k = 1}^n k^2p(k)\frac{\beta^k_{v_1}}{k\lambda_{v}(V_{2}^{*})+(\mu+d_{v}+\gamma_{v})}-1. $

    As $ \lambda_{v} $ increases $ f $ decreases. $ \lim_{\lambda_{v} \to \infty} f(\lambda_{v}) = -1 $. But $ f(0) = \mathcal{R}_{v} -1 > 0 $. Then if $ \mathcal{R}_{v} > 1 $, $ f(\lambda_{v}) $ has a unique zero, giving us a unique boundary equilibrium for the system.

    To find the invasion number of HIV and stability of $ E_{1}^{*} $ we first linearize the system (2.2) around $ E_{1}^{*} $. We set $ S_{k}(t) = x_{k}(t)+S_{k1}^{*} $, $ U_{k}(t) = u_{k}(t)+U_{k1}^{*} $, $ V_{k}(t) = v_{k}(t) $, $ i_{k}(t, \tau) = y_{k}(t, \tau) $ and $ N_{k}(t) = n_{k}(t)+N_{k1}^{*} $, the system for the perturbations becomes,

    $ \left\{ dxk(t)dt=Sk1λu(u,y)xk(t)βunj=1Uj1Nj1+Sk1βunj=1Uj1nkN2j1μxk(t)+δuk(t)kSk1λv(v,y),duk(t)dt=Sk1λu(u,y)+xk(t)βunj=1Uj1Nj1Sk1βunj=1Uj1nkN2j1kqvUk1λv(v,y)(μ+du+δ)uk(t),dvk(t)dt=kSk1λv(v,y)quvk(t)βunj=1Uj1Nj1(μ+dv+γv)vk(t)+δ0σ(τ)yk(t,τ)dτ,yk(t,τ)dt+ksyk(t,τ)dτ=(μ+di(τ)+γi(τ)+δσ(τ))yk(t,τ),ksyk(t,0)=kqvUk1λv(v,y)+quvk(t)βunj=1Uj1Nj1,
    \right. $
    (4.10)

    where,

    $ \lambda_{u}(u, y) = \beta_{u}\sum\limits_{j = 1}^{n}\frac{u_{j}+\int_{0}^{\infty}y_{j}(t, \tau)}{N_{j1}^{*}}, $
    $ \lambda_{v}(v, y) = \frac 1{ < k > }\sum\limits_{j = 1}^n jp(j) \frac{\beta^j_{v_1}v_j(t)+ \int_0^\infty \beta^j_{v_2}(\tau)y_j(t, \tau)d\tau}{N_{j1}^{*}}. $

    We look for solutions of the form $ x_{k}(t) = x_{k}e^{\lambda t} $, $ u_{k}(t) = u_{k}e^{\lambda t} $, $ v_{k}(t) = v_{k}e^{\lambda t} $, $ y_{k}(\tau, t) = y_{k}(\tau)e^{\lambda t} $ and obtain the following eigenvalue problem,

    $ {λxk=Sk1λu(u,y)xkβunj=1Uj1Nj1+Sk1βunj=1Uj1nkN2j1μxk+δukkSk1λv(v,y),λuk=Sk1λu(u,y)+xkβunj=1Uj1Nj1Sk1βunj=1Uj1nkN2j1kqvUk1λv(v,y)(μ+du+δ)uk,λvk=kSk1λv(v,y)quvkβunj=1Uj1Nj1(μ+dv+γv)vk+δ0σ(τ)yk(τ)dτ,ksyk(τ)dτ+λyk=(μ+di(τ)+γi(τ)+δσ(τ))yk(τ),ksyk(0)=kqvUk1λv(v,y)+quvkβunk=1Uk1Nk1,
    $
    (4.11)

    where,

    $ \lambda_{u}(u, y) = \beta_{u}\sum\limits_{j = 1}^{n}\frac{u_{j}+\int_{0}^{\infty}y_{j}(\tau)d\tau}{N_{j1}^{*}}, $

    and

    $ λv(v,y)=1<k>nj=1jp(j)βjv1vj+0βjv2(τ)yj(τ)dτNj1.
    $
    (4.12)

    Now, using the third, fourth and fifth equation of (4.11) we will compute the invasion number of HIV.

    $ {(λ+μ+dv+γv)vk=kSk1λv(v,y)quvkβunj=1Uj1Nj1+δ0σ(τ)yk(τ)dτ,ksyk(τ)dτ+λyk=(μ+di(τ)+γi(τ)+δσ(τ))yk(τ),ksyk(0)=kqvUk1λv(v,y)+quvkβunj=1Uj1Nj1.
    $

    From the second equation of (4.2) we get $ y_{k}(\tau) = y_k(0)\pi(\tau)e^{-\frac{\lambda\tau}{k_s}} $, where $ \pi(\tau) $ is as defined before. Suppose $ K = \beta_{u}q_{u}\sum_{k = 1}^{n}\frac{U_{k1}^{*}}{N_{k1}^{*}} $ and $ Q(\lambda) = \delta[\int_0^{\infty} \sigma(\tau)\pi(\tau) e^{-\frac{\lambda\tau}{k_s}} d\tau] $. Then from the first and third equations of (4.2) we get,

    $ {(λ+μ+dv+γv+K)vkQ(λ)yk(0)=kSk1λv(v,y),Kvk+ksyk(0)=kqvUk1λv(v,y).
    $
    (4.13)

    Solving for $ v_k $ and $ y_k(0) $ we obtain,

    $  vk=kskSk1+kqvUk1Q(λ)ks(λ+μ+dv+γv+K)KQ(λ)λv(v,y), yk(0)=kqv(λ+μ+dv+γv+K)Uk1+KkSk1ks(λ+μ+dv+γv+K)KQ(λ)λv(v,y).
    $

    Supplying these values in Eq (4.12), and cancelling $ \lambda_v(v, y) $ from both sides of the equation, we obtain,

    $ 1=1<k>nj=1j2p(j)[βjv1Nj1ksjSj1+jqvUj1Q(λ)ks(λ+μ+dv+γv+K)KQ(λ)+0βjv2(τ)π(τ)dτNj1jqv(λ+μ+dv+γv+K)Uj1+KjSj1ks(λ+μ+dv+γv+K)KQ(λ)].
    $
    (4.14)

    We define

    $ R1vi=1<k>nj=1j2p(j)[βjv1Nj1ksjSj1+jqvUj1δ[0σ(τ)π(τ)dτ]ks(μ+dv+γv+K)Kδ[0σ(τ)π(τ)dτ]+0βjv2(τ)π(τ)dτNj1jqv(μ+dv+γv+K)Uj1+KjSj1ks(μ+dv+γv+K)Kδ[0σ(τ)π(τ)dτ]].
    $
    (4.15)

    We call $ \mathcal{R}^{1}_{v_i} $ the invasion reproduction number of HIV infection. Now suppose,

    $ Gvi(λ)=1<k>nj=1j2p(j)[βjv1Nj1ksSj1+qvUj1Q(λ)ks(λ+μ+dv+γv+K)KQ(λ)+0βjv2(τ)eλτπ(τ)dτNj1qv(λ+μ+dv+γv+K)Uj1+KSj1ks(λ+μ+dv+γv+K)KQ(λ)]
    $

    $ \int_{0}^{\infty} \beta^{j}_{v_2}(\tau) e^{-\lambda \tau} \pi(\tau) d\tau = \beta_{j}(\lambda) $. $ \beta_{j}(\lambda) $ is bounded above by $ \beta_{j}(0) $ and $ Q(\lambda) $ is bounded above by $ Q(0) $. Then $ \mathcal{G}_{vi}(0) = \mathcal{R}^{1}_{v_i} $ and $ \lim_{\lambda \to \infty} \mathcal{G}_{vi}(\lambda) = 0 $. Suppose (4.14) has a solution $ \lambda = x+iy $ with $ \Re({\lambda}) = x\ge0 $ and $ \mathcal{R}_{vi}^{1} < 1 $. First we prove the following result.

    $ |qv(λ+μ+dv+γv+K)Uj1Nj1|+|KSj1Nj1||ks(λ+μ+dv+γv+K)||KQ(0)||qv(μ+dv+γv+K)Uj1Nj1|+|KSj1Nj1||ks(μ+dv+γv+K)||KQ(0)|
    $
    (4.16)

    Proof. To prove (4.16) we write down the left hand side of the inequality,

    $ \frac{\left|q_v(\lambda+\mu+d_v+\gamma_v+K)\frac{U_{j1}^{*}}{N_{j1}^{*}}\right|+\left|K\frac{S_{j1}^{*}}{N_{j1}^{*}}\right|}{\left|k_s(\lambda+\mu+d_v+\gamma_v+K)\right|-\left|KQ(0)\right|} = \frac{q_v C_1 z +K C_2}{k_s z-KQ(0)} = f(z), $

    where, $ C_1 = \frac{U_{j1}^{*}}{N_{j1}^{*}} $, $ C_2 = \frac{S_{j1}^{*}}{N_{j1}^{*}} $ and $ z = \sqrt{(x+\mu+d_v+\gamma_v+K)^{2}+y^{2}} $. Since $ f'(z) < 0 $, $ f(z) $ is a decreasing function. That is when $ z(0, 0)\leq z(x, y) $, $ f(z(0, 0))\geq f(z(x, y)) $. But $ f(z(0, 0)) $ is just the right hand side of (4.16). Using (4.16) we can now state the following,

    $ 1=|Gvi(λ)|=|1<k>nj=1j2p(j)[βjv1Nj1ksSj1+qvUj1Q(λ)ks(λ+μ+dv+γv+K)KQ(λ)+βj(λ)Nj1qv(λ+μ+dv+γv+K)Uj1+KSj1ks(λ+μ+dv+γv+K)KQ(λ)]|1<k>nj=1j2p(j)[|βjv1Nj1ksSj1+qvUj1Q(0)ks(μ+dv+γv+K)KQ(0)|+|βj(0)qv(λ+μ+dv+γv+K)Uj1Nj1+KSj1Nj1ks(λ+μ+dv+γv+K)KQ(0)|]1<k>nj=1j2p(j)[|βjv1Nj1ksSj1+qvUj1Q(0)ks(μ+dv+γv+K)KQ(0)|+βj(0)|qv(λ+μ+dv+γv+K)Uj1Nj1|+|KSj1Nj1||ks(λ+μ+dv+γv+K)||KQ(0)|]1<k>nj=1j2p(j)[|βjv1Nj1ksSj1+qvUj1Q(0)ks(μ+dv+γv+K)KQ(0)|+|qv(μ+dv+γv+K)Uj1Nj1|+|KSj1Nj1||ks(μ+dv+γv+K)||KQ(0)|]=|Gvi(0)|=R1vi<1.
    $
    (4.17)

    This is a contradiction. So (4.14) only has solutions with non-negative real parts when $ \mathcal{R}_{vi}^{1} < 1 $.

    Now let us suppose, $ \mathcal{R}^{1}_{vi} > 1 $. It can be shown that $ \mathcal{G}_{vi}(\lambda) $ is decreasing. Then since $ \mathcal{G}_{vi}(0) = \mathcal{R}^{1}_{v_i} > 1 $ and $ \lim_{\lambda \to \infty} \mathcal{G}_{vi}(\lambda) = 0 $ for real and positive $ \lambda $, (4.14) must have at least one positive root when $ \mathcal{R}_{vi} > 1 $. Now, from (4.1)–(4.3) we get,

    $ S_{k1}^* = \frac{n}{\mathcal{R}_{u}}\frac{U_{k1}^{*}}{n\left(1-\frac{1}{\mathcal{R}_{u}}\right)} = \frac{U_{k1}^{*}}{\mathcal{R}_{u}-1}, $

    i.e.,

    $ \frac{\frac{U_{k1}^{*}}{\mathcal{R}_{u}-1}+U_{k1}^{*}}{N_{k1}^*} = \frac{U_{k1}^{*}}{N_{k1}^*}\left(\frac{1}{\mathcal{R}_{u}-1}+1\right) = 1. $

    Solving the equation we get,

    $  Uk1Nk1=11Ru, Sk1Nk1=1Ru.
    $
    (4.18)

    To find the remaining eigenvalues, satisfying the third, fourth and fifth equation of (4.11), $ y_{j}(\tau) = 0 $ and $ v_j = 0 $ for all $ j = 1, 2, \cdots, n $. The first two equations then just reduce to

    $  λxk=Sk1βunj=1ujNjxkβunj=1Uj1Nj1+Sk1βunj=1Uj1nkN2j1μxk+δuk, λuk=Sk1βunj=1ujNj+xkβunj=1Uj1Nj1Sk1βunj=1Uj1nkN2j1(μ+du+δ)uk,
    $
    (4.19)

    Adding the two equations and solving for $ n_{k} $ we get

    $ n_{k} = -\frac{d_{u}u_{k}}{(\lambda+\mu), } $

    i.e.,

    $ x_{k} = -\frac{\lambda+\mu+d_{u}}{\lambda+\mu}u_{k}. $

    Replacing $ x_{k} $ and $ n_{k} $ in the second equation of (4.19) we get,

    $  (λ+μ+du+δ)uk=Sk1βunj=1ujNj1λ+μ+duλ+μukβunj=1Uj1Nj1+Sk1βunj=1Uj1N2j1duuj(λ+μ) (λ+μ+du+δ+βun(11Ru)λ+μ+duλ+μ)uk=Sk1βunj=1ujNj1(1+duUj1Nj1λ+μ).
    $
    (4.20)

    Multiplying both sides of the equation $ \frac{1}{N_{k1}^{*}} $ and summing over $ 1 $ to $ n $,

    $ (λ+μ+du+δ+βun(11Ru)λ+μ+duλ+μ)nk=1ukNk1=nk=1Sk1Nk1βunj=1ujNj1(1+duUj1Nj1λ+μ).
    $
    (4.21)

    $ \sum_{j = 1}^{n}\frac{u_{j}}{N_{j1}^{*}} = 0 $ implies from (4.20) all $ u_k $ would be zero, which would not be of interest. $ \sum_{j = 1}^{n}\frac{u_{j}}{N_{j1}^{*}}\neq 0 $ for non-equilibrium points, and cancelling the expression on both sides, then the characteristic equation becomes,

    $ (λ+μ+du+δ)(λ+μ)+βu(λ+μ+du)n(11Ru)=βu(λ+μ+du(11Ru))nRu.
    $
    (4.22)

    Rewriting this equation as a quadratic equation, we get

    $ \lambda^{2}+(2\mu+d_u+\delta+\beta_u n -\beta_{u}\frac{n}{\mathcal{R}_{u}})\lambda+\left(\mu+d_{u}+\delta\right)\mu+\beta_{u}n\left(1-\frac{1}{\mathcal{R}_{u}}\right)(\mu+d_{u})-\beta_{u}n\frac{1}{\mathcal{R}_{u}}\left(\mu+d_{u}\left(1-\frac{1}{\mathcal{R}_{u}}\right)\right) = 0 $ (4.23)

    Simplifying the equation, we have $ \lambda^2 +b \lambda +c = 0 $ where $ b = \mu+\beta_{u}n > 0 $ and $ c = \left(1-\frac{1}{\mathcal{R}_{u}}\right)\left(\beta_{u}n(\mu+d_u)-(\mu+d_u+\delta)d_{u}\right) > 0 $ since, $ \beta_{u}n > (\mu+d_{u}+\delta) $ when $ \mathcal{R}_{u} > 1 $. Hence this quadratic equation has only roots with negative real parts. Combining the work above we can conclude,

    Theorem 2. The unique boundary equilibrium $ E_{1}^{*} $ is locally asymptotically stable if $ \mathcal{R}^{1}_{v_i} < 1 $, and it is unstable if $ \mathcal{R}^{1}_{v_i} > 1 $.

    To find the invasion number of opioid addiction and stability of $ E_{2}^{*} $ we first linearize the system (2.2) around $ E_{2}^{*} $. We set $ S_{k}(t) = x_{k}(t)+S_{k2}^{*} $, $ U_{k}(t) = u_{k}(t) $, $ V_{k}(t) = v_{k}(t)+V_{k2}^{*} $, $ i_{k}(t, \tau) = y_{k}(t, \tau) $ and $ N_{k}(t) = n_{k}(t)+N_{k2}^{*} $. The system for the perturbations becomes,

    $ \left\{ dxk(t)dt=Sk2λ2u(u,y)μxk(t)+δuk(t)kSk2λ2v(v,y)C1kxk+kSk21<k>nj=1jp(j)βjv1Vj2  njN j22,duk(t)dt=Sk2λ2u(u,y)kqvukC1(μ+du+δ)uk(t),dvk(t)dt=kxkC1+kSk2λ2v(v,y)kSk21<k>nj=1jp(j)βjv1Vj2njN2j2quVk2λ2u(u,y)(μ+dv+γv)vk(t)+δ0σ(τ)yk(t,τ)dτ,yk(t,τ)dt+ksyk(t,τ)dτ=(μ+di(τ)+γi(τ)+δσ(τ))yk(t,τ),ksyk(t,0)=kqvC1uk+quVk2λ2u(u,y),
    \right. $
    (4.24)

    where,

    $ \lambda_{u}^{2}(u, y) = \beta_{u}\sum\limits_{j = 1}^{n}\frac{u_{j}+\int_{0}^{\infty}y_{j}(t, \tau)}{N_{j2}^{*}} $
    $ \lambda_{v}^{2}(v, y) = \frac 1{ < k > }\sum\limits_{j = 1}^n jp(j) \frac{\beta^j_{v_1}v_j(t)+ \int_0^\infty \beta^j_{v_2}(\tau)y_j(t, \tau)d\tau}{N_{j2}^{*}} $

    and

    $ C_{1} = \frac 1{ < k > }\sum\limits_{j = 1}^n jp(j) \frac{\beta^j_{v_1}V_{j2}^{*}}{N_{j2}^{*}}. $

    We look for solutions of the form $ x_{k}(t) = x_{k}e^{\lambda t} $, $ u_{k}(t) = u_{k}e^{\lambda t} $, $ v_{k}(t) = v_{k}e^{\lambda t} $, $ y_{k}(t, \tau) = y_{k}(\tau)e^{\lambda t} $ and obtain the following eigenvalue problem,

    $ {λxk=Sk2λ2u(u,y)μxk+δukkSk2λ2v(v,y)C1kxk+kSk21<k>nj=1jp(j)βkv1Vj2njN2j2,λuk=Sk2λ2u(u,y)kqvukC1(μ+du+δ)uk,λvk=kxkC1+kSk2λ2v(v,y)kSk21<k>nj=1jp(j)βkv1Vj2njN2j2quVk2λ2u(u,y)(μ+dv+γv)vk+δ0σ(τ)yk(τ)dτ,ksyk(τ)dτ+λyk=(μ+di(τ)+γi(τ)+δσ(τ))yk(τ),ksyk(0)=kqvC1uk+quVk2λ2u(u,y).
    $
    (4.25)

    From the fourth equation of (4.25) we get

    $ yk(τ)=yk(0)π(τ)eλksτ.
    $
    (4.26)

    Let $ Q(\lambda) = \int_{0}^{\infty}\pi(\tau)e^{\frac{-\lambda}{k_s}\tau}d\tau. $ From the second equation of (4.25) we get

    $ u_{k} = \frac{S_{k2}^{*}}{\lambda+\mu+d_{u}+\delta+kq_v C_1}\beta_{u}\sum\limits_{j = 1}^{n}\frac{u_{j}+y_{j}(0)Q(\lambda)}{N_{j2}^{*}}. $

    Multiplying both sides of this equation with $ \frac{1}{N_{k2}^{*}} $ we get,

    $ ukNk2=βuSk2(λ+μ+du+δ+kqvC1)Nk2nj=1uj+yj(0)Q(λ)Nj2.
    $
    (4.27)

    From the fifth equation of (4.25) we get,

    $ k_sy_{k}(0) = kq_vC_{1}u_{k}+q_uV_{k2}^{*}\beta_{u}\sum\limits_{j = 1}^{n}\frac{u_{j}+y_{j}(0)Q(\lambda)}{N_{j2}^{*}}. $

    Multiplying both sides of this equation with $ \frac{Q(\lambda)}{N_{k2}^{*}} $ we get,

    $ yk(0)Q(λ)Nk2=kqvC1Q(λ)βuSk2ks(λ+μ+du+δ+kqvC1)Nk2nj=1uj+yj(0)Q(λ)Nj2+Q(λ)quVk2βuNk2ksnj=1uj+yj(0)Q(λ)Nj2.
    $
    (4.28)

    Summing both side of (4.27) and (4.28) from $ 1 $ to $ n $ and adding together we get,

    $ S(λ)=S(λ)[nk=1(1+kqvC1Q(λ))βuSk2ks(λ+μ+du+δ+kqvC1)Nk2+Q(λ)quβuksnk=1Vk2Nk2]
    $
    (4.29)

    where, $ S(\lambda) = \sum_{j = 1}^{n}\frac{u_{j}+y_{j}(0)Q(\lambda)}{N_{j2}^{*}}. $ Since $ S(\lambda) = 0 $ implies from (4.27) $ u_k = 0 $ for all $ k $, $ S(\lambda)\neq0 $. We cancel $ S(\lambda) $ from both sides and get,

    $ 1=nk=1(1+kqvC1Q(λ))βuSk2ks(λ+μ+du+δ+kqvC1)Nk2+Q(λ)quβuksnk=1Vk2Nk2.
    $
    (4.30)

    Let $ \Pi = \int_{0}^{\infty}\pi(\tau)d\tau $. We define

    $ R2ui=nk=1(1+kqvC1Π)βuSk2ks(μ+du+δ+kqvC1)Nk2+Πquβuksnk=1Vk2Nk2.
    $
    (4.31)

    We call $ \mathcal{R}^{2}_{u_i} $ the invasion reproduction number of opioid addiction. We claim that when $ \mathcal{R}^{2}_{u_i} < 1 $ the boundary equilibrium $ E_{2}^{*} $ is locally asymptotically stable, that is all the roots of (4.30) have negative real parts. Suppose

    $ Gui(λ)=nk=1(1+kqvC1Q(λ))βuSk2ks(λ+μ+du+δ+kqvC1)Nk2+Q(λ)quβuksnk=1Vk2Nk2.
    $

    Then $ \mathcal{G}_{ui}(0) = \mathcal{R}^{2}_{u_i} $ and $ \lim_{\lambda \to \infty}\mathcal{G}_{ui}(\lambda) = 0 $. Assume the Eq (4.30) has roots with non-negative real part $ \Re (\lambda) > 0 $. The Eq (4.30) satisfies,

    $  1=|nk=1(1+kqvC1Q(λ))βuSk2ks(λ+μ+du+δ+kqvC1)Nk2+Q(λ)quβuksnk=1Vk2Nk2| |nk=1(1+kqvC1Q(λ))βuSk2ks(λ+μ+du+δ+kqvC1)Nk2|+|Q(λ)quβuksnk=1Vk2Nk2| |nk=1(1+kqvC1Q((λ)))βuSk2ks(λ+μ+du+δ+kqvC1)Nk2|+|Q((λ))quβuksnk=1Vk2Nk2| |nk=1(1+kqvC1Π)βuSk2ks(μ+du+δ+kqvC1)Nk2|+|Πquβuksnk=1Vk2Nk2| R2ui<1
    $
    (4.32)

    This is a contradiction. Hence all roots of (4.30) have negative real parts when $ \mathcal{R}^{2}_{u_i} < 1 $. Now let us suppose, $ \mathcal{R}^{2}_{u_i} > 1 $. Then since $ \mathcal{G}'_{ui}(\lambda) < 0 $ when $ \lambda > 0 $, $ \mathcal{G}_{ui}(\lambda) $ is decreasing when $ \lambda > 0 $. But we have, $ \mathcal{G}_{ui}(0) = \mathcal{R}^{2}_{u_i} > 1 $ and $ \lim_{\lambda \to \infty}\mathcal{G}_{ui}(\lambda) = 0 $. Then (4.30) has at least one positive root when $ \mathcal{R}^{2}_{u_i} > 1 $. If $ \lambda $ is not a solution of characteristic equation (4.30), we have $ u_j = 0 $, $ y_j(0) = 0 $, the remaining two equations of (4.25) then just reduce to

    $  λxk=μxkkSk21<k>nj=1jp(j)βjv1vjNj2C1kxk+kSk21<k>nj=1jp(j)βjv1Vj2njN2j2,λvk=kxkC1+kSk21<k>nj=1jp(j)βjv1vjNj2kSk21<k>nj=1jp(j)βjv1Vj2njN2j2(μ+dv+γv)vk.
    $
    (4.33)

    Adding the two equations and solving for $ n_{k} $ we get

    $ n_{k} = -\frac{(d_{v}+\gamma_v)v_{k}}{(\lambda+\mu)}, $

    i.e.,

    $ x_{k} = -\frac{\lambda+\mu+d_{v}+\gamma_v}{\lambda+\mu}v_{k}. $

    Using these values for $ n_{k} $ and $ x_{k} $ in the second equation of (4.33) we get

    $ (λ+μ+dv+γv+kC1(λ+μ+dv+γv)λ+μ)vk=kSk21<k>nj=1jp(j)βjv1Vj2Nj2(λ+μ+(dv+γv)Vj2Nj2λ+μ).
    $
    (4.34)

    Dividing both sides by $ N_{k2}^{*} $ and readjusting we obtain,

    $ vkNk2=kSk2Nk21(λ+μ+dv+γv)(λ+μ+kC1)1<k>nj=1jp(j)βjv1vjNj2(λ+μ+(dv+γv)Vj2Nj2).
    $
    (4.35)

    Multiplying both sides of this equation with $ \frac 1{ < k > }\sum_{k = 1}^n kp(k) \beta^k_{v_1}(\lambda+\mu+(d_v+\gamma_v)\frac{V_{k2}^{*}}{N_{k2}^{*}}) $, we get,

    $ T(λ)=1<k>nk=1k2p(k)βkv1Sk2Nk2(λ+μ+(dv+γv)Vk2Nk2)(λ+μ+dv+γv)(λ+μ+kC1)T(λ),
    $
    (4.36)

    where $ T(\lambda) = \frac 1{ < k > }\sum_{j = 1}^n jp(j) \frac{\beta^j_{v_1}v_{j}}{N_{j2}^{*}} \left(\lambda+\mu+(d_v+\gamma_v)\frac{V_{j2}^{*}}{N_{j2}^{*}}\right) $. Since $ T(\lambda) = 0 $ implies from (4.35), $ v_k = 0 $ for all $ k $, $ T(\lambda) \neq 0 $ and we get the following characteristic equation,

    $ 1=1<k>nk=1k2p(k)βkv1Sk2Nk2(λ+μ+(dv+γv)Vk2Nk2)(λ+μ+dv+γv)(λ+μ+kC1).
    $
    (4.37)

    From (4.8) we obtain,

    $ \frac{V_{k2}^{*}}{N_{k2}^{*}} = \frac{kC_1}{kC_1+\mu+d_v+\gamma_v}, $
    $ \frac{S_{k2}^{*}}{N_{k2}^{*}} = \frac{\mu+d_v+\gamma_v}{kC_1+\mu+d_v+\gamma_v}. $

    Assume the Eq (4.37) has roots with non-negative real part. Using (4.9), $ \Re \lambda \ge 0 $ the Eq (4.37) satisfies,

    $  1=|1<k>nk=1k2p(k)βkv1Sk2Nk2(λ+μ+(dv+γv)Vk2Nk2)(λ+μ+dv+γv)(λ+μ+kC1)| =|1<k>nk=1k2p(k)βkv1μ+dv+γvkC1+μ+dv+γv(λ+μ+(dv+γv)kC1kC1+μ+dv+γv)(λ+μ+dv+γv)(λ+μ+kC1)| <|1<k>nk=1k2p(k)βkv1μ+dv+γvkC1+μ+dv+γv1(λ+μ+dv+γv)| 1<k>nk=1k2p(k)βkv1μ+dv+γvkC1+μ+dv+γv1(μ+dv+γv)=1
    $
    (4.38)

    This is a contradiction. So we can state the following theorem,

    Theorem 3. The unique boundary equilibrium $ E_{2}^{*} $ is locally asymptotically stable if $ \mathcal{R}^{2}_{u_i} < 1 $, and is unstable if $ \mathcal{R}^{2}_{u_i} > 1 $.

    We present a numerical scheme for the immuno-epidemiological models (2.1) and (2.2). The within-host model, consisting of ordinary differential equations can be solved by a stiff ODE solver in MATLAB.

    For the between-host model we introduce a finite-difference method. We discretize the domain

    $ \mathcal{D} = \{ (t, \tau): \ 0 \leq t \leq T, \ 0\leq\tau\leq A\} $

    where A is a maximal infection age and time $ T < \infty $, a maximal time. We take $ \Delta t = \Delta \tau $, with $ k_s = 1 $, and so the points in age and line direction can be computed as,

    $ \tau_m = m\Delta t \quad t_j = j\Delta t. $

    Setting $ M = \left[\frac{A}{\Delta t}\right] $ and $ N = \left[\frac{T}{\Delta t}\right] $, we obtain $ A = M\Delta t $, $ T = N\Delta t $. The numerical method computes approximations to the solution at the mesh points. We assume $ S_k(t_j) = S_k^j $, $ U_k(t_j) = U_k^j $, $ V_k(t_j) = V_k^j $ and $ i_k(t_j, \tau_m) = i^j_{m, k} $. We summarize the numerical method below,

    $ { λjv=1knk=1βv1Vjk+Mm=1Δtβv2mij+1m,kNjk,j=0,,N1; λju=nk=1βuUjk+Mm=1Δtij+1m,kNjk,j=0,,N1; Sj+1k=Sjk+ΛΔt+δUjkΔt1+kλjvΔt+λjuΔt+μΔt,j=0,,N1; Uj+1k=Ujk+λjuSj+1kΔt1+kqvλjvΔt+(μ+du+δ)Δt,j=0,,N1; Vj+1k=Vjk+kλjvSj+1kΔt+δΔtMm=1σmij+1m,kΔt1+quλjuΔt+(μ+dv+γv)Δt,j=0,,N1; ij+1m+1,k=ijm,k1+μΔt+(dm+γm+δσm)Δt,j=0,,N1; m=0,,M1, ij+10,k=kqvUj+1kλjv+quVj+1kλju,j=0,,N1; Nj+1k=Njk+Δt(ΛduUj+1k(dv+γv)Vj+1kMm=1(dm+γm)ijm,k)1+μΔt,j=0,,N1.
    $
    (5.1)

    To study the coexistence equilibrium analytically is not feasible for this model. So we take the help of simulations to predict the existence of coexistence equilibrium in a scale free network scenario. We consider specific parameter values for which $ \mathcal{R}^{2}_{u_i} $ and $ \mathcal{R}^{1}_{v_i} $ are greater than 1. The simulations suggest that the coexistence equilibrium exists and is stable. Given parameter values are constant, the invasion number of HIV, $ \mathcal{R}^{1}_{v_i} $, seem to show dependence on the size of the network used. With the same parameter values given in Table 3, when the network contains 200 nodes, $ \mathcal{R}^{1}_{v_i} $ is close to $ 1.4 $, while with 300 nodes, $ \mathcal{R}^{1}_{v_i} $ increases to $ 2.9 $. The invasion number of opioid addiction $ \mathcal{R}^{1}_{v_i} $ remains stable near the same value $ 1.2 $ when size of the network is increased from 200 nodes to 300 nodes. This is to be expected since the spread of opioid has been considered homogeneous over the network. Simulations suggesting a coexistence equilibrium are shown in Figure 1.

    Table 1.  List of parameters of the within-host model.
    Notation Meaning
    $ s $ $ \text{Production rate of healthy T-cells} $
    $ d_T $ $ \text{Clearance rate of healthy T-cells} $
    $ T $ $ \text{Number of healthy cells} $
    $ T_{i} $ $ \text{Number of infected cells} $
    $ V_i $ $ \text{Number of virions} $
    $ \delta_i $ $ \text{Death rate of infected cells} $
    $ c $ $ \text{Clearance rate of virions} $
    $ k_0 $ $ \text{Transmission coeffiecient in absence of opioid} $
    $ k_1 $ $ \text{Maximal increase in infection rate due to opioid} $
    $ C_0 $ $ \text{Half saturation constant} $
    $ d_c $ $ \text{Rate of degradation of opioid} $

     | Show Table
    DownLoad: CSV
    Table 2.  Definitions of parameters and dependent variables of the between-host model.
    Parameter/Variable Description
    $ S_{k}(t) $ $ \text{Number of susceptible individuals at time $t$ with degree $k$ } $
    $ V_{k}(t) $ $ \text{Number of HIV infected individuals at time $t$ with degree $k$ } $
    $ U_{k}(t) $ $ \text{Number of opioid addicted individuals at time $t$ with degree $k$ } $
    $ i_{k}(t, \tau) $ $ \text{Density of co-affected individuals with coinfection age $\tau$ at time $t$ with degree $k$ } $
    $ A_{k}(t) $ $ \text{Number of individuals with AIDS at time $t$ with degree $k$ } $
    $ \Lambda_k $ $ \text{Constant recruitment rate for nodes with degree $k$ } $
    $ \beta_u $ $ \text{Transmission rate of opioid addiction} $
    $ \beta_{v_1} $ $ \text{Transmission rate of HIV infection of HIV infected only individuals} $
    $ \beta_{v_2}(\tau) $ $ \text{Transmission rate of HIV infection of coinfected individuals} $
    $ \mu $ $ \text{Natural death rate} $
    $ d_u $ $ \text{Death rate induced by opioid addiction} $
    $ d_v $ $ \text{Death rate induced by HIV infection} $
    $ d_i(\tau) $ $ \text{Death rate induced by coaffection at coaffection age $\tau$ } $
    $ d_a $ $ \text{Death rate induced by AIDS} $
    $ \delta $ $ \text{Recovery rate from opioid addiction} $
    $ \sigma(\tau) $ $ \text{Recovery rate from coinfection to HIV only} $
    $ q_u $ $ \text{Increase coefficient of HIV infection due to opioid usage} $
    $ q_v $ $ \text{Increase coefficient of opioid usage due to HIV infection} $
    $ \gamma_v $ $ \text{Rate of transition from HIV to AIDS} $
    $ \gamma_i(\tau) $ $ \text{Rate of transition from coinfection to AIDS} $

     | Show Table
    DownLoad: CSV
    Table 3.  Parameter estimation results from [38].
    Parameter Estimated Value Units
    $ \beta_u $ 0.385676 1/time
    $ \beta_{v_1} $ 0.0551 1/time
    $ k_0 $ 0.00011046 1/time
    $ B $ 15318.9 vRNA/ml
    $ \delta $ 0.118227 1/time
    $ q_u $ 0.867138 Unitless
    $ q_v $ 30.6189 Unitless
    $ d_u $ 0.00817752 1/time
    $ d_v $ 0.0144092 1/time
    $ d_a $ 1.2766e+11 1/time
    $ d_0 $ 2.72895e-07 ml/($ \text{time} \times \text{cells} $)
    $ d_1 $ 3.4671e-06 1/time
    $ \gamma_v $ 0.0223488 1/time
    $ \gamma_0 $ 1.63927e-12 1/time
    $ \sigma $ 0.000270006 Unitless
    $ s $ 22843.6 CD4 count/($ \text{time} \times \text{ml} $)
    $ d $ 0.0766824 1/time
    $ k_1 $ 2.02785e-05 vRNA/($ \text{CD4 count} \times \text{time} $)
    $ \delta_i $ 0.725266 1/time
    $ N_v \delta_i $ 8465.63 vRNA/($ \text{CD4 count} \times \text{time} $)

     | Show Table
    DownLoad: CSV
    Figure 1.  Simulations with the network model. In this simulation the average degree is 7.63 and $ \mathcal{R}^{1}_{vi} = 1.2173 $, $ \mathcal{R}^{2}_{ui} = 1.1607 $. The number of nodes is 100. The maximal degree is 28 but there are no occupied nodes with degree 1, 2, 3.

    We define U(t) = $ \sum_{k = 1}^{n}U_k(t) $ as the total opioid addicted population. Similarly we define $ V(t) $, $ I(t) $ and $ S(t) $ as the total HIV infected, total co-infected and total susceptible population respectively. The network utilized had $ 200 $ nodes. The parameters $ \beta_u $ and $ \beta^{k}_{v_1} $ are estimated by the following formulas, $ \beta_{u} = \mathcal{R}_u(\mu + d_u + \delta); $ and $ \beta^{k}_{v1} = \mathcal{R}_{v}(\mu+d_v+\gamma_v) $. Since the model we consider does not include treatment for HIV, we consider $ \mathcal{R}_v = 5.5 $. This estimate is an average value collected from [39], which gives an estimation of basic reproduction number of HIV in Rural South West Uganda. Estimated value of $ \mathcal{R}_u $ according to [27] would be close to $ 1.1 $. Given the fact that a high percentage of US citizens mentioned it would be easy for someone to access opioids for illicit purposes, according to a poll conducted in 2018 (46 percent) and people who misuse opioids often get them from a family member or friend who has a prescription [40], we considered the $ \mathcal{R}_u $ to be higher in range, around $ 3.25 $. The maximal degree $ n $ for the following simulations (except Figure 1) is $ 43 $, with number of nodes being $ 200 $.

    The two parameters in (2.2) that are of particular interest are $ q_u $ and $ q_v $, which determine how much one epidemics impacts the other. That is $ q_v $ determines how likely opioid users are to get infected by HIV compared to non-users, and $ q_u $ determines how likely HIV infected people are to become opioid addicted compared to non-infected people. In [27] the estimated value of the $ q_v $ term equivalent was 94.5, and in [38], the estimated value of $ q_v $ was 30.62. Both estimates suggest a high dependence effect of opioid usage on HIV infection. In our simulations we take the lesser estimate of $ q_v = 30.62 $. The total co-affected population due to varying values of $ q_v $ is simulated, such as $ 0.5 q_v, q_v, 1.5q_v, 2q_v, 2.5q_v, 3q_v. $ The estimated value for $ q_u $ in [38] was below 1, but HIV-infected persons are more likely to have chronic pain, receive opioid analgesic treatment, receive higher doses of opioids, and to have substance use disorders and mental illness compared with the general population, putting them at increased risk for opioid use disorder [41]. So the simulations were done with the fitted value for $ q_u $ along with double, five times and ten times the fitted value. While the total number of co-affected varies according to the network, the trend seems to be similar, with $ q_v $ increasing, the total number of co-affected increases (Figure 2). A definite situation of interest is when $ q_u = 4.3 $, the maximum value seems to be achieved when $ q_v $ is close to $ 60 $, not at the highest value of approximately $ 90 $. The simulation was repeated with these values for different network sizes and provided the same result.

    Figure 2.  Figure shows total co-affected individuals for six different values of $ q_v $. Top Left: $ q_u = 0.86 $, Top Right: $ q_u = 1.72 $, Bottom Left: $ q_u = 4.3 $, Bottom Right: $ q_u = 8.6 $. The other parameters used are given in Table 3.

    A second set of simulations were performed, taking the base value of $ q_u = 1.72 $. The four individual cases have all other parameters and network values same, only $ q_v $ is varied, from $ 15, 31, 62 $ and $ 93 $ respectively. While some of the cases do have permutations, over all the trend is similar, and opposite to the simulations in Figure 2. That is $ q_u $ increasing causes the total number of co-affected to decrease (Figure 3). Again we notice a situation of interest, when $ q_v = 60 $, the maximum value seems to be achieved when $ q_u $ is close to $ 1.72 $, not at the lowest value of approximately $ 0.86 $. Repeated simulation with differing sized networks did not show changes in the results.

    Figure 3.  Figure shows total co-affected individuals for six different values of $ q_u $, Top Left: $ q_v = 15 $, Top Right: $ q_v = 30 $, Bottom Left: $ q_v = 60 $, Bottom Right: $ q_v = 90 $. The other parameters used are given in Table 3.

    Another parameter of interest is $ \delta $, which denotes the rate of recovery from addiction in the epidemiological model. The estimate for that in [27] is close to $ 0.033 $, while in [38] the estimate is approximately $ 0.11 $. If successful treatment is considered without subtracting the relapses, the rate would probably be close to $ 0.05 $ [42]. To simulate for differing values of $ \delta $, we consider $ \beta_u $ and $ \beta^{k}_{v1} $ as constants, directly taking the values from Table 3. All the other parameter values are as mentioned in Table 3. We consider four differing situations for $ \delta $, with the value being $ 0.02 $, the fitted value $ 0.11 $, and target high values of $ 0.5 $ and $ \delta = 1 $. We also investigated different scenarios with differing network sizes, with number of nodes being $ 100 $, $ 300 $ and $ 500 $ respectively. Interestingly, irrespective of network size, the total number of co-affected people appeared to be higher with the value of $ \delta $ increasing (Figure 4). One quite plausible explanation would be the higher recovery would decrease the number of opioid overdose deaths significantly, thereby increasing the co-affected prevalence.

    Figure 4.  Figure shows total co-affected individuals for four different values of $ \delta $, Top Left: Network Size 100, Top Right: Network Size 300, Bottom: Network size 500. The other parameters used are given in Table 3.

    We formulate a within-host model linked with a dynamic network HIV/opioid coinfection epidemiological model with demography, through epidemiological parameters. The system is described by ordinary differential equations coupled with partial differential equations in a nested fashion. The network multi-scale model here is an extension of the multi-scale model considered in [38]. The disease free equilibrium of the system always exists and is locally asymptotically stable when both the basic reproduction numbers of opioid and HIV, $ \mathcal{R}_u $ and $ \mathcal{R}_v $ are less than 1.

    The boundary equilibrium $ E_{1}^* $ exists when $ \mathcal{R}_u $ is more than 1 and $ E_{2}^* $ exist when $ \mathcal{R}_v $ is more than 1. We define the invasion reproduction numbers $ \mathcal{R}^{1}_{v_i} $ and $ \mathcal{R}^{2}_{u_i} $. The invasion reproduction number $ \mathcal{R}^{2}_{u_i} $ gives the reproduction of the opioid users when the population is at the equilibrium $ E_2^{*} $, that is, when HIV infection alone is at equilibrium in the single population. The invasion reproduction number $ \mathcal{R}^{1}_{v_i} $ gives the reproduction of the HIV infection at the equilibrium $ E_1^{*} $, that is when the opioid transmission alone is at equilibrium in the single population. When $ \mathcal{R}^{1}_{v_i} < 1 $, $ E_{1}^* $ is locally stable and when $ \mathcal{R}^{1}_{v_i} > 1 $, $ E_{1}^* $ is unstable. When $ \mathcal{R}^{2}_{u_i} < 1 $, $ E_{2}^* $ is locally stable and when $ \mathcal{R}^{2}_{u_i} > 1 $, $ E_{2}^* $ is unstable. The model is too complicated to compute or consider the stability of an endemic equilibrium, analytically, but simulations suggest that there is an interior equilibrium potentially under the condition that both invasion numbers are larger than one.

    We use fitted parameters from [38], to perform simulations, to explore the effect of the change of the parameters $ q_u $, $ q_v $ and $ \delta $. The parameters $ q_u $ and $ q_v $ represent the effect of one epidemic on the other, and we simulated for plausible values of $ q_u $ and $ q_v $, to get estimates of co-affected prevalence. The estimate of $ q_v $, the likelihood of a heroin user to be infected with HIV, in [27] was approximately 94, and since a large number of opioid users progress to becoming heroin users, (Data from 2011 showed that an estimated 4 to 6 percent who misuse prescription opioids switch to heroin and about 80 percent of people who used heroin first misused prescription opioids) [43], chances of the $ q_v $ estimate for all illicit opioid usage being close to the "heroin only" estimate is high. We notice that increase in $ q_v $ causes the co-affected population to rise significantly, which indicates that control strategies focusing on reducing HIV infection among the opioid addicted population would be effective in decoupling the epidemics, corroborating with the conclusions in [27].

    To the best of our knowledge, there have not been previous models with a parameter similar to $ q_u $, that provides an estimate for the effect of HIV infection on opioid use. The simulations from our model show that increase of $ q_u $ in general causes the number of co-affected to decline. That is the prevalence of co-affected people declines sharply with the increase to more HIV infected people being addicted to opioids. This does corroborate real life data, since deaths due to overdose in the US per year were aproximately 120,000 in the year 2020, and the number increased 15 percent by 2021 [44], compared to the number of deaths due to HIV being around 18,500 in 2020 [45]. We noticed significant increase in the prevalence of co-affected individuals, when the parameter $ \delta $, representing the recovery from opioid usage was increased significantly. The sharp increase points to the idea that control measures should focus more on treatment of opioid use disorder, and that uncoupling the two epidemics is a priority to prevent loss of human lives. Counseling about the dangers of opioid addiction for HIV infected people must be provided, and similarly counseling about the dangers of getting infected with HIV should be provided to reported opioid addicts.

    In summary, we have developed a novel multi-scale network model of HIV and opioid epidemics. We have analized the model and obtained conditions for HIV-only to persist or opioid-only to persist. Simulations suggest that the two epidemics can co-exist for some parameter values. Simulations further suggest that decreasing $ q_v $ decreases the number of co-affected and may lead to decopuling the epidemics. Thus control measures targeted at reducing $ q_v $ should be coupled with treatment of opioid affected individuals which is consistent with our previois findings.

    Maia Martcheva was partially supported by NSF DMS-1951595. Necibe Tuncer was partially supported by NSF DMS-1951626.

    The authors declare there is no conflict of interest.

    The code was written by one of the authors and is available on request by contacting Churni Gupta.

    [1] Herbst RS, Soria JC, Kowanetz M, et al. (2014) Predictive correlates of response to the anti-PD-L1 antibody MPDL3280A in cancer patients. Nature 515: 563-567. doi: 10.1038/nature14011
    [2] Tumeh PC, Harview CL, Yearley JH, et al. (2014) PD-1 blockade induces responses by inhibiting adaptive immune resistance. Nature 515: 568-571. doi: 10.1038/nature13954
    [3] Ansell SM, Lesokhin AM, Borrello I, et al. (2015) PD-1 blockade with nivolumab in relapsed or refractory Hodgkin's lymphoma. N Engl J Med 372: 311-319. doi: 10.1056/NEJMoa1411087
    [4] Coley WB (1907) Sarcoma of the Long Bones: The Diagnosis, Treatment and Prognosis, with a Report of Sixty-Nine Cases. Ann Surg 45: 321-368.
    [5] Topalian SL, Weiner GJ, Pardoll DM (2011) Cancer immunotherapy comes of age. J Clin Oncol 29: 4828-4836. doi: 10.1200/JCO.2011.38.0899
    [6] Grosso JF, Jure-Kunkel MN (2013) CTLA-4 blockade in tumor models: an overview of preclinical and translational research. Cancer Immun 13: 5.
    [7] Pardoll DM (2012) The blockade of immune checkpoints in cancer immunotherapy. Nat Rev Cancer 12: 252-264. doi: 10.1038/nrc3239
    [8] de Pasquale MD, Castellano A, de Sio L, et al. (2011) Bevacizumab in pediatric patients: how safe is it? Anticancer Res 31: 3953-3957.
    [9] Okada K, Yamasaki K, Tanaka C, et al. (2013) Phase I study of bevacizumab plus irinotecan in pediatric patients with recurrent/refractory solid tumors. Jpn J Clin Oncol 43: 1073-1079. doi: 10.1093/jjco/hyt124
    [10] Ebb D, Meyers P, Grier H, et al. (2012) Phase II trial of trastuzumab in combination with cytotoxic chemotherapy for treatment of metastatic osteosarcoma with human epidermal growth factor receptor 2 overexpression: a report from the children's oncology group. J Clin Oncol 30: 2545-2551. doi: 10.1200/JCO.2011.37.4546
    [11] Fry TJ, Lankester AC (2010) Cancer immunotherapy: will expanding knowledge lead to success in pediatric oncology? Hematol Oncol Clin North Am 24: 109-127. doi: 10.1016/j.hoc.2009.11.010
    [12] Dalziel M, Crispin M, Scanlan CN, et al. (2014) Emerging principles for the therapeutic exploitation of glycosylation. Science 343: 1235681. doi: 10.1126/science.1235681
    [13] Kanda Y, Yamada T, Mori K, et al. (2007) Comparison of biological activity among nonfucosylated therapeutic IgG1 antibodies with three different N-linked Fc oligosaccharides: the high-mannose, hybrid, and complex types. Glycobiology 17: 104-118.
    [14] Morris MJ, Divgi CR, Pandit-Taskar N, et al. (2005) Pilot trial of unlabeled and indium-111-labeled anti-prostate-specific membrane antigen antibody J591 for castrate metastatic prostate cancer. Clin Cancer Res 11: 7454-7461. doi: 10.1158/1078-0432.CCR-05-0826
    [15] Linden O, Hindorf C, Cavallin-Stahl E, et al. (2005) Dose-fractionated radioimmunotherapy in non-Hodgkin's lymphoma using DOTA-conjugated, 90Y-radiolabeled, humanized anti-CD22 monoclonal antibody, epratuzumab. Clin Cancer Res 11: 5215-5222. doi: 10.1158/1078-0432.CCR-05-0172
    [16] Sharkey RM, Brenner A, Burton J, et al. (2003) Radioimmunotherapy of non-Hodgkin's lymphoma with 90Y-DOTA humanized anti-CD22 IgG (90Y-Epratuzumab): do tumor targeting and dosimetry predict therapeutic response? J Nucl Med 44: 2000-2018.
    [17] Boerman OC, Koppe MJ, Postema EJ, et al. (2007) Radionuclide therapy of cancer with radiolabeled antibodies. Anticancer Agents Med Chem 7: 335-343. doi: 10.2174/187152007780618126
    [18] Spigel DR, Ervin TJ, Ramlau RA, et al. (2013) Randomized phase II trial of Onartuzumab in combination with erlotinib in patients with advanced non-small-cell lung cancer. J Clin Oncol 31: 4105-4114. doi: 10.1200/JCO.2012.47.4189
    [19] Kim ES, Neubauer M, Cohn A, et al. (2013) Docetaxel or pemetrexed with or without cetuximab in recurrent or progressive non-small-cell lung cancer after platinum-based therapy: a phase 3, open-label, randomised trial. Lancet Oncol 14: 1326-1336. doi: 10.1016/S1470-2045(13)70473-X
    [20] Attias D, Weitzman S (2008) The efficacy of rituximab in high-grade pediatric B-cell lymphoma/leukemia: a review of available evidence. Curr Opin Pediatr 20: 17-22. doi: 10.1097/MOP.0b013e3282f424b0
    [21] Barth MJ, Goldman S, Smith L, et al. (2013) Rituximab pharmacokinetics in children and adolescents with de novo intermediate and advanced mature B-cell lymphoma/leukaemia: a Children's Oncology Group report. Br J Haematol 162: 678-683. doi: 10.1111/bjh.12434
    [22] Moreno L, Vaidya SJ, Pinkerton CR, et al. (2013) Long-term follow-up of children with high-risk neuroblastoma: the ENSG5 trial experience. Pediatr Blood Cancer 60: 1135-1140. doi: 10.1002/pbc.24452
    [23] Yalcin B, Kremer LC, Caron HN, et al. (2013) High-dose chemotherapy and autologous haematopoietic stem cell rescue for children with high-risk neuroblastoma. Cochrane Database Syst Rev 8: CD006301.
    [24] Yu AL, Gilman AL, Ozkaynak MF, et al. (2010) Anti-GD2 antibody with GM-CSF, interleukin-2, and isotretinoin for neuroblastoma. N Engl J Med 363: 1324-1334. doi: 10.1056/NEJMoa0911123
    [25] Ozkaynak MF, Sondel PM, Krailo MD, et al.(2000) Phase I study of chimeric human/murine anti-ganglioside G(D2) monoclonal antibody (ch14.18) with granulocyte-macrophage colony-stimulating factor in children with neuroblastoma immediately after hematopoietic stem-cell transplantation: a Children's Cancer Group Study. J clin oncol 18: 4077-4085.
    [26] Chen X, Soma LA, Fromm JR (2013) Targeted therapy for Hodgkin lymphoma and systemic anaplastic large cell lymphoma: focus on brentuximab vedotin. Onco Targets Ther 7: 45-56.
    [27] Brotelle T, Lemal R, Molucon-Chabrot C, et al. (2014) [Gemtuzumab ozogamicin for treatment of acute myeloid leukemia]. Bull Cancer 101: 211-218.
    [28] Daver N, O'Brien S (2013) Novel therapeutic strategies in adult acute lymphoblastic leukemia--a focus on emerging monoclonal antibodies. Curr Hematol Malig Rep 8: 123-131. doi: 10.1007/s11899-013-0160-7
    [29] Kreitman RJ, Pastan I (2011) Antibody fusion proteins: anti-CD22 recombinant immunotoxin moxetumomab pasudotox. Clin Cancer Res 17: 6398-6405. doi: 10.1158/1078-0432.CCR-11-0487
    [30] Mussai F, Campana D, Bhojwani D, et al. (2010) Cytotoxicity of the anti-CD22 immunotoxin HA22 (CAT-8015) against paediatric acute lymphoblastic leukaemia. Br J Haematol 150: 352-358. doi: 10.1111/j.1365-2141.2010.08251.x
    [31] Thomas X (2014) Blinatumomab: a new era of treatment for adult ALL? Lancet Oncol 16: 6-7
    [32] Topp MS, Kufer P, Gokbuget N, et al. (2011) Targeted therapy with the T-cell-engaging antibody blinatumomab of chemotherapy-refractory minimal residual disease in B-lineage acute lymphoblastic leukemia patients results in high response rate and prolonged leukemia-free survival. J Clin Oncol 29: 2493-2498. doi: 10.1200/JCO.2010.32.7270
    [33] Shah AH, Bregy A, Heros DO, et al. (2013) Dendritic cell vaccine for recurrent high-grade gliomas in pediatric and adult subjects: clinical trial protocol. Neurosurgery 73: 863-867. doi: 10.1227/NEU.0000000000000107
    [34] Ciocca DR, Cayado-Gutierrez N, Maccioni M, et al. (2012) Heat shock proteins (HSPs) based anti-cancer vaccines. Curr Mol Med 12: 1183-1197. doi: 10.2174/156652412803306684
    [35] Yang I, Fang S, Parsa AT (2010) Heat shock proteins in glioblastomas. Neurosurg Clin N Am 21: 111-123. doi: 10.1016/j.nec.2009.09.002
    [36] Graner MW, Bigner DD (2006) Therapeutic aspects of chaperones/heat-shock proteins in neuro-oncology. Expert Rev Anticancer Ther 6: 679-695. doi: 10.1586/14737140.6.5.679
    [37] Hodi FS, O'Day SJ, McDermott DF, et al. (2010) Improved survival with ipilimumab in patients with metastatic melanoma. N Engl J Med 363: 711-723. doi: 10.1056/NEJMoa1003466
    [38] Kantoff PW, Schuetz TJ, Blumenstein BA, et al. (2010) Overall survival analysis of a phase II randomized controlled trial of a Poxviral-based PSA-targeted immunotherapy in metastatic castration-resistant prostate cancer. J Clin Oncol 28: 1099-1105. doi: 10.1200/JCO.2009.25.0597
    [39] Heimberger AB, Sampson JH (2009) The PEPvIII-KLH (CDX-110) vaccine in glioblastoma multiforme patients. Expert Opin Biol Ther 9: 1087-1098. doi: 10.1517/14712590903124346
    [40] Fest S, Huebener N, Bleeke M, et al. (2009) Survivin minigene DNA vaccination is effective against neuroblastoma. Int J Cancer 125: 104-114. doi: 10.1002/ijc.24291
    [41] Gilboa E (2007) DC-based cancer vaccines. J Clin Invest 117: 1195-1203. doi: 10.1172/JCI31205
    [42] Napolitani G, Rinaldi A, Bertoni F, et al. (2005) Selected Toll-like receptor agonist combinations synergistically trigger a T helper type 1-polarizing program in dendritic cells. Nat Immunol 6: 769-776. doi: 10.1038/ni1223
    [43] Shen L, Evel-Kabler K, Strube R, et al. (2004) Silencing of SOCS1 enhances antigen presentation by dendritic cells and antigen-specific anti-tumor immunity. Nat biotechnol 22: 1546-1553. doi: 10.1038/nbt1035
    [44] Gilboa E (2004) Knocking the SOCS1 off dendritic cells. Nature biotechnology 22: 1521-1522. doi: 10.1038/nbt1204-1521
    [45] Cohen N, Mouly E, Hamdi H, et al. (2006) GILZ expression in human dendritic cells redirects their maturation and prevents antigen-specific T lymphocyte response. Blood 107: 2037-2044. doi: 10.1182/blood-2005-07-2760
    [46] Konkankit VV, Kim W, Koya RC, et al. (2011) Decitabine immunosensitizes human gliomas to NY-ESO-1 specific T lymphocyte targeting through the Fas/Fas ligand pathway. J Transl Med 9: 192. doi: 10.1186/1479-5876-9-192
    [47] Chou J, Voong LN, Mortales CL, et al. (2012) Epigenetic modulation to enable antigen-specific T-cell therapy of colorectal cancer. J Immunother 35: 131-141. doi: 10.1097/CJI.0b013e31824300c7
    [48] Chang CN, Huang YC, Yang DM, et al. (2011) A phase I/II clinical trial investigating the adverse and therapeutic effects of a postoperative autologous dendritic cell tumor vaccine in patients with malignant glioma. J Clin Neurosci 18: 1048-1054. doi: 10.1016/j.jocn.2010.11.034
    [49] Cho DY, Yang WK, Lee HC, et al. (2012) Adjuvant immunotherapy with whole-cell lysate dendritic cells vaccine for glioblastoma multiforme: a phase II clinical trial. World Neurosurg 77: 736-744. doi: 10.1016/j.wneu.2011.08.020
    [50] Liau LM, Prins RM, Kiertscher SM, et al. (2005) Dendritic cell vaccination in glioblastoma patients induces systemic and intracranial T-cell responses modulated by the local central nervous system tumor microenvironment. Clin Cancer Res 11: 5515-5525. doi: 10.1158/1078-0432.CCR-05-0464
    [51] Phuphanich S, Wheeler CJ, Rudnick JD, et al. (2013) Phase I trial of a multi-epitope-pulsed dendritic cell vaccine for patients with newly diagnosed glioblastoma. Cancer Immunol Immunother 62: 125-135. doi: 10.1007/s00262-012-1319-0
    [52] Wheeler CJ, Black KL, Liu G, et al. (2008) Vaccination elicits correlated immune and clinical responses in glioblastoma multiforme patients. Cancer Res 68: 5955-5964. doi: 10.1158/0008-5472.CAN-07-5973
    [53] Yu JS, Wheeler CJ, Zeltzer PM, et al. (2001) Vaccination of malignant glioma patients with peptide-pulsed dendritic cells elicits systemic cytotoxicity and intracranial T-cell infiltration. Cancer Res 61: 842-847.
    [54] Dannull J, Haley NR, Archer G, et al. (2013) Melanoma immunotherapy using mature DCs expressing the constitutive proteasome. J Clin Invest 123: 3135-3145.
    [55] Slingluff CL, Jr., Lee S, Zhao F, et al. (2013) A randomized phase II trial of multiepitope vaccination with melanoma peptides for cytotoxic T cells and helper T cells for patients with metastatic melanoma (E1602). Clin Cancer Res 19: 4228-4238. doi: 10.1158/1078-0432.CCR-13-0002
    [56] Dillman R, Barth N, Selvan S, et al. (2004) Phase I/II trial of autologous tumor cell line-derived vaccines for recurrent or metastatic sarcomas. Cancer Biother Radiopharm 19: 581-588. doi: 10.1089/1084978042484812
    [57] Perroud MW, Jr., Honma HN, Barbeiro AS, et al. (2011) Mature autologous dendritic cell vaccines in advanced non-small cell lung cancer: a phase I pilot study. J Exp Clin Cancer Res 30: 65. doi: 10.1186/1756-9966-30-65
    [58] Navada SC, Steinmann J, Lubbert M, et al. (2014) Clinical development of demethylating agents in hematology. J Clin Invest 124: 40-46. doi: 10.1172/JCI69739
    [59] Burdach S, van Kaick B, Laws HJ, et al. (2000) Allogeneic and autologous stem-cell transplantation in advanced Ewing tumors. An update after long-term follow-up from two centers of the European Intergroup study EICESS. Stem-Cell Transplant Programs at Dusseldorf University Medical Center, Germany and St. Anna Kinderspital, Vienna, Austria. Ann Oncol 11: 1451-1462.
    [60] Wu R, Forget MA, Chacon J, et al. (2012) Adoptive T-cell therapy using autologous tumor-infiltrating lymphocytes for metastatic melanoma: current status and future outlook. Cancer J 18: 160-175. doi: 10.1097/PPO.0b013e31824d4465
    [61] Bridgeman JS, Hawkins RE, Hombach AA, et al. (2010) Building better chimeric antigen receptors for adoptive T cell therapy. Curr Gene Ther 10: 77-90. doi: 10.2174/156652310791111001
    [62] Grupp SA, Prak EL, Boyer J, et al. (2012) Adoptive transfer of autologous T cells improves T-cell repertoire diversity and long-term B-cell function in pediatric patients with neuroblastoma. Clin Cancer Res 18: 6732-6741. doi: 10.1158/1078-0432.CCR-12-1432
    [63] Peres E, Wood GW, Poulik J, et al. (2008) High-dose chemotherapy and adoptive immunotherapy in the treatment of recurrent pediatric brain tumors. Neuropediatrics 39: 151-156. doi: 10.1055/s-0028-1093333
    [64] Bao L, Cowan MJ, Dunham K, et al. (2012) Adoptive immunotherapy with CMV-specific cytotoxic T lymphocytes for stem cell transplant patients with refractory CMV infections. J Immunother 35: 293-298. doi: 10.1097/CJI.0b013e31824300a2
    [65] Park JR, Digiusto DL, Slovak M, et al. (2007) Adoptive transfer of chimeric antigen receptor re-directed cytolytic T lymphocyte clones in patients with neuroblastoma. Mol Ther 15: 825-833.
    [66] Ramos CA, Savoldo B, Dotti G (2014) CD19-CAR Trials. Cancer J 20: 112-118. doi: 10.1097/PPO.0000000000000031
    [67] Maude SL, Frey N, Shaw PA, et al. (2014) Chimeric antigen receptor T cells for sustained remissions in leukemia. N Engl J Med 371: 1507-1517. doi: 10.1056/NEJMoa1407222
    [68] Tasian SK, Pollard JA, Aplenc R (2014) Molecular Therapeutic Approaches for Pediatric Acute Myeloid Leukemia. Front Oncol 4: 55.
    [69] Hoffman LM, Gore L (2014) Blinatumomab, a Bi-Specific Anti-CD19/CD3 BiTE((R)) Antibody for the Treatment of Acute Lymphoblastic Leukemia: Perspectives and Current Pediatric Applications. Front Oncol 4: 63.
    [70] Hamanishi J, Mandai M, Iwasaki M, et al. (2007) Programmed cell death 1 ligand 1 and tumor-infiltrating CD8+ T lymphocytes are prognostic factors of human ovarian cancer. Proc Natl Acad Sci U S A 104: 3360-3365. doi: 10.1073/pnas.0611533104
    [71] Dong H, Zhu G, Tamada K, et al. (1999) B7-H1, a third member of the B7 family, co-stimulates T-cell proliferation and interleukin-10 secretion. Nat Med 5: 1365-1369. doi: 10.1038/70932
    [72] Toomer KH, Chen Z (2014) Autoimmunity as a Double Agent in Tumor Killing and Cancer Promotion. Front Immunol 5: 116.
    [73] Wolchok JD, Chan TA (2014) Cancer: Antitumour immunity gets a boost. Nature 515: 496-498. doi: 10.1038/515496a
    [74] Parsa AT, Waldron JS, Panner A, et al. (2007) Loss of tumor suppressor PTEN function increases B7-H1 expression and immunoresistance in glioma. Nat Med 13: 84-88. doi: 10.1038/nm1517
    [75] Shen JK, Cote GM, Choy E, et al. (2014) Programmed cell death ligand 1 expression in osteosarcoma. Cancer Immunol Res 2: 690-698. doi: 10.1158/2326-6066.CIR-13-0224
    [76] Gilboa E, McNamara J, 2nd, Pastor F (2013) Use of oligonucleotide aptamer ligands to modulate the function of immune receptors. Clin Cancer Res 19: 1054-1062. doi: 10.1158/1078-0432.CCR-12-2067
    [77] Chatterton Z, Burke D, Emslie KR, et al. (2014) Validation of DNA methylation biomarkers for diagnosis of acute lymphoblastic leukemia. Clin Chem 60: 995-1003. doi: 10.1373/clinchem.2013.219956
    [78] Ahsan S, Raabe EH, Haffner MC, et al. (2014) Increased 5-hydroxymethylcytosine and decreased 5-methylcytosine are indicators of global epigenetic dysregulation in diffuse intrinsic pontine glioma. Acta Neuropathol Commun 2: 59. doi: 10.1186/2051-5960-2-59
  • This article has been cited by:

    1. Chelsea Spence, Mary E. Kurz, Thomas C. Sharkey, Bryan Lee Miller, Scoping Literature Review of Disease Modeling of the Opioid Crisis, 2024, 0279-1072, 1, 10.1080/02791072.2024.2367617
    2. Eric Numfor, Necibe Tuncer, Maia Martcheva, Optimal control of a multi-scale HIV-opioid model, 2024, 18, 1751-3758, 10.1080/17513758.2024.2317245
  • Reader Comments
  • © 2015 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5468) PDF downloads(1105) Cited by(1)

Figures and Tables

Tables(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog