On modeling Maze solving ability of slime mold via a hyperbolic model of chemotaxis

ON MODELING MAZE SOLVING ABILITY OF SLIME MOLD VIA A HYPERBOLIC MODEL OF CHEMOTAXIS arXiv:1601.01137v1 [math.NA] 6 Jan 2016 G. BRETTI AND R. NATALINI Abstract. Many studies have shown that Physarum polycephalum slime mold is able to find the shortest path in a maze. In this paper we study this behavior in a network, using a hyperbolic model of chemotaxis. Suitable transmission and boundary conditions at each node are considered to mimic the behavior of such an organism in the feeding process. Several numerical tests are presented for special network geometries to show the qualitative agreement between our model and the observed behavior of the mold. Keywords: Chemotaxis, networks, finite difference schemes, shortest path problem, physarum polycephalum, hyperbolic equations. AMS subject classification: 92C17, 92C42, 65M06, 05C38, 05C85. 1. Introduction In recent years, the behavior of Physarum polycephalum, a “many-headed” slime mold with multiple nuclei, has been studied in depth and considered as a model organism representing amoeboid movement and cell motility. The body of the plasmodium, representing the main vegetative phase of Physarum polycephalum, contains a network of tubes (pseudopodia), by means of which nutrients and chemical signals circulate throughout the organism [21]; the flow in the tube is bidirectional, as it can be observed to switch back and forth. The movement is characterized by cytoplasmic streaming, driven by rhythmic contractions of organism, that sustains and reorganizes tubes: a larger flow leads to a wider tube. An example of the behavior of starved plasmodium searching for food has been shown in an experiment performed by Nakagaki, Yamada, and Tóth in [21]. First, they built a maze on an agar surface. Therefore, they put over the agar surface some parts of the plasmodium, which started to spread and coalesced to form a single organism that filled the maze. After that, they placed agar blocks containing oatmeal at the exits and, after some hours, they observed the dead ends of the plasmodium shrank (the tubes with a smaller flux ended to disappear) and the colony, after exploring all possible connections, produced the formation of a single thick pseudopodium spanning the minimum distance between the food sources. See for instance some popular videos freely available on the web, as for instance this one [29]. Istituto per le Applicazioni del Calcolo “M.Picone”, Rome, Italy, g.bretti@iac.cnr.it, roberto.natalini@cnr.it. 1 2 G. BRETTI AND R. NATALINI Then, to summarize, the main features in the evolution of the slime mold system are the following two steps: dead end cutting, and the selection of the solution path among the competitive paths. This behavior can be applied to both path-finding in a maze and path selection in a transport network. The ability of solving this kind of problems of an unicellular organism has been deeply investigated using various mathematical models. Some kinetic models were proposed and studied in [18, 19, 25]. The hydrodynamic properties of the transportation of flux by the tubes of plasmodium and the mechanisms under the rhythmic oscillation of Physarum were investigated in [24], and biologically inspired models for adaptive network construction were developed [26] and [28]. A mathematical model based on ordinary differential equations to describe the dynamics of Physarum was introduced in [25]. It is based on a regulation mechanism, which balances the system between the thickness of a tube and the flux through it. In that paper, it was shown numerically that, in the asymptotic steady state of the model, the solution converges to the minimum-length solution between a source and a sink on any input graph. To show this property rigorously, it is necessary to prove that the globally asymptotically stable equilibrium point corresponds to the shortest path connecting two exits of a network and the system has neither the oscillating nor chaotic solution. In [18] there was proof of convergence for two parallel links, while the convergence of this model for any network to the shortest path connecting the source and the sink was analytically proved in [1]. Another approach is represented by modeling by using partial differential equations. In fact, the movement of individuals under the effect of a chemical stimulus of chemoattractant has been widely studied in mathematics in the last two decades, see [15, 20, 23], and various models involving partial differential equations have been proposed. This process is known as chemotaxis: the organism will move with higher probability towards food sources - positive chemotaxis and, while eating, it will produce a chemical substance, usually a chemokine. The opposite situation is that it moves away (negative chemotaxis), if it experience the presence of a harmful substance. The basic unknowns in these chemotactic models are the density of individuals and the concentrations of some chemical attractants. One of the most considered models of this type is the Patlak-Keller-Segel system [16], where the evolution of the density of cells is described by a parabolic equation, and the concentration of a chemoattractant is generally given by a parabolic or elliptic equation, depending on the different regimes to be described and on authors’ choices. The analytical study of the KS model on networks was presented in the recent paper [5], where the existence of a time global and spatially continuous solution for both the doubly parabolic and the parabolic-elliptic systems was proved. By contrast, models based on hyperbolic/kinetic equations for the evolution of the density of individuals, are characterized by a finite speed of propagation [8, 23, 7, 6, 13]. MAZE SOLVING ABILITY OF SLIME MOLD 3 In this paper, we deal with a hyperbolic-parabolic model, originally considered in [27], and later reconsidered in [12], which arises as a simple model for chemotaxis on a line: (1)   ut + vx = 0, vt + λ2 ux = χφx u − v,  φt − D φxx = au − bφ. The function u is the density of individuals in the considered medium, v is their averaged flux and φ denotes the density of chemoattractant. The individuals move at a constant speed λ ≥ 0, changing their direction along the axis during the time. The positive constant D is the diffusion coefficient of the chemoattractant; the positive coefficients a and b, are respectively its production and degradation rates, and χ is the chemotactic sensitivity. Let us underline that the flux v in model (1) corresponds to v = −λ2 ux + χφx u for the parabolic KS system. This system was analytically studied on the whole line and on bounded intervals in [13], while an effective numerical approximation, the Asymptotic High Order (AHO) scheme, was proposed in [22], see also [10, 11]. The AHO scheme was introduced in order to balance correctly the source term with the differential part and avoid an incorrect approximation of the density flux at equilibrium. These schemes are based on standard finite differences methods, modified by a suitable treatment of the source terms, and they take into account, using a Taylor expansion of the solution, of the behavior of the solutions near non constant stationary states. Model (1) can also be considered on a network, where solutions on each arc are coupled through transmission conditions ensuring the total continuity of the fluxes at each node, while the densities can have jumps. The first analytical study of this model on a network was given in [14], where a global existence result of solution for suitably small initial data was established. Concerning the numerical approximation, in [4], we extended the AHO scheme of [22] to the case of networks. In this case, a particular attention was paid to the proper setting of conditions at internal and external nodes in order to guarantee the conservation of the total mass. Here, inspired by the study on a doubly parabolic chemotaxis system presented in [3, 5], where the one-dimensional KS model was extended to networks, we present a numerical study of amoeboid movements, using the hyperbolic chemotaxis model (1) on networks connecting two or more exits. The main difference with previous works is not only the hyperbolicity of the equations describing the cell movement, but the transmission conditions, which in our case, unlike what was assumed in [3, 5], are designed to allow for discontinuous densities at each node. Actually, this condition is believed to be more appropriate when considering a flow of individuals on a network, see for instance [9] and references therein and also the recent discussion in [2]. 4 G. BRETTI AND R. NATALINI As in [4], the numerical approximation of the hyperbolic part of the system is based on the AHO scheme, while the parabolic part is approximated with the Crank-Nicolson scheme. However, since here we include the modeling of the inflow of the mass and the introduction of food at the network exits, we assume inflow boundary conditions at external nodes both for the density of cells and the chemoattractant. Therefore we need to deal with more general Neumann boundary conditions with respect to [4], both for the hyperbolic and the parabolic part of the system. After proposing a modified scheme, this scheme will be used to obtain numerical solutions for networks with special geometries, motivated by experimental observations. Let us notice that the properties of dead-end cutting and shortest path selection in our model are represented by different mass distribution on the edges of the assigned network described as an oriented graph. The paper is organized as follows. In Section 2 we describe the main setting of problem (1). Section 3 is devoted to the description of the numerical approximation of the problem based on the AHO scheme of second order with a suitable discretization of the transmission and boundary conditions for the case of non-null Neumann boundary conditions. In particular, a strong emphasis will be given to formulate the correct boundary conditions for the AHO schemes for the inflow and outflow cases. Finally, in Section 4, we report some numerical experiments for the distribution of densities within networks of different topologies, having in mind the laboratory experiments made with Physarum, in order to test the correspondence between our simulations and the real behavior of such an organism. In this framework, a comparison between the results obtained with the two models, the one in [3] and the present one, is also provided, showing substantially the same asymptotic behavior of the solutions, even if the transmission conditions at each node are not the same, and so the transitory profiles. 2. Analytical framework First, let us represent a network as a connected graph. We define a connected graph G = (N , A) as formed by two finite sets, a set of P nodes (or vertices) N and a set of N arcs (or edges) A, such that an arc connects a pair of nodes. Since arcs are bidirectional the graph is non-oriented, but we need to fix an artificial orientation in order to fix a sign to the velocities. The network is therefore composed of ”oriented” arcs and there are two different types of intervals at a node p ∈ N : incoming ones – the set of these intervals is denoted by Ip – and outgoing ones – whose set is denoted by Op . For example, on the network depicted in Figure 1, 1 ∈ I and 2, 3 ∈ O. We will also denote in the following by Iout and Oout the set of the arcs incoming or outgoing from the outer boundaries. The N arcs of the MAZE SOLVING ABILITY OF SLIME MOLD 5 network are parametrized as intervals [0, Li ], i = 1, . . . , N , and for an incoming arc, Li is the abscissa of the node, whereas it is 0 for an outgoing arc. We consider system (1) on each arc and rewrite it in diagonal variables for its hyperbolic part by setting v 1 u± . (2) u± = 2 λ Here u+ and u− are the Riemann invariants of the system and u+ (resp. u− ) denotes the density of cells following the orientation of the arc (resp. the density of cells going in the opposite direction). This transformation is inverted by u = u+ + u− and v = λ(u+ − u− ), and yields: 
1 + + − +    ut + λux = 2λ (χφx − λ)u + (χφx + λ)u ,  
1 (3) (χφx − λ)u+ + (χφx + λ)u− , u− =− − λu−  x t  2λ    φt − Dφxx = a(u+ + u− ) − bφ. We complement this system by initial conditions at t = 0 on each arc u+ (x, 0) = u+ 0 (x), u− (x, 0) = u− 0 (x), φ(x, 0) = φ0 (x), for x ∈ [0, L], with u+ u− φ0 some C 3 functions. 0, 0, Up to now, we omitted the indexes related to the arc number since no confusion was possible. From now on, however, we need to distinguish the quantities on different arcs and we denote by u± i , ui , vi , φi etc., the values of the corresponding variables on the i-th arc. On the outer boundaries, we consider the more general boundary conditions: ( − u+ i (0, t) = gi (t, ui (0, t)), if i ∈ Iout , (4) + u− i (Li , t) = gi (t, ui (Li , t)), if i ∈ Oout , for some nonlinear smooth functions gi (t, u), which correspond to general boundary conditions on the flux function v: ( vi (0, t) = ψi (t, u(0, t)), if i ∈ Iout , (5) vi (Li , t) = ψi (t, u(Li , t)), if i ∈ Oout , where ψi (t, u) are the same functions after the change of variables. We can also assume non-vanishing conditions for the flux of chemoattractant, to reproduce the flow of nutrient substances at the outer boundaries (6) where φ̄i are given functions. ∂x φi (., t) = φ̄i (t), 6 G. BRETTI AND R. NATALINI 2.1. Transmission conditions. Let us describe how to define the conditions at each node; this is an important point, since the behavior of the solution will be very different according to the conditions we choose. Moreover, let us recall that the coupling between the densities on the arcs are obtained through these conditions. At node p ∈ N , we have to give values to the components such that the corresponding characteristics are going out of the node. Therefore, we consider the following transmission conditions at each node:  X X  ξi,j u− ξi,j u+ (Lj , t) + u− (Li , t) =  j (0, t), if i ∈ Ip , j i   j∈Op j∈Ip (7) X X + +   ui (0, t) = ξi,j uj (Lj , t) + ξi,j u−  j (0, t), if i ∈ Op ,  j∈Ip j∈Op where the constant ξi,j ∈ [0, 1] are the transmission coefficients: they represent the probability that a cell at a node decides to move from the i−th to the j−th arc of the network, also including the turnabout on the same arc. Let us notice that the condition differs when the arc is an incoming or an outgoing arc. Indeed, for an incoming (resp. outgoing) arc, the value of the − function u+ i (resp. ui ) at the node is obtained through the system and we + need only to define u− i (resp. ui ) at the boundary. These transmission conditions ensure the continuity of fluxes at each node, but not the continuity of the densities; however, the continuity of the fluxes is enough since it implies that we cannot loose nor gain any cells during the passage through a node. This is obtained using a condition mixing the transmission coefficients ξi,j and the velocities of the arcs connected at each node p. Fixing a node and denoting the velocities of the arcs by λi ≥ 0, i ∈ Ip ∪ Op , in order to have the flux conservation at node p, which is given by: X X − − (8) λi (u+ λi (u+ i (Li , t) − ui (Li , t)) = i (0, t) − ui (0, t)), i∈Ip i∈Op it is enough to impose the following condition: X λi ξi,j = λj , j ∈ Ip ∪ Op . (9) i∈Ip ∪Op The coefficients ξi,j belong to [0, 1], and at every node p ∈ N , we have: X ξi,j = 1 for all i ∈ Ip ∪ Op . (10) j∈Ip ∪Op In the case of no-flux at the outer boundaries, condition (9) ensures that the global mass µ(t) of the system is conserved along the time, namely: N Z Li N Z Li X X (11) µ(t) = ui (x, t)dx = µ0 := ui (x, 0)dx, for all t > 0. i=1 0 i=1 0 MAZE SOLVING ABILITY OF SLIME MOLD 7 In the case of non-vanishing Neumann conditions at the outer boundaries the conservation of the global mass of the system obviously does not hold. If we consider the case treated in the numerical tests in Section 4, we have two arcs connected to the outer boundaries (a source and a sink), say arcs l and m. In such a case we have that the derivative in time of the total mass µ is given by: (12) µ′ (t) = vl (0, t) − vm (Lm , t), for all t > 0. Now let us consider the transmission conditions for φ in system (1). We complement conditions (4) and (7) with a transmission condition for φ. As previously, we do not impose the continuity of the density of chemoattractant φ, but only the continuity of the flux at node p ∈ N . Therefore, we use the Kedem-Katchalsky permeability condition [17], which has been first proposed in the case of flux through a membrane. For some positive coefficients κi,j , we impose at each node X (13) Di ∂n φi = κi,j (φj − φi ), i ∈ Ip ∪ Op , j∈Ip ∪Op that, under the condition κi,j = κj,i , i, j = 1, . . . , N the conservation of the fluxes at node p holds, that is to say: X Di ∂n φi = 0. i∈Ip ∪Op Note that setting κi,i = 0, i = 1, . . . , N , does not change condition (13), and the positivity of the transmission coefficients κi,j , guarantees the energy dissipation for the equation for φ in (1), when the term in u is absent. 3. A numerical approximation for system (1) Here we recall the numerical scheme presented in [4] and adapt this scheme to the case of non-vanishing Neumann boundary conditions. 3.1. Approximation of the hyperbolic equations of the system in the case of a network. Let us consider a network as previously defined. Each arc ai ∈ A, 1 ≤ i ≤ N , is parametrized as an interval ai = [0, Li ] and is discretized with a space step hi and discretization points xji for j = 0, . . . , Mi +1. We still denote by k the time step, which is the same for all the arcs of the network. In this subsection, we denote by win,j the discretization on the grid at time tn and at point xji of a function wi , i = 1, . . . , N on the i-th arc for j = 0, . . . , Mi + 1 and n ≥ 0. Here we describe the discretization of system (1) on each arc, denoting by f = χφx u and omitting the parabolic equation for φ. Since we also work with Neumann boundary conditions for the φ function, the function f will satisfy the following conditions on the boundary : f (0, t) = f (L, t) = 0. We therefore consider the following system  ut + vx = 0, (14) vt + λ2 ux = f − v, 8 G. BRETTI AND R. NATALINI and rewrite it in a diagonal form, using the usual change of variables (2),  1 + 1 − −   u− f, t − λux = (u − u ) − 2 2λ (15)   u+ + λu+ = 1 (u− − u+ ) + 1 f. x t 2 2λ To have a reliable scheme, with a correct resolution of fluxes at equilibrium, we have to deal with Asymptotically High Order schemes on each arc i, see [22] for the general conditions on such coefficients in order to have a consistent, monotone scheme, that is of higher order on the stationary solutions. In particular, the second order scheme in diagonal variables reads as:     n,j+1 n,j+1 n+1,j kλi k − + k4 (un,j u−,i = 1 − λi hki − k4 un,j +,i + u+,i ) −,i + hi 4 u−,i
(16) − 4λk i fin,j+1 + fin,j , j = 0, . . . , Mi ,     n+1,j n,j−1 n,j kλi k k k u+,i = 1 − λi hi − 4 u+,i + hi − 4 un,j−1 + k4 (un,j −,i + u−,i ) +,i
(17) + 4λk i fin,j−1 + fin,j , j = 1, . . . , Mi + 1. Note that we have a second–order AHO scheme on every stationary solutions, which is enough to balance the flux of the system at equilibrium. It is possible to see that component by component monotonicity conditions (neglecting the source term f ) are satisfied if 4h , h + 4λ see [22] for more details. For simplicity here we only deal with second order schemes, but third order schemes could also be considered (these schemes require a fourth–order AHO scheme for the parabolic equation with a fivepoints discretization for φx ). Finally, we rewrite the second order scheme (16)-(17) in the variables u and v as :  λk k  n,j+1 i n,j−1 (un,j+1 − 2un,j ) vi − vin,j−1 + uin+1,j = un,j i i + ui i − 2hi 2hi k k + (−vin,j−1 + vin,j+1 ) + (fin,j−1 − fin,j+1 ), 4λi 4λi (19)  λk 2k  λ i n,j−1 vin+1,j = vin,j − i (v n,j+1 − 2vin,j + vin,j−1 ) un,j+1 − u + i i 2hi 2hi i k k (f n,j−1 + 2fin,j + fin,j+1 ). − (vin,j−1 + 2vin,j + vin,j+1 ) + 4 4λi i (18) h ≤ 4λ and k ≤ 3.1.1. Boundary conditions. Now, we define the numerical boundary conditions associated to this scheme. As shown in [4], to implement this scheme on each interval we need two boundary conditions at the endpoints linked to external nodes and two transmission conditions at the endpoints linked to internal nodes. Considering an arc and its initial and end nodes, there MAZE SOLVING ABILITY OF SLIME MOLD 9 are two possibilities: either they are external nodes, namely nodes from the outer boundaries linked to only one arc, or they are internal nodes connecting several arcs together. The boundary and transmission conditions will therefore depend on this feature. Below, we will impose the boundary conditions for v: (21) or (28) and (23) at outer nodes, and two transmission conditions (22)–(24) at inner nodes. The first type of boundary conditions for will come from no-flux condition at outer nodes : ( vin+1,0 = 0, if i ∈ Iout , (20) vin+1,Mi +1 = 0, if i ∈ Oout , where Iout (resp. Oout ) means that the arc is incoming from (resp. outgoing to) the outer boundary. In the u± -variables, these conditions become : ( un+1,0 = un+1,0 +,i −,i , if i ∈ Iout , (21) i +1 i +1 , if i ∈ Oout . un+1,M = un+1,M +,i −,i The boundary condition at a node p will come from a discretization of the transmission condition (7), that is to say  X X n,M +1 n,Mi +1  ξi,j un,0 ξi,j u+,j j + u−,i =  −,j , if i ∈ Ip ,   j∈Op j∈Ip (22) X X n,M +1 n,0 j   ξi,j un,0 ξi,j u+,j + u+,i =  −,j , if i ∈ Op .  j∈Ip j∈Op However, these relations link all the unknowns together and they cannot be used alone. An effective way to compute all these quantities will be presented after equation (24) below. We still have two missing conditions per arc, which can be recovered by using the boundary conditions on v. In particular, in the case of null Neumann boundary conditions (21), from the exact mass conservation between two successive computational steps,we obtain for the (16)-(17) the following boundary conditions:  k λi k n+1,0  − 1)un,1 un+1,0 = u−,i = (1 − λi )un,0  −,i + ( −,i +,i   hi 4 hi     k k n,1   + un,1 f , if i ∈ Iout ,  +,i − 4 4λi i (23) k k λi  n+1,M +1 n+1,M +1 +1   u+,i = u−,i = (1 − λi )un,M + ( + 1)un,M  +,i +,i  hi 4 hi     k n,Mi k   f , if i ∈ Oout , − un,M −,i − 4 4λi i where Iout and Oout have the same meaning as previously. i +1 Using the transmission conditions (22) for un+1,M if i ∈ Ip and for −,i n+1,0 u+,i if i ∈ Op , we obtain the following numerical boundary conditions, 10 G. BRETTI AND R. NATALINI  −1 P with δi = hi hi + j∈Ip ∪Op hj ξj,i :  k λi k
n,Mi +1 n,Mi +1 i +1 un+1,M = δ (1 − ) + u−,i 1 − 2k + i u+,i +,i 2 hi 2 (24) 
1 k 1 1 k λ i n,Mi n,Mi +1 n,Mi i − ( − 2 + kun,M u + f + f ) , if i ∈ Ip , +,i hi 2 2 −,i λi 2 i 2 i  k
k n,1 λi k n,0 n,0 1 − + ) + u + u+,i un+1,0 = δ u (1 − 2k i −,i −,i +,i hi 2 2 2 λi 1
k 1 n,1 1 n,0  n,1 + ku−,i 2 − − ( fi + fi ) , if i ∈ Op . hi 2 λi 2 2 Once these quantities are computed, we can use equations (22) at time tn+1 , n+1,0 i +1 to obtain un+1,M if i ∈ Ip and u+,i if i ∈ Op . −,i As proved in [4], the consistency of (24) at a node holds under the condition on each arc (25) hi = 2kλi . Please notice that using this condition in combination with (18), implies that k ≤ 2. In conclusion, we have imposed four boundary conditions (21), (22), (23), and (24) on each interval. Conditions (21) and (23) deal with the outer boundary, whereas conditions (22) and (24) deal with the node. Under these conditions, the total numerical mass is conserved at each step. Let us now consider the more general case of non vanishing Neumann boundary conditions, which is finally the main goal of this section. In particular, we consider, as in [3], the inflow condition at the nodes p1 and p2 , respectively: 2 , l ∈ O p1 1 + ul (0, t) (26) vl (0, t) = (27) vm (Lm , t) = − 2 m ∈ Ip2 . 1 + um (Lm , t) The discretization of conditions (26)-(27) is:  2   vln+1,0 =  n+1,0 n+1,0 , l ∈ Op1  1 + u+,l + u−,l (28) 2  n+1,Mm +1   =− , m ∈ Ip2 .  vm n+1,Mm +1 n+1,Mm +1 1 + u+,m + u−,m Since in the general case of non-null Neumann conditions the conservation of the total mass of the system does not hold, we need to compute the numerical approximation scheme at the outer boundaries in such a case, taking into account the relation (12). The discrete total mass is given by MAZE SOLVING ABILITY OF SLIME MOLD n Itot = N X 11 I ni , where the mass corresponding to the arc i is defined as: i=1   Mi n,0 n,Mi +1 X u u i  I ni = hi  i + un,j i + 2 2 j=1   n,0 n,Mi +1 n,Mi +1 Mi X un,0 + u u + u +,i −,i +,i −,i n,j . + = hi  (un,j +,i + u−,i ) + 2 2 j=1 n+1 n , we find: Computing Itot − Itot N X hi k 1 n+1,0 1 n+1,0 λi 1 n,0 (u+,i − un,0 − un,0 − )u +,i ) + (u−,i −,i ) + (2 2 k k hi 2 +,i i=1  ! λi 1 n,1 1 1 n,1 1 n,0 1 n,0 1 n,1 f + fi + u−,i − u+,i − (2 − )u−,i + 2 2 hi 2 λi 2 i 2 n+1 Itot + − hi k 2 n Itot = 1 n+1,Mi +1 1 n+1,Mi +1 1 n,Mi +1 n,Mi +1 i +1 (u − u+,i ) + (u−,i − un,M ) + u+,i −,i k +,i k 2 λi 1 n,Mi +1 λi 1 1 n,Mi i − (2 − )un,M − )u−,i +,i − u−,i hi 2 hi 2 2  ! 1 1 n,Mi +1 1 n,Mi . f + fi − λi 2 i 2 + (2 We are going to impose boundary conditions such that the right-hand side in the previous difference is equal to the difference between the fluxes at the outer boundaries. In particular, in the case of the networks considered in Section 4, at the outer boundaries we have the inflow conditions (26) for an arc l exiting from the source node and (27) for an arc m entering a sink node. Then, we discretize the relation (12), and, using the trapezoidal rule it can be written as: k n+1 n,Mm +1 n+1,Mm +1 n + vln+1,0 − vm ) Itot = (vln,0 − vm − Itot 2 and then, splitting the terms, for the left boundary (on arc l) and right boundary (on arc m), we get, respectively, ! 2 k 2 n,0 + 2 1 + un,0 1 + un+1,0 + un+1,0 +,l + u−,l +,l −,l (29) 1 n+1,0 1 λl 1 (u+,l − un,0 ) + (un+1,0 − un,0 ) + (2 − )un,0 +,l −,l −,l k k hl 2 +,l  ! λl 1 n,1 1 1 n,1 1 n,0 1 n,0 1 n,1 f + fl + u−,l − u+,l − (2 − )u−,l + 2 2 hl 2 λl 2 l 2 = h1 k 2 12 G. BRETTI AND R. NATALINI and (30) 2 k 2 = m +1 m +1 1 + un,M + un,M +,m −,m hm k 2 + 2 m +1 m +1 1 + un+1,M + un+1,M +,m −,m ! 1 n+1,Mm +1 1 1 n,Mm +1 n,Mm +1 n,Mm +1 m +1 (u − u+,m ) + (un+1,M − u−,m ) + u+,m k +,m k −,m 2 λm 1 n,Mm +1 λm 1 n,Mm 1 n,Mm − u−,m − )u−,m − (2 − )u hm 2 hm 2 +,m 2 !   1 n,Mm +1 1 n,Mm 1 f + fm . − λm 2 m 2 + (2 Let us now focus on the condition for arc l. Setting 2 λl 1 n,0 1 n,0 + hl − (un,0 +,l + u−,l ) + (2 h − 2 )u+,l k 1+ + l  ! λl 1 n,1 1 1 n,1 1 n,0 1 n,0 1 n,1 , f + fl + u−,l − u+,l − (2 − )u−,l + 2 2 hl 2 λl 2 l 2 An = − (31) n,0 u+,l un,0 −,l and the other quantities hl , β1n = α1 + An , γ1n = An − 2, k we can write the equation (29) as (33) n+1,0 n+1,0 n+1,0 2 n n α1 (un+1,0 + γ1n + un+1,0 +,l ) + (β1 + 2α1 u−,l )u+,l −,l (α1 (1 + u−,l ) + A ) = 0. (32) α1 = Then, choosing the positive root of equation (33), we get the formula: (34) n+1,0 u+,l " 1 −β1n − 2α1 un+1,0 := g1n (un+1,0 −,l −,l ) = 2α 1 r   + β1n + 2α1 un+1,0 −,l 2 − 4α1 γ1n + n+1,0 u−,l (An + α1 ) + 2 α1 (un+1,0 −,l )  # We compute (34) using the value of un+1,0 obtained from the numerical −,l scheme (16) at the boundary j = 0, with n,0 fln,0 = χun,0 l φx,l , and fln,1 = n,1 χun,1 l φx,l = n,2 − φn,0 n,1 φl l χul . 2hl . MAZE SOLVING ABILITY OF SLIME MOLD 13 n,0 n,0 Then, plugging the expression above into (16) and denoting u+,l = g1n−1 (u−,l ), we have: (35)     λl k k k k k n+1,0 n,0 n,0 n,0 n−1 n,0 u−,l = u−,l 1 − + g1 (u−,l ) − − χφ χφ − hl 4 4λl x,l 4 4λl x,l ! ! n,2 n,0 n,2 n,0 kχ φl − φl kχ φl − φl λl k k k n,1 n,1 + u−,l − − − + u+,l . hl 4 4λl 2hl 4 4λl 2hl Let us now analyze the sign of the coefficients in (35) in order to ensure the monotonicity of the scheme. Under the condition (25), we want to satisfy the following relations: (36)  n,2 n,0  kχ φl − φl λl k k   − − ≥ 0,  hl 4 4λl 2hl n,2 n,0  k kχ φl − φl    − ≥ 0, 4 4λl 2hl and we also have to guarantee the monotonicity with respect to un,0 −,l of the first two terms on the right-hand side of equation (35):     k k λl k k k n,0 n,0 n,0 n−1 n,0 + g1 (u−,l ) . − − χφ χφ − (37) u−,l 1 − hl 4 4λl x,l 4 4λl x,l For (36) we find the condition:   1 1 2λl n,1 (38) φx,l ≤ − , with k ≤ 1. k 2 χ It is not obvious to have that the condition (38) is respected, since, in general, the gradient of φ can increase during chemotactic process if the mass grows, thus causing blow-up of solutions. Then at each time step we need to check if the condition is satisfied, in order to have a finite solution. Setting un,0 −,l = u in (37) and passing to the derivative respect to u, we study the expression     k k k λl k k n,0 n,1 n−1 ′ + (g1 ) (u) ≥ 0. − − χφ χφ − (39) 1− hl 4 4λl x,l 4 4λl x,l Since we have dg1n−1 (u) = du   1  4α1 (β1n−1 − α1 + An−1 )  = −1, −2α1 + q 2α1 2 (β1n−1 − 2α1 u)2 − 4α1 (γ1n−1 + u(An−1 + α1 ) + α1 u2 ) then (39) under the condition (25) reduces to k ≤ 1. 14 G. BRETTI AND R. NATALINI Reasoning as above, for arc m we obtain: m +1 un+1,M +,m (40)  k λm k k n,Mm +1 n,Mm +1 u+,m − − χφ = 1− hm 4 4λm x,m   k k n−1 n,Mm +1 n,Mm +1 χφ + + g2 (u+,m ) 4 4λm x,m ! n,Mm +1 n,Mm −1 − φm kχ φm k n,Mm + + u−,m 4 4λm 2hm  + m un,M +,m n,Mm +1 n,Mm −1 λm k k kχ φm − φm − + hm 4 4λm 2hm ! . As above, the condition (25) reduces to k ≤ 1. 3.2. Approximation of the parabolic equation for φ. We solve the parabolic equation, using a finite differences scheme in space and a CrankNicolson method in time, namely an explicit-implicit method in time. Therefore, we will have the following equation for φin+1,j , 1 ≤ j ≤ Mi ,  Di k  n,j+1 n,j n,j−1 −φ + 2φ − φ i i i 2h2i  Di k  n+1,j n+1,j−1 + 2φ − φ − 2 −φn+1,j+1 i i i 2hi ai k n+1,j bi k n+1,j + (ui + un,j (φ + φn,j i )− i ). 2 2 i φin+1,j = φn,j i − (41) Now, let us find the two boundary conditions needed on each interval. As in [4], the boundary conditions will be given in the case of an outer node and in the case of an inner node. On the outer boundary we discretize φx using a second order approximation. Then, the numerical condition, corresponding to (6), for a general φ̄, not necessarily depending on φn+1 , reads as: i (42)  4 1 2hi   φn+1,0 = φin+1,1 − φin+1,2 − φ̄, if i ∈ Iout , i 3 3 3   φn+1,Mi +1 = 4 φn+1,Mi − 1 φn+1,Mi −1 + 2hi φ̄, if i ∈ O . out i 3 i 3 i 3 Let us now describe our numerical approximation for the transmission condition (13) which, as the transmission condition for the hyperbolic part (7), couples the φ functions of arcs having a node in common. MAZE SOLVING ABILITY OF SLIME MOLD 15 At node p condition (13) is discretized using the second-order discretization as described in [4]: 1 2 hi X 4 n+1,Mj +1 i − φin+1,Mi −1 + ηip φin+1,Mi +1 = φn+1,M κi,j φj i 3 3 3 Di j∈Ip + (43) 2 hi X κi,j φn+1,0 , if i ∈ Ip , j 3 Di j∈Op 1 n+1,2 2 hi X 4 n+1,Mj +1 − φi + κi,j φj ηip φin+1,0 = φn+1,1 i 3 3 3 Di j∈Ip 2 hi X + κi,j φn+1,0 , if i ∈ Op , j 3 Di j∈Op hi with ηip = 1+ 23 D j∈Ip ∪Op κi,j . Let us remark that the previous discretizai tions are compatible with relations (42) with φ̄ = 0 (homogeneous Neumann boundary condtions) considering that for outer boundaries the coefficients κi,j are null. Therefore, in this case, the value of ηiout is just equal to 1. Since equations (43) are coupling of the unknowns of all arcs altogether, we have to solve a large system which contains all the equations of type (41) and also the discretizations of transmission conditions (43). Note that for the computational resolution of the mentioned system, characterized by a sparse banded matrix, we used the LAPACK-Linear Algebra PACKage routine DGBSV designed for banded matrix. Once the values of φin+1,j are known, we can compute a second-order discretization of the derivatives of φ which gives the values of the f function, namely :   1  n+1,j+1 n+1,j−1   φ − φ , 1 ≤ j ≤ Mi ,  i i  2 hi    1   n+1,1 n+1,0 n+1,j −φn+1,2 + 4φ − 3φ , j = 0, φx,i = i i i  2 hi      1  n+1,Mi −1   φi − 4φin+1,Mi + 3φin+1,Mi +1 , j = Mi + 1. 2 hi The discretization of f needed at equations (19),(23), and (24) is therefore n+1,j n+1,j given by fin+1,j = χφx,i ui . P 4. Tests on different network geometries Here we present some computational tests on different network topologies to investigate the behaviour of individuals moving through them. We firstly consider the basic network composed of one incoming and two outgoing arcs (T-shaped network) in order to show that the model is able to reproduce the main features of the plasmodium. Then we consider a network of 7 arcs and 6 nodes (diamond-like graph) and a more complex network of 26 arcs and 18 nodes, both of them connecting two exits (a source and a sink), where inflow conditions are implemented to mimic the presence of food sources. 16 G. BRETTI AND R. NATALINI Finally, we consider a network of 21 arcs and 15 nodes with multiple exits (two sources and three sinks). For all the networks we deal with we assume to have equal velocities λi = λ and equally distributed transmission coefficients at each node. If, for instance, we have a node p connected with the total number of N arcs (no matter how many are incoming or outgoing), we assume to have: 1 (44) ξij = , ∀i ∈ Ip , ∀j ∈ Op , N in such a way to satisfy the conditions (9) and (10). 4.1. T-shaped graph: four nodes and three edges. Here we consider the network of three arcs and four nodes (1 internal node and three external nodes) shown in Figure 1. We assume to have √ arcs of length Li = 1 and we set parameters as ai = 1, bi = 0.1, λi = 0.33, Di = 1, with χi = 1 representing positive chemotaxis, for i = 1, 2, 3. Furthermore, for the transmission conditions for u at the internal node we set dissipative coefficients ξi,j = 13 for i, j = 1, 2, 3 and for the transmission conditions for φ we assume κi,j = 1 for i 6= j and κi,i = 0, for i, j = 1, 2, 3. Figure 1. Network of three edges: one incoming and two outgoing arcs (T-shaped network). 4.1.1. Example 1: the zero flux case. In this first example we set homogeneous Neumann boundary conditions at the outer boundaries for both u and φ, in order to have zero flux conditions as in [3]. The initial values for ui are randomly equally distributed in [0.25, 0.35] for each i = 1, 2, 3. Moreover, we set φ1 (x, 0) = φ3 (x, 0) = 0 and φ2 (x, 0) = 2, see Fig. 2. MAZE SOLVING ABILITY OF SLIME MOLD 17 Figure 2. Initial data for the test in Example 1. As shown in Fig. 1, the non-zero initial condition for chemoattractant on arc 2, leads the cells to move towards it. Since there the concentration of chemoattractant is largest, they accumulate at the right boundary of the arc, see the left picture of Fig. 3. Moreover, due to the diffusion flux, chemoattractant distributes on the other arcs, thus determining a fast decrease in the quantity of the chemoattractant on arc 2 as shown in the right picture of Fig. 3. The mentioned results are in accordance to the ones reported in [3]. We also computed the asymptotic solution achieved at time T = 43.5 and in Figg. 5-6 we reported the profiles of the final densities and fluxes. Figure 3. Example 1. The distribution of the density ui (x, t) (on the left) and of the chemoattractant φi (x, t) (on the right), for i = 1, 2, 3, inside the T-shaped graph at time t = 0.5. 18 G. BRETTI AND R. NATALINI Figure 4. Example 1. The distribution of the density ui (x, t), (on the left) and of the chemoattractant φi (x, t) (on the right), i = 1, 2, 3, inside the T-shaped graph at time t = 1. Figure 5. Example 1. The distribution of the density ui (x, t), (on the left) and of the chemoattractant φi (x, t) (on the right), i = 1, 2, 3, inside the T-shaped graph at time t = 43.5. 4.1.2. Example 2: non-homogeneous Neumann boundary conditions. Now, in order to reproduce the experiment with the T-shaped plasmodium tubes of Physarum presented in [25] to show the property of dead-end cutting, we set non-homogeneous Neumann boundary conditions at the outer boundaries for φ, to mimic the presence of food at the left and lower ends. For u we impose homogeneous Neumann boundary conditions at the outer boundaries as above. The initial values for ui are randomly equally distributed in [0.25, 0.35] for each i = 1, 2, 3 to mimic the plasmodium spread in the network, while we set φi (x, 0) = 0, i = 1, 2, 3. MAZE SOLVING ABILITY OF SLIME MOLD 19 Figure 6. Example 1. The asymptotic fluxes of the arcs of the T-shaped graph at time t = 43.5. Figure 7. Example 2. The distribution of the density ui (x, t) on each arc of the T-shaped network at time t = 0.5 (on the left) and t = 1 (on the right). Let us recall that in our modeling the dead-end cutting property corresponds to have a mass equal to zero on the egdes where the flux is null. In accordance to the cited experiment, as the organism starts moving to feed itself, we observe that the mass concentrate on the left and lower edges (arcs 1 and 2) since food sources are located at the outer nodes 0 and 2, and decreases quickly until it becomes null on the upper edge (arc 3), whose endpoint has no food source. 20 G. BRETTI AND R. NATALINI Figure 8. Example 2. The distribution of the density ui (x, t) on each arc of the T-shaped network at time t = 3.5 (on the left) and t = 7 (on the right). 4.2. Study on the path followed by amoeboid-like organism. Here we are interested in the identification of the shortest path followed by individuals in the case of positive (attractive) chemotaxis. 4.2.1. Wheatstone bridge-shaped network with a source and a sink (two exits). Let us consider the Wheatstone bridge-shaped network where, for our convenience, we inserted one incoming and one outgoing arc, respectively, exiting or entering an external node, as shown in Fig. 9. Since this modification does not change the structure of the shortest paths in the network, the analysis on the Wheatstone bridge network can be found in [18], where it was proven that the globally asymptotically stable equilibrium point corresponds to the shortest path connecting the exits. If, for instance, we want Figure 9. Wheatstone bridge-shaped network: a network with 7 arcs and 4 internal nodes. the path 1 → 2 → 6 → 4 → 7 to be shortest one, we need to set: (45) L2 + L6 < L5 and L4 + L6 < L3 . MAZE SOLVING ABILITY OF SLIME MOLD 21 In particular, we assume L1 = L7 = 0.2, L2 = L4 = 0.3, L3 = L5 = 2, L6 = 0.3. Note that in Fig. 9 we depicted the arcs composing the shortest path with the dashed line. √ We set parameters as ai = 1, bi = 0.1, λi = 0.33, χi = 1, Di = 1, for i = 1, . . . , 7. As before, for the transmission conditions for u, at each node p we set dissipative coefficients ξi,j = 31 for i, j = 1, . . . , 7 and for the transmission conditions for φ we assume κi,j = 1 for i 6= j and κi,i = 0, for i, j = 1, . . . , 7. The initial values for ui are randomly equally distributed in [0.45, 0.55] and we set φi (x, 0) = 0 for each i = 1, . . . , 7. We set non-homogeneous Neumann boundary conditions for the chemottractant φ at the node 0 and at the node 5 (inflow conditions) : (46) ∂x φ1 (0, t) = −1, ∂x φ7 (L7 , t) = 1. Firstly, we assume for u homogeneous Neumann boundary conditions ux (·, t) = 0 corresponding to v(·, t) = 0 at the node 0 and at the node 5. Figure 10. Wheatstone bridge-shaped network. The distribution of the density ui (x, t) on each arc of the diamond graph-like network at time 0.5 (on the left) and t = 1.5 (on the right). The effect on the movement of the organism is that it is able to find the shortest path connecting node 0 with node 5 in such a way to minimize the sum of the length of the arcs composing the path. As a result, see Figg. 10-11, the mass concentration is higher on the arcs composing the shortest path, as underlined by the colourscale and thickness of the line connecting nodes. Note that the solution obtained at time T = 6 is a stationary state. 22 G. BRETTI AND R. NATALINI Figure 11. Wheatstone bridge-shaped network. The distribution of the density ui (x, t) on each arc of the diamond graph-like network at time t = 2.5 (on the left) and t = 6 (on the right). Then, to have a certain quantity of cells enters into the network by the node 0 and the node 5, we set non-zero flux conditions at the outer boundaries (26) at the node 0 and (27) at the node 5: 2 2 , v7 (L7 , t) = − . v1 (0, t) = 1 + u1 (0, t) 1 + u7 (L7 , t) Figure 12. Wheatstone bridge-shaped network. The distribution of the density ui (x, t) on each arc of the diamond graph-like network at time 10 (on the left) and t = 50 (on the right). As in the previous case, cells migrate on the arcs on the shortest path connecting node 0 with node 5 in such a way to minimize the sum of the MAZE SOLVING ABILITY OF SLIME MOLD 23 Figure 13. Wheatstone bridge-shaped network. The distribution of the density ui (x, t) on each arc of the diamond graph-like network at time t = 130 (on the left) and t = 200 (on the right). length of the arcs composing the path. As shown in the following Figures 1213, the mass concentration is higher on the path of minimum total length. In this case the connection between the two exits in more evident, as underlined by the range of the density in the colourscale. 4.2.2. A network with two exits: the maze. As in [3], the network we consider here is a maze composed of 26 arcs and 18 nodes, with the exits placed at the node 0 (source) and at the node 17 (sink). We set the length Li = 0.5 on arcs i = 1, 5, 9, 10, 14, 21, 25, 26 and Li = 10 on the others. Note that in Figure 14 we represent such network, with the shortest arcs depicted with the straight line and the longest arcs with the√curve line. We set parameters as ai = 1, bi = 0.1, λi = 0.33, Di = 1, χi = 1, for all i, considering again positive chemotaxis. Furthermore, for the transmission conditions for u we set dissipative coefficients ξi,j , see Table 1, and for the transmission conditions on φ we assume again κi,j = 1 for i 6= j and κi,i = 0, for i, j = 1, . . . , 26. The initial values for ui are randomly equally distributed in [0.45, 0.55] and we assume φi (x, 0) = 0, ∀i. node p= 1, 2, 3, 5, 8, 9, 12, 14, 15, 16 ξi,j = 13 , i, j ∈ Ip ∪ Op node p= 4, 13 ξi,j = 12 , i, j ∈ Ip ∪ Op node p= 6, 7, 10, 11 ξi,j = 14 , i, j ∈ Ip ∪ Op Table 1. Transmission coefficients used for the numerical simulations of the maze in Fig. 14 given node by node. Here we set non-homogeneous Neumann boundary conditions at the outer boundaries for φ: inflow conditions leading to a high concentration of the 24 G. BRETTI AND R. NATALINI Figure 14. The maze: a network with 26 arcs and 18 nodes. chemoattractant at the specified point are assumed at the node 0 and at the node 17: (47) ∂x φ1 (0, t) = −1, ∂x φ26 (L26 , t). Furthermore, to have a certain quantity of cells enters into the network by the node 0 and the node 17, we set non-zero flux conditions at the outer boundaries: 2 2 v1 (0, t) = , v26 (L26 , t) = − . 1 + u1 (0, t) 1 + u26 (L26 , t) The effect on the movement of the cells is that they try to find the minimumlength way connecting node 0 with node 17 in order to minimize the sum of the length of the arcs composing the path. Let us recall that the deadend cutting of the zones with no-flux and the shortest path selection is here represented, as shown in the following Figures 15-16, by the higher mass concentration on the edges composing the shortest path. We also observe that at a certain time blow up of solutions verifies, due to the growth of the total mass. We remark that, though the evolution is slower respect to that presented in [3], the results are asymptotically the same obtained in [3]. The simulation of the network above was performed by a personal computer, processor Intel Centrino core 2 Duo 2 Ghz, RAM 3 Gb and the CPU time for 1333 time iterations was 19.71 s. 4.2.3. A network with more than two exits. Here we consider a network with more than two exits (NM2E): in fact in this case we have two sources and three sinks. Such network, which is composed of 21 arcs and 15 nodes, with the sources placed at the nodes 0 and 1, while the sinks are at the nodes 9, 10 and 11, we set the length Li = 0.5 on arcs i = 1, 2, 3, 4, 8, 13, 14, 15, 19, 20, 21 and Li = 10 on the others. The network we consider is shown in Figure 14, where the shortest arcs are depicted with the dashed line. MAZE SOLVING ABILITY OF SLIME MOLD 25 Figure 15. The maze. The distribution of the density ui (x, t) on each arc of the network in Fig. 14 at time t = 3 (on the left) and t = 12 (on the right). Figure 16. The maze. The distribution of the density ui (x, t) on each arc of the network in Fig. 14 at time t = 30 (on the left) and t = 60 (on the right). √ As before, we set parameters as ai = 1, bi = 0.1, λi = 0.33, Di = 1, χi = 1, for all i. Furthermore, for the transmission conditions for u, we set dissipative coefficients ξi,j , see Table 2, and for the transmission conditions on φ we assume κi,j = 1 for i 6= j and κi,i = 0, for i, j = 1, . . . , 21. We assume the inflow boundary conditions: (48) ∂x φ1 (0, t) = −1 = ∂x φ2 (0, t), ∂x φ19 (L19 , t) = 1 = ∂x φ20 (L20 , t) = ∂x φ21 (L21 , t). Moreover, we assume for u homogeneous Neumann boundary conditions ux (·, t) = 0 at the nodes 0, 1 and at the nodes 9, 10, 11. 26 G. BRETTI AND R. NATALINI Figure 17. Network of 21 arcs and 15 nodes (NM2E): the dashed line represents the shortest arcs. The initial values for ui are randomly equally distributed in [0.45, 0.55] and we assume φi (x, 0) = 0, ∀i. node p= 2, 3, 4, 6, 9, 10, 11 ξi,j = 13 , i, j ∈ Ip ∪ Op node p= 7, 8 ξi,j = 16 , i, j ∈ Ip ∪ Op node p= 5 ξi,j = 14 , i, j ∈ Ip ∪ Op Table 2. Transmission coefficients used for the numerical simulations of network NM2E in Fig. 14 given node by node. Figure 18. The distribution of the density ui (x, t) on each arc of the network NM2E in Fig. 17 at time t = 10 (on the left) and t = 30 (on the right). Even in this case, the organism traces the path connecting nodes 0 and 1 with nodes 19, 20, 21 in such a way to minimize the sum of the total length of the arcs. As a result, see the following Figures 18-19, besides the arcs connected to the exits (1, 2, 19, 20, 21), the path includes the arcs with minimum length 3, 4, 8, 13, 14, 15. MAZE SOLVING ABILITY OF SLIME MOLD 27 Figure 19. The distribution of the density ui (x, t) on each arc of the network NM2E in Fig. 17 at time t = 60 (on the left) and t = 180 (on the right). References [1] V. Bonifaci, K. Mehlhorn and G. Varma. Physarum can compute the shortest path. Journal of Theoretical Biology, 309 (2012), 121-133. [2] S. Borsche, A. Klar and T.N.H. Pham Kinetic and related macroscopic models for chemotaxis on networks preprint ArXiv: 1512.07001v1 (2015). [3] S. Borsche, S. G ottlich, A. Klar and P. Schillen. The scalar Keller-Segel model on networks. Mathematical Models and Methods in Applied Sciences, 24(2), (2014), 221– 247. [4] G. Bretti, R. Natalini and M.Ribot. A hyperbolic model of chemotaxis on a network: a numerical study. Mathematical Modelling and Numerical Analysis, 48(1), (2014), 231–258. [5] F. Camilli and L. Corrias. Parabolic models for chemotaxis on weighted networks. preprint arXiv:1511.07279 (2015). [6] Y. Dolak and T. Hillen. Cattaneo models for chemosensitive movement. Numerical solution and pattern formation, J. Math. Biol., 46 (2003), 153–170; corrected version after misprinted p.160 in J. Math. Biol., 46 (2003), 461–478. [7] F. Filbet, P. Laurençot, and B. Perthame. Derivation of hyperbolic models for chemosensitive movement. J. Math. Biol., 50(2) (2005), 189–207. [8] A. Gamba, D. Ambrosi, A. Coniglio, A. de Candia, S. Di Talia, E. Giraudo, G. Serini, L. Preziosi, and F. Bussolino. Percolation, morphogenesis, and Burgers dynamics in blood vessels formation, Phys. Rev. Letters, 90 (2003), 118101.1–118101.4. [9] M.Garavello,B.Piccoli, Traffic flow on networks. Conservation laws models. AIMS Series on Applied Mathematics, 1. American Institute of Mathematical Sciences (AIMS), Springfield, MO, 2006. xvi+243 pp. ISBN: 978-1-60133-000-0; 1-60133-000-6 [10] L. Gosse. Asymptotic-preserving and well-balanced schemes for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes, J. Math. Anal. Appl., 388(2) (2012), 964–983. [11] L. Gosse. Well-balanced numerical approximations display asymptotic decay toward Maxwellian distributions for a model of chemotaxis in a bounded interval, SIAM J. Sci. Comput., 34(1) (2012), A520–A545. 28 G. BRETTI AND R. NATALINI [12] J.M. Greenberg and W. Alt. Stability results for a diffusion equation with functional drift approximating a chemotaxis model, Trans. Amer. Math. Soc., 300 (1987), 235– 258. [13] F.R.Guarguaglini, C.Mascia, R.Natalini, M.Ribot, Stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis, Discrete Contin.Dyn.Syst.Ser.B, 12 (2009), 39-76. [14] F.R. Guarguaglini and R. Natalini. Global smooth solutions for a hyperbolic chemotaxis model on a network. SIAM J. Math. Anal., in press. [15] D. Horstmann. From 1970 until present: the Keller-Segel model in chemotaxis and its consequences, I. Jahresber. Deutsch. Math.-Verein., 105 (2003), 103–165. [16] E. F. Keller and L.A. Segel. Initiation of slime mold aggregation viewed as an instability. J. Theor. Biol., 26 (1970), 399–415. [17] O. Kedem and A. Katchalsky. Thermodynamic analysis of the permeability of biological membrane to non-electrolytes. Biochimica et Biophysica Acta, 27 (1958), 229–246. [18] T. Miyaji and I. Ohnishi. Mathematical analysis to an adaptive network of the Plasmodium system, Hokkaido Math. J., 36, (2007), 445–465. [19] T. Miyaji and I. Ohnishi. Physarum can solve the shortest path problem on Riemannian surface mathematically rigorously. Int. J. Pure Appl. Math., 46, (2008), 353369. [20] J. D. Murray. Mathematical biology. I. An introduction. Third edition. Interdisciplinary Applied Mathematics, 17. Springer-Verlag, New York, 2002; II. Spatial models and biomedical applications. Third edition. Interdisciplinary Applied Mathematics, 18. Springer-Verlag, New York, 2003. [21] T. Nakagaki, H. Yamada and A. Tóth. Maze-solving by an amoeboid organism. Nature, 407, (2000), 470. [22] R. Natalini and M.Ribot. An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis. SIAM Journal of Numerical Analysis, 50 (2012), 883–905. [23] B. Perthame. Transport equations in biology, Frontiers in Mathematics, Birkhäuser, 2007. [24] A. Tero, R. Kobayashi, and T. Nakagaki. A coupled-oscilator model with a conservation law for the rhythmic amoeboid movement of plasmodial slime molds, Phys. D., 205, (2005), 125-135. [25] A. Tero, R. Kobayashi, and T. Nakagaki. A mathematical model for adaptive transport network in path finding by true slime mold, J. Theoret. Biol., 244, (2007), 553-564. [26] A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker, K. Yumiki, R. Kobayashi, T. Nakagaki. Rules for Biologically Inspired Adaptive Network Design, Science, 327, (2010), 439-442. [27] L. A. Segel. A theoretical study of receptor mechanisms in bacterial chemotaxis. SIAM J. Appl. Math. 32 (1977), 653–665. [28] S. Watanabe, A. Tero, A. Takamatsu and T. Nakagaki. Traffic optimization in railroad networks using an algorithm mimicking an amoeba-like organism, Physarum plasmodium, ByoSystems, 105, (2011), 225-232. [29] https://youtu.be/czk4xgdhdY4.