Comparative analysis of copper and zinc based agrichemical biocide products: materials characteristics, phytotoxicity and in vitro antimicrobial efficacy

Parthiban Rajasekaran; Harikishan Kannan; Smruti Das; Mikaeel Young; Swadeshmukul Santra; Parthiban Rajasekaran; Harikishan Kannan; Smruti Das; Mikaeel Young; Swadeshmukul Santra

doi:10.3934/environsci.2016.3.439

AIMS Environmental Science

2016, Volume 3, Issue 3: 439-455. doi: 10.3934/environsci.2016.3.439

Previous Article Next Article

Research article

Comparative analysis of copper and zinc based agrichemical biocide products: materials characteristics, phytotoxicity and in vitro antimicrobial efficacy

1.
NanoScience Technology Center, University of Central Florida, 12424 Research Parkway, Suite 400, Orlando, FL 32826, United States
2.
Department of Chemistry, University of Central Florida, 12424 Research Parkway, Suite 400, Orlando, FL 32826, United States
3.
Department of Material Science and Engineering, University of Central Florida, 12424 Research Parkway, Suite 400, Orlando, FL 32826, United States
4.
Burnett School of Biomedical Sciences, University of Central Florida, 12424 Research Parkway, Suite 400, Orlando, FL 32826, United States.

Received: 18 May 2016 Accepted: 20 July 2016 Published: 26 July 2016

In the past few decades, copper based biocides have been extensively used in food crop protection including citrus, small fruits and in all garden vegetable production facilities. Continuous and rampant use of copper based biocides over decades has led to accumulation of this metal in the soil and the surrounding ecosystem. Toxic levels of copper and its derivatives in both the soil and in the run off pose serious environmental and public health concerns. Alternatives to copper are in great need for the agriculture industry to produce food crops with minimal environmental risks. A combination of copper and zinc metal containing biocide such as Nordox 30/30 or an improved version of zinc-only containing biocide would be a good alternative to copper-only products if the efficacy can be maintained. As of yet there is no published literature on the comparative study of the materials characteristics and phyto-compatibility properties of copper and zinc-based commercial products that would allow us to evaluate the advantages and disadvantages of both versions of pesticides. In this report, we compared copper hydroxide and zinc oxide based commercially available biocides along with suitable control materials to assess their efficacy as biocides. We present a detailed material characterization of the biocides including morphological studies involving electron microscopy, molecular structure studies involving X-ray diffraction, phytotoxicity studies in model plant (tomato) and antimicrobial studies involving surrogate plant pathogens (Xanthomonas alfalfae subsp. citrumelonis, Pseudomonas syringae pv. syringae and Clavibacter michiganensis subsp. michiganensis). Zinc based compounds were found to possess comparable to superior antimicrobial properties while exhibiting significantly lower phytotoxicity when compared to copper based products thus suggesting their potential as an alternative.

Keywords:

Citation: Parthiban Rajasekaran, Harikishan Kannan, Smruti Das, Mikaeel Young, Swadeshmukul Santra. Comparative analysis of copper and zinc based agrichemical biocide products: materials characteristics, phytotoxicity and in vitro antimicrobial efficacy[J]. AIMS Environmental Science, 2016, 3(3): 439-455. doi: 10.3934/environsci.2016.3.439

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

Abstract

1. Introduction

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.

2. Kinetic and macroscopic models for chemotaxis

2.1. Kinetic equation

Denote with $f=f(x,t,v)$ the cell distribution at position $x\in\mathbb{R}$ , moving with velocity $v\in V\equiv [-1,\space 1]$ at the time $t\in [0,\space T]$ , $\rho(x,t)=\int_Vf(x,t,v)dv$ is the macroscopic density of cells, $m=m(x,t)$ is the density of the chemoattractant. Let $\alpha \in \mathbb{R}^+$ denote the rate of turning into the favourable direction given by the chemoattractant, let $\lambda \in \mathbb{R}^+$ be the turning rate of the cells and $D, \gamma_\rho, \gamma_m \in$ $\mathbb{R}^+$ be the diffusivity, production and decay rate of the chemoattractant, respectively. As in [6], we consider the kinetic equation

$\left\{ \begin{array}{lcl} \partial_t f + v\partial_x f&=&-\lambda \left( f - \frac{\rho}{2} \right)+\frac{1}{2}\alpha v\overline{\partial_x m} \rho \\ \partial_t m -D(\partial_{xx}m)&=&\gamma_{\rho}\rho -\gamma_m m\ . \end{array} \right.$

(1)

Here, we have defined

$\overline{\partial_x m}=\frac{\partial_x m}{\sqrt{1+\vert\partial_xm\vert^2}}.$

We assume

$\label{condprev} \lambda\geq \alpha$

(2)

in order to guarantee positivity of the turning kernel

$\lambda+\alpha v\overline{\partial_x m}\geq 0$

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

$\left\{ \begin{array}{lcl} \partial _t f + \frac{1}{\epsilon} v\partial_x f &=& -\frac{\lambda}{\epsilon^2}\left( f - \frac{\rho}{2}\right)+\frac{1}{2\epsilon}\alpha v \overline{\partial_x m}\rho\\ \partial _t m - D(\partial_{xx}) m&=& \gamma_\rho \rho-\gamma_m m .\\ \end{array} \right.$

(3)

Condition (2) turns into $\lambda\geq \epsilon \alpha$ or

$\label{cond} \epsilon \le \frac{\lambda}{\alpha}.$

(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]

$\begin{align}\label{eq:evenodd} r(x,t,v)=\frac{1}{2}\left(f(x,t,v)+f(x,t,-v)\right) ,&& j(x,t,v)=\frac{1}{2\epsilon}\left(f(x,t,v)-f(x,t,-v)\right) \end{align}$

(5)

for positive velocities $(v\geq 0)$ . This leads to the system

$\partial_t r + v\partial_x j =-\frac{\lambda}{\epsilon^2}\left( r-\frac{\rho}{2}\right)$

(6)

$\begin{align} \partial_t j + \frac{1}{\epsilon^2}v\partial_x r&=-\frac{1}{\epsilon^2}\left(\lambda j -\frac{1}{2}\alpha v \overline{\partial_x m} \rho\right)\ . \end{align}$

(7)

When $\epsilon\to 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 $v \ge 0$ we obtain

$\begin{align} \label{cont} \partial_t \rho + 2 \partial_x \int_{0}^1 v j dv =0. \end{align}$

(8)

Moreover we observe $r= \frac{\rho}{2} + \mathcal{O}(\epsilon^2)$ and $j = \frac{\alpha \overline{\partial_x m} }{2 \lambda} v \rho - \frac{v}{\lambda}\partial_x r +\mathcal{O}(\epsilon^2)$ . Plugging the expansions into (8) gives the diffusion limit of the kinetic equation, the modified or flux-limited Keller-Segel equation

$\left\{ \begin{array}{lcl} \partial_t\rho - \partial_x \left(\frac{1}{3\lambda}\partial_x\rho - \frac{\alpha}{3\lambda}\overline{\partial_x m} \rho\right)&=&0\\ \partial_t m -D\partial_{xx}m&=&\gamma_{\rho}\rho -\gamma_m m.\\ \end{array} \right.$

(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.

2.2. Full-moment non-linear models

Consider the following averaged quantities

$\begin{align*} \rho(x,t)=\int_{-1}^1f(x,t,v)dv\ ,&& q(x,t)=\int_{-1}^1 vf(x,t,v)dv\ , \end{align*}$

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

$\begin{align}\label{equ:fm_realizability_con} &&\rho\geq 0,\ &&P\geq 0,\ &&\frac{q^2}{\rho^2}\leq\frac{P}{\rho}\leq 1.\ \end{align}$

(10)

This leads to the equations

$\left\{ \begin{array}{lcl} \partial_t\rho +\frac{1}{\epsilon}\partial_x q&=&0 \\[0.3 cm] \partial_t q +\frac{1}{\epsilon}\partial_xP& =&-\frac{1}{\epsilon^2}\left( \lambda q -\epsilon\frac{\alpha}{3}\overline{\partial_xm}\rho\right)\\[0.3 cm] \partial_tm - D\partial_{xx}m &=& \gamma_{\rho}\rho-\gamma_mm\ . \end{array} \right.$

(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 $\rho$ and $q$ to obtain expressions for $P$ . A linear closure with $f(x,t,v) = \frac{1}{2}\rho(x,t) + \frac{3}{2} vq(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) = a\exp (vb)\$

(12)

with $a \in \mathbb{R}^+$ and $b \in \mathbb{R}$ . This closure function is chosen using an entropy minimization principle, see for example [23]. This leads to

$P =\int_Vv^2fdv = \rho h(u), h(u)=\left(1-\frac{2}{b}u\right), u:=\frac{q}{\rho}=coth(b)-\frac1b\ ,$

where

$\begin{align*} &\lim\limits_{u\to 0} h(u)=\frac{1}{3},\ &u^2\leq h(u)\leq 1 \ \end{align*}$

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)=\frac{1}{3}+\frac{2}{3}u^2$ or the Levermore-Lorentz Eddington factor $h(u)=\frac{1}{3}+\frac{2u^2}{2+\sqrt{4-3u^2}}$ . The resulting full-moment macroscopic model for chemotaxis derived by the nonlinear closure (12) reads

$\left\{ \begin{array}{lcl} \partial_t\rho + \frac{1}{\epsilon}\partial_x q&=&0 \\[0.3 cm] \partial_t q +\frac{1}{\epsilon}\partial_x\left(\rho h\left(\frac{q}{\rho}\right)\right)&=&-\frac{1}{\epsilon}\left( \lambda \frac{1}{\epsilon} q - \frac{\alpha}{3}\overline{\partial_xm}\rho\right)\\[0.3 cm] \partial_tm - D\partial_{xx}m &=& \gamma_{\rho}\rho-\gamma_mm\ . \end{array} \right.$

(13)

Scaling and denoting $q_{\epsilon}:=\frac{q}{\epsilon}$ it can be rewritten as

$\left\{ \begin{array}{lcl} \partial_t\rho + \partial_x q_{\epsilon}&=&0 \\[0.3 cm] \partial_t q_{\epsilon} +\frac{1}{\epsilon^2}\partial_x\left(\rho h\left(\epsilon\frac{q_{\epsilon}}{\rho}\right)\right)&=&-\frac{1}{\epsilon^2}\left( \lambda q_{\epsilon} - \frac{\alpha}{3}\overline{\partial_xm}\rho\right)\\[0.3 cm] \partial_tm - D\partial_{xx}m &=& \gamma_{\rho}\rho-\gamma_mm\ . \end{array} \right.$

(14)

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

2.3. Half-moment non-linear models

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

$\begin{aligned} \rho ^- &=\int_{-1}^0 f(v)dv\ , & \rho^+ &=\int_{0}^1f(v)dv\ , \\ q ^- &=\int_{-1}^0 vf(v)dv \ , & q^+ &=\int_{0}^1vf(v)dv \ , \\ P ^- &=\int_{-1}^0 v^2f(v)dv \ , & P^+ &=\int_{0}^1v^2f(v)dv \ . \end{aligned}$

(15)

Using $f(x,t,v)\geq 0$ , one obtains the conditions $\rho^{\pm}\geq 0$ , $\pm q^{\pm}\geq 0$ , $P^{\pm}\geq 0$ and

$\left(\frac{q^{\pm}}{\rho^{\pm}}\right)^2\leq\frac{P^{\pm}}{\rho^{\pm}}\leq\frac{\pm q^{\pm}}{\rho^{\pm}}\leq 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

$\left\{ \begin{array}{lcl} \epsilon\partial_t\rho^{\pm} + \partial_x q^{\pm}&=&-\frac{1}{\epsilon}\lambda\left(\rho^{\pm}-\frac{\rho^+ +\rho^-}{2}\right)\pm \frac{1}{4}\alpha\overline{\partial_xm}\left(\rho^++\rho^-\right) \\ \epsilon\partial_tq^{\pm} + \partial_x P^{\pm}&=&-\frac{1}{\epsilon}\lambda\left(q^{\pm}\mp \frac{\rho^+ +\rho^-}{4}\right) + \frac{1}{6}\alpha\overline{\partial_xm}\left(\rho^++\rho^-\right) \\ \partial_tm - D\partial_{xx}m &=& \gamma_{\rho}\rho-\gamma_mm\ . \end{array} \right.$

(17)

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

$\begin{align}\label{eq:halfmoment_nonlinear_closure} f(v < 0) = a^- \exp{(vb^-)}\ , && f(v > 0) = a^+ \exp{(vb^+)}\ . \end{align}$

(18)

Plugging (18) in (15) leads to

$\begin{align*} \begin{aligned} &P^-=\rho^-h^-(u^-) \ , & &P^+ = \rho^+h^+(u^+) \ ,\\ &h^-(u^-) = u^-\left(-1-\frac{2}{b^-}\right)-\frac{1}{b^-}\ , & &h^+(u^+) = u^+\left(1-\frac{2}{b^+}\right)+\frac{1}{b^+}\ , \\ &u^-:=\frac{q^-}{\rho ^-}=\frac{1}{\exp{(b^-)}-1}-\frac{1}{(b^-)}\ , & &u^+:=\frac{q^+}{\rho ^+} = \frac{\exp{(b^+)}}{\exp{(b^+)}-1}-\frac{1}{(b^+)}\ , \end{aligned} \end{align*}$

such that

$\begin{align}\label{equ:eddington_half} \lim\limits_{u^\pm\to \pm\frac12}h^\pm(u^\pm) =\frac{1}{3} \ , && (u^{\pm})^2\leq h^{\pm}(u^{\pm})\leq \pm u^{\pm}\ . \end{align}$

(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]

$\begin{align}\label{equ:h_pm} h^-(u^-)=\frac23 (u^-)^2 - \frac13 u^- \ , && h^+(u^+)=\frac23 (u^+)^2 + \frac13 u^+\ . \end{align}$

(20)

Finally, the half-moment model for chemotaxis reads

$\left\{ \begin{array}{lcl} \partial_t\rho^\pm + \frac{1}{\epsilon}\partial_xq^\pm &=& -\frac{1}{\epsilon^2}\left[ \lambda\left(\rho^\pm -\frac{\rho^++\rho^-}{2}\right)\mp\frac{\epsilon\alpha}{4}\overline{\partial_x m}\left(\rho^++\rho^-\right) \right]\\ \partial_tq^\pm +\frac{1}{\epsilon} \partial_x\left(\rho^\pm h^\pm\left(\frac{q^\pm}{\rho^\pm}\right)\right) &=& -\frac{1}{\epsilon^2} \left[ \lambda\left(q^\pm \mp\frac{\rho^++\rho^-}{4}\right) - \frac{\epsilon\alpha}{6}\overline{\partial_x m}\left(\rho^++\rho^-\right)\right]\\ \partial_tm - D\partial_{xx}m &=& \gamma_{\rho}\rho-\gamma_mm\ . \end{array} \right.$

(21)

Moreover, by introducing new quantities

$\begin{align}\label{equ:new_cha_hat} \begin{aligned} \rho &= \rho^+ + \rho^-\ ,&&\hat{\rho}=\frac{\rho^+ - \rho^-}{\epsilon}\ ,&\\ q_{\epsilon} &=\frac{q^++q^-}{\epsilon}\ ,& &\hat{q}=q^+-q^-\ ,&\end{aligned} \end{align}$

(22)

the half-moment model can be rewritten as

$\left\{ \begin{array}{lcl} \partial_t\rho + \partial_xq_{\epsilon} &=& 0\\ \partial_tq_{\epsilon} + \frac{1}{\epsilon^2}\partial_xP(\rho,q_{\epsilon},\hat{\rho},\hat{q}) &=&-\frac{1}{\epsilon^2}\left(\lambda q_{\epsilon} - \frac{\alpha}{3}\overline{\partial_x m}\rho\right)\\ \partial_t\hat{\rho} + \frac{1}{\epsilon^2}\partial_x\hat{q} &=& -\frac{1}{\epsilon^2}\left(\lambda\hat{\rho} -\frac{1}{2}\alpha \overline{\partial_x m}\rho\right)\\ \partial_t\hat{q} + \partial_x\hat{P}(\rho,q_{\epsilon},\hat{\rho},\hat{q}) &=& -\frac{1}{\epsilon^2}\lambda\left(\hat{q}-\frac{\rho}{2}\right)\\ \partial_tm - D\partial_{xx}m &=& \gamma_{\rho}\rho-\gamma_mm\ \end{array} \right.$

(23)

with

$\begin{align}\label{equ:HM_M1_P} \begin{array}{lclcr} P(\rho,q_{\epsilon},\hat{\rho},\hat{q})&:=& P^+\left(\frac{q^+}{\rho^+}\right)&+&P^-\left(\frac{q^-}{\rho^-}\right)\\ &=& \frac{\rho+\epsilon\hat{\rho}}{2}h^+\left(\frac{\epsilon q_{\epsilon}+\hat{q}}{\rho+\epsilon\hat{\rho}}\right) &+& \frac{\rho-\epsilon\hat{\rho}}{2}h^-\left(\frac{\epsilon q_{\epsilon}-\hat{q}}{\rho-\epsilon\hat{\rho}}\right),\\ \hat{P}(\rho,q_{\epsilon},\hat{\rho},\hat{q})&:=& \frac{1}{\epsilon}P^+\left(\frac{q^+}{\rho^+}\right)&-&\frac{1}{\epsilon}P^-\left(\frac{q^-}{\rho^-}\right)\\ &=&\frac{\rho+\epsilon\hat{\rho}}{2\epsilon}h^+\left(\frac{\epsilon q_{\epsilon}+\hat{q}}{\rho+\epsilon\hat{\rho}}\right) &-& \frac{\rho-\epsilon\hat{\rho}}{2\epsilon}h^-\left(\frac{\epsilon q_{\epsilon}-\hat{q}}{\rho-\epsilon\hat{\rho}}\right).\\ \end{array} \end{align}$

(24)

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

3. Asymptotic preserving well-balanced (APWB) schemes for the nonlinear models on an interval

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.

3.1. Framework for APWB scheme

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

$\left\{ \begin{array}{lcl} \partial_t \rho +\partial_x q_{\epsilon} &=&0\\ \partial_t q_{\epsilon} +\frac{a}{\epsilon^2}\partial_x \rho &=& -\frac{1}{\epsilon^2}\left(\lambda q_{\epsilon}- c(\partial_xm)\rho\right)\\ \end{array} \right.$

(25)

with $c(x)=\frac{\alpha}{3}\frac{x}{\sqrt{1+x^2}}$ , $a,\lambda,\alpha>0$ , $0<\epsilon\leq\frac{\lambda}{\alpha}$ . The limit problem as $\epsilon \rightarrow 0$ is

$\partial_t \rho + \partial_x\left(\frac{c(\partial_xm)}{\lambda}\rho - \frac{a}{\lambda}\partial_x\rho\right)= 0.$

By introducing the diagonal variables $U = \frac12\rho + \frac{\epsilon}{2\sqrt{a}} q_{\epsilon}, V = \frac12\rho -\frac{\epsilon}{2\sqrt{a}} q_{\epsilon}$ , the system (25) is rewritten as

$\begin{align}\label{equ:exam_model_re} \begin{array}{lcl} \partial_tU +\frac{\sqrt{a}}{\epsilon}\partial_xU &=& \left(\frac{c(\partial_xm)}{2\sqrt{a}\epsilon}-\frac{\lambda}{2\epsilon^2}\right)U + \left(\frac{c(\partial_xm)}{2\sqrt{a}\epsilon}+\frac{\lambda}{2\epsilon^2}\right)V\\ \partial_tV -\frac{\sqrt{a}}{\epsilon}\partial_xV &=& -\left(\frac{c(\partial_xm)}{2\sqrt{a}\epsilon}-\frac{\lambda}{2\epsilon^2}\right)U - \left(\frac{c(\partial_xm)}{2\sqrt{a}\epsilon}+\frac{\lambda}{2\epsilon^2}\right)V. \end{array} \end{align}$

(26)

Using

$\left\{ \begin{array}{lcl} U^*(\Delta x) &=& V^*(\Delta x) +\frac{A+B}{1+B}(U^*(0)-V^*(\Delta x)) + \frac{A+B-2}{1+B}V^*(\Delta x)\\ V^*(0) &=& U^*(0) + \frac{A+B}{1+B}(V^*(\Delta x)-U^*(0)) - \frac{A+B-2}{1+B}V^*(\Delta x)\\ && A = \left(1-\frac{\lambda \Delta x}{\epsilon\sqrt{a}}\right)\exp\left(\frac{\Delta x }{a}c\left(\partial_x m\right)\right)\\ && B= \left(1+\frac{\lambda \Delta x}{\epsilon\sqrt{a}}\right)\exp\left(\frac{\Delta x }{a}c\left(\partial_x m\right)\right)\\ \end{array} \right.$

(27)

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

$\left\{ \begin{array}{lcl} U^{n+1}_i &=& U^n_i-l\left(U^{n+1}_i-U^*_{i-\frac12}\right) = U^n_i - l\left(U^{n+1}_i -V^{n+1}_i\right) \\&& - l\frac{A_{i-\frac12}^n+B_{i-\frac12}^n}{1+B_{i-\frac12}^n}\left(V_i^n-U^n_{i-1}\right)+l\frac{A_{i-\frac12}^n+B_{i-\frac12}^n-2}{1+B_{i-\frac12}^n}V_i^n\\ V^{n+1}_i &=& V^n_i+l\left(V^*_{i+\frac12}-V_i^{n+1}\right) = V^n_i + l\left(U^{n+1}_i -V^{n+1}_i\right)\\ &&+ l\frac{A_{i+\frac12}^n+B_{i+\frac12}^n}{1+B_{i+\frac12}^n}\left(V_{i+1}^n-U^n_{i}\right)-l\frac{A_{i+\frac12}^n+B_{i+\frac12}^n-2}{1+B_{i+\frac12}^n}V_{i+1}^n\\ \end{array} \right.$

(28)

with $l: = \frac{\sqrt{a}}{\epsilon}\frac{\Delta t}{\Delta x}$ and $\partial_xm^n_{i+\frac12}: = \frac{m^n_{i+1}-m^n_i}{\Delta x}$ . The APWB numerical solution (28) can be rewritten explicitly as

$\begin{align}\label{equ:ex_wb_exam_model} \left\{ \begin{array}{lcrcr} U^{n+1}_i &=& \frac{2l^2}{2l+1}\frac{1}{1+B^n_{i+\frac12}}V_{i+1}^n &+& \left(\frac{l+1}{2l+1}-\frac{l^2}{2l+1}\frac{A_{i+\frac12}^n+B_{i+\frac12}^n}{1+B_{i+\frac12}^n}\right)U^n_i\\ &+& \left(\frac{l}{2l+1}-\frac{2l(l+1)}{2l+1}\frac{1}{1+B_{i-\frac12}^n}\right)V_i^n &+& \frac{l(l+1)}{2l+1}\frac{A_{i-\frac12}^n+B_{i-\frac12}^n}{1+B_{i-\frac12}^n}U_{i-1}^n \\ V^{n+1}_i &=& \frac{l(l+1)}{2l+1}\frac{2}{1+B^n_{i+\frac12}}V_{i+1}^n &+& \left(\frac{l}{2l+1}-\frac{l(l+1)}{2l+1}\frac{A_{i+\frac12}^n+B_{i+\frac12}^n}{1+B_{i+\frac12}^n}\right)U^n_i\\ &+& \left(\frac{l+1}{2l+1}-\frac{2l^2}{2l+1}\frac{1}{1+B_{i-\frac12}^n}\right)V_i^n &+& \frac{l^2}{2l+1}\frac{A_{i-\frac12}^n+B_{i-\frac12}^n}{1+B_{i-\frac12}^n}U_{i-1}^n . \end{array} \right.\ \end{align}$

(29)

It is easy to check that for small enough values of $\epsilon\leq\sqrt{a}\frac{\Delta t}{\Delta x}\leq 1$ , the APWB scheme (29) is positivity preserving for $U$ and $V$ , if the following parabolic CFL restriction holds

$\Delta t \leq \frac{\lambda\Delta x^2}{4\left(a+\Delta x\right)}.$

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

3.2. APWB scheme for the nonlinear full-moment model

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

$\left\{ \begin{array}{lcl} \partial_t\rho + \partial_x z&=&0 \\ \partial_t z + \frac{a}{\epsilon^2}\partial_x\rho + \frac{\lambda}{\epsilon^2}z-\frac{\alpha}{3\epsilon^2}\overline{\partial_xm}\rho &=& \frac{1}{\beta}(q_{\epsilon}-z)\\ \partial_t q_{\epsilon} +\frac{1}{\epsilon^2}\partial_xw +\frac{\lambda}{\epsilon^2}q - \frac{\alpha}{3\epsilon^2}\overline{\partial_xm}\underbrace{\left(\frac{1}{h\left(\epsilon\frac{q_{\epsilon}}{\rho}\right)}w\right)}_{\rho}&=&0\\ \partial_t w +a\partial_x q &=& \frac{1}{\beta}\left(\rho h\left(\epsilon\frac{q_{\epsilon}}{\rho}\right)-w\right) \ . \end{array} \right.$

(30)

For $\beta\to 0$ the system converges to (14). Note that for any value of $\beta$ 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 $\beta=0$

$\begin{align}\label{equ:full_m1_relaxation} \rho^{(1)}=\rho^n, q_{\epsilon}^{(1)}= z^{(1)} = q_{\epsilon}^n, w^{(1)}=\rho^n h\left(\epsilon\frac{q_{\epsilon}^n}{\rho^n}\right)\ , \end{align}$

(31)

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

$\begin{align}\label{equ:new variable_m1} \begin{aligned} U=\frac12\rho+\frac{\epsilon}{2\sqrt{a}}z\ ,&V=\frac12\rho-\frac{\epsilon}{2\sqrt{a}}z\ , \\ \bar{U}=\frac{\epsilon}{2} q_{\epsilon}+\frac{1}{2\sqrt{a}}w\ ,&\bar{V}=-\frac{\epsilon}{2} q_{\epsilon}+\frac{1}{2\sqrt{a}}w\ , \end{aligned} \end{align}$

(32)

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

$\begin{align}\label{equ:telegrap_m1} \begin{aligned} \left\{ \begin{array}{lcl} \partial_tU +\frac{\sqrt{a}}{\epsilon}\partial_xU &=& \left(\frac{c_1}{2\sqrt{a}\epsilon}-\frac{\lambda}{2\epsilon^2}\right)U + \left(\frac{c_1}{2\sqrt{a}\epsilon}+\frac{\lambda}{2\epsilon^2}\right)V\\ \partial_tV -\frac{\sqrt{a}}{\epsilon}\partial_xV &=& -\left(\frac{c_1}{2\sqrt{a}\epsilon}-\frac{\lambda}{2\epsilon^2}\right)U - \left(\frac{c_1}{2\sqrt{a}\epsilon}+\frac{\lambda}{2\epsilon^2}\right)V\\ c_1(\partial_xm) &=& \frac{\alpha}{3}\overline{\partial_xm}\\ \end{array} \right. \\ \left\{ \begin{array}{lcl} \partial_t\bar{U} +\frac{\sqrt{a}}{\epsilon}\partial_x\bar{U} &=& \left(\frac{c_2}{2\sqrt{a}\epsilon}-\frac{\lambda}{2\epsilon^2}\right)\bar{U} + \left(\frac{c_2}{2\sqrt{a}\epsilon}+\frac{\lambda}{2\epsilon^2}\right)\bar{V}\\ \partial_t\bar{V} -\frac{\sqrt{a}}{\epsilon}\partial_x\bar{V} &=&- \left(\frac{c_2}{2\sqrt{a}\epsilon}-\frac{\lambda}{2\epsilon^2}\right)\bar{U} -\left(\frac{c_2}{2\sqrt{a}\epsilon}+\frac{\lambda}{2\epsilon^2}\right)\bar{V}\\ c_2(\rho,q,\partial_xm) &=& \frac{a}{h\left(\epsilon\frac{q_{\epsilon}}{\rho}\right)}c_1\ . \end{array} \right. \end{aligned} \end{align}$

(33)

Now, we can apply the APWB schemes presented in 3.1 for each subsystem in (33). Remark that for $a:=\max\Big\vert h\left(\epsilon\frac{q_{\epsilon}}{\rho}\right)\Big\vert$ the quantities in (32) are nonnegative and the scheme converges formally to the Keller-Segel equations as $\epsilon \rightarrow 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 $\rho$ remains positive. In the linear case, this is not an issue to be investigated, since the density $\rho$ in the Cattaneo model is not necessarily positive. In contrast, in the nonlinear case, it is an important open issue to be investigated.

3.3. APWB scheme for half-moment M1 model

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)

$\begin{align} \begin{aligned} &\left\{ \begin{array}{lcl} \partial_t\rho + \partial_xz &=& 0\\ \partial_t z +\frac{a}{\epsilon^2}\partial_x\rho + \frac{\lambda}{\epsilon^2}z-\frac{\alpha}{3\epsilon^2}\overline{\partial_xm}\rho &=& \frac{1}{\beta}(q_{\epsilon}-z)\\ \partial_tq_{\epsilon} +\frac{1}{\epsilon^2} \partial_x w +\frac{\lambda}{\epsilon^2}q_{\epsilon} -\frac{\alpha}{3\epsilon^2}\overline{\partial_x m}\underbrace{\left(\frac{2w-\epsilon (h^+-h^-)\hat{\rho}}{h^++h^-}\right)}_{\rho}&=& 0\\ \partial_tw + a\partial_x q_{\epsilon} &=&\frac{1}{\beta}\left(P-w\right)\\ \end{array} \right.\ \\ &\left\{ \begin{array}{lcl} \partial_t\hat{\rho} +\frac{1}{\epsilon^2}\partial_x\hat{z} +\frac{\lambda}{\epsilon^2}\hat{\rho} - \frac{\alpha}{2\epsilon^2}\overline{\partial_x m}\rho&=&0\\ \partial_t \hat{z}+ a\partial_x\hat{\rho} + \frac{\lambda}{\epsilon^2}\hat{z}&=& \frac{1}{\beta}\left(\hat{q}-\hat{z}\right)\\ \partial_t\hat{q} + \partial_x \hat{w} +\frac{\lambda}{\epsilon^2}\left(\hat{q}-\frac{\rho}{2}\right)&=&0\\ \partial_t\hat{w} +\frac{a}{\epsilon^2}\partial_x\hat{q} &=&\frac{1}{\beta}\left(\hat{P}-\hat{w}\right) \end{array} \right.\end{aligned} \end{align}$

(34)

with $h^{\pm}:= h^{\pm}\left(\frac{\epsilon q_{\epsilon}\pm\hat{q}}{\rho\pm\epsilon\hat{\rho}}\right)$ and $P=P\left(\rho,q_{\epsilon},\hat{\rho},\hat{q}\right),\hat{P} = {P}\left(\rho,q_{\epsilon},\hat{\rho},\hat{q}\right)$ defined in section 2.3. To obtain the correct asymptotic behaviour using a splitting scheme, we introduce a parameter $\sigma < \lambda$ , and rewrite the system (34) as

$\begin{align}\label{equ:relaxed_m1_half_splitting} \begin{aligned} &\left\{ \begin{array}{lcl} \partial_t\rho + \partial_xz &=& 0\\ \partial_t z +\frac{a}{\epsilon^2}\partial_x\rho + \frac{\lambda}{\epsilon^2}z-\frac{\alpha}{3\epsilon^2}\overline{\partial_xm}\rho &=& \frac{1}{\beta}(q_{\epsilon}-z)\\ \partial_tq_{\epsilon} +\frac{1}{\epsilon^2} \partial_x w +\frac{\lambda-\sigma}{\epsilon^2}q_{\epsilon} -\frac{2 w \alpha \overline{\partial_x m}}{3\epsilon^2(h^+ + h^-)} &=& -\frac{\sigma}{\epsilon^2}q_{\epsilon} -\frac{\alpha\overline{\partial_x m}(h^+ - h^-)}{3\epsilon(h^+ + h^-)}\hat{\rho}\\ \partial_tw +a\partial_x q_{\epsilon} &=&\frac{1}{\beta}(P-w)\\ \end{array} \right.\ \\ &\left\{ \begin{array}{lcl} \partial_t\hat{\rho} +\frac{1}{\epsilon^2}\partial_x\hat{z} &=&-\frac{\lambda}{\epsilon^2}\hat{\rho} + \frac{\alpha}{2\epsilon^2}\overline{\partial_x m}\rho\\ \partial_t \hat{z}+ a\partial_x\hat{\rho} + \frac{\lambda}{\epsilon^2}\hat{z}&=& \frac{1}{\beta}\left(\hat{q}-\hat{z}\right)\\ \partial_t\hat{q} + \partial_x \hat{w} + \frac{\sigma}{\epsilon}\hat{q} &=& -\frac{\lambda-\sigma}{\epsilon^2}\hat{q} + \frac{\lambda}{2\epsilon^2}\rho\\ \partial_t\hat{w} +\frac{a}{\epsilon^2}\partial_x\hat{q} &=& \frac{1}{\beta}\left(\hat{P}-\hat{w}\right)\ . \end{array} \right.\ \end{aligned} \end{align}$

(35)

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

$\begin{align}\label{equ:relaxation_half} \begin{aligned} &\left\{ \begin{array}{lcl} \rho^{(1)} &=& \rho^n\\ q_{\epsilon}^{(1)}&=&\frac{\epsilon^2}{\epsilon^2+\sigma\Delta t}q_{\epsilon}^n - \frac{\epsilon\Delta t\alpha}{3(\epsilon^2+\sigma\Delta t)}\overline{\partial_x m^n}\left(\frac{h^{+,n}-h^{-,n}}{h^{+,n}+h^{-,n}}\right)\hat{\rho}^{(1)}\\ z^{(1)} &=& q_{\epsilon}^{(1)}\\ w^{(1)} &=& P\left(\rho^{(1)},q_{\epsilon}^{(1)},\hat{\rho}^{(1)},\hat{q}^{(1)}\right)\\ \end{array} \right.\ \\ &\left\{ \begin{array}{lcl} \hat{\rho}^{(1)} &=& \frac{\epsilon^2}{\epsilon^2+\lambda\Delta t}\hat{\rho}^n + \frac{\Delta t\alpha}{2(\epsilon^2+\lambda\Delta t)}\overline{\partial_x m^n}\rho^n\\ \hat{q}^{(1)} &=& \frac{\epsilon^2}{\epsilon^2+(\lambda-\sigma)\Delta t}\hat{q}^n + \frac{\Delta t\lambda}{2(\epsilon^2+(\lambda-\sigma)\Delta t)}\rho^n\\ \hat{z}^{(1)} &=& \hat{q}^{(1)}\\ \hat{w}^{(1)} &=& \hat{P}\left(\rho^{(1)},q_{\epsilon}^{(1)},\hat{\rho}^{(1)},\hat{q}^{(1)}\right)\ . \end{array} \right.\ \end{aligned} \end{align}$

(36)

This relaxation scheme ensures that system (35) has the correct asymptotic limit as $\epsilon\to 0$ . Thus (35) converges to the Keller-Segel equation (9). Moreover, the solution in the relaxation part (36) preserves the realizability condition i.e., $\vert\epsilon q_{\epsilon}^{(1)}\vert\leq \rho^{(1)}$ .

The transport part in (35) includes four independent subsystems with (36) as the initial data. Introducing the characteristic variables $(U_1,V_1)$ , $(\bar{U}_1,\bar{V}_1)$ , $(U_2,V_2)$ , $(\bar{U}_2,\bar{V}_2)$ ,

$\begin{align}\label{equ:c_va_hm1} \begin{aligned} &U_1=\frac12\rho+\frac{\epsilon}{2\sqrt{a}}z\ ,&&\bar{U}_1=\frac{\epsilon}{2} q_{\epsilon}+\frac{1}{2\sqrt{a}}w\ , \\ &U_2=\frac{\epsilon}{2}\hat{\rho} + \frac{1}{2\sqrt{a}}\hat{z}\ ,&&\bar{U}_2=\frac{1}{2} \hat{q} + \frac{\epsilon}{2\sqrt{a}}\hat{w}\ , \\ &V_1=\frac12\rho-\frac{\epsilon}{2\sqrt{a}}z\ ,&&\bar{V}_1=-\frac{\epsilon}{2} q_{\epsilon}+\frac{1}{2\sqrt{a}}w\ , \\ &V_2=\frac{\epsilon}{2}\hat{\rho} - \frac{1}{2\sqrt{a}}\hat{z}\ ,&&\bar{V}_2=-\frac{1}{2} \hat{q} + \frac{\epsilon}{2\sqrt{a}}\hat{w}\, &\\ \end{aligned} \end{align}$

(37)

we use the APWB scheme for the left hand side of (35). With $a:=\max\left(\Big\vert\frac{w}{\rho}\Big\vert\right)$ 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 $\beta\to 0$ and $\epsilon\to 0$ .

4. Coupling conditions for chemotaxis models on networks

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.

4.1. Coupling condition for kinetic equations [6]

Consider the kinetic equation

$\partial_t f +\frac{1}{\epsilon} v\partial_x f = -\frac{\lambda}{\epsilon^2}\left( f - \frac{\rho}{2} \right)+\frac{1}{2 \epsilon}\alpha v \overline{\partial_x m} \rho\ .$

(38)

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

$\begin{align*} f^+ = A f^-, \end{align*}$

Figure 1. Sketch of a tripod network.

DownLoad: Full-Size Img PowerPoint

where $f_i^+=f_i(v)$ and $f_i^-=f_i(-v)$ for $v\in[0,1]$ and $i = 1,\dots,N$ . In order to conserve the total mass in the system the matrix $A\in \mathbb{R}^{N\times N}$ has to fulfill

$\begin{align*} \sum\limits_{i=1}^{N}a_{i,j} = 1 \ \ \ \ \ \ \ \ \forall j = 1,\dots,N\ . \end{align*}$

Moreover, positivity requires $a_{i,j}\geq 0$ , $i,j=1,\dots,N$ . If all edges are treated equally, the coupling conditions for three edges are

$\left[ \begin{array}{c} f_1^+\\ f_2^+\\ f_3^+\\ \end{array} \right] = \left[ \begin{array}{ccc} 0&1/2&1/2\\ 1/2&0&1/2\\ 1/2&1/2&0\\ \end{array}\right]\left[ \begin{array}{c} f_1^-\\ f_2^-\\ f_3^-\\ \end{array} \right]\ .$

(39)

This can be generalized to the cases $N>3$ by choosing $a_{i,j}=\frac{1}{N-1}$ $i\neq j$ and $a_{i,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.

4.2. Coupling condition for the nonlinear half-moment model

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)

$\begin{align}\label{equ:coup_half_moment} \left\{ \begin{array}{lcl} \left[ \begin{array}{c} \rho_1^+\\ \rho_2^+\\ \rho_3^+\\ \end{array} \right] &=& \left[ \begin{array}{ccc} 0&1/2&1/2\\ 1/2&0&1/2\\ 1/2&1/2&0\\ \end{array}\right]\left[ \begin{array}{c} \rho_1^-\\ \rho_2^-\\ \rho_3^-\\ \end{array} \right] \\ \left[ \begin{array}{c} q_1^+\\ q_2^+\\ q_3^+\\ \end{array} \right] &=& - \left[ \begin{array}{ccc} 0&1/2&1/2\\ 1/2&0&1/2\\ 1/2&1/2&0\\ \end{array}\right]\left[ \begin{array}{c} q_1^-\\ q_2^-\\ q_3^-\\ \end{array} \right].\\ \end{array} \right. \end{align}$

(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

$\sum\limits_{i=1}^3q_{\epsilon,i} = \frac{1}{\epsilon}\sum\limits_{i=1}^3(q_i^+ + q_i^-)=0\ .$

In particular, the coupling conditions for the half-moment model yield the positivity of the densities, i.e. $\rho^+_i\geq 0\ \forall i=1,\dots,N$ if $\rho^-_j\geq 0\ \forall j=1,\dots,N$ .

4.3. Numerical coupling of the nonlinear half-moment model

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

$\begin{align}\label{equ:transport_m1_half} \begin{aligned} &\left\{ \begin{array}{lcl} \partial_t\rho + \partial_xz &=& 0\\ \partial_t z +\frac{a}{\epsilon^2}\partial_x\rho + \frac{\lambda}{\epsilon^2}z-\frac{\alpha}{3\epsilon^2}\overline{\partial_xm}\rho &=&0\\ \end{array} \right.\ \\ &\left\{ \begin{array}{lcl} \partial_tq_{\epsilon} +\frac{1}{\epsilon^2} \partial_x w +\frac{\lambda-\sigma}{\epsilon^2}q_{\epsilon} -\frac{\alpha}{3\epsilon^2}\overline{\partial_x m}\frac{2}{(h^+ + h^-)}w &=&0\\ \partial_tw + a\partial_x q_{\epsilon} &=&0\\ \end{array} \right. \\ &\left\{ \begin{array}{lcl} \partial_t\hat{\rho} +\frac{1}{\epsilon^2}\partial_x\hat{z} &=&0\\ \partial_t \hat{z}+ a\partial_x\hat{\rho}+\frac{\lambda}{\epsilon^2}\hat{z} &=&0\\ \end{array} \right. \\ &\left\{ \begin{array}{lcl} \partial_t\hat{q} + \partial_x \hat{w} +\frac{\sigma}{\epsilon^2}\hat{q} &=& 0\\ \partial_t\hat{w} +\frac{a}{\epsilon^2}\partial_x\hat{q}&=& 0\ . \end{array} \right.\ \end{aligned} \end{align}$

(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: $z_i = q_i,\hat{z}_i = \hat{q}_i$ and $w_i = P_i,\hat{w}_i = \hat{P}_i$ : Thus, conditions for the states $z$ and $\hat z$ are found from those prescribed for $q^\mp$ . To obtain conditions for $w$ and $\hat w$ we have to find conditions for $P_i^\pm$ . These are easily derived from the kinetic coupling conditions (39). Multiply (39) by $v^2$ and integrate over positive velocities to obtain

$\begin{align}\label{equ:cc_P_halfmoments} \left[ \begin{array}{c} P_1^+\\ P_2^+\\ P_3^+\\ \end{array} \right] &= \left[ \begin{array}{ccc} 0&1/2&1/2\\ 1/2&0&1/2\\ 1/2&1/2&0\\ \end{array}\right]\left[ \begin{array}{c} P_1^-\\ P_2^-\\ P_3^-\\ \end{array} \right]. \end{align}$

(42)

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

$\begin{align}\label{equ:cc_rho_q} \begin{array}{lcl} \left[ \begin{array}{c} \rho_1+ \epsilon\hat\rho_1\\ \rho_2+ \epsilon\hat\rho_2\\ \rho_3+ \epsilon\hat\rho_3 \end{array} \right] &=& \left[ \begin{array}{ccc} 0&1/2&1/2\\ 1/2&0&1/2\\ 1/2&1/2&0\\ \end{array}\right]\left[ \begin{array}{c} \rho_1- \epsilon\hat\rho_1\\ \rho_2- \epsilon\hat\rho_2\\ \rho_3- \epsilon\hat\rho_3 \end{array} \right] \\ \left[ \begin{array}{c} \epsilon q_1+ \hat q_1\\ \epsilon q_2+ \hat q_2\\ \epsilon q_3+ \hat q_3 \end{array} \right] &=& - \left[ \begin{array}{ccc} 0&1/2&1/2\\ 1/2&0&1/2\\ 1/2&1/2&0\\ \end{array}\right]\left[ \begin{array}{c} \epsilon q_1- \hat q_1\\ \epsilon q_2- \hat q_2\\ \epsilon q_3- \hat q_3 \end{array} \right] \end{array} \end{align}$

(43)

and with (42) we obtain

$\left[ \begin{array}{c} P_1+ \epsilon \hat P_1\\ P_2+ \epsilon \hat P_2\\ P_3+\epsilon \hat P_3 \end{array} \right] = \left[ \begin{array}{ccc} 0&1/2&1/2\\ 1/2&0&1/2\\ 1/2&1/2&0\\ \end{array}\right]\left[ \begin{array}{c} P_1-\epsilon \hat P_1\\ P_2-\epsilon \hat P_2\\ P_3-\epsilon \hat P_3 \end{array} \right]\ .$

(44)

This leads to

$\begin{align}\label{equ:cc_z_w} \begin{aligned} \left[ \begin{array}{c} \epsilon z_1+ \hat z_1\\ \epsilon z_2+ \hat z_2\\ \epsilon z_3+ \hat z_3 \end{array} \right] &= - \left[ \begin{array}{ccc} 0&1/2&1/2\\ 1/2&0&1/2\\ 1/2&1/2&0\\ \end{array}\right]\left[ \begin{array}{c} \epsilon z_1- \hat z_1\\ \epsilon z_2- \hat z_2\\ \epsilon z_3- \hat z_3 \end{array} \right]\\ \left[ \begin{array}{c} w_1+ \epsilon\hat w_1\\ w_2+ \epsilon\hat w_2\\ w_3+ \epsilon\hat w_3 \end{array} \right] &= \left[ \begin{array}{ccc} 0&1/2&1/2\\ 1/2&0&1/2\\ 1/2&1/2&0\\ \end{array}\right]\left[ \begin{array}{c} w_1- \epsilon\hat w_1\\ w_2- \epsilon\hat w_2\\ w_3- \epsilon\hat w_3 \end{array} \right]\ . \end{aligned} \end{align}$

(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 $\rho$ , 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.

5. Numerical tests

In this section we investigate the proposed models in several numerical test cases for different values of $\epsilon$ . For simplicity we use for the nonlinear half-moment model the explicit Kershaw closure, i.e. $h^{\pm}$ given by (20). In all the considered examples we will use the following values for the parameters $\lambda =\alpha =1$ , $D = 1$ , $\gamma_{\rho}=1$ and $\gamma_m=0.1$ . The spatial resolution is $\Delta 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 $N_v=50$ cells.

5.1. Numerical solutions on an interval

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

$f(x,v,0)=F(v)\rho(x,0)=\frac{1}{2}\rho(x,0)\ ,$

with

$\rho(x,0)= \left\{ \begin{array}{lr} 1&\text{ if } 0\leq x\leq 1 \\ 0&\text{ if } 1 \leq x\leq 2\ . \end{array} \right.$

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 $\epsilon=1,\ 0.5,\ 0.1$ and $\epsilon=\ 10^{-6}$ are shown. For the kinetic model we observe that diffusion depends strongly on the value of $\epsilon$ 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 $P_1$ model for $\epsilon=1$ and $\epsilon=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 $\epsilon$ . 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 $\epsilon$ 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

$\epsilon = 1$ (top) and

$\epsilon = 0.5$ (bottom).

DownLoad: Full-Size Img PowerPoint

Figure 3. Numerical solutions of the four models on an interval at time

$t=0.2$ with

$\epsilon = 0.1$ (top) and

$\epsilon = 10^{-6}$ (bottom).

DownLoad: Full-Size Img PowerPoint

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)= \left\{ \begin{array}{lcl} a\exp(vb)&if &-1\leq v\leq 0\\ 0& if &0\leq v\leq 1\\ \end{array} \right.$

(46)

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

$\begin{align}\label{equ:initial_half_moment} \begin{aligned} &\rho^-(x,0)= \int_{-1}^0 f dv = \rho_0\ ,&& q^-(x,0)=\int_{-1}^0 vf dv = -0.9\times\rho_0\ ,&\\ &\rho^+(x,0)= \int_{0}^1 f dv = 0\ ,&& q^+(x,0)=\int_{0}^1 vf dv = 0\ ,&\end{aligned} \end{align}$

(47)

with

$\rho_0(x)= \left\{\begin{array}{lcl} 1 &\text{ if }& 0\leq x\leq 1\\ 10^{-6} &\text{ if }& 1 \leq x\leq 2\\ \end{array}\right., m(x,0)=0,\ \forall x\in[0,2] \ .$

The initial conditions for the full-moment models are

$\rho(x,0):=\rho^+ + \rho^- = \rho_0, q_{\epsilon}(x,0):=\frac{q^++q^-}{\epsilon}=-\frac{0.9}{\epsilon}\rho_0\ .$

The solutions are plotted in Figure 4, for $\Delta x = 0.02$ at $t=0.5$ . Even though the initial conditions satisfy the realizability conditions $\rho^{\pm}\geq\vert q^{\pm}\vert$ (for half-moment models) and $\rho\geq\vert q\vert$ (for full-moment models), only the nonlinear full-and half-moment models preserve the positivity of the cell density.

Figure 4. Linear models (

$P_1$ ) with negative full or half range densities and positivity preserving nonlinear models (

$M_1$ ).

DownLoad: Full-Size Img PowerPoint

5.2. Numerical solutions on a tripod network

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

$\begin{align*} \rho_1(x,0) = 1\ ,&& \rho_2(x,0) = 4\ ,&& \rho_3(x,0) = 3\ , \end{align*}$

which is consistent with the following values for the kinetic equation

$f_i(x,v,0)=F(v)\rho_i(x,0)=\frac{1}{2}\rho_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 $\epsilon$ 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 $\epsilon$ decreases all solutions approach the solution of the Keller-Segel equation.

Figure 5. Numerical solutions on a tripod network at time

$t=0.2$ ,

$\Delta x=0.02$ ,

$\epsilon = 1$ .

DownLoad: Full-Size Img PowerPoint

Figure 6. Numerical solutions on a tripod network at time

$t=0.2$ ,

$\Delta x=0.02$ ,

$\epsilon = 0.5$ .

DownLoad: Full-Size Img PowerPoint

Figure 7. Numerical solutions on a tripod network at time

$t=0.2$ ,

$\Delta x=0.02$ ,

$\epsilon = 0.1$ .

DownLoad: Full-Size Img PowerPoint

Figure 8. Numerical solutions on a tripod network at time

$t=0.2$ ,

$\Delta x=0.02$ ,

$\epsilon = 10^-6$ .

DownLoad: Full-Size Img PowerPoint

5.3. Numerical solutions on a larger network

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 $\sqrt{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)=\frac12$ 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 $\epsilon=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$ .

DownLoad: Full-Size Img PowerPoint

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 $P_1$ model with approximately $15 \%$ of the computation time of the kinetic model. The Keller-Segel model needs $18\%$ , the half moment $P_1$ model needs $22\%$ and the half moment $M_1$ 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.

DownLoad: Full-Size Img PowerPoint

6. Conclusions

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 $\epsilon$ . 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 $\rho$ in the nonlinear half moment model on the graph.

References

[1]	Wen A, Balogh B, Momol MT, et al. (2009) Management of bacterial spot of tomato with phosphorous acid salts. Crop prot 28: 859-863.
[2]	Graham J, Gruber B, Bock C (2014) Research progress for integrated canker management. Citrus Industry 95: 20-24.
[3]	Arias-Estévez M, López-Periago E, Martínez-Carballo E, et al. (2008) The mobility and degradation of pesticides in soils and the pollution of groundwater resources. Agr Ecosyst Environ 123: 247-260.
[4]	Hoang TC, Rogevich EC, Rand GM, et al. (2008) Copper desorption in flooded agricultural soils and toxicity to the Florida apple snail (Pomacea paludosa): Implications in Everglades restoration. Environ Poll 154: 338-347.
[5]	Behlau F, Belasque J, Graham J, et al. (2010) Effect of frequency of copper applications on control of citrus canker and the yield of young bearing sweet orange trees. Crop Prot 29: 300-305. doi: 10.1016/j.cropro.2009.12.010
[6]	Schuler LJ, Hoang TC, Rand GM (2008) Aquatic risk assessment of copper in freshwater and saltwater ecosystems of South Florida. Ecotoxicol 17: 642-659. doi: 10.1007/s10646-008-0236-7
[7]	Behlau F, Canteros BI, Jones JB, et al. (2012) Copper resistance genes from different xanthomonads and citrus epiphytic bacteria confer resistance to Xanthomonas citri subsp. citri. Eur j plant pathol 133: 949-963. doi: 10.1007/s10658-012-9966-8
[8]	Behlau F, Canteros BI, Minsavage GV, et al. (2011) Molecular characterization of copper resistance genes from Xanthomonas citri subsp. citri and Xanthomonas alfalfae subsp. citrumelonis. Appl environ microb 77: 4089-4096.
[9]	Ed Etxeberria PG, Priyanka Bhattacharya, Parvesh Sharma, et al. (2016) Determining the size exclusion for nanoparticles in citrus leaves. Hortic Sci 51: 1-6.
[10]	Hendricks KE, Donahoo RS, Roberts PD, et al. (2013) Effect of copper on growth characteristics and disease control of the recently introduced Guignardia citricarpa on Citrus in Florida. Am J Plant Sci 4: 282.
[11]	Jorgensen JH, Turnidge JD (2015) Susceptibility Test Methods: Dilution and Disk Diffusion Methods*.
[12]	Young M, Santra S (2014) Copper (Cu)–Silica Nanocomposite Containing Valence-Engineered Cu: A New Strategy for Improving the Antimicrobial Efficacy of Cu Biocides. J agr food chem 62: 6043-6052. doi: 10.1021/jf502350w
[13]	Hans M, Erbe A, Mathews S, et al. (2013) Role of copper oxides in contact killing of bacteria. Langmuir 29: 16160-16166. doi: 10.1021/la404091z
[14]	Balogh B, Jones J, Momol M, et al. (2003) Improved efficacy of newly formulated bacteriophages for management of bacterial spot on tomato. Plant Dis 87: 949-954. doi: 10.1094/PDIS.2003.87.8.949
[15]	Baron M, Arellano JB, Gorge JL (1995) Copper and photosysten II-A controversial relationship. Physiol Plant 94: 174-180. doi: 10.1111/j.1399-3054.1995.tb00799.x
[16]	Kuepper H, Goetz B, Mijovilovich A, et al. (2009) Complexation and Toxicity of Copper in Higher Plants. I. Characterization of Copper Accumulation, Speciation, and Toxicity in Crassula helmsii as a New Copper Accumulator. Plant Physiol 151: 702-714.
[17]	Gupta N, Ram H, Kumar B (2016) Mechanism of Zinc absorption in plants: uptake, transport, translocation and accumulation. Rev Environ Sci Bio-Technol 15: 89-109. doi: 10.1007/s11157-016-9390-1
[18]	Haslett BS, Reid RJ, Rengel Z (2001) Zinc mobility in wheat: Uptake and distribution of zinc applied to leaves or roots. Ann Bot 87: 379-386. doi: 10.1006/anbo.2000.1349
[19]	Aslani F, Bagheri S, Julkapli NM, et al. (2014) Effects of Engineered Nanomaterials on Plants Growth: An Overview. Scientific World J 2014: 641759.
[20]	Chichiricco G, Poma A (2015) Penetration and Toxicity of Nanomaterials in Higher Plants. Nanomater 5: 851-873. doi: 10.3390/nano5020851

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

Your name:*

Email:*
© 2016 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)