Review Special Issues

Monoclonal Antibodies as Treatment Modalities in Head and Neck Cancers

  • Received: 29 August 2015 Accepted: 28 October 2015 Published: 03 November 2015
  • The standard treatments of surgery, radiation, and chemotherapy in head and neck squamous cell carcinomas (HNSCC) causes disturbance to normal surrounding tissues, systemic toxicities and functional problems with eating, speaking, and breathing. With early detection, many of these cancers can be effectively treated, but treatment should also focus on retaining the function of the proximal nerves, tissues and vasculature surrounding the tumor. With current research focused on understanding pathogenesis of these cancers in a molecular level, targeted therapy using monoclonal antibodies (MoAbs), can be modified and directed towards tumor genes, proteins and signal pathways with the potential to reduce unfavorable side effects of current treatments. This review will highlight the current MoAb therapies used in HNSCC, and discuss ongoing research efforts to develop novel treatment agents with potential to improve efficacy, increase overall survival (OS) rates and reduce toxicities.

    Citation: Vivek Radhakrishnan, Mark S. Swanson, Uttam K. Sinha. Monoclonal Antibodies as Treatment Modalities in Head and Neck Cancers[J]. AIMS Medical Science, 2015, 2(4): 347-359. doi: 10.3934/medsci.2015.4.347

    Related Papers:

    [1] Raimund Bürger, Christophe Chalons, Rafael Ordoñez, Luis Miguel Villada . A multiclass Lighthill-Whitham-Richards traffic model with a discontinuous velocity function. Networks and Heterogeneous Media, 2021, 16(2): 187-219. doi: 10.3934/nhm.2021004
    [2] Maya Briani, Emiliano Cristiani . An easy-to-use algorithm for simulating traffic flow on networks: Theoretical study. Networks and Heterogeneous Media, 2014, 9(3): 519-552. doi: 10.3934/nhm.2014.9.519
    [3] Mauro Garavello, Roberto Natalini, Benedetto Piccoli, Andrea Terracina . Conservation laws with discontinuous flux. Networks and Heterogeneous Media, 2007, 2(1): 159-179. doi: 10.3934/nhm.2007.2.159
    [4] 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
    [5] Helge Holden, Nils Henrik Risebro . Follow-the-Leader models can be viewed as a numerical approximation to the Lighthill-Whitham-Richards model for traffic flow. Networks and Heterogeneous Media, 2018, 13(3): 409-421. doi: 10.3934/nhm.2018018
    [6] Adriano Festa, Simone Göttlich, Marion Pfirsching . A model for a network of conveyor belts with discontinuous speed and capacity. Networks and Heterogeneous Media, 2019, 14(2): 389-410. doi: 10.3934/nhm.2019016
    [7] Christophe Chalons, Paola Goatin, Nicolas Seguin . General constrained conservation laws. Application to pedestrian flow modeling. Networks and Heterogeneous Media, 2013, 8(2): 433-463. doi: 10.3934/nhm.2013.8.433
    [8] Raimund Bürger, Stefan Diehl, M. Carmen Martí, Yolanda Vásquez . A difference scheme for a triangular system of conservation laws with discontinuous flux modeling three-phase flows. Networks and Heterogeneous Media, 2023, 18(1): 140-190. doi: 10.3934/nhm.2023006
    [9] Giuseppe Maria Coclite, Lorenzo di Ruvo, Jan Ernest, Siddhartha Mishra . Convergence of vanishing capillarity approximations for scalar conservation laws with discontinuous fluxes. Networks and Heterogeneous Media, 2013, 8(4): 969-984. doi: 10.3934/nhm.2013.8.969
    [10] Wen Shen . Traveling wave profiles for a Follow-the-Leader model for traffic flow with rough road condition. Networks and Heterogeneous Media, 2018, 13(3): 449-478. doi: 10.3934/nhm.2018020
  • The standard treatments of surgery, radiation, and chemotherapy in head and neck squamous cell carcinomas (HNSCC) causes disturbance to normal surrounding tissues, systemic toxicities and functional problems with eating, speaking, and breathing. With early detection, many of these cancers can be effectively treated, but treatment should also focus on retaining the function of the proximal nerves, tissues and vasculature surrounding the tumor. With current research focused on understanding pathogenesis of these cancers in a molecular level, targeted therapy using monoclonal antibodies (MoAbs), can be modified and directed towards tumor genes, proteins and signal pathways with the potential to reduce unfavorable side effects of current treatments. This review will highlight the current MoAb therapies used in HNSCC, and discuss ongoing research efforts to develop novel treatment agents with potential to improve efficacy, increase overall survival (OS) rates and reduce toxicities.



    Traffic flow models based on scalar conservation laws with continuous flux functions are widely used in the literature. For a general presentation of the models, we refer to the books [11,12,23] and the references therein. Extensions to road traffic networks have been also established. We mention in particular the contributions [6,15], where the authors introduce the coupled network problem and show the existence of solutions. Within this article, we are concerned with the special case of scalar conservation laws with discontinuous flux in the unknown that are motivated in the traffic flow theory by the observation of a gap between the free flow and the congested flow regime [4,5,8]. This phenomenon generates an interesting dynamical behavior called zero waves, i.e., waves with infinite (negative) speed but zero wave strength, and has been investigated in recent years either from a theoretical or numerical point of view, see for instance [2,19,20,22,24] or more generally [1,3,7,13].

    To the best of our knowledge, the study of scalar conservation laws with discontinuous flux functions on networks is still missing in the traffic flow literature. However, in the context of supply chains with discontinuous flux such considerations have been already done [10,14]. We remark that supply chain models differ essentially from traffic flow models due to simpler dynamics and different coupling conditions.

    In this work, we aim to derive a traffic network model, where the dynamics on each road are governed by a scalar conservation law with discontinuous flux function in the unknown. For simplicity, we restrict to piecewise linear flux functions. Special emphasis is put on the coupling at junction points to ensure a unique admissible weak solution. In particular, we focus on dispersing junctions where the number of incoming roads does not exceed the number of outgoing roads and merging with two incoming and one outgoing road. The latter type of junction can be extended to the case of multiple incoming roads and a single outgoing one. In order to construct a suitable numerical scheme that is not based on regularization techniques we adapt the splitting algorithm originally introduced in [22]. Therein, the discontinuous flux is decomposed into a Lipschitz continuous flux and a Heaviside flux such that a two-point monotone flux scheme, e.g., Godunov, can be employed in an appropriate manner. This algorithm has been studied in [22] for the case of a single road only. However, in the network case, multiple roads with possibly disjunctive flux functions need to be considered at a junction point to ensure mass conservation. Hence, the key challenge is to determine the correct flux through the junction in an appropriate manner. Therefore, a detailed case distinction in accordance with the theoretical investigations is provided for the different types of junctions. The numerical results validate the proposed algorithm for some relevant network problems.

    The paper is organized as follows: in Section 2 we discuss the basic model and Riemann problems which permit to derive an exact solution. We extend the modeling framework to networks in Section 3 and focus on the coupling conditions. In Section 4, we introduce how the splitting algorithm [22] can be extended to also deal with the different types of junctions. Finally, we present a suitable discretization and numerical simulations in Section 5.

    In this section, we briefly recall the case of the Lighthill-Whitham-Richards (LWR) model [18,21] on a single edge with a flux function having a single decreasing jump at $ u^* \in [0, u_{{\rm{max}}}] $ that has been intensively studied in e.g., [19,22,24]. The shape of this flux function is inspired by empirical research which suggest that the fundamental diagram is of a reverse $ \lambda $ shape and consist of a free flow area and a drop to a congested area, see [4,5,8].

    Following [22], we consider the scalar conservation law Eq (2.1),

    $ {ut+f(u)x=0,(x,t)(a,b)×(0,T)=:ΠT,u(x,0)=u0(x)[0,umax],x(a,b),u(a,t)=r(t)[0,umax],t(0,T),u(b,t)=s(t)[0,umax],F(t)˜f(s(t)),t(0,T).
    $
    (2.1)

    More precisely the flux function is defined as follows Eq (2.2),

    $ f(u)={f1(u) ˆ=Flux1,ifu[0,u],f2(u) ˆ=Flux2,ifu(u,umax].
    $
    (2.2)

    We denote $ f_1(u^*) = f(u^*-) $ and $ f_2(u^*) = f(u^*+) $, where $ f(u^*-) > f(u^*+) $. In order to assign a single numerical value to $ f(u^*) $ we arbitrarily define $ f(u^*) = f(u^*-) $. One example is shown in Figure 1. The jump has magnitude

    $ α:=f(u)f(u+).
    $
    (2.3)
    Figure 1.  An example for a piecewise linear discontinuous flux function with discontinuity at $ u^* $, cf. [19].

    As usual we require for the flux function $ f : [0, u_{\rm{max}}] \rightarrow \mathbb{R}^+ $ with $ f(0) = f(u_{{\rm{max}}}) = 0 $ and in addition, we assume $ f_1 $ to be linear increasing and $ f_2 $ linear decreasing.

    As in [22], the multivalued version of $ f $ is defined by Eq (2.4),

    $ ˜f(u)={f(u),u[0,u),[f(u+),f(u)],u=u,f(u),u(u,umax].
    $
    (2.4)

    Finally, we have to discuss the imposed boundary conditions at $ x = b $ in Eq (2.1). Due to the flux discontinuity we have a non-standard boundary condition $ \mathcal{F}(t) \in \tilde f (s(t)) $. If the boundary condition at the right is $ s(t) \neq u^* $, then the additional boundary condition is redundant, i.e., $ \mathcal{F}(t) = f(s(t)) $. Otherwise, if $ s(t) = u^* $, we have to choose between two possible flux values $ \tilde f(u^*) \in \{f(u^*+), f(u^*-)\} $. Then, the function $ \mathcal{F} $ assigns a flux value to $ u(b,t) = s(t) $ according to [22]

    $ F(t)={f(u),ifthetrafficaheadofx=bisfreeflowing,f(u+),ifthetrafficaheadofx=biscongested.
    $
    (2.5)

    The state of traffic ahead of $ x = b $ is additional information that is determined independently from the other data of the problem.

    Remark 2.1. [22,Remark 1.3] We note that for the boundary condition at the left end the state of traffic ahead of $ x = a $ is known from the start values. Hence, no additional boundary conditions are necessary.

    The following assumptions are important for the proof of existence and uniqueness of solutions.

    Assumption 2.2. [22,Assumption 1.1] The initial data satisfies $ u_0(x) \in [0, u_{\mathit{{\rm{max}}}}] $ for $ x \in (a,b) $ and $ u_0 \in \mathit{{\rm{BV}}}([a,b]). $ The boundary data satisfies $ r(t),s(t) \in [0, u_{\mathit{{\rm{max}}}}] $ for $ t \in [0,T] $, and $ r,s \in \mathit{{\rm{BV}}}([0,T]) $.

    A weak solution is intended in the following sense:

    Definition 2.3 [22,Definition 1.1] A function $ u \in L^{\infty} (\Pi_T) $ is called a weak solution to the initial-boundary value problem (2.1) if there exists a function $ v \in L^{\infty}(\Pi_T) $ satisfying

    $ v(x,t) \in \tilde f(u(x,t)) \rm{ a.e.} $

    such that for each $ \psi \in C_0^1((a,b) \times[0,T)), $

    $ \int_0^T \int_a^b (u\psi_t + v\psi_x){\rm{d}}x{\rm{d}}t + \int_a^bu_0(x)\psi(x,0){\rm{d}}x = 0. $

    As usual, weak solutions do not lead to a unique solution and additional criteria are necessary to rule out physically incorrect solutions. In particular, the discontinuity of the flux prohibits from directly using the classical approaches. Note that in [22] an adapted version of Oleinik's entropy condition [9] is used to single out the correct solution, while in [24] the convex hull construction [17] is used to construct solutions to Riemann problems.

    Here, we will concentrate on the convex hull construction. For completeness we will shortly recall the solutions to Riemann problems considered in [24] as they are essential in order to construct a Riemann solver at a junction.

    We consider a Riemann problem with initial data $ (u_L,u_R) $ and its exact solution. For simplicity we will focus on a piecewise linear flux function as in [24] and in Figure 1. The exact solution is then derived with a regularization of the flux and the convex hull construction. As regularizations are not unique we choose a simpler regularization than the one in [24]. In fact, we use the regularization introduced in [19] for a piecewise quadratic flux function and apply it to the piecewise linear flux function. In the following we simply recall the results from [19] in case of a linear flux function for the sake of completeness and in order to understand better the upcoming analysis. Note that the the results also coincide with the ones of [24].

    We consider the following flux function by Eq (2, 6),

    $ f(u)={d1u+d0,xu(Flux1),e1u+e0,x>u(Flux2)
    $
    (2.6)

    with the regularized flux function given by Eq (2.7),

    $ fϵ(u)={f1(u)=d1u+d0,0uu,fϵmid(u)=1ϵ(f1(u)fϵ2(u+ϵ))(uu)+f1(u),u<u<u+ϵ,fϵ2(u)=e1u+e0,u+ϵ<uumax.
    $
    (2.7)

    We define $ u_\epsilon : = u^* + \epsilon $. An example of the regularization is shown in Figure 2. Even though the regularization differs slightly from the one used in [24], the different cases to consider and the approach to solve the Riemann problems are completely analogous. We refer to [24,Section 3] for further details. We define $ h: [0, u_{\rm{max}}] \rightarrow \mathbb{R} $ to be a function connecting the points $ (u^*, f(u^*+)) $ and the initial value lying on Flux 2. In addition, $ \gamma $ denotes the intersection of $ h $ and Flux 1, see also Figure 3. The point $ \gamma $ distinguishes the cases 3 and 4 in the following discussion.

    Figure 2.  Example of the regularization of a piecewise linear flux function.
    Figure 3.  The dotted line shows the function $ h $, connecting $ (u^*,f(u^*+)) $ and $ u_R $. The value $ \gamma $ is the intersection point.

    Case 1: Either $ u_L < u^*, u_R < u^* $ or $ u_L > u^*, u_R > u^* $

    This case corresponds to the classical case of solving Riemann problems, where the solution consists of a single rarefaction wave or shock, see [17].

    Case 2: $ u_R < u^* < u_L $

    By using the smallest convex hull approach the solution consists of a contact line following $ f(u) $, connecting $ u_R $ and $ u^* $ and a shock, connecting $ u^* $ and $ u_L $. The speed of the contact discontinuity is given by $ d_1>0 $ and the speed of the shock can be calculated with the Rankine-Hugoniot condition, i.e.,

    $ s=f(uL)f(u)uLu<0,
    $

    where we recall that $ f(u^*) = f(u^*-) $ holds. The exact solution is then given by Eq (2.8),

    $ u(x,t)={uL,ifx<st,u,ifstxd1t,uR,ifx>d1t.
    $
    (2.8)

    Case 3: $ \gamma < u_L < u^* < u_R $

    Here, the solution is given by a shock connecting $ u_L $ and $ u^* $ and contact discontinuity reconnecting it to $ u_R $. The speed of the shock equals:

    $ s = \frac{f(u^*+)-f(u_L)}{u^*-u_L} < 0. $

    Note that due to $ u_L<u_R $ the convex hull construction chooses the flux $ f(u^*+) $ at $ u^* $. The contact discontinuity connecting $ u^* $ and $ u_R $ moves at the speed $ e_1<0 $. Hence,

    $ u(x,t)={uL,ifx<st,u,ifstxe1t,uR,ifx>e1t.
    $
    (2.9)

    Case 4: $ u_L < u^* < u_R $ and $ \gamma \geq u_L $

    In this case, we get only one shock connecting $ u_L $ and $ u_R $ as a solution due to the condition $ \gamma \geq u_L $. The speed of the shock can be calculated by the Rankin-Hugoniot condition. The exact solution is given by Eq (2.10),

    $ u(x,t)={uL,ifx<st,uR,ifxst.
    $
    (2.10)

    Remark 2.4. As aforementioned in [19] Riemann solutions for piecewise quadratic discontinuous flux functions are derived. They also cover the case of a piecewise linear flux function if the quadratic terms are zero. For a general quadratic discontinuous flux, the solutions are more involved since no contact discontinuities occur.

    Up to now, we have not addressed the case, where one of the boundary conditions equals the critical density $ u^* $. In this situation so called zero waves can occur, see [24]. Those waves travel with speed $ \mathcal{O}(\frac{1}{\epsilon}) $ from right to left and have strength $ \mathcal{O}(\epsilon) $, depending on the regularization parameter $ \epsilon $. Although they have no strength, information can be exchanged between Riemann problems lying next to each other influencing the solution. For more details on the origin of zero waves in the context of traffic flow we refer to [24]. In particular, they arise when double Riemann problems are considered. For a detailed discussion of these double Riemann problems we refer to [24,Section 4].

    Next, we focus on networks where we allow for discontinuous flux functions. The key idea is to consider the regularized flux function $ f^\epsilon $ and take the limit $ \epsilon \rightarrow 0 $ to obtain results for the original problem at intersections. In doing so, the standard approaches can be directly applied to the regularized flux function $ f^\epsilon $.

    We start with a short introduction to the network setting. For more details on traffic flow network models we refer the reader to [11] and the references therein.

    Let $ G = (\mathcal{V}, \mathcal{E}) $ be a directed finite graph consisting of a nonempty set $ \mathcal{V} $ of vertices and a nonempty set of edges $ \mathcal{E} $ representing the roads. We denote by $ \delta_v^- $ and $ \delta_v^+ $ the set of all incoming and outgoing edges for every vertex $ v \in \mathcal{V} $. Every edge $ e \in \mathcal{E} $ is interpreted as an interval $ I_e = (a_e,b_e) \subset \mathbb{R} $ representing the spatial extension of the road $ e $. For all edges $ e $, which do not discharge into a vertex, i.e. $ e \notin \cup_{v \in \mathcal{V}} \delta_v^- $, we set an infinite road length by $ b_e = +\infty $. Analogously, for all edges $ e $, which do not originate from a vertex, i.e. $ e \notin \cup_{v \in \mathcal{V}} \delta_v^+ $, we set an infinite road length by $ a_e = -\infty $.

    We call the couple $ (\mathcal{I}, \mathcal{V}) $ with $ \mathcal{I} = \{I_e: e\in \mathcal{E}\} $ a road network.

    In order to derive the network solution, we restrict to the description of a single junction $ v\in\mathcal{V} $. Hence, from now on we consider a fixed junction with $ n $ incoming and $ m $ outgoing roads. Then, we need additional information how vehicles are distributed from one road to another. Hence, we need to define a distribution matrix $ A \in [0,1]^{n \times m} $ as in [11]:

    $ A = \left( {β1,1β1,mβn,1βn,m
    } \right). $

    To conserve the mass we assume $ \sum_{j = 1}^m \beta_{i,j} = 1 $ for every $ i \in \{1, \dots, n\} $. Furthermore, some technical assumptions on the distribution matrix are necessary which can be found in ref. [11] (Eq (4.2.4)). These assumptions are necessary to obtain unique solutions at the junction for the case $ n\leq m $. Following ref. [11] (definition 4.2.4) and definition 2.3, weak solutions at a junction are intended in the following sense

    Definition 3.1 (Weak solution at a junction). Let $ v \in \mathcal{V} $ be a vertex with $ n $ incoming and $ m $ outgoing roads, represented by the intervals $ I_1^{{\rm{in}}}, \dots I_n^{{\rm{in}}} $ and $ I_1^{{\rm{out}}}, \dots I_n^{{\rm{out}}} $. A weak solution at the junction $ v $ is a collection $ u $ and $ w $ of functions $ u_i^{{\rm{in}}}: I_i^{{\rm{in}}} \times \mathbb{R}_{\geq 0} \rightarrow \mathbb{R} $, $ w_i^{{\rm{in}}}: I_i^{{\rm{in}}} \times \mathbb{R}_{\geq 0} \rightarrow \mathbb{R} $ with $ w_i^{{\rm{in}}}(x,t)\in \tilde f(u_i^{{\rm{in}}}(x,t)) $ a.e., where $ i = 1, \dots,n $, $ u_j^{{\rm{out}}} $ : $ I_j^{{\rm{out}}} \times \mathbb{R}_{\geq 0} \rightarrow \mathbb{R} $ and $ w_j^{{\rm{out}}}: I_j^{{\rm{out}}} \times \mathbb{R}_{\geq 0} \rightarrow \mathbb{R} $ with $ w_j^{{\rm{out}}}(x,t)\in \tilde f(u_j^{{\rm{out}}}(x,t)) $ a.e., where $ j = 1, \dots,m $, respectively, such that

    $ ni=1(T0Iini(uini(x,t)tϕini(x,t)+wini(x,t)xϕini(x,t))dxdt)+mj=1(T0Ioutj(uoutj(x,t)tϕoutj(x,t)+woutj(x,t)xϕoutj(x,t))dxdt)=0,
    $

    for every collection of test functions $ \phi_i^{{\rm{in}}} \in C_0^1((a_i^{{\rm{in}}},b_i^{{\rm{in}}}] \times [0,T]), i = 1, \dots,n $ and $ \phi_j^{{\rm{out}}} \in C_0^1((a_j^{{\rm{out}}},b_j^{{\rm{out}}}] \times [0,T]), j = 1, \dots,m $, satisfying

    $ \phi_i^{{\rm{in}}}(b_i^{{\rm{in}}}, \cdot) = \phi_j^{{\rm{out}}}(a_j^{{\rm{out}}}, \cdot), \quad \partial_x \phi_i^{{\rm{in}}}(b_i^{{\rm{in}}}, \cdot) = \partial_x \phi_j^{{\rm{out}}}(a_j^{{\rm{out}}}, \cdot), $

    for $ i \in {1,\dots n} \rm{ and } j \in {1, \dots,m}. $

    Additionally, in order to get unique solutions, we will consider the following concept of admissible solutions, which adapts ref. [11] (rule (A) and (B), p. 81) to the discontinuous setting:

    Definition 3.2 (Admissbile Weak Solution). We call $ u $ and $ w $ (as defined in definition 3.1) an admissible weak solution at a junction $ v \in \mathcal{V} $ with the corresponding distribution matrix $ A $ if

    1. $ u_i^{{\rm{in}}}(\cdot,t) \in {\rm{BV}}(I_i^{{\rm{in}}}), u_j^{{\rm{out}}}(\cdot,t) \in {\rm{BV}}(I_j^{{\rm{out}}}) \;{\rm{ for\; all }}\; i,j, $

    2. $ u $ and $ w $ constitute a weak solution at the junction $ v $,

    3. $ {w_j^{{\rm{out}}}}(a_j^{{\rm{out}}} +,t)) = \sum_{i = 1}^n \beta_{i,j} {w_i^{{\rm{in}}}}(b_i^{{\rm{in}}} -,t)) \;{\rm{ for\; all }}\; j = 1, \dots, m \rm{ and } t \geq 0, $

    4. $ \sum_{i = 1}^n {w_i^{{\rm{in}}}}(b_i^{{\rm{in}}} -,t)) $ is maximal subject to 2. and 3.

    In particular, the maximization of the inflow with respect to the distributions parameters and the technical assumption [11] guarantee the uniqueness of solutions for a continuous flux.

    If $ n > m $, that means more incoming than outgoing roads, we need so-called right of way parameters which prescribe in which order the vehicles are going to drive. We consider the special case $ n = 2 $ and $ m = 1 $. Let $ C \in { \mathbb{R}} $ denote the amount of vehicles that can move on the road. The right of way parameter $ q \in (0,1) $ describes the following situation: If not all vehicles can move to the outgoing road, then $ qC $ is the amount of vehicles that enter from the first and $ (1-q)C $ the amount from the second road. If there are more roads, we need to define several right of way parameters.

    For solving the maximization problems imposed by the definition 3.2 so-called supply and demand functions can be used, see [16]. The demand describes the maximal flux the incoming road wants to send. In contrast, the supply describes the maximal flux the outgoing road is able to absorb. The definition of the supply and demand functions of the regularized function is straightforward. As $ f(u^*) \in \{ f(u^*+), f(u^*-)\} $ is not unique, we have to consider the limit for $ \epsilon \rightarrow 0 $ to determine the correct values of the supply and demand for the discontinuous flux function. Nevertheless, it becomes apparent that the maximal flux in both cases is given by $ f(u^*-) $. This leads to the following definition:

    Definition 3.3. For a network with flux function $ f $ having a discontinuity at $ u^* $, the demand is given by

    $ D(u)={f(u),u[0,u),f(u),u[u,umax].
    $
    (3.1)

    On the contrary, the supply reads as

    $ S(u)={f(u),u[0,u),f(u),u(u,umax]
    $
    (3.2)

    and

    $ S(u)={f(u),freeflowing,f(u+),congested.
    $
    (3.3)

    Remark 3.4. If we consider the regularized flux function $ f^\epsilon $, the definition of supply and demand function is completely analogous by replacing $ f $ with $ f^\epsilon $.

    In order to show existence and uniqueness in the discontinuous case we need to define an additional function. For a regularized flux function we notice that for every flux value, we get two different density values, see left picture in Figure 4. As different density values lead to different solutions, we need to be able to distinguish them and choose the correct solution. In the continuous or regularized case a mapping usually called $ \tau $ is used to determine the density values for which the flux values coincide, see also left picture in Figure 4 and ref. [11] (definition 3.2.6) for the precise definition of $ \tau $. This situation becomes more complicated for the discontinuous flux functions. For $ \epsilon \rightarrow 0 $ and $ f(u) > f(u^*+) $, the corresponding density value given by $ \tau $ can converge to the discontinuity area. This situation is pictured in Figure 4 on the right. Therefore, we need to adapt the definition of the mapping $ \tau $ to also cover this case and determine the corresponding density value. This mapping is defined in the following and is denoted by $ \eta $. For readability, we also define as an abbreviation $ f^{-1}_+: = f^{-1}(f(u^*+)) $.

    Figure 4.  The left panel shows the regularized version. The right panel shows the discontinuous flux $ f $. The grey area displays the interval, where $ f(\eta(u)) \leq f(u) $ and $ \eta(u) = u^* $.

    Definition 3.5. Let the function $ \eta : [0, u_{\rm{max}}] \rightarrow [0, u_{\rm{max}}] $ satisfy:

    1. $ f(\eta(u)) = f(u) $ for every $ u \in \bigl([0, { f^{-1}_+}] \cup [u^*,u_{\rm{max}}]\bigr) $ and $ \eta(u) \neq u $ for every $ u \in \bigl([0, { f^{-1}_+}] \cup (u^*,u_{\rm{max}}]\bigr) $.

    2. For $ u \in \bigl({ f^{-1}_+}, u^*\bigr) $ it holds, $ \eta(u) = u^* $ and $ f(\eta(u)) = f(u) $.

    Note that if $ u \in \bigl({ f^{-1}_+}, u^*\bigr) $, the flux value is given by $ f(u)\in(f(u^*+),f(u^*-)) $.

    Remark 3.6. We note that the mapping $ \tau $ from ref. [11] (definition 3.2.6) is in principle given by the first statement in definition 3.5 with the separated interval $ [0, { f^{-1}_+}] \cup [u^*,u_{\rm{max}}] $ being replaced by $ [0,u_{\rm{max}}] $. Hence, the only difference between $ \tau $ and $ \eta $ lies in the treatment of $ u \in \bigl({ f^{-1}_+}, u^*\bigr) $.

    Now, we present a Riemann solver for two types of junctions. First, we consider a junction with $ n $ incoming and $ m $ outgoing roads with $ n\leq m $, second, a junction with $ n = 2 $ incoming and $ m = 1 $ outgoing roads. For the continuous case the existence of solutions satisfying in particular the definition 3.2 has been shown in [11]. Since the continuous flux is not multivalued, the definition in 3.2 coincides with ref. [11] (Definition 4.2.4). As we have seen in Section 2.1, the solution to a Riemann problem may produce more than one wave. In order to have a well-defined network problem it is important that the waves produced at the junction have negative speed on the incoming and positive speed on the outgoing roads. This is ensured by the following theorem:

    Theorem 3.7. Let $ w \in \mathcal{V} $ be a junction with $ n $ incoming and $ m $ outgoing roads with $ n\leq m $. For every constant initial datum $ u_{i,0}^{\mathit{{\rm{in}}}} , u_{j,0}^{\mathit{{\rm{out}}}} \in [0, u_{\mathit{{\rm{max}}}}],\ i = 1,\ldots,n,\ j = 1,\ldots,m $, there exists a unique admissible weak solution $ u $ in the sense of definition 3.2. Moreover, for the incoming road the solution is given by the wave $ (u_{i,0}^{\mathit{{\rm{in}}}} , u_{i}^{\mathit{{\rm{in}}}}),\ i = 1,\ldots,n $ and for every outgoing road the solution is the wave $ (u_{j,0}^{\mathit{{\rm{out}}}},u_{j}^{\mathit{{\rm{out}}}}), j = 1,\ldots,m $, where

    $ uini{{uini,0}(η(uini,0),umax],ifuini,0[0,f1+],{uini,0}[u,umax],ifuini,0(f1+,u),[u,umax],ifuini,0[u,umax],
    $
    (3.4)

    and

    $ uoutj{[0,u],ifuoutj,0[0,u),[0,u],ifuoutj,0=uandfreeflowing,{u}[0,f1+),ifuoutj,0=uandcongested,{uoutj,0}[0,η(uoutj,0)),ifuoutj,0(u,umax].
    $
    (3.5)

    Proof. Using the definition of the supply and demand functions in definition 3.3 and the results from [12,section 5.2.3] we can follow the proof of [12,theorem 5.1.2] and uniquely determine the inflows which maximize the flux through junctions subject to the distribution parameters. It remains to show that by the choice of the density values the correct waves are induced.

    We start with considering the outgoing roads. If $ u_{j,0}^{{\rm{out}}}\in [0,u^*) $, the choice $ u_j^{{\rm{out}}}\in[0,u^*] $ produces the classical case of a shock or contact discontinuity, both with the positive speed $ d_1 $, since the flux is a linear increasing function. If $ u_{j,0}^{{\rm{out}}} = u^* $, we need to know if the traffic ahead is congested or not in order to determine the correct flux value of $ u^* $ and the corresponding density. In particular, we are now concerned with a double Riemann problem by $ u_j^{{\rm{out}}} $, $ u_{j,0}^{{\rm{out}}} = u^* $ and the information about the traffic ahead. By the analysis done in [24,Section 4] the choices of $ u_j^{{\rm{out}}} $ give rise to waves with positive velocities, even though the velocities differ depending on the traffic situation. Note that the free flowing case corresponds to case 1 in [24,Section 4] and the congested case to case 2. In particular, in the latter case $ f_j^{{\rm{out}}}\leq f(u^*+) = f(u^*) $ holds such that the solution is a shock with velocity always greater or equal to zero, compare also case 4 of Section 2.1.

    For $ u_{j,0}^{{\rm{out}}}\in (u^*,u_{{\rm{max}}}] $ we need to choose the density $ u_{j}^{{\rm{out}}} \in \{ u_{j,0}^{{\rm{out}}} \} \cup [0, \eta(u_{j,0}^{{\rm{out}}})) $ which corresponds to the incoming flux. Here, we are in case 4 of Section 2.1, as $ u_{j,0}^{{\rm{out}}}<\gamma $, such that the solution is a shock with positive velocity.

    On the contrary, considering the incoming roads and $ u_{i,0}^{\rm{in}} \in [0, { f^{-1}_+}] $, the choice $ u_i^{{\rm{in}}}\in \{ u_{i,0}^{{\rm{in}}} \} \cup (\eta(u_{i,0}^{{\rm{in}}}),u_{{\rm{max}}}] $ also corresponds to case 4, but this time the shock has a negative velocity.

    Now, let $ u_{i,0}^{\rm{in}} \in (u^*, u_{\rm{max}}] $ and $ f_i^{{\rm{in}}} $ denote the uniquely determined inflow. If $ f_i^{{\rm{in}}}<f(u^*+) $, the choice $ u_i^{{\rm{in}}}\in (u^*, u_{\rm{max}}] $ corresponds to the solution of a Riemann problem with a linear decreasing function and hence either a shock or contact discontinuity with negative velocity. On the contrary if $ f_i^{{\rm{in}}}\geq f(u^*+) $, we choose $ u_i^{{\rm{in}}} = u^* $. Here, the convex hull reconstruction induces a shock as a solution with velocity

    $ s=finif(uini,0)uuini,0<0.
    $

    Further, if $ u_{i,0}^{\rm{in}} = u^* $, we would normally also need to consider a double Riemann problem for which we would need the information of boundary data of the incoming roads. Nevertheless, the solutions at the junction are the same. As can be seen in [24,section 4] cases 2 and 4, if $ f_i^{{\rm{in}}}<f(u^*+) $ the choice $ u_i^{{\rm{in}}}\in (u^*, u_{\rm{max}}] $ induces a contact discontinuity with negative velocity as the solution at the junction. If $ f_i^{{\rm{in}}}\geq f(u^*+) $ we choose $ u_i^{{\rm{in}}} = u^* $ which gives a constant solution.

    Now, let us turn to the remaining case of $ u_{i,0}^{\rm{in}} \in ({ f^{-1}_+}, u^*) $. Again the solution depends on the inflow. If $ f_i^{{\rm{in}}}<f(u^*+) $, the choice of $ u_{i}^{{\rm{in}}}\in\{u_{i,0}^{{\rm{in}}}\} \cup (u^*,u_{{\rm{max}}}] $ induces negative waves. Nevertheless, depending on $ u_{i,0}^{{\rm{in}}}> \gamma $ or $ u_{i,0}^{{\rm{in}}}\leq \gamma $ we are either in case 3 or case 4 of Section 2.1. Even though case 3 induces two waves, all have a negative velocity. If $ f_i^{{\rm{in}}}\geq f(u^*+) $, choosing $ u_i^{{\rm{in}}} = u^* $ gives a shock with velocity

    $ s=finif(uini,0)uuini,00,
    $

    as $ f_i^{{\rm{in}}} $ is bounded from above by $ f(u_{i,0}^{{\rm{in}}}) $ due to the definition of the demand.

    Hence, the choices of the densities induce the correct waves.

    Now, we consider the case of more incoming than outgoing roads. Exemplary, we study the 2–to–1 situation, even though the results can be easily extended to the $ n $–to–1 case. We assume that a right of way parameter $ q $ is given and the unique density of the discontinuous flux can be determined according to the following theorem:

    Theorem 3.8. Let $ w \in \mathcal{V} $ be a junction with $ n = 2 $ incoming and $ m = 1 $ outgoing roads with a right way parameter $ q \in (0, 1) $. For every constant initial datum $ u_{i,0}^{\mathit{{\rm{in}}}} , u_{1,0}^{\mathit{{\rm{out}}}} \in [0, u_{\mathit{{\rm{max}}}}] $ with $ i \in \{1,2\} $, there exists a unique admissible weak solution $ u $ in the sense of definition 3.2. Moreover, for every incoming road $ i $ the solution is given by the wave $ (u_{i,0}^{\mathit{{\rm{in}}}} , u_{i}^{\mathit{{\rm{in}}}}) $ and for the outgoing edge the solution is the wave $ (u_{1,0}^{\mathit{{\rm{out}}}},u_{1}^{\mathit{{\rm{out}}}}) $, where

    $ uini{{uini,0}(η(uini,0),umax],ifuini,0[0,f1+],{uini,0}[u,umax],ifuini,0(f1+,u),[u,umax],ifuini,0[u,umax],
    $
    (3.6)

    and

    $ uout1{[0,u],ifuout1,0[0,u),[0,u],ifuout1,0=uandfreeflowing,{u}[0,f1+),ifuout1,0=uandcongested,{uout1,0}[0,η(uout1,0)),ifuout1,0(u,umax].
    $
    (3.7)

    Proof. Following [11,Section 3.2.2], the flux values at the junction can be calculated with the following steps:

    1. Calculate the maximal possible flux $ f_{\rm{max}} = \min\{ D(u_{1,0}^{\rm{in}}) + D(u_{2,0}^{\rm{in}}), S(u_{1,0}^{\rm{out}})\} $.

    2. Consider the right of way parameter and the flux maximization and calculate the intersection $ P = (qf_{\rm{max}}, (1-q)f_{\rm{max}}) $.

    3. If $ P $ is in the feasible area $ \Omega = \{ z\in [0, D(u_{1,0}^{in})] \times [0, D(u_{2,0}^{\rm{in}})]: z_1 + z_2 \leq f_{\rm{max}}\} $, set $ (z_1,z_2) = P $. If not, determine the point in $ \Omega \cap \{(z_1,z_2): z_1+z_2 = f_{{\rm{max}}}\} $ which is the closest to the line $ z_2 = (1-q)/q z_1 $.

    4. The flux values are given by $ f^{\rm{in}}_1 = z_1,\ f^{\rm{in}}_2 = z_2,\ f^{\rm{out}}_1 = z_1+z_2 $.

    Completely analogous to theorem 3.7 we can show that the choice of the densities admits the correct wave speeds.

    The Riemann solutions proposed in theorem 3.7 and theorem 3.8 are the key ingredients for the splitting algorithm on networks in the next section.

    Different problems might occur when designing a numerical scheme for a conservation law with discontinuous flux. However, the main difficulties are induced by the zero waves. Since these waves have infinite speed, the regular CFL condition is scaled by the regularization parameter $ \epsilon $. This leads to very inefficient step sizes and high computational effort. Possible ways to avoid using the regular CFL condition are implicit methods [20] or the splitting algorithm introduced in [22], see also [2]. In this section we recall the main ideas of the splitting algorithm and also present its extension to networks.

    We consider a flux function $ f $ with discontinuity at $ u^* $, which satisfies assumption 2.2. The jump has magnitude $ \alpha = f(u^*-)-f(u^*+) $. The main idea of the splitting algorithm is to split the discontinuous function $ f $ into two parts, namely $ p $ and $ g $, such that $ f(u) = p(u) + g(u) $. Here, $ g $ is a piecewise constant function defined by

    $ g(u)=αH(uu),
    $

    where $ H $ denotes the Heaviside function. This function extracts the jump from the function $ f $. In addition, we define $ p(u) = f(u) - g(u). $ Hence, $ p $ is continuous and can be used for numerical algorithms. An example can be seen in Figure 5.

    Figure 5.  Based on [22,Figure 1]. The left panel shows the discontinuous flux function, the right panel shows the splitted version. The dotted line represents the piecewise constant function $ g $, the solid line shows the continuous function $ p $.

    This case has been already treated in [22] and will be the basis for the splitting algorithm on networks. When solving the scalar conservation law (2.1) on a single road, the boundary value in the case $ s(t) = u^* $ is not unique. In the continuous case, the flux value is determined by $ \mathcal{F}(t) $ in Eq (2.5). When splitting the function $ f $, the discontinuity is shifted to the function $ g $. So analogously, we define the multivalued version of $ g $ in order to be able to assign a flux value in the critical case, i.e.,

    $ ˜g(u)={0,u[0,u),[α,0],u=u,α,u(u,umax].
    $
    (4.1)

    Furthermore, we define $ \mathcal{G}(t) \in \tilde g(s(t)). $ It holds,

    $ F(t)=f(u)G(t)=0,F(t)=f(u+)G(t)=α.
    $
    (4.2)

    Additionally to the assumptions 2.2, we assume:

    Assumption 4.1. [22,Assumption 1.1] The initial data satisfies $ u_0(x) \in [0, u_{\mathit{{\rm{max}}}}] $ for $ x \in (a,b),\ u_0 \in \mathit{{\rm{BV}}}([a,b]) $ and $ g(u_0) \in \mathit{{\rm{BV}}}([a,b]). $ We also assume that $ \mathcal{G}(t) \in [-\alpha,0] $ for $ t \in [0,T] $, and $ \mathcal{G} \in \mathit{{\rm{BV}}}([0,T]) $.

    Remark 4.2. We emphasize that the original splitting algorithm for a single road [22] is not limited to piecewise linear discontinuous flux functions. Another prominent example might be concave piecewise quadratic flux functions with discontinuity again at $ u^* $, cf. [19,22,24].

    Then, we are able to handle the flux function $ p $ and the jump $ g $ separately. As usual, we discretize the domain $ \Pi_T $ by a spatial mesh and a temporal mesh. Let $ \Delta x $ denote the spatial and $ \Delta t $ the temporal step size. We define the grid constant via

    $ \lambda = \frac{\Delta t}{\Delta x}. $

    For an integer $ {K} \in \mathbb{N} $ the grid size is given by $ \Delta x = \frac{b-a}{{K}+1} $ and the grid points by $ x_{{k}} = a+{k} \Delta x $. Further, we define the $ {K} $ disjunct intervals $ I_{{k}} = [x_{{k}} - \frac{\Delta x}{2}, x_{{k}} + \frac{\Delta x}{2}) $ and $ \mathcal{{K}} = \{ 1, \cdots {K}\} $ and $ \mathcal{{K}}^+ = \{ 0, \cdots {K}+1\} $. The time steps are defined by $ t^n = n \Delta t $, for $ n = 0,1, \dots,N $. The value $ N \in \mathbb{N} $ is chosen such that the time horizon fulfills $ T \in [t^N, t^N+\Delta t). $

    As the algorithm splits the function $ f $ into two parts, $ f(u) = p(u) + g(u) $, every step of the algorithm consists of two half steps. Here, the first half step corresponds to approximately solving $ u_t + g(u)_x = 0 $ denoted by $ U_{{k}}^{n+\frac{1}{2}} \approx u(x_{{k}},t^n). $ Consequently, the second step approximates $ u_t + p(u)_x = 0 $ and is denoted by $ U_{{k}}^n \approx u(x_{{k}},t^n). $

    We denote the backward spatial difference by $ \Delta_- U_{{k}}^n = U_{{k}}^n - U_{{k}-1}^n $ and the forward spatial difference by $ \Delta_+ U_{{k}}^n = U_{{k}+1}^n - U_{{k}}^n $. The initial values are denoted by $ U_{{k}}^0 = u_0(x_{{k}}) $, the boundary values by

    $ rn=rn+12=r(tn),Un+120=Un0=rnsn=sn+12=s(tn)Un+12K+1=UnK+1=sn.
    $

    The function $ g $ covers the issue of $ s(t^n) = u^* $. We define

    $ gn+12K+1=gnK+1=G(tn)={0,ifs(tn)<u,0,ifs(tn)=u,trafficaheadofx=bisfreeflowing,α,ifs(tn)=u,trafficaheadofx=biscongested,α,ifs(tn)>u.
    $
    (4.3)

    That means, we can describe the boundary value $ \mathcal{F}(t) $ via $ \mathcal{G}(t) $. For the splitting algorithm in its original fashion [22] we additionally need the following function.

    Definition 4.3 [22,Eq (3.7)] Let $ G(z) : = z - \lambda g(z) $ with grid constant $ \lambda $ and $ \tilde G $ its multivalued version. The function $ \tilde G $ is strictly increasing and has a unique inverse, $ z \rightarrow \tilde G^{-1}(z) $, which is a single valued function. It is Lipschitz continuous and nondecreasing. It holds, $ 0 \leq \tilde G^{-1}(z) \leq z $. The functions are given by

    $ ˜G(u)={u,u[0,u),[u,u+λα],u=u,u+λα,u(u,umax], ˜G1(u)={u,u[0,u),u,u[u,u+λα),uλα,u[u+λα,umax+λα].
    $

    The splitting algorithm [22] can be then expressed as

    $ {{Un+1/2k=˜G1(Unkλgn+1/2k+1),k=K,K1,,1,gn+1/2k=(Un+1/2kUnk+λgn+1/2k+1)/λ,k=K,K1,,1,Un+1k=Un+1/2kλΔpg(Un+1/2k+1,Un+1/2k), kK.
    $
    (4.4)

    Note that the first half step, which includes the first two equations, is implicit. Nevertheless, instead of solving a nonlinear system of equations, the equation can be solved backwards in space starting with $ {k} = {K} $. The last equation in (4.4) uses the Godunov scheme for the second half step which flux function is denoted by $ p^g $.

    We note that for the implicit equation a CFL condition is not needed, but it is required for the third step. As $ p $ is continuous, we can use the standard CFL condition of the Godunov scheme. Note that similar to [22] also other two-point monotone schemes, e.g., the Lax-Friedrichs method, can be used.

    As shown in [22,Theorem 5.1] the splitting algorithm (4.4) converges to a weak solution of Eq (2.1). However, obtaining a similar statement about weak entropy solutions is still an open problem.

    The key idea to numerically solve such discontinuous conservation laws on networks is to use the splitting algorithm only on the roads and determine the correct in- and outflows at the boundaries by the help of the Riemann solver established in, e.g., theorem 3.7. As the splitting algorithm works with flux values, there is no need to compute the exact densities at the junction. Instead we need to know how the solution at the junction influences the flux values. The algorithm for a single junction is depicted in algorithm 1. The general description of the algorithm allows for either junctions with given distribution or right of way parameters. For simplicity, we assume that each road is represented by the same interval $ (a,b) $, such that we have $ I_1^{{\rm{in}}} = \dots = I_n^{{\rm{in}}} = I_1^{\rm{out}} = \dots = I_m^{\rm{out}} = (a,b) $.

    Remark 4.4. Note that this simplification enables the use of the same grid points on each road which spares further sub- or superindices. However, the algorithm can be easily adapted to different road lengths.

    We assume in the following that the space and time grid is the same as in the previous subsection. The approximate solutions are denoted by $ U_{i,k}^{n,{\rm{in}}}\approx u_{i,k}^{n,{\rm{in}}}(x_k,t^n) $ for $ i = 1,\dots,n $ and $ U_{j,k}^{n,{\rm{out}}}\approx u_{j,k}^{n,{\rm{out}}}(x_{{k}},t^n) $ for $ j = 1,\dots,m $. The half steps are denoted accordingly. We assume that the initial values $ u^{\rm{in}}_{i,0}(x) $ and $ u^{\rm{out}}_{j,0}(x) $ and the boundary values $ u^{\rm{in}}_{i,0}(a,t) $ and $ u^{\rm{out}}_{j,0}(b,t) $ are given. Then, we can directly set the starting values as $ U_{i,k}^{0,{\rm{in}}} = u^{\rm{in}}_{i,0}(x_k) $ for $ i = 1,\dots,n $ and $ U_{j,k}^{0,{\rm{out}}} = u^{\rm{out}}_{j,0}(x_k) $ for $ j = 1,\dots,m $. The left boundary values for the incoming roads are given by $ U_{i,0}^{n,{\rm{in}}} = u^{\rm{in}}_{i,0}(a,t^n) $ and the right boundary values for the outgoing roads by $ U_{j,0}^{n,{\rm{out}}} = u^{\rm{out}}_{j,0}(b,t^n) $. The missing boundary data, i.e. the right one for the incoming road and the left one for the outgoing road, are determined by the flux values at the junction.

    The overall strategy of the splitting algorithm on a network consists of three important steps:

    1. Solve the optimization problem induced by definition 3.2 (in particular item 4) at the junction to calculate the flux values $ f_i^{\rm{in}}, f_j^{\rm{out}} $ with the discontinuous flux $ f $.

    Here, it is crucial to use the discontinuous flux function $ f $ when calculating demand and supply at the junction, which are necessary for the optimization problem, since otherwise not only the flux values can be different but also the decision whether supply or demand is active. Therefore, in line 6, we need to compute first the fluxes at the junction with the original discontinuous flux function $ f $.

    These flux values bring us now to:

    $ 2. $ Consider the corresponding density values at the junction applying either theorem 3.7 or theorem 3.8 depending on the type of junction. If the density value is greater or equal than $ u^* $, i.e., the road is congested, the flux value needs to be adjusted. Hence, we determine adjusted flux values $ f_\rm{adj.}^\rm{in/out} $.

    Table Algorithm 1.  Splitting algorithm for a fixed junction of dispersing or merging type.
    Require: number of incoming roads $ n $, number of outgoing roads $ m $, distribution matrix $ A $, (right of way parameters), flux $ f $, discontinuity $ u^* $, jump magnitude $ \alpha $, interval boundaries $ a,\ b $, number of grid points $ K $, time step size $ \Delta t $, initial values $ u^{\rm{in}}_{i,0}(x) $ and $ u^{\rm{out}}_{j,0}(x) $, boundary values $ u^{\rm{in}}_{i,0}(a,t) $ and $ u^{\rm{out}}_{j,0}(b,t) $
    Ensure: approximate solutions $ U_{i,k}^{n,{\rm{in}}} $ and $ U_{j,k}^{n,{\rm{out}}} $
    1: Initilization:
    2: $ \Delta x=(b-a)/(K+1) $ and $ \lambda=\Delta t/\Delta x $,
    3: $ U_{i,k}^{0,{\rm{in}}}=u^{\rm{in}}_{i,0}(x_k) $ for $ i=1,\dots,n $, $ U_{j,k}^{0,{\rm{out}}}=u^{\rm{out}}_{j,0}(x_k) $ for $ j=1,\dots,m $
    4: $ U_{i,0}^{n,{\rm{in}}}=u^{\rm{in}}_{i,0}(a,t^n) $, $ U_{j,0}^{n,{\rm{out}}}=u^{\rm{out}}_{j,0}(b,t^n) $
    5: for $ n=0,\dots,N-1 $ do
    6:    Solve the by definition 3.2 induced optimization problem at the junction based on the flux $ f $, which gives the fluxes $ f_i^{\rm{in}},\ f_j^{{\rm{out}}} $
    7:    Compute the densities at the junction with an appropriate Riemann solver
    8:    Compute the adjusted flux values for the incoming roads $ f_{i,\rm{adj.}}^{\rm{in}} $
    9:    for $ i=1,\dots,n $ do
    10:      $ g_{i,K+1}^{n+1/2,{\rm{in}}} = g_{i,K+1}^{n,{\rm{in}}} =f_i^{{\rm{in}}}-f_{i,\rm{adj.}}^{\rm{in}} $
    11:     for $ k=K,K-1,\dots,1 $ do
    12:        $ U_{i,k}^{n+1/2,{\rm{in}}} = \tilde G^{-1}\bigl(U_{i,k}^{n,{\rm{in}}} - \lambda g_{i,k+1}^{n+1/2,{\rm{in}}}\bigr) $
    13:       $ g_{i,k}^{n+1/2,{\rm{in}}} = \bigl(U_{i,k}^{n+1/2,{\rm{in}}} - U_{i,k}^{n,{\rm{in}}} + \lambda g_{i,k+1}^{n+1/2,{\rm{in}}}\bigl)/\lambda $
    14:      end for
    15:     $ p^g\bigl(U_{i,K+1}^{n+1/2,{\rm{in}}},U_{i,K}^{n+1/2,{\rm{in}}}\bigr)=f_{i,\rm{adj.}}^{\rm{in}} $
    16:     for $ k\in \mathcal{K} $ do
    17:        $ U_{i,k}^{n+1,{\rm{in}}} = U_{i,k}^{n+1/2,{\rm{in}}} - \lambda \Delta_{-}p^g\bigl(U_{i,k+1}^{n+1/2,{\rm{in}}},U_{i,k}^{n+1/2,{\rm{in}}}\bigr) $
    18:      end for
    19:    end for
    20:    for $ j=1,\dots,m $ do
    21:      Compute $ g_{j,K+1}^{n+1/2,{\rm{out}}} = g_{j,K+1}^{n,{\rm{out}}} $ as in (4.3)
    22:      for $ k=K,K-1,\dots,1 $ do
    23:        $ U_{j,k}^{n+1/2,{\rm{out}}} = \tilde G^{-1}\bigl(U_{j,k}^{n,{\rm{in}}} - \lambda g_{j,k+1}^{n+1/2,{\rm{out}}}\bigr) $
    24:        $ g_{j,k}^{n+1/2,{\rm{out}}} = \bigl(U_{j,k}^{n+1/2,{\rm{out}}} - U_{j,k}^{n,{\rm{out}}} + \lambda g_{j,k+1}^{n+1/2,{\rm{out}}}\bigl)/\lambda $
    25:     end for
    26:      $ p^g\bigl(U_{j,1}^{n+1/2,{\rm{in}}},U_{j,0}^{n+1/2,{\rm{out}}}\bigr)=f_j^{\rm{out}} $
    27:      for $ k\in \mathcal{K} $ do
    28:        $ U_{j,k}^{n+1,{\rm{out}}} = U_{j,k}^{n+1/2,{\rm{out}}} - \lambda \Delta_{-}p^g\bigl(U_{j,k+1}^{n+1/2,{\rm{out}}},U_{j,k}^{n+1/2,{\rm{out}}}\bigr) $
    29:      end for
    30:    end for
    31: end for

     | Show Table
    DownLoad: CSV

    Using the calculated (unadjusted) flux values from step one, we can determine the densities at the junction with the help of the appropriate Riemann solver (theorem 3.7 and Eqs (3.4)–(3.5) or theorem 3.8 and Eqs (3.6)–(3.7)) at the junction (line 7). Then, these density values can be used to calculate the corresponding flux value of $ p $ to get the boundary values and decide about an adjustment. This means that if the density given by the Riemann solver suggests that the traffic ahead is free flowing, no adjustment needs to be made - as the flux value already lies on $ p $. If the traffic ahead is congested we take the flux value on the curve $ p $ for the corresponding density value and not its original flux value on the curve $ f $. Both the first step, i.e solving the optimization problem at the junction, and the second step, i.e., how to adjust the flux, strongly depend on the junction type. Hence, we postpone a detailed discussion to the following subsections, in which we will consider the 1–to–1, 1–to–2 and 2–to–1 case in more detail.

    In addition, the first steps give us all the ingredients for the final step:

    $ 3. $ Away from the junction, the splitting algorithm can be used. At the junction itself, the adjusted fluxes are used to determine the missing boundary data.

    The adjusted flux values from the previous step are important for the second half step (line 17 or 28) of the splitting algorithm which uses a Godunov type scheme based on $ p $. Here, the missing fluxes in and out of the junction the Godunov scheme needs are simply given by these adjusted flux values. The flux values can be used directly as boundary conditions, see lines 15 and 26, and the exact density values are not needed for the scheme.

    Further boundary data is needed in the first half step of the algorithm, lines 12–13 and 23–24. Here, we start with $ k = K $ to avoid solving a nonlinear equation system. For example, the value $ g^{n,{\rm{out}}}_{j,K+1} $ is determined by the density value at the boundary. Nevertheless, for every outgoing road the boundary values are known and the value $ g^{n,{\rm{out}}}_{j,K+1} $ can be determined as before in Eq (4.3). However, for the incoming road, i.e., for $ g^{n,{\rm{in}}}_{i,K+1} $, the density value at the junction is necessary. From the adjustment made in step two we know whether the traffic ahead is congested or not. In addition, instead of working with the exact density values at the junction directly, we can also use the adjusted flows computed in the second step (or line 8). Here, we can simply compute the boundary value $ g^{n,{\rm{in}}}_{i,K+1} $ by the difference between the fluxes calculated with $ f $ (before adjustment) and the ones with $ p $ (after adjustment), i.e.,

    $ gn+1/2,ini,K+1=gn,ini,K+1=finifini,adj.
    $
    (4.5)

    Note that the definition of $ g^{n,{\rm{in}}}_{i,K+1} $ in Eq (4.5) is slightly different in comparison to Eq (4.3) or for the outgoing roads. Nevertheless, if the density value on the incoming road is smaller than $ u^* $, $ g^{n,{\rm{in}}}_{i,K+1} $ equals zero (no adjustment needs to be made) and if it is greater $ u^* $, it equals $ -\alpha $. Only in the case the density value equals $ u^* $ the situation is more involved, as the traffic on the incoming road can be congested and $ f(u^*-)>f^{{\rm{in}}}>f(u^*+) $ can occur such that we need $ g^{n,{\rm{in}}}_{i,K+1}\in (-\alpha,0) $.

    Furthermore, we can decrease the computational costs of the algorithm: In the second step (or in line 7 of algorithm 1) the density at the junction is computed. This can be very expensive and hence we aim to avoid this. In the third step we have seen that for the missing boundary data only the adjusted flux values are necessary and not the densities at the junction themselves. Therefore, studying first each junction type in detail allows to determine the corresponding flux values based on the density values and the supply and demand functions and the intermediate expensive step for the computation of the exact densities can be skipped. Hence, we combine the first and second step of the strategy in one single step. In the following, we will study a 1–to–1 junction in detail and present the tailored algorithm. As the strategy is completely analogous for the 1–to–2 junction and 2–to–1 junction, we will only present the algorithms and discuss important properties. The algorithms can then be used to replace the lines 6 to 8 in algorithm 1.

    Remark 4.5. The extension to $ n $–to–1 or 1–to–$ m $ junctions follows immediately from the upcoming discussion. For more complex junctions, it is possible to adapt the proposed strategy by studying the different solutions of the optimization problem at the junction.

    Remark 4.6. Further note that the presented strategy and also algorithm 1 can be used for arbitrary junctions and nonlinear discontinuous flux functions once an appropriate Riemann solver similar to theorem 3.7 and 3.8 is established.

    First, we consider a junction with one incoming and one outgoing road in detail. Let $ u_{1,0}^{\rm{in}}, u_{1,0}^{\rm{out}} $ be constant initial values. The flux values are given by $ f^{{\rm{in}}}_1 = f^{{\rm{out}}}_1 = \min\{D(u_{1,0}^{\rm{in}}),S(u_{1,0}^{\rm{out}})\} $. We distinguish the following cases:

    Case A: demand and supply are equal

    There are two different situations depending on the density value of the incoming road where demand and supply can be equal.

    1. $ \bf{u_{1,0}^{\rm{in}}} \bf{\geq u^*} $: Here, the traffic on the outgoing road needs to be free-flowing or equal to $ u^* $ and the traffic ahead of the outgoing road is free-flowing. Supply and demand equal the maximal possible flux $ f(u^*-) $. For these flux values the solution of the Riemann solver on the incoming and outgoing road is given by $ u^* $. As corresponding flux values already lie on $ p $, we do not need to adjust them.

    2. $ \bf{u_{1,0}^{\rm{in}}} \bf{< u^*} $: Demand and supply can also be equal, if the number of vehicles the demand wants to send is equal to the number of vehicles the supply can take. Therefore, $ u_{1,0}^{\rm{in}} < \gamma $ and $ u_{1,0}^{{\rm{out}}}>u^* $ need to hold. The solution proposed by the Riemann solver at the junction is the initial condition for each road. Hence, the traffic on the outgoing road stays congested and we need to adjust the flux value on the outgoing road, i.e., $ f_{1,\rm{adj}}^{\rm{out}} = f_1^{\rm{out}} + \alpha $. The flux on the incoming road stays free-flowing.

    Case B: supply is restrictive

    If the supply is restrictive, i.e. $ f^{{\rm{in}}}_1 = f^{{\rm{out}}}_1 = S(u_{1,0}^{\rm{out}})<D(u_{1,0}^{\rm{in}}) $, the outgoing road cannot take the whole amount of cars, the demand wants to send. From this follows that the traffic on the outgoing road is congested, i.e. $ u_{1,0}^{\rm{out}} > u^* $ or $ u_{1,0}^{\rm{out}} = u^* $ and the traffic ahead being congested. This has an effect on the incoming road. Again, we face two situations:

    1. $ \bf{u_{1,0}}^{\rm{in}} \bf{< u^*} $: Initially, the traffic on the incoming road is free-flowing and the traffic on the outgoing road is congested. As the supply is restrictive, the traffic on the incoming road congests as well, the solution of the Riemann solver is on both roads given by $ u_{1,0}^{\rm{out}} $. Therefore, we have to adjust the flux value, $ f_{1,\rm{adj}}^{\rm{in}} = f_1^{\rm{in}}+\alpha. $ In addition, we need to set $ {g_{1,K+1}^{n,{\rm{in}}}} = -\alpha $. The outgoing flux needs to be also adjusted, i.e., $ f_{1,\rm{adj}}^{\rm{out}} = f_1^{\rm{out}}+\alpha. $

    2. $ \bf{u_{1,0}}^{\rm{in}} \bf{> u^*} $: If the traffic on the incoming road is congested as well, the solution of the Riemann solver on both roads is again given by $ u_{1,0}^{\rm{out}} $. This leads to the same adjustments as in the previous case.

    Case C: demand is restrictive

    If the demand is restrictive, the outgoing road is able to take the whole amount of vehicles the demand wants to send. Here, we only have one possible situation for the initial condition on the incoming road:

    1. $ \bf{u_{1,0}}^{\rm{in}} \bf{\leq u^*} $: The solution of the Riemann problem on each road is now given by $ u_{1,0}^{\rm{in}} $, such that no adjustment is necessary.

    Note that we have the same flux function on each road. Hence, in the case $ {u_{1,0}}^{\rm{in}} {> u^*} $ either the supply is restrictive or demand and supply are equal.

    The whole procedure is summarized in algorithm 2.

    Table Algorithm 2.  1–to–1 junction.
    Require: Demand $ D_1 $, Supply $ S_1 $, flux $ f $, discontinuity $ u^* $, jump magnitude $ \alpha $
    Ensure: Flux values $ f_1^{\rm{in}}, f_1^{\rm{out}} $
    1: $ f_1^{\rm{in}} = \min \bigl\{D_1, S_1\bigr \} $
    2: $ f_1^{\rm{out}} = f_1^{\rm{in}} $
    3: if $ S_1 = D_1 $ and $ D_1 \neq f(u^*-) $ then
    4:    $ f_1^{\rm{out}} = f_1^{\rm{out}} + \alpha $
    5: else if $ S_1<D_1 $ then
    6:     $ f_1^{\rm{out}} = f_1^{\rm{out}} + \alpha $
    7:     $ f_1^{\rm{in}} = f_1^{\rm{in}} + \alpha $
    8: end if

     | Show Table
    DownLoad: CSV

    Remark 4.7. Theoretically, the Riemann solver in theorem 3.7 coincides for a 1–to–1 junction with the Riemann solver on a single road. Hence, the procedure described in algorithm 2 leads to the same solution. In contrast to that on the numerical level, the splitting algorithm for a 1–to–1 junction does not exactly coincide with the splitting algorithm used for a single road [22]. The reason for the computational difference is that the splitting algorithm for the 1–to–1 junction considers the exact solution of the Riemann problem at the junction point and hence uses exact values for $ g_{1,K+1}^{n,{\rm{in}}} $ while for a single road this value is only approximated using $ g_{1,K+1}^{n,{\rm{out}}} $.

    The difference to the 1–to–1 junction is now that we have to consider two supplies. So the case distinctions to determine the flux values are a bit more complex. However, the procedure itself does not change. The details can be seen in algorithm 3.

    We remark that if the inflow on the first road equals the demand and the restriction given by at least one of the supplies, the latter road needs to be congested. The Riemann solver states congestion such that the flux needs to be adjusted. Then, the incoming road either stays free flowing or the solution is given by $ u^* $. Here, no adjustment is necessary. Note that the conditions $ S_i\neq f(u^*-) $ are necessary for the specific case of $ \beta_i = 1 $.

    Table Algorithm 3.  1–to–2 junction.
    Require: Demand $ D_1 $, Supply $ S_1, S_2 $, distribution parameter $ \beta_1, \beta_2 $, flux $ f $, discontinuity $ u^* $, jump magnitude $ \alpha $
    Ensure: Flux values $ f_1^{\rm{in}}, f_1^{\rm{out}}, f_2^{\rm{out}} $
    1: $ f_1^{\rm{in}} = \min \bigl\{D_1, S_1/\beta_1, S_2/\beta_2 \bigr \} $
    2: $ f_1^{\rm{out}} = \beta_1 f_1^{\rm{in}}, f_2^{\rm{out}} = \beta_2 f_1^{\rm{in}} $
    3: if $ f_1^{\rm{in}} = D_1 $ then
    4:    if $ S_1 = D_1 $ and $ S_1 \neq f(u^*-) $ then
    5:      $ f_1^{\rm{out}} = f_1^{\rm{out}} + \alpha $
    6:   else if $ S_2 = D_1 $ and $ S_2 \neq f(u^*-) $ then
    7:      $ f_2^{\rm{out}} = f_2^{\rm{out}} + \alpha $
    8:    end if
    9: else if $ f_1^{\rm{in}}=S_1/\beta_1 $ or $ f_1^{\rm{in}}=S_2/\beta_2 $ then
    10:    if $ f_1^{\rm{in}}>f(u^*+) $ then
    11:     $ f_1^{\rm{in}}=f(u^*-) $
    12:    else
    13:     $ f_1^{\rm{in}}=f_1^{\rm{in}}+\alpha $
    14:    end if
    15:    if $ S_1/\beta_1=S_2\beta_2 $ and $ S_1\beta_1<D_1 $ then
    16:      $ f_1^{\rm{out}} = f_1^{\rm{out}} + \alpha $
    17:      $ f_2^{\rm{out}} = f_2^{\rm{out}} + \alpha $
    18:   else if $ S_1/\beta_1<S_2\beta_2 $ and $ S_1\beta_1<D_1 $ then
    19:      $ f_1^{\rm{out}} = f_1^{\rm{out}} + \alpha $
    20:   else if $ S_2/\beta_2<S_1\beta_1 $ and $ S_2\beta_2<D_1 $ then
    21:      $ f_2^{\rm{out}} = f_2^{\rm{out}} + \alpha $
    22:    end if
    23: end if

     | Show Table
    DownLoad: CSV

    If at least one of the supply restrictions is active, we need to adjust the corresponding flux values as in the 1–to–1 case and the incoming road. Nevertheless, here an interesting case (which is not possible in the 1–to–1 situation) can occur. The solution of the Riemann problem at the incoming road can be given by $ u^* $ and the inflow is in $ f^{{\rm{in}}}_1\in(f(u^*+),f(u^*-)) $, see also the proof of theorem 3.7. We need to adjust the flux to $ p(u^*) = f(u^*-) $. Hence, in this case if we use Eq (4.5) the adjustment is between $ 0 $ and $ -\alpha $. Finally, as in the 1–to–1 case there is no need for an adjustment if the demand is restrictive.

    Recall that for two incoming and one outgoing roads, the maximal possible flux on the outgoing road is given by $ f_{\rm{max}} = \min\{ D(u_{1,0}^{\rm{in}}) + D(u_{2,0}^{\rm{in}}), S(u_{1,0}^{\rm{out}})\} $. The corresponding flux values on the incoming roads are always smaller or equal to the actual demand. The algorithm for the corresponding adjusted flux values is listed in algorithm 4.

    Table Algorithm 4.  2–to–1 junction.
    Require: Demand $ D_1 $, $ D_2 $, Supply $ S_1 $, right of way parameter $ q_1, q_2 $, flux $ f $, discontinuity $ u^* $, jump magnitude $ \alpha $
    Ensure: Flux values $ f_1^{\rm{in}}, f_2^{\rm{in}}, f_1^{\rm{out}} $
    1: $ f_{\rm{max}} = \min \bigl\{D_1+D_2, S_1 \bigr \} $
    2: $ f_1^{\rm{out}}= f_{\rm{max}}, z_1= q_1 f_{\rm{max}}, z_2=q_2 f_{\rm{max}} $
    3: if $ z_1>D_1 $ then
    4:   $ z_1=D_1, z_2=f_{\rm{max}}-z_2 $ 5: else if $ z_2>D_2 $ then
    6:   $ z_2=D_2, z_1=f_{\rm{max}}-z_1 $
    7: end if
    8: $ f_1^{\rm{in}}=z_1, f_2^{\rm{in}}=z_2 $
    9: if $ D_1+D_2 = S_1 $ and $ S_1 \neq f(u^*-) $ then
    10:    $ f_1^{\rm{out}} = f_1^{\rm{out}} + \alpha $
    11:else if $ D_1+D_2<S_1 $ then
    12:    if $ S_1 \neq f(u^*-) $ then
    13:      $ f_1^{\rm{out}} = f_1^{\rm{out}} + \alpha $
    14:    end if
    15:    if $ f_1^{\rm{in}}<D_1 $ then
    16:      if $ f_1^{\rm{in}}>f(u^*+) $ then
    17:       $ f_1^{\rm{in}}=f(u^*-) $
    18:      else
    19:       $ f_1^{\rm{in}}=f_1^{\rm{in}}+\alpha $
    20:     end if
    21:    end if
    22:    if $ f_2^{\rm{in}}<D_2 $ then
    23:      if $ f_2^{\rm{in}}>f(u^*+) $ then
    24:       $ f_2^{\rm{in}}=f(u^*-) $
    25:      else
    26:       $ f_2^{\rm{in}}=f_2^{\rm{in}}+\alpha $
    27:      end if
    28:    end if
    29: end if

     | Show Table
    DownLoad: CSV

    As before, no adjustment is needed when the demand on both roads is smaller than the supply. If the supply is active we might need to adjust the outgoing road and in most cases at least one incoming road. As in the 1–to–2 case, $ f_i^{{\rm{in}}} $ can be greater than $ f(u^*+) $ such that the Riemann solver of theorem 3.8 gives $ u^* $ as a solution with the corresponding flux $ p(u^*) = f(u^*-) $ leading to an adjustment smaller than $ \alpha $.

    In this section, we present some numerical examples to compare the splitting algorithm on networks with the Riemann solver. Further, we compute the solution using a regularized flux and the Godunov scheme. We consider the following flux function:

    $ f(u)={u,u[0,0.5),0.5(1u),u[0.5,1].
    $
    (5.1)

    The corresponding regularization is given by Eq (2.7).

    We consider in particular the 1–to–2 and 2–to–1 situations. For our comparison, we choose constant initial data on each road. The junction is always located at $ x = 0 $, such that the incoming roads have negative $ x $-coordinates and the outgoing ones positive. All roads are assumed to have the same length of 2. Furthermore, we compute the sum of the $ L^1 $ errors on each road at $ t = 1 $ to compare the numerical approaches. As the splitting algorithm requires $ \lambda\leq 1 $ as CFL condition, but the regularization approach $ \lambda\leq \epsilon $, we compare different CFL constants and regularization parameters with each other.

    In both scenarios the supply of the first outgoing road is restrictive. The parameter settings are as follows:

    $ (1) $ $ u_{1,0}^{\rm{in}} = 0.4 $, $ u_{1,0}^{\rm{out}} = 0.9 $, $ u_{2,0}^{\rm{out}} = 0.7 $, $ \beta = (0.75, 0.25) $.

    The exact solution is given by

    $ uin1(x,t)={0.4,x32t,0.5,32t<x<12t,1315,12t<x<0,uout1(x,t)=0.9,uout2(x,t)={115β20x841t,0.7,841t<x.
    $

    Apparently, the solution induces two waves on the incoming road.

    $ (2) $ $ u_{1,0}^{\rm{in}} = 0.4 $, $ u_{1,0}^{\rm{out}} = 0.7 $, $ u_{2,0}^{\rm{out}} = 0.2 $, $ \beta = (0.5, 0.5) $.

    The exact solution is given by

    $ uin1(x,t)={0.4,xt,0.5,t<x<0,uout1(x,t)=0.7,1x3,uout2(x,t)={0.3β21xt,0.2,t<x.
    $

    This example generates $ f_1^{{\rm{in}}} = 0.3 $ such that we are in the specific case that the flux has to be adjusted, but only to $ 0.5 $.

    In Table 1, the $ L^1 $-errors are compared for different $ \Delta x $ and the numerical convergence rate (CR) is determined by a least square fitting method. For the splitting algorithm we are allowed to choose a rather large CFL constant of 0.75 while for a direct comparison to the regularized approach $ \lambda = 0.1 $. is chosen. Furthermore, the regularization approach is also computed with a smaller regularization parameter. For the regularization approach we choose $ \lambda = \epsilon $.

    Table 1.  The table shows the $ L^1 $ error for the splitting algorithm (Splitt.) and the regularized problem (Reg.) for two different $ \epsilon $ values for the 1–to–2 examples.
    Example 1
    Splitt. Reg.
    $ \Delta x $ $ \lambda=0.75 $ $ \lambda=0.1 $ $ \epsilon = 0.1 $ $ \epsilon = 0.01 $
    $ 0.04 $ 33.44e-03 46.77e-03 82.26e-03 51.82e-03
    $ 0.02 $ 24.17e-03 29.05e-03 65.59e-03 30.69e-03
    $ 0.01 $ 14.16e-03 20.12e-03 60.86e-03 24.45e-03
    $ 0.005 $ 8.97e-03 12.49e-03 58.37e-03 20.44e-03
    CR 0.64695 0.62453 0.1593 0.4353
    Example 2
    Splitt. Reg.
    $ \Delta x $ $ \lambda=0.75 $ $ \lambda=0.1 $ $ \epsilon = 0.1 $ $ \epsilon = 0.01 $
    $ 0.04 $ 4.58e-03 7.41e-03 44.70e-03 14.57e-03
    $ 0.02 $ 2.97e-03 4.24e-03 43.47e-03 11.28e-03
    $ 0.01 $ 2.03e-03 2.89e-03 42.50e-03 10.06e-03
    $ 0.005 $ 1.24e-03 1.99e-03 41.80e-03 9.21e-03
    CR 0.61911 0.62327 0.0322 0.2150

     | Show Table
    DownLoad: CSV

    We can see that the error terms obtained by the splitting algorithm are the lowest and so the computational costs. Obviously, the error terms increase with a lower CFL due to numerical diffusion. For a direct comparison with the regularized approach, the CFL condition should be the same. Meaning that the second column in Table 1 for the splitting algorithm should be compared with the first one of the regularized approach. In this case, we can see that the splitting algorithm also performs better in both examples. By choosing a smaller regularization parameter the performance of the regularized approach increases, but also the computational costs. To obtain similar error terms as for the splitting algorithm the regularization parameter needs to be further reduced at very high computational costs.

    For $ \lambda = 0.1 $ and $ \Delta x = 0.01 $ the approximate solutions and the exact solution are displayed in figure 6 for the first example. On the incoming road the splitting algorithm approximates the solution much better, while for the regularized approach an even smaller regularization parameter would be needed for a correct approximation. On the outgoing roads both schemes capture the right waves. Note that we zoomed in on the second outgoing road to make the small differences visible.

    Figure 6.  First example with $ \lambda = 0.1 $ and $ \Delta x = 0.01 $.

    Here, we consider two scenarios, where in the first scenario the demand is restrictive while in the second one the supply. The parameter settings are as follows:

    $ (1) $ $ u_{1,0}^{\rm{in}} = 0.2 $, $ u_{2,0}^{\rm{in}} = 0.25 $, $ u_{1,0}^{\rm{out}} = 0.3 $, $ q = 0.75 $.

    The exact solution is given by

    $ uin1(x,t)=0.2,uout1(x,t)={0.450xt,0.3,t<xuin2(x,t)=0.25.
    $

    $ (2) $ $ u_{1,0}^{\rm{in}} = 0.6 $, $ u_{2,0}^{\rm{in}} = 0.7 $, $ u_{1,0}^{\rm{out}} = 0.4 $, $ q = 0.8 $.

    The exact solution is given by

    $ uin1(x,t)={0.5x2t,0.6,2t<x<0,uout1(x,t)={0.50xt,0.4,t<x,uin2(x,t)={0.8x0.5t,0.7,0.5t<x<0.
    $

    In particular, the flux value for the first incoming road needs to be adjusted from $ 0.4 $ to $ 0.5 $ in the splitting algorithm, while on the second road it is adjusted as usual from $ 0.1 $ to $ 0.35 = 0.1+\alpha $. Note that due to the high velocity on the first incoming road we evaluate the error at $ t = 0.5 $.

    Considering the $ L^1 $ errors in Table 2 similar observations as in the 1–to–2 case can be made. The splitting algorithm performs best and has the lowest computational costs. In the first example the solution of the splitting algorithm coincides with the regularized approach as the flows are equal at the junction due to the active demand, see second column of the splitting algorithm and first one of the regularized approach in Table 2. Nevertheless, in the second example the splitting algorithm performs much better for the same CFL constant. In addition, the regularization parameter needs to be even smaller than 0.01 to reach the precision of the faster splitting algorithm.

    Table 2.  The table shows the $ L^1 $ error for the splitting algorithm (Splitt.) and the regularized problem (Reg.) for two different $ \epsilon $ values for the 2–to–1 examples.
    Example 1
    Splitt. Reg.
    $ \Delta x $ $ \lambda=0.75 $ $ \lambda=0.1 $ $ \epsilon = 0.1 $ $ \epsilon = 0.01 $
    $ 0.04 $ 9.25e-03 16.22e-03 16.22e-03 17.01e-03
    $ 0.02 $ 5.90e-03 11.63e-03 11.63e-03 12.19e-03
    $ 0.01 $ 2.98e-03 8.13e-03 8.13e-03 8.52e-03
    $ 0.005 $ 8.97e-03 5.71e-03 5.71e-03 5.99e-03
    CR 0.53838 0.50353 0.50353 0.50347
    Example 2
    Splitt. Reg.
    $ \Delta x $ $ \lambda=0.75 $ $ \lambda=0.1 $ $ \epsilon = 0.1 $ $ \epsilon = 0.01 $
    $ 0.04 $ 14.12e-03 20.10e-03 85.69e-03 27.55e-03
    $ 0.02 $ 9.65e-03 13.86e-03 79.98e-03 21.65e-03
    $ 0.01 $ 6.41e-03 9.57e-03 75.96e-03 17.49e-03
    $ 0.005 $ 4.51e-03 6.69e-03 73.22e-03 14.65e-03
    CR 0.55295 0.52959 0.07551 0.30432

     | Show Table
    DownLoad: CSV

    For $ \lambda = 0.1 $ and $ \Delta x = 0.01 $ the approximate solutions and the exact solution of the second example are displayed in Figure 7. In particular, on the first incoming road the splitting algorithm approximates the solution much better. For the considered resolution the difference between the exact solution and its approximation by the splitting algorithm is not visible on the first road. On the other roads the approximate solutions are very similar. Here, we zoomed in to display the differences to the exact solution.

    Figure 7.  Second example with $ \lambda = 0.1 $ and $ \Delta x = 0.01 $.

    We have presented a Riemann solver at a junction for conservation laws with discontinuous flux. We have adapted the splitting algorithm of [22] to networks and demonstrated its validity in comparison with the exact solution. We have also pointed out that the splitting algorithm on networks is faster and more accurate than the approach with a regularized flux. Future research includes the investigation of other network models, where the flux is discontinuous.

    J.F. was supported by the German Research Foundation (DFG) under grant HE 5386/18-1 and S.G. under grant GO 1920/10-1.

    The authors declare there is no conflict of interest.

    [1] American Cancer Society. Cancer Facts & Figures 2015. Atlanta: American Cancer Society; 2015.
    [2] Rousseau A, Badoual C (2011) Squamous cell carcinoma: an overvie Atlas Genet Cytogenet Oncol Haematol. Head and Neck, in press.
    [3] Schantz SP, Harrison LB, Forastiere A (2001) Tumors of the nasal cavity and paranasal sinuses, nasopharynx, oral cavity, and oropharynx. In: DeVita VT, Hellman SA, Rosenberg SA, eds. Cancer: principles and practice of oncology. 6th ed. Philadelphia: Lippincott Williams & Wilkins: 797-860.
    [4] Rubin Grandis J, Melhem MF, Gooding WE, et al. (1988) Levels of TGF-alpha and EGFR protein in head and neck squamous cell carcinoma and patient survival. J Natl Cancer Inst 90: 824-832.
    [5] Chung CH, Ely K, McGavran L, et al. (2006) Increased epidermal growth factor receptor gene copy number is associated with poor prognosis in head and neck squamous cell carcinomas. J Clin Oncol 24: 4170-4176.
    [6] Temam S, Kawaguchi H, El-Naggar AK, et al. (2007) Epidermal growth factor receptor copy number alterations correlate with poor clinical outcome in patients with head and neck squamous cancer. J Clin Oncol 25: 2164-2170. doi: 10.1200/JCO.2006.06.6605
    [7] Bonner JA, Harari PM, Giralt J (2006) Radiotherapy plus cetuximab for squamous-cell carcinoma of the head and neck. N Engl J Med 354: 567-578. doi: 10.1056/NEJMoa053422
    [8] Goldstein NI, Prewett M, Zuklys K, et al. (1995) Biological efficacy of a chimeric antibody to the epidermal growth factor receptor in a human tumor xenograft model. Clin Cancer Res 1: 1311-1318.
    [9] Li S, Schmitz KR, Jeffrey PD, et al. (2005) Structural basis for inhibition of the epidermal growth factor receptor by cetuximab. Cancer Cell 7: 301-311. doi: 10.1016/j.ccr.2005.03.003
    [10] Sato JD, Kawamoto T, Le AD, et al. (1983) Biological effects in vitro of monoclonal antibodies to human epidermal growth factor receptors. Mol Biol Med 1: 511-529.
    [11] Kang X, Patel D, Ng S, et al. (2007) High affinity Fc receptor binding and potent induction of antibody-dependent cellular cytotoxicity (ADCC) in vitro by anti-epidermal growth factor receptor antibody cetuximab. J Clin Oncol 25: 128s.
    [12] Kimura H, Sakai K, Arao T, et al. (2007) Antibody-dependent cellular cytotoxicity of cetuximab against tumor cells with wild-type or mutant epidermal growth factor receptor. Cancer Sci 98: 1275-1280.
    [13] Zhang W, Gordon M, Schultheis AM, et al. (2007) FCGR2A and FCGR3A polymorphisms associated with clinical outcome of epidermal growth factor receptor expressing metastatic colorectal cancer patients treated with single-agent cetuximab. J Clin Oncol 25: 3712-3718. doi: 10.1200/JCO.2006.08.8021
    [14] Fan Z, Baselga J, Masui H, et al. (1993) Antitumor effect of anti-epidermal growth factor receptor monoclonal antibodies plus cis-diamminedichloroplatinum on well established A431 cell xenografts. Cancer Res 53: 4637-4642.
    [15] Burtness B, Goldwasser MA, Flood W, et al. (2005) Phase III randomized trial of cisplatin plus placebo compared with cisplatin plus cetuximab in metastatic/recurrent head and neck cancer: an Eastern Cooperative Oncology Group study. J Clin Oncol 23: 8646-8654 doi: 10.1200/JCO.2005.02.4646
    [16] Thomas KH, Patrick JS (2013) Antigen-specific immunotherapy in head and neck cancer. Adv Cell Mol Otolaryngol 1.
    [17] Ira Mellman, George Coukos, Glenn Dranoff (2011) Cancer immunotherapy comes of age. Nature 480: 480-489.
    [18] Andrew MS, James PA, Jedd DW (2012) Monoclonal antibodies in cancer therapy. Cancer Immun 12: 14.
    [19] Brahmer JR, Drake CG, Wollner I, et al. (2010) Phase I study of single-agent anti-programmed death-1 (MDX-1106) in refractory solid tumors: safety, clinical activity, pharmacodynamics and immunologic correlates. J Clin Oncol 28: 3167-3175. doi: 10.1200/JCO.2009.26.7609
    [20] 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
    [21] Suntharalingam, G, Perry MR, Ward S, et al. (2006) Cytokine storm in a phase 1 trial of the anti-CD28 monoclonal antibody TGN1412. N Engl J Med 355: 1018-1028. doi: 10.1056/NEJMoa063842
    [22] Brahmer JR, Drake CG, Wollner I, et al. (2010). Phase I study of single-agent anti-programmed death-1 (MDX-1106) in refractory solid tumors: safety, clinical activity, pharmacodynamics and immunologic correlates. J Clin Oncol 28: 3167-3175. doi: 10.1200/JCO.2009.26.7609
    [23] Keir ME, Butte MJ, Freeman GJ, et al. (2008) PD-1 and its ligands in tolerance and immunity. Annu Rev Immunol 26: 677-704 doi: 10.1146/annurev.immunol.26.021607.090331
    [24] Parsa AT, Waldron JS, Panner A, et al. (2006) Loss of tumor suppressor PTEN function increases B7-H1 expression and immunoresistance in glioma. Nature Med 13: 84-88.
    [25] Gadiot J, Hooijkaas AI, Kaiser AD, et al. (2011) Overall survival and PD-L1 expression in metastasized malignant melanoma. Cancer 117: 2192-2201. doi: 10.1002/cncr.25747
    [26] Gao Q, Wang XY, Qiu SJ, et al. (2009) Overexpression of PD-L1 significantly associates with tumor aggressiveness and postoperative recurrence in human hepatocellular carcinoma. Clin Cancer Res 15: 971-979.
    [27] Leach DR, Krummel MF, Allison JP (1996) Enhancement of antitumor immunity by CTLA-4. Science 271: 1734-1736. doi: 10.1126/science.271.5256.1734
    [28] Van Cutsem E, Köhne CH, Hitre E, et al. (2009) Cetuximab and chemotherapy as initial treatment for metastatic colorectal cancer. N Engl J Med 360: 1408-1417.
    [29] Schliemann C, Neri D (2010) Antibody-based vascular tumor targeting. Cancer Res 180: 201-216.
    [30] Pillay V, Gan HK, Scott AM (2011) Antibodies in oncology. N Biotechnol 28: 518-529. doi: 10.1016/j.nbt.2011.03.021
    [31] Divgi CR, Welt S, Kris M, et al. (1991) Phase I and imaging trial of indium 111-labeled anti-epidermal growth factor receptor monoclonal antibody 225 in patients with squamous cell lung carcinoma. J Natl Cancer Inst 83: 97-104. doi: 10.1093/jnci/83.2.97
    [32] Ellis LM, Hicklin DJ (2008) VEGF-targeted therapy: mechanisms of anti-tumor activity. Nat Rev Cancer 8: 579-591. doi: 10.1038/nrc2403
    [33] Friedman HS, Prados MD, Wen PY, et al. (2009) Bevacizumab alone and in combination with irinotecan in recurrent glioblastoma. J Clin Oncol 27: 4733-4740. doi: 10.1200/JCO.2008.19.8721
    [34] Heiss MM, Murawa P, Koralewski P, et al. (2010) The trifunctional antibody catumaxomab for the treatment of malignant ascites due to epithelial cancer: results of a prospective randomized phase II/III trial. Int J Cancer 127: 2209-2221. doi: 10.1002/ijc.25423
    [35] Boland WK, Bebb G. (2009) Nimotuzumab: a novel anti-EGFR monoclonal antibody that retains anti-EGFR activity while minimizing skin toxicity. Expert Opin Biol Ther 9:1199-1206.
    [36] Curran D, Giralt J, Harari PM, et al. (2007) Quality of life in head and neck cancer patients after treatment with high-dose radiotherapy alone or in combination with cetuximab. J Clin Oncol 25: 2191-2197. doi: 10.1200/JCO.2006.08.8005
    [37] Bonner JA, Harari PM, Giralt J, et al. (2010) Radiotherapy plus cetuximab for locoregionally advanced head and neck cancer: 5-year survival data from a phase 3 randomized trial, and relation between cetuximab-induced rash and survival. Lancet Oncol 11: 21-28. doi: 10.1016/S1470-2045(09)70311-0
    [38] Koutcher L, Sherman E, Fury M, et al. (2011) Concurrent cisplatin and radiation versus cetuximab and radiation for locally advanced head-and-neck cancer. Int J Radiat Oncol Biol Phys 81: 915-922. doi: 10.1016/j.ijrobp.2010.07.008
    [39] Chew A, Hay J, Laskin JJ, et al. (2011) Toxicity in combined modality treatment of HNSCC: Cisplatin versus cetuximab. J Clin Oncol 29: 5526.
    [40] Shapiro LQ, Sherman EJ, Koutcher L, et al. (2012) Efficacy of concurrent cetuximab (C225) versus 5-fluorouracil/carboplatin (5FU/CBDCA) or cisplatin (CDDP) with intensity-modulated radiation therapy (IMRT) for locally advanced head and neck cancer (LAHNSCC). J Clin Oncol 30: 5537.
    [41] Pryor DI, Porceddu SV, Burmeister BH, et al. (2009) Enhanced toxicity with concurrent cetuximab and radiotherapy in head and neck cancer. Radiother Oncol 90: 172-176. doi: 10.1016/j.radonc.2008.09.018
    [42] Vermorken JB, Mesia R, Rivera F, et al. (2008) Platinum-based chemotherapy plus cetuximab in head and neck cancer. N Engl J Med 359: 1116-1127. doi: 10.1056/NEJMoa0802656
    [43] Hitt R, Irigoyen A, Cortes-Funes H, et al. (2012). Spanish Head and Neck Cancer Cooperative Group (TTCC) Phase II study of the combination of cetuximab and weekly paclitaxel in the first-line treatment of patients with recurrent and/or metastatic squamous cell carcinoma of head and neck. Ann Oncol 23: 1016-1022. doi: 10.1093/annonc/mdr367
    [44] Ang KK, Zhang QE, Rosenthal DI, et al. (2011) A randomized phase III trial (RTOG 0522) of concurrent accelerated radiation plus cisplatin with or without cetuximab for stage III-IV head and neck squamous cell carcinomas (HNC). J Clin Oncol 30: 360.
    [45] Ley J, Mehan P, Wildes TM, et al. (2012) Concurrent cisplatin vs. cetuximab with definitive radiation therapy (RT) for head and neck squamous cell carcinoma (HNSCC): A retrospective comparison. Multidisciplinary Head and Neck Cancer Symposium (Phoenix, AZ)
    [46] Balz V, Scheckenbach K, Gotte K, et al. (2003) Is the p53 inactivation frequency in squamous cell carcinomas of the head and neck underestimated? Analysis of p53 exons 2-11 and human papillomavirus 16/18 E6 transcripts in 123 unselected tumor specimens. Cancer Res 63: 1188-1191.
    [47] Deleo AB (1998) p53-based immunotherapy of cancer. Crit Rev Immunol 18: 29-35. doi: 10.1615/CritRevImmunol.v18.i1-2.40
    [48] Hoffmann TK, Nakano K, Elder EM, et al. (2000) Generation of T cells specific for the wild-type sequence p53(264-272) peptide in cancer patients: Implications for immunoselection of epitope loss variants. J Immunol 165: 5938-5944. doi: 10.4049/jimmunol.165.10.5938
    [49] Andrade FPA, Ito D, Deleo AB, et al. (2010) CD8+ T cell recognition of polymorphic wild-type sequence p53(65-73) peptides in squamous cell carcinoma of the head and neck. Cancer Immunol Immunother 59: 1561-1568. doi: 10.1007/s00262-010-0886-1
    [50] Chikamatsu K, Albers A, Stanson J, et al. (2003) P53(110-124)-specific human CD4+ T-helper cells enhance in vitro generation and antitumor function of tumor-reactive CD8+ T cells. Cancer Res 63: 3675-3681.
    [51] Albers AE, Ferris RL, Kim GG, et al. (2005) Immune responses to p53 in patients with cancer: Enrichment in tetramer+p53 peptide-specific T cells and regulatory T cells at tumor sites. Cancer Immunol Immunother 54: 1072-1081. doi: 10.1007/s00262-005-0670-9
    [52] Zhang Y, Sturgis EM, Huang Z, et al. (2012) Genetic variants of the p53 and p73 genes jointly increase risk of second primary malignancies in patients after index squamous cell carcinoma of the head and neck. Cancer 118: 485-492. doi: 10.1002/cncr.26222
    [53] Clayman GL, El-Naggar AK, Lippman SM, et al. (1998) Adenovirus-mediated p53 gene transfer in patients with advanced recurrent head and neck squamous cell carcinoma. J Clin Oncol 16: 2221-2232.
    [54] Heusinkveld M, Goedemans R, Briet RJ, et al. (2012) Systemic and local human papillomavirus 16-specific T-cell immunity in patients with head and neck cancer. Int J Cancer 131: E74-85. doi: 10.1002/ijc.26497
    [55] Wansom D, Light E, Worden F, et al. (2010) Correlation of cellular immunity with human papillomavirus 16 status and outcome in patients with advanced oropharyngeal cancer. Arch Otolaryngol Head Neck Surg 136: 1267-1273. doi: 10.1001/archoto.2010.211
    [56] The Cancer Genome Atlas Network, 2015. Comprehensive Genomic Characterization of Head and Neck Squamous Cell Carcinomas. Nature 517: 576-582.
    [57] Wei G, John ZHL, Jimmy YWC, et al. (2012) mTOR pathway and mTOR inhibitors in head and neck cancer department of surgery. Otolaryngology.
    [58] Hay N, Sonenberg N (2004) Upstream and downstream of mTOR. Genes Dev 18: 1926-1945. doi: 10.1101/gad.1212704
    [59] Lionello SM, Loreggian BL (2012) High mTOR expression is associated with a worse oncological outcome in laryngeal carcinoma treated with postoperative radiotherapy: a pilot study. J Oral Pathol Med 41: 136-140. doi: 10.1111/j.1600-0714.2011.01083.x
  • 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(6026) PDF downloads(1250) Cited by(2)

Figures and Tables

Tables(2)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog