Multiplayer games and HIV transmission via casual encounters

  • Received: 27 January 2015 Revised: 27 June 2016 Published: 01 April 2017
  • MSC : Primary: 91A80; Secondary: 91C99

  • Population transmission models have been helpful in studying the spread of HIV. They assess changes made at the population level for different intervention strategies. To further understand how individual changes affect the population as a whole, game-theoretical models are used to quantify the decision-making process. Investigating multiplayer nonlinear games that model HIV transmission represents a unique approach in epidemiological research. We present here 2-player and multiplayer noncooperative games where players are defined by HIV status and age and may engage in casual (sexual) encounters. The games are modelled as generalized Nash games with shared constraints, which is completely novel in the context of our applied problem. Each player's HIV status is known to potential partners, and players have personal preferences ranked via utility values of unprotected and protected sex outcomes. We model a player's strategy as their probability of being engaged in a casual unprotected sex encounter (USE), which may lead to HIV transmission; however, we do not incorporate a transmission model here. We study the sensitivity of Nash strategies with respect to varying preference rankings, and the impact of a prophylactic vaccine introduced in players of youngest age groups. We also study the effect of these changes on the overall increase in infection level, as well as the effects that a potential prophylactic treatment may have on age-stratified groups of players. We conclude that the biggest impacts on increasing the infection levels in the overall population are given by the variation in the utilities assigned to individuals for unprotected sex with others of opposite HIV status, while the introduction of a prophylactic vaccine in youngest age group (15-20 yr olds) slows down the increase in HIV infection.

    Citation: Stephen Tully, Monica-Gabriela Cojocaru, Chris T. Bauch. Multiplayer games and HIV transmission via casual encounters[J]. Mathematical Biosciences and Engineering, 2017, 14(2): 359-376. doi: 10.3934/mbe.2017023

    Related Papers:

    [1] Raul Borsche, Axel Klar, T. N. Ha Pham . Nonlinear flux-limited models for chemotaxis on networks. Networks and Heterogeneous Media, 2017, 12(3): 381-401. doi: 10.3934/nhm.2017017
    [2] Markus Musch, Ulrik Skre Fjordholm, Nils Henrik Risebro . Well-posedness theory for nonlinear scalar conservation laws on networks. Networks and Heterogeneous Media, 2022, 17(1): 101-128. doi: 10.3934/nhm.2021025
    [3] Michael Herty, Niklas Kolbe, Siegfried Müller . Central schemes for networked scalar conservation laws. Networks and Heterogeneous Media, 2023, 18(1): 310-340. doi: 10.3934/nhm.2023012
    [4] Darko Mitrovic . Existence and stability of a multidimensional scalar conservation law with discontinuous flux. Networks and Heterogeneous Media, 2010, 5(1): 163-188. doi: 10.3934/nhm.2010.5.163
    [5] Raimund Bürger, Harold Deivi Contreras, Luis Miguel Villada . A Hilliges-Weidlich-type scheme for a one-dimensional scalar conservation law with nonlocal flux. Networks and Heterogeneous Media, 2023, 18(2): 664-693. doi: 10.3934/nhm.2023029
    [6] José Antonio Carrillo, Yingping Peng, Aneta Wróblewska-Kamińska . Relative entropy method for the relaxation limit of hydrodynamic models. Networks and Heterogeneous Media, 2020, 15(3): 369-387. doi: 10.3934/nhm.2020023
    [7] Qing Li, Steinar Evje . Learning the nonlinear flux function of a hidden scalar conservation law from data. Networks and Heterogeneous Media, 2023, 18(1): 48-79. doi: 10.3934/nhm.2023003
    [8] Panpan Xu, Yongbin Ge, Lin Zhang . High-order finite difference approximation of the Keller-Segel model with additional self- and cross-diffusion terms and a logistic source. Networks and Heterogeneous Media, 2023, 18(4): 1471-1492. doi: 10.3934/nhm.2023065
    [9] Raimund Bürger, Kenneth H. Karlsen, John D. Towers . On some difference schemes and entropy conditions for a class of multi-species kinematic flow models with discontinuous flux. Networks and Heterogeneous Media, 2010, 5(3): 461-485. doi: 10.3934/nhm.2010.5.461
    [10] Yannick Holle, Michael Herty, Michael Westdickenberg . New coupling conditions for isentropic flow on networks. Networks and Heterogeneous Media, 2020, 15(4): 605-631. doi: 10.3934/nhm.2020016
  • Population transmission models have been helpful in studying the spread of HIV. They assess changes made at the population level for different intervention strategies. To further understand how individual changes affect the population as a whole, game-theoretical models are used to quantify the decision-making process. Investigating multiplayer nonlinear games that model HIV transmission represents a unique approach in epidemiological research. We present here 2-player and multiplayer noncooperative games where players are defined by HIV status and age and may engage in casual (sexual) encounters. The games are modelled as generalized Nash games with shared constraints, which is completely novel in the context of our applied problem. Each player's HIV status is known to potential partners, and players have personal preferences ranked via utility values of unprotected and protected sex outcomes. We model a player's strategy as their probability of being engaged in a casual unprotected sex encounter (USE), which may lead to HIV transmission; however, we do not incorporate a transmission model here. We study the sensitivity of Nash strategies with respect to varying preference rankings, and the impact of a prophylactic vaccine introduced in players of youngest age groups. We also study the effect of these changes on the overall increase in infection level, as well as the effects that a potential prophylactic treatment may have on age-stratified groups of players. We conclude that the biggest impacts on increasing the infection levels in the overall population are given by the variation in the utilities assigned to individuals for unprotected sex with others of opposite HIV status, while the introduction of a prophylactic vaccine in youngest age group (15-20 yr olds) slows down the increase in HIV infection.


    The Keller-Segel equations have been widely used to model chemotaxis. These equations and in particular, the properties of their solutions have been intensively investigated, see for example [9,26,32,36]. For adapted and expanded versions of the Keller-Segel model we refer to [7,13,14,29,32,33,12]. Moreover, in recent literature improved flux-limited Keller-Segel models, taking into account the finite speed of propagation, have been developed in [1]. For surveys and extended reference lists, see, for example [2,3].

    Our starting point is the classical kinetic chemotaxis equation [11]. Scaling it with the so called diffusive scaling leads to the Keller-Segel equation, see [11]. In general, the derivation of Keller-Segel type models, including flux-limited diffusion models and Fokker-Planck type models, from underlying kinetic or microscopic models is discussed for example in [2,12,11]. In particular, using moment closure approaches one may obtain macroscopic equations intermediate between kinetic and Keller-Segel equations, see the above mentioned references or [19,28,6]. Linear moment closure models called the P1-model (in the radiative transfer literature) or the Cattaneo equations are in many situations resonable approximations of the kinetic equation [6]. However, a major drawback of these approximations is that they do not guarantee the positivity of the cell density. Thus, in this paper we consider a non-linear closure method to deal with this problem. Additionally, we also introduce half-moment non-linear closures for the kinetic equations and obtain associated hydrodynamic macroscopic equations similar to [17].

    To model chemotaxis on a network, the crucial point is to define suitable coupling conditions. In previous work [5,6] coupling conditions for the scalar Keller-Segel equation and kinetic equation have been discussed. An analytical investigation can be found in [10]. These investigations deal with the case of a network problem for a parabolic system of equations. Conservation of mass and continuity of density and chemo-attractant at the nodes are used as coupling conditions. Moreover, in [8,25] the hyperbolic-parabolic Cattaneo equations on a graph are studied. For general work on kinetic equations on networks, for example for traffic flow equations, we refer to [18] or [27].

    In the present work we will derive the coupling conditions for nonlinear half-moment models from the coupling condition for the kinetic model. These coupling conditions guarantee on the one hand the conservation of mass through nodes, and on the other hand, they satisfy a positivity condition for the density. In the present case, as in [8,25], the hyperbolic equations require another approach to the coupling conditions as in the Keller-Segel context. However, in the diffusive limit, when the scaling parameter goes to 0, all coupling conditions converge to those of the Keller-Segel model, i.e. the conservation of mass through nodes and a continuity condition, compare for example [5].

    Numerical methods for these nonlinear models have to deal with issues like the numerical moment realizability and (as in the linear case) with a large range of the scaling parameter describing different regimes. Thus, it is desirable to use an asymptotic preserving well-balanced (APWB) scheme. These numerical methods work uniformly with respect to the scaling parameter and one obtains in particular a scheme suitable for computations near the diffusive limit. We refer to [22] for the basic idea of such schemes. Combining this with a numerical scheme to solve the coupling problem we obtain results for the full network problem.

    The paper is organized as follows. In section 2 we discuss the kinetic model and its nonlinear hydrodynamic full-and half-moment approximations which are derived using non-linear closure function. Section 3 considers the coupling conditions for the kinetic equations and derived conditions for the macroscopic models. Section 4 and 5 contain details about the numerical schemes on an interval and the numerical treatment of the coupling conditions. Finally, the numerical results for tripod and more general networks are shown in Section 6.

    Denote with f=f(x,t,v) the cell distribution at position xR, moving with velocity vV[1, 1] at the time t[0, T], ρ(x,t)=Vf(x,t,v)dv is the macroscopic density of cells, m=m(x,t) is the density of the chemoattractant. Let αR+ denote the rate of turning into the favourable direction given by the chemoattractant, let λR+ be the turning rate of the cells and D,γρ,γm R+ be the diffusivity, production and decay rate of the chemoattractant, respectively. As in [6], we consider the kinetic equation

    {tf+vxf=λ(fρ2)+12αv¯xmρtmD(xxm)=γρργmm . (1)

    Here, we have defined

    ¯xm=xm1+|xm|2.

    We assume

    λα (2)

    in order to guarantee positivity of the turning kernel

    λ+αv¯xm0

    for all v, compare [11]. We introduce a parameter ϵ[0,) and scale the kinetic equation (1) in the following way

    {tf+1ϵvxf=λϵ2(fρ2)+12ϵαv¯xmρtmD(xx)m=γρργmm. (3)

    Condition (2) turns into λϵα or

    ϵλα. (4)

    This equation will be used as the starting point for the macroscopic models considered in the next section.

    The kinetic equation (3) can be transformed using the even-and odd-parities [30]

    r(x,t,v)=12(f(x,t,v)+f(x,t,v)),j(x,t,v)=12ϵ(f(x,t,v)f(x,t,v)) (5)

    for positive velocities (v0). This leads to the system

    tr+vxj=λϵ2(rρ2) (6)
    tj+1ϵ2vxr=1ϵ2(λj12αv¯xmρ) . (7)

    When ϵ0, the limit of the kinetic equation (3) is the Keller-Segel equation with flux-limited chemoattractant.

    The procedure is shortly revisited. Start with the scaled version of the kinetic equation (6). Integrating the equation for r over v0 we obtain

    tρ+2x10vjdv=0. (8)

    Moreover we observe r=ρ2+O(ϵ2) and j=α¯xm2λvρvλxr+O(ϵ2). Plugging the expansions into (8) gives the diffusion limit of the kinetic equation, the modified or flux-limited Keller-Segel equation

    {tρx(13λxρα3λ¯xmρ)=0tmDxxm=γρργmm. (9)

    In addition, from the scaled version for kinetic chemotaxis equation (3), using moment closure approaches one may obtain macroscopic, hyperbolic equations intermediate between kinetic and Keller-Segel equations, see [19,29]. There are two kinds of moment models for chemotaxis, full-moment models and half-moment models, that are corresponding to full-and half-moment closure methods respectively.

    Consider the following averaged quantities

    ρ(x,t)=11f(x,t,v)dv ,q(x,t)=11vf(x,t,v)dv ,

    where ρ represents the density and q the flow of the cells. The fact that f(x,t,v)0, it leads to the realizability conditions for full-moment models, see [16,37]:

    ρ0, P0, q2ρ2Pρ1.  (10)

    This leads to the equations

    {tρ+1ϵxq=0tq+1ϵxP=1ϵ2(λqϵα3¯xmρ)tmDxxm=γρργmm . (11)

    In order to close the model (11), a moment closure is specified by using a suitable approximation of f in terms of a function of ρ and q to obtain expressions for P. A linear closure with f(x,t,v)=12ρ(x,t)+32vq(x,t) leads to the P1 model for chemotaxis as in [8,25,6]. However this approach does not guarantee the realizability condition (10) and the cell density f can be negative.

    Therefore, one may assume that the cell density f(x,t,v) is approximated by an exponential function

    f(x,t,v)=aexp(vb)  (12)

    with aR+ and bR. This closure function is chosen using an entropy minimization principle, see for example [23]. This leads to

    P=Vv2fdv=ρh(u),h(u)=(12bu),u:=qρ=coth(b)1b ,

    where

    limu0h(u)=13, u2h(u)1 

    guaranteeing the realizability condition. The function h(u) is called the Eddington factor. Without refering to an underlying closure function the function h(u) can also be chosen directly such that the realizability condition is fulfilled. In [35,23] models with different forms for this function have been proposed, for example, the Kershaw closure h(u)=13+23u2 or the Levermore-Lorentz Eddington factor h(u)=13+2u22+43u2. The resulting full-moment macroscopic model for chemotaxis derived by the nonlinear closure (12) reads

    {tρ+1ϵxq=0tq+1ϵx(ρh(qρ))=1ϵ(λ1ϵqα3¯xmρ)tmDxxm=γρργmm . (13)

    Scaling and denoting qϵ:=qϵ it can be rewritten as

    {tρ+xqϵ=0tqϵ+1ϵ2x(ρh(ϵqϵρ))=1ϵ2(λqϵα3¯xmρ)tmDxxm=γρργmm . (14)

    Under the realizability condition, the model (14) is hyperbolic and converges to the Keller-Segel equations when ϵ0.

    We will construct a macroscopic, half-moment model for (3) by applying a half-moment nonlinear closure, which was introduced in [17], see [6] for a linear version. Consider the following quantities

    ρ=01f(v)dv ,ρ+=10f(v)dv ,q=01vf(v)dv ,q+=10vf(v)dv ,P=01v2f(v)dv ,P+=10v2f(v)dv . (15)

    Using f(x,t,v)0, one obtains the conditions ρ±0, ±q±0, P±0 and

    (q±ρ±)2P±ρ±±q±ρ±1 . (16)

    These conditions again guarantee the existence of a non-negative density for these equations [16,37].

    Integrating the kinetic equation (3) leads to the half-moment model

    {ϵtρ±+xq±=1ϵλ(ρ±ρ++ρ2)±14α¯xm(ρ++ρ)ϵtq±+xP±=1ϵλ(q±ρ++ρ4)+16α¯xm(ρ++ρ)tmDxxm=γρργmm . (17)

    The nonlinear closure for the half-moment models is derived similarly as for the full-moment case. One uses the ansatz

    f(v<0)=aexp(vb) ,f(v>0)=a+exp(vb+) . (18)

    Plugging (18) in (15) leads to

    P=ρh(u) ,P+=ρ+h+(u+) ,h(u)=u(12b)1b ,h+(u+)=u+(12b+)+1b+ ,u:=qρ=1exp(b)11(b) ,u+:=q+ρ+=exp(b+)exp(b+)11(b+) ,

    such that

    limu±±12h±(u±)=13 ,(u±)2h±(u±)±u± . (19)

    Condition (19) guarantees the realizability condition (16). We note that an explicit maximum-entropy Eddington factor approximating the above closure is obtained by using the Kershaw closure [34]

    h(u)=23(u)213u ,h+(u+)=23(u+)2+13u+ . (20)

    Finally, the half-moment model for chemotaxis reads

    {tρ±+1ϵxq±=1ϵ2[λ(ρ±ρ++ρ2)ϵα4¯xm(ρ++ρ)]tq±+1ϵx(ρ±h±(q±ρ±))=1ϵ2[λ(q±ρ++ρ4)ϵα6¯xm(ρ++ρ)]tmDxxm=γρργmm . (21)

    Moreover, by introducing new quantities

    ρ=ρ++ρ ,ˆρ=ρ+ρϵ ,qϵ=q++qϵ ,ˆq=q+q , (22)

    the half-moment model can be rewritten as

    {tρ+xqϵ=0tqϵ+1ϵ2xP(ρ,qϵ,ˆρ,ˆq)=1ϵ2(λqϵα3¯xmρ)tˆρ+1ϵ2xˆq=1ϵ2(λˆρ12α¯xmρ)tˆq+xˆP(ρ,qϵ,ˆρ,ˆq)=1ϵ2λ(ˆqρ2)tmDxxm=γρργmm  (23)

    with

    P(ρ,qϵ,ˆρ,ˆq):=P+(q+ρ+)+P(qρ)=ρ+ϵˆρ2h+(ϵqϵ+ˆqρ+ϵˆρ)+ρϵˆρ2h(ϵqϵˆqρϵˆρ),ˆP(ρ,qϵ,ˆρ,ˆq):=1ϵP+(q+ρ+)1ϵP(qρ)=ρ+ϵˆρ2ϵh+(ϵqϵ+ˆqρ+ϵˆρ)ρϵˆρ2ϵh(ϵqϵˆqρϵˆρ). (24)

    The half-moment model (23) is hyperbolic and has again the Keller-Segel equations as macroscopic diffusive limit when ϵ goes to 0.

    We concentrate in the following on the discretization of the equations for the movements of the cell density. The equation for the chemoattractant m is always solved with a forward central difference scheme as in [5]. We construct an APWB numerical scheme for the nonlinear chemotaxis model following closely the strategies developed in [21] for the linear Cattaneo model and [4,15] for radiation hydrodynamics. Moreover, we refer to [24,20,31] for further details on AP and well-balanced schemes.

    First, let us recall the framework of the APWB scheme in [24,20,21] for the Cattaneo model

    {tρ+xqϵ=0tqϵ+aϵ2xρ=1ϵ2(λqϵc(xm)ρ) (25)

    with c(x)=α3x1+x2, a,λ,α>0, 0<ϵλα. The limit problem as ϵ0 is

    tρ+x(c(xm)λρaλxρ)=0.

    By introducing the diagonal variables U=12ρ+ϵ2aqϵ,V=12ρϵ2aqϵ, the system (25) is rewritten as

    tU+aϵxU=(c(xm)2aϵλ2ϵ2)U+(c(xm)2aϵ+λ2ϵ2)VtVaϵxV=(c(xm)2aϵλ2ϵ2)U(c(xm)2aϵ+λ2ϵ2)V. (26)

    Using

    {U(Δx)=V(Δx)+A+B1+B(U(0)V(Δx))+A+B21+BV(Δx)V(0)=U(0)+A+B1+B(V(Δx)U(0))A+B21+BV(Δx)A=(1λΔxϵa)exp(Δxac(xm))B=(1+λΔxϵa)exp(Δxac(xm)) (27)

    the well-balanced Godunov scheme for (26) is given as

    {Un+1i=Unil(Un+1iUi12)=Unil(Un+1iVn+1i)lAni12+Bni121+Bni12(VniUni1)+lAni12+Bni1221+Bni12VniVn+1i=Vni+l(Vi+12Vn+1i)=Vni+l(Un+1iVn+1i)+lAni+12+Bni+121+Bni+12(Vni+1Uni)lAni+12+Bni+1221+Bni+12Vni+1 (28)

    with l:=aϵΔtΔx and xmni+12:=mni+1mniΔx. The APWB numerical solution (28) can be rewritten explicitly as

    {Un+1i=2l22l+111+Bni+12Vni+1+(l+12l+1l22l+1Ani+12+Bni+121+Bni+12)Uni+(l2l+12l(l+1)2l+111+Bni12)Vni+l(l+1)2l+1Ani12+Bni121+Bni12Uni1Vn+1i=l(l+1)2l+121+Bni+12Vni+1+(l2l+1l(l+1)2l+1Ani+12+Bni+121+Bni+12)Uni+(l+12l+12l22l+111+Bni12)Vni+l22l+1Ani12+Bni121+Bni12Uni1.  (29)

    It is easy to check that for small enough values of ϵaΔtΔx1, the APWB scheme (29) is positivity preserving for U and V, if the following parabolic CFL restriction holds

    ΔtλΔx24(a+Δx).

    Furthermore, as mentioned in [21], (28) converges to a centered discretization of the Keller-Segel equation, when ϵ0.

    Based on the above considerations and the schemes developed in the above cited references we consider the following relaxation system for (14)

    {tρ+xz=0tz+aϵ2xρ+λϵ2zα3ϵ2¯xmρ=1β(qϵz)tqϵ+1ϵ2xw+λϵ2qα3ϵ2¯xm(1h(ϵqϵρ)w)ρ=0tw+axq=1β(ρh(ϵqϵρ)w) . (30)

    For β0 the system converges to (14). Note that for any value of β the system (30) converges to the Keller-Segel equation with a diffusion coefficient equal to a.

    System (30) is split into two into separate sub systems. The right hand side treated with a simple implicit numerical scheme and we obtain for β=0

    ρ(1)=ρn,q(1)ϵ=z(1)=qnϵ,w(1)=ρnh(ϵqnϵρn) , (31)

    which is used as the initial condition for the system given by the left hand side. By introducing the new variables

    U=12ρ+ϵ2az ,V=12ρϵ2az ,ˉU=ϵ2qϵ+12aw ,ˉV=ϵ2qϵ+12aw , (32)

    the left hand side of system (30) can be rewritten as

    {tU+aϵxU=(c12aϵλ2ϵ2)U+(c12aϵ+λ2ϵ2)VtVaϵxV=(c12aϵλ2ϵ2)U(c12aϵ+λ2ϵ2)Vc1(xm)=α3¯xm{tˉU+aϵxˉU=(c22aϵλ2ϵ2)ˉU+(c22aϵ+λ2ϵ2)ˉVtˉVaϵxˉV=(c22aϵλ2ϵ2)ˉU(c22aϵ+λ2ϵ2)ˉVc2(ρ,q,xm)=ah(ϵqϵρ)c1 . (33)

    Now, we can apply the APWB schemes presented in 3.1 for each subsystem in (33). Remark that for a:=max|h(ϵqϵρ)| the quantities in (32) are nonnegative and the scheme converges formally to the Keller-Segel equations as ϵ0.

    Remark 1. We note that under similar conditions as in the linear case, it can be shown, that the scheme preserves positivity in U and V. However, this does not mean that we can prove that the density ρ remains positive. In the linear case, this is not an issue to be investigated, since the density ρ in the Cattaneo model is not necessarily positive. In contrast, in the nonlinear case, it is an important open issue to be investigated.

    Following the same procedure as for the full-moment model in the previous section, we consider the following relaxation system for the half-moment model (23)

    {tρ+xz=0tz+aϵ2xρ+λϵ2zα3ϵ2¯xmρ=1β(qϵz)tqϵ+1ϵ2xw+λϵ2qϵα3ϵ2¯xm(2wϵ(h+h)ˆρh++h)ρ=0tw+axqϵ=1β(Pw) {tˆρ+1ϵ2xˆz+λϵ2ˆρα2ϵ2¯xmρ=0tˆz+axˆρ+λϵ2ˆz=1β(ˆqˆz)tˆq+xˆw+λϵ2(ˆqρ2)=0tˆw+aϵ2xˆq=1β(ˆPˆw) (34)

    with h±:=h±(ϵqϵ±ˆqρ±ϵˆρ) and P=P(ρ,qϵ,ˆρ,ˆq),ˆP=P(ρ,qϵ,ˆρ,ˆq) defined in section 2.3. To obtain the correct asymptotic behaviour using a splitting scheme, we introduce a parameter σ<λ, and rewrite the system (34) as

    {tρ+xz=0tz+aϵ2xρ+λϵ2zα3ϵ2¯xmρ=1β(qϵz)tqϵ+1ϵ2xw+λσϵ2qϵ2wα¯xm3ϵ2(h++h)=σϵ2qϵα¯xm(h+h)3ϵ(h++h)ˆρtw+axqϵ=1β(Pw) {tˆρ+1ϵ2xˆz=λϵ2ˆρ+α2ϵ2¯xmρtˆz+axˆρ+λϵ2ˆz=1β(ˆqˆz)tˆq+xˆw+σϵˆq=λσϵ2ˆq+λ2ϵ2ρtˆw+aϵ2xˆq=1β(ˆPˆw) .  (35)

    Using again an implicit scheme for the right hand side of (35) and setting β=0, one obtains for the first step of the splitting scheme

    {ρ(1)=ρnq(1)ϵ=ϵ2ϵ2+σΔtqnϵϵΔtα3(ϵ2+σΔt)¯xmn(h+,nh,nh+,n+h,n)ˆρ(1)z(1)=q(1)ϵw(1)=P(ρ(1),q(1)ϵ,ˆρ(1),ˆq(1)) {ˆρ(1)=ϵ2ϵ2+λΔtˆρn+Δtα2(ϵ2+λΔt)¯xmnρnˆq(1)=ϵ2ϵ2+(λσ)Δtˆqn+Δtλ2(ϵ2+(λσ)Δt)ρnˆz(1)=ˆq(1)ˆw(1)=ˆP(ρ(1),q(1)ϵ,ˆρ(1),ˆq(1)) .  (36)

    This relaxation scheme ensures that system (35) has the correct asymptotic limit as ϵ0. Thus (35) converges to the Keller-Segel equation (9). Moreover, the solution in the relaxation part (36) preserves the realizability condition i.e., |ϵq(1)ϵ|ρ(1).

    The transport part in (35) includes four independent subsystems with (36) as the initial data. Introducing the characteristic variables (U1,V1), (ˉU1,ˉV1), (U2,V2), (ˉU2,ˉV2),

    U1=12ρ+ϵ2az ,ˉU1=ϵ2qϵ+12aw ,U2=ϵ2ˆρ+12aˆz ,ˉU2=12ˆq+ϵ2aˆw ,V1=12ρϵ2az ,ˉV1=ϵ2qϵ+12aw ,V2=ϵ2ˆρ12aˆz ,ˉV2=12ˆq+ϵ2aˆw (37)

    we use the APWB scheme for the left hand side of (35). With a:=max(|wρ|) we have again the positivity for the characteristic variables (37). Moreover, the half-moment model (34) converges formally to the Keller-Segel equation (9) when β0 and ϵ0.

    In this section we state a hierarchy of coupling conditions for the kinetic and the nonlinear half-moment model proposed in the previous section. In [6] coupling conditions have been derived for the kinetic model and linear macroscopic approximations. These conditions can be used as well for the nonlinear half moment model discussed in the previous section. Finding conditions for the nonlinear full moment model is much harder since the signs of the eigen values can change with the states of the system. First the coupling conditions for the kinetic model are recalled to derive thereof the nonlinear half moment coupling conditions.

    Consider the kinetic equation

    tf+1ϵvxf=λϵ2(fρ2)+12ϵαv¯xmρ . (38)

    As stated in [6], the coupling conditions for (38) should assign a value to all f(v) for v[0,1]. A possible choice is

    f+=Af,
    Figure 1.  Sketch of a tripod network.

    where f+i=fi(v) and fi=fi(v) for v[0,1] and i=1,,N. In order to conserve the total mass in the system the matrix ARN×N has to fulfill

    Ni=1ai,j=1        j=1,,N .

    Moreover, positivity requires ai,j0, i,j=1,,N. If all edges are treated equally, the coupling conditions for three edges are

    [f+1f+2f+3]=[01/21/21/201/21/21/20][f1f2f3] . (39)

    This can be generalized to the cases N>3 by choosing ai,j=1N1 ij and ai,i=0.

    Remark 2. Proving existence, uniqueness or stability of the kinetic equation with the above coupling conditions on the graph are interesting further issues, which would deserve a deeper analysis.

    Since the kinetic coupling conditions (39) are linear and independent of v we obtain directly by integrating the coupling conditions for the half moment model (21)

    {[ρ+1ρ+2ρ+3]=[01/21/21/201/21/21/20][ρ1ρ2ρ3][q+1q+2q+3]=[01/21/21/201/21/21/20][q1q2q3]. (40)

    Recall that all properties of (39) are inherited by the above coupling conditions. For example the total mass in the system is conserved since

    3i=1qϵ,i=1ϵ3i=1(q+i+qi)=0 .

    In particular, the coupling conditions for the half-moment model yield the positivity of the densities, i.e. ρ+i0 i=1,,N if ρj0 j=1,,N.

    Here only the case of a junction coupling three edges will be treated, but all procedures easily extend to arbitrary junctions. Consider the left hand side of (35) which includes four independent subsystems of eight variables

    {tρ+xz=0tz+aϵ2xρ+λϵ2zα3ϵ2¯xmρ=0 {tqϵ+1ϵ2xw+λσϵ2qϵα3ϵ2¯xm2(h++h)w=0tw+axqϵ=0{tˆρ+1ϵ2xˆz=0tˆz+axˆρ+λϵ2ˆz=0{tˆq+xˆw+σϵ2ˆq=0tˆw+aϵ2xˆq=0 .  (41)

    For the additional numerical quantities we need further coupling conditions. Note that for i=1,...,N from the relaxation part of (35) we have: zi=qi,ˆzi=ˆqi and wi=Pi,ˆwi=ˆPi: Thus, conditions for the states z and ˆz are found from those prescribed for q. To obtain conditions for w and ˆw we have to find conditions for P±i. These are easily derived from the kinetic coupling conditions (39). Multiply (39) by v2 and integrate over positive velocities to obtain

    [P+1P+2P+3]=[01/21/21/201/21/21/20][P1P2P3]. (42)

    Inserting (22) into (40), we express the coupling conditions in the variables of (41)

    [ρ1+ϵˆρ1ρ2+ϵˆρ2ρ3+ϵˆρ3]=[01/21/21/201/21/21/20][ρ1ϵˆρ1ρ2ϵˆρ2ρ3ϵˆρ3][ϵq1+ˆq1ϵq2+ˆq2ϵq3+ˆq3]=[01/21/21/201/21/21/20][ϵq1ˆq1ϵq2ˆq2ϵq3ˆq3] (43)

    and with (42) we obtain

    [P1+ϵˆP1P2+ϵˆP2P3+ϵˆP3]=[01/21/21/201/21/21/20][P1ϵˆP1P2ϵˆP2P3ϵˆP3] . (44)

    This leads to

    [ϵz1+ˆz1ϵz2+ˆz2ϵz3+ˆz3]=[01/21/21/201/21/21/20][ϵz1ˆz1ϵz2ˆz2ϵz3ˆz3][w1+ϵˆw1w2+ϵˆw2w3+ϵˆw3]=[01/21/21/201/21/21/20][w1ϵˆw1w2ϵˆw2w3ϵˆw3] . (45)

    Since on each edge there are four waves with positive speed, we need exactly twelve equations (43), (45) for three edges, as provided.

    Remark 3. We have not investigated the positivity preserving property of the scheme for the density ρ, see Remark 1 for the case of a single line. It would be even more interesting to investigate this property on the graph taking into account the discretization of the coupling conditions.

    In this section we investigate the proposed models in several numerical test cases for different values of ϵ. For simplicity we use for the nonlinear half-moment model the explicit Kershaw closure, i.e. h± given by (20). In all the considered examples we will use the following values for the parameters λ=α=1, D=1, γρ=1 and γm=0.1. The spatial resolution is Δx=0.02 and the time step is chosen according to the CFL condition. At end points where no coupling conditions are imposed zero Neumann boundary conditions are applied. For the kinetic model we discretize the velocity space V=[1,1] with Nv=50 cells.

    Test 1. We consider the interval [0,2] and Δx=0.005. As initial conditions for the kinetic equation we use

    f(x,v,0)=F(v)ρ(x,0)=12ρ(x,0) ,

    with

    ρ(x,0)={1 if 0x10 if 1x2 .

    The remaining initial values for the hydrodynamic equations can be derived from the kinetic initial condition. At t=0 no chemoattractant m(x,0)=0 is present.

    In the Figures 2 and 3 the densities at t=0.2 for the values ϵ=1, 0.5, 0.1 and ϵ= 106 are shown. For the kinetic model we observe that diffusion depends strongly on the value of ϵ as expected. The closest approximation to the kinetic equations are the half moment models, in particular, the nonlinear half-moment model. For the half moment P1 model for ϵ=1 and ϵ=0.5 one can clearly observe the four waves generated by the advective part. In the full moment P1 model only two waves are used to approximate the kinetic solution. The flux-limited Keller-Segel equation is evolving too fast for large ϵ. In general, the half-moment models are clearly superior to the full moment models and the nonlinear models provide better approximations than the linear models. For small ϵ all models converge to the solution of the Keller-Segel equation.

    Figure 2.  Numerical solutions of the four models on an interval at time t=0.2 with ϵ=1 (top) and ϵ=0.5 (bottom).
    Figure 3.  Numerical solutions of the four models on an interval at time t=0.2 with ϵ=0.1 (top) and ϵ=106 (bottom).

    Test 2. In the second test we investigate the positivity preserving of the nonlinear moment models which is not guaranteed in the case of the linear models. We use non-equilibrium initial conditions

    f(x,0,v)={aexp(vb)if1v00if0v1 (46)

    for the kinetic equation (3). a and b are chosen such that the initial data of the half-moment models are given by

    ρ(x,0)=01fdv=ρ0 ,q(x,0)=01vfdv=0.9×ρ0 ,ρ+(x,0)=10fdv=0 ,q+(x,0)=10vfdv=0 , (47)

    with

    ρ0(x)={1 if 0x1106 if 1x2,m(x,0)=0, x[0,2] .

    The initial conditions for the full-moment models are

    ρ(x,0):=ρ++ρ=ρ0,qϵ(x,0):=q++qϵ=0.9ϵρ0 .

    The solutions are plotted in Figure 4, for Δx=0.02 at t=0.5. Even though the initial conditions satisfy the realizability conditions ρ±|q±| (for half-moment models) and ρ|q| (for full-moment models), only the nonlinear full-and half-moment models preserve the positivity of the cell density.

    Figure 4.  Linear models (P1) with negative full or half range densities and positivity preserving nonlinear models (M1).

    This test case focuses on the coupling conditions for the non-linear half moment model proposed in section 4. We study a junction connecting three outgoing edges, as depicted in Figure 1. The results of the nonlinear half-moment model with the coupling conditions derived in section 4 are compared with the results of the kinetic, linear half and full moment and Keller-Segel model with coupling conditions described in [6]. On each of the edges we consider the interval [0,1]. As initial conditions we choose

    ρ1(x,0)=1 ,ρ2(x,0)=4 ,ρ3(x,0)=3 ,

    which is consistent with the following values for the kinetic equation

    fi(x,v,0)=F(v)ρi(x,0)=12ρi(x,0)i=1,2,3 .

    All other quantities are initially zero.

    In Figures 5-8 the densities at t=0.2 for different values of ϵ are shown. The results are similar to those on a single interval. The nonlinear half moment model is again the best approximation to the kinetic result. As the value of ϵ decreases all solutions approach the solution of the Keller-Segel equation.

    Figure 5.  Numerical solutions on a tripod network at time t=0.2, Δx=0.02, ϵ=1.
    Figure 6.  Numerical solutions on a tripod network at time t=0.2, Δx=0.02, ϵ=0.5.
    Figure 7.  Numerical solutions on a tripod network at time t=0.2, Δx=0.02, ϵ=0.1.
    Figure 8.  Numerical solutions on a tripod network at time t=0.2, Δx=0.02, ϵ=106.

    In this last test case we consider a larger network of 31 edges and 23 nodes as shown in Figure 9. The length of the short edges is 0.5, for the longer ones we have 1 and 2 respectively. Note that this network does not only contain three way junctions, but also nodes connecting up to five edges. At the open ends of the network Dirichlet boundary conditions are imposed. For the kinetic model we prescribe f(0,v,t)=12 for v>0. The boundary values for the other models can be derived thereof. As initial conditions all values are set to zero except the density in the edges at the outer boundaries, which is set to 1. For the numerical computation we used 15, 30 and 42 cells for the spatial discretization of the edges, respectively. The velocity space in the kinetic model is resolved with 30 points. The scaling parameter is chosen as ϵ=1.

    In Figure 9 the density at time t=5 for all four models is shown. There is an inflow from both sides of the network. As observed in the previous tests, the states of the Keller-Segel model propagate faster, such that the network is filled at an earlier time. The solution of the kinetic model and the one of the half moment models almost coincide.

    Figure 9.  Comparison of the numerical solutions on a larger network at t=5.

    In Figure 10 the evolution of the total mass in the network up to T=30 is shown. As before, we observe that the Keller-Segel model fills the network much faster, than the other three. The values of the half moment and of the kinetic model almost coincide. Concerning the computation times, the least expensive model is the full moment P1 model with approximately 15% of the computation time of the kinetic model. The Keller-Segel model needs 18%, the half moment P1 model needs 22% and the half moment M1 model needs 23% of the kinetic computation time. The overall CPU time for the kinetic computation with 30 spatial and 30 velocity cells per edge has been approximately 210 min on an Intel Core i5-3230M with 2.6 GHz.

    Figure 10.  Total mass over time in the large network.

    In this paper we developed nonlinear half-moment models for chemotaxis models which approximate very efficiently the solution of a network problem governed by kinetic equations. An APWB scheme is investigated for these nonlinear moment models. We derive from coupling conditions for the kinetic models corresponding conditions for the macroscopic models. Note that this procedure also might be applied to other kinetic models. In the numerical tests we investigated the dynamics for different coupling conditions and for varying values of ϵ. A simulation on a larger network showed the applicability to more complicated structures and differences in behavior of macroscopic and kinetic models. Finally, we remark that there are several open issues. Most important would be analytical investigations of the kinetic network problem and a thorough investigation of the positivity preserving property of the APWB scheme for the density ρ in the nonlinear half moment model on the graph.

    [1] [ K. J. Arrow,G. Debreu, Existence of an equilibrium for a competitive economy, Econometrica, 22 (1954): 265-290.
    [2] [ J. P. Aubin and A. Cellina, Differential Inclusions: Set-Valued Maps and Viability Theory, Springer-Verlag New York, 1984.
    [3] [ R. M. Axelrod, The complexity of cooperation: Agent-based models of competition and collaboration, Princeton University Press, 1997.
    [4] [ T. Basar and G. J. Olsder, Dynamic Noncooperative Game Theory, SIAM, 1995.
    [5] [ A. Bensoussan, Points de Nash dans le cas de fonctionnelles quadratiques et jeux differentiels lineaires a N personnes, SIAM J. Control, 12 (1974): 460-499.
    [6] [ Center for Disease Control and Prevention, Diagnoses of HIV Infection in the United States and Dependent Areas, HIV Surveillance Report, 2011.
    [7] [ B. Coburn,S. Blower, A major HIV risk factor for young men who have sex with men is sex with older partners, J. Acquir. Immune Defic. Syndr., 54 (2010): 113-114.
    [8] [ M. G. Cojocaru,L. B. Jonker, Existence of solutions to projected differential equations on Hilbert spaces, Proceedings of the American Mathematical Society, 132 (2004): 183-193.
    [9] [ Cojocaru, M. -G., Wild, E., On Describing the Solution Sets of Generalized Nash Games with Shared Constraints, submitted to: Optimization and Engineering, 2016.
    [10] [ M. G. Cojocaru,C. T. Bauch,M. D. Johnston, Dynamics of vaccination strategies via projected dynamical systems, Bulletin of mathematical biology, 69 (2007): 1453-1476.
    [11] [ M. G. Cojocaru, Dynamic equilibria of group vaccination strategies in a heterogeneous population, Journal of Global Optimization, 40 (2008): 51-63.
    [12] [ J. P. Dodds,A. Nardone,D. E. Mercey,A. M. Johnson, Increase in high risk sexual behaviour among homosexual men, London 1996-8: Cross sectional, questionnaire study, BMI, 320 (2000): 1510-1511.
    [13] [ F. Facchinei,A. Fischer,V. Piccialli, On generalized Nash games and variational inequalities, Operations Research Letters, 35 (2007): 159-164.
    [14] [ F. Fachinei,C. Kanzow, Generalized Nash Equilibrium Problems, Ann Oper Res, 175 (2010): 177-211.
    [15] [ J. W. Friedman, Game theory with applications to economics, Oxford University Press New York, 1990.
    [16] [ D. Gabay,H. Moulin, On the uniqueness and stability of Nash-equilibria in noncooperative games, Applied Stochastic Control in Econometrics and Management Science, 130 (1980): 271-293.
    [17] [ R. H. Gray,M. J. Wawer,R. Brookmeyer,N. K. Sewankambo,D. Serwadda,F. Wabwire-Mangen,T. Lutalo,X. Li,T. VanCott,T. C. Quinn, Probability of HIV-1 transmission per coital act in monogamous, heterosexual, HIV-1-discordant couples in Rakai, Uganda, The Lancet., 357 (2001): 1149-1153.
    [18] [ D. Greenhalgh,F. Lewis, Stochastic models for the spread of HIV amongst intravenous drug users, Stochastic Models, 17 (2001): 491-512.
    [19] [ T. B. Hallett,S. Gregson,O. Mugurungi,E. Gonese,G. P. Garnett, Assessing evidence for behaviour change affecting the course of HIV epidemics: a new mathematical modelling approach and application to data from Zimbabwe, Epidemics, 1 (2009): 108-117.
    [20] [ Y. H. Hsieh,C. H. Chen, Modelling the social dynamics of a sex industry: Its implications for spread of HIV/AIDS, Bulletin of Mathematical Biology, 66 (2004): 143-166.
    [21] [ J. Hofbauer and K. Sigmund, Evolutionary Games and Population Dynamics, Cambridge University Press, 1998.
    [22] [ D. M. Huebner,M. A. Gerend, The relation between beliefs about drug treatments for HIV and sexual risk behavior in gay and bisexual men, Annals of Behavioral Medicine, 4 (2001): 304-312.
    [23] [ T. Ichiishi, Game Theory for Economic Analysis, New York: Academic Press, New York, 1983.
    [24] [ D. Kinderlehrer,D. Stampacchia, null, An Introduction to Variational Inequalities and their Application, Academic Press, New York, 1980.
    [25] [ K. Nabetani, P. Tseng and M. Fukushima, Parametrized Variational Inequality Approaches to Generalized Nash Equilibrium Problems with Shared Constraints, (Technical Report 2008), Department of Applied Mathematics and Physics, Graduate School of Informatics, Kyoto University, Japan.
    [26] [ J. F. Nash, Equilibrium points in n-person games, CMS, 2 (2005): 21-56.
    [27] [ J.-S. Pang,M. Fukushima, Quasi-variational inequalities, generalized Nash equilibria, and multi-leader-follower games, CMS, 2 (2005): 21-56.
    [28] [ R. H. Remien,P. N. Halkitis,A. O'Leary,R. J. Wolitski,C. A. Gómez, Risk perception and sexual risk behaviors among HIV-positive men on antiretroviral therapy, AIDS and Behavior, 9 (2005): 167-176.
    [29] [ J. B. Rosen, Existence and uniqueness of equilibrium points for concave n-person games, Econometrica, 33 (1965): 520-534.
    [30] [ K. D. Schroeder,F. G. Rojas, A game theoretical analysis of sexually transmitted disease epidemics, Rationality and Society, 14 (2002): 353-383.
    [31] [ R. J. Smith,S. M. Blower, Could disease-modifying HIV vaccines cause population-level perversity?, The Lancet Infectious Diseases, 4 (2004): 636-639.
    [32] [ J. M. Stephenson,J. Imrie,M. M. D. Davis,C. Mercer,S. Black,A. J. Copas,G. J. Hart,O. R. Davidson,I. G. Williams, Is use of antiretroviral therapy among homosexual men associated with increased risk of transmission of HIV infection?, Sexually Transmitted Infections, 79 (2003): 7-10.
    [33] [ S. Tully,M. G. Cojocaru,C. T. Bauch, Coevolution of risk perception, sexual behaviour, and HIV transmission in an agent-based model, Journal of theoretical biology, 337 (2013): 125-132.
    [34] [ P. T. Harker, Generalized Nash games and quasi-variational inequalities, European Journal of Operational Research, 54 (1991): 81-94.
    [35] [ S. Tully, M. G. Cojocaru and C. T. Bauch, Sexual behaviour, risk perception, and HIV transmission can respond to HIV antiviral drugs and vaccines through multiple pathways, Scientific Reports, 5(2015).
    [36] [ J. X. Velasco-Hernandez,Y. H. Hsieh, Modelling the effect of treatment and behavioral change in HIV transmission dynamics, Journal of mathematical biology, 32 (1994): 233-249.
    [37] [ J. Von Neumann,O. Morgenstern, Theory of games and economic behavior, Bulletin of American Mathematical Society, 51 (1945): 498-504.
    [38] [ U. S. Census Bureau, Current Population Survey: Annual Social and Economic Supplement, 2012. Available at: https://www.census.gov/topics/population.html, (Accessed: Aug 2013).
    [39] [ University of Western News, HIV vaccine produces no adverse effects in trials, 2013, Available at: http://communications.uwo.ca/western_news/stories/2013/\September/hiv_vaccine_produces_no_adverse_effects_in_trials.html, (Accessed: Sept 2014).
    [40] [ Global Health Observatory, World Health Organization, 2014, Available at: http://www.who.int/gho/hiv/en/, (Accessed: Sept 2014).
    [41] [ Index Mundi, Zimbabwe Age structure, 2007, Available at: http://www.indexmundi.com/zimbabwe/age_structure.html, (Accessed: Aug 2014).
  • This article has been cited by:

    1. Raul Borsche, Axel Klar, Florian Schneider, 2019, Chapter 1, 978-3-030-20296-5, 1, 10.1007/978-3-030-20297-2_1
    2. G. Corbin, A. Hunt, A. Klar, F. Schneider, C. Surulescu, Higher-order models for glioma invasion: From a two-scale description to effective equations for mass density and momentum, 2018, 28, 0218-2025, 1771, 10.1142/S0218202518400055
    3. Alberto Tesei, 2024, Chapter 12, 978-3-031-60772-1, 241, 10.1007/978-3-031-60773-8_12
    4. Sergei A. Avdonin, Alexander S. Mikhaylov, Victor S. Mikhaylov, Abdon E. Choque‐Rivero, Discretization of the Wave Equation on a Metric Graph, 2024, 0170-4214, 10.1002/mma.10630
  • Reader Comments
  • © 2017 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(3941) PDF downloads(552) Cited by(4)

Figures and Tables

Figures(8)  /  Tables(2)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog