DEGREE-BIASED ADVECTION–DIFFUSION ON UNDIRECTED GRAPHS/NETWORKS

. There are several phenomena in nature governed by simultaneous or intermittent diﬀusion and advection processes. Many of these systems are networked by their own nature. Here we propose a degree-biased advection processes to undirected networks. For that purpose we deﬁne and study the degree-biased advection operator. We then develop a degree-biased advection–diﬀusion equation on networks and study its general properties. We give computational evidence of the utility of this new model by studying artiﬁcial graphs as well as a real-life patched landscape network in southern Madagascar. In the last case we show that the foraging movement of the species L. catta in this environment occurs mainly in a diﬀusive way with important contributions of advective motions in agreement with previous empirical observations


Introduction
Advection/diffusion processes appear very frequently in many geological, biological and man-made systems [26].This kind of processes are characterized by having a quantity that tends to expand evenly in all directions (diffusion) while it is affected by an external force, that drifts the expansion of that quantity in some preferred directions (advection).In the case of ecological systems, for instance, such advective phenomenon describes motion of species towards good resources like nutrients or against noxious stimuli like predators.This combination of advective and diffusive motions could end up with species moving across the whole system with higher concentrations around the nutrients and/or lower concentrations near the predators.
In geology, for instance, when there is a layered permeable system with flow predominately parallel to the bedding, the transport of solute through the more permeable layers is dominated by advection, while in less permeable layers molecular diffusion controls the transport [20].Similar processes governs also other geological and man-made systems such as transport in rivers and drainage networks [1].
In biological systems, including fungi, organs, vascular networks and plants, the transport of resources for metabolism and growth is frequently carried by advective and diffusive processes through different kinds of transport networks [2,4,21,22,28,30,34,39]. Advection-diffusion is also fundamental for understanding the growing processes of tumours where nutrients and oxygen in the microenvironment are transported by these mechanisms to support tumour cells and stimulates their proliferation and migration [29,35,36].
Another scenario in which advection-diffusion plays a fundamental role is in animal foraging and dispersal [10].Many animals use combined local diffusion with longer range advective motions to change from unfavorable parts of their environment to more suitable areas [24].The advective motion can emerge either from the behavior of individuals or from physical transport processes, such as winds, currents in rivers, and so on.The combined advection-diffusion strategy of foraging has been shown to be advantageous for species [17].That is, foraging strategy is improved when animals change their behaviour from diffusion to advection-diffusion and back depending on their surrounding space [6].Moreover, it has been shown that populations using strategies with no advection or high advection can be invaded by other populations [5].Therefore, according to these results, strategies with no advection or high advection would be predicted to be unstable relative to evolution.
A common characteristic of many of the complex systems where advection-diffusion is ubiquitous is that they have a networked nature.That is, in many of these systems there are discrete entities which are interconnected to each other forming a complex network/graph.They range from molecular and cellular networks to landscape networks (see [12] for many examples).In the case of networks, the studies about advection-diffusion have mainly focused on directed graphs.In these studies the advection is created by the directionality of the edges in the graph, which is accounted for a divergence operator giving rise to an advection Laplacian operator [7-9, 27, 33].The problem is that many complex networks are undirected.For instance, in a landscape network formed by regions of a landscape connected by corridors, there is no sign indicating prohibited directions.However, in many real-life situations the own characteristics of the nodes may represent a preference of movement towards these points.The simplest scenario is to consider that the number of connections that a node has-its degree-is a proxy for attracting or repelling an object towards a given node.Such degree-biased motion paradigm has been under development by our group in recent works [14,15,18].Here we extend this approach to account for degreebiased advection-diffusion on undirected graphs.Our approach is mainly a mathematical one, proving some of the general and more important properties of the degree-biased advection operator as well as of the degreebiased advection-diffusion equation.Additionally, we provide some computational evidence on the applicability of this model to the study of foraging of L. catta in a patched landscape in southern Madagascar.

Definition and set up
Let G = (V, E) be a simple undirected graph with |E| = m edges and |V | = n vertices, none of them isolated.Let A be the adjacency matrix of the graph G and let k v denote the degree of the node v ∈ V , so the diagonal matrix of node degrees is defined as K = diag A • 1 .Denote N v as the set of neighbours of the node v, which means the set of nodes w for which the edge (v, w) ∈ E exists.Finally, let 2 (V ) be the Hilbert space of square summable functions on V .Definition 2.1.Let α ∈ {−1, 0, 1} and define the degree-biased Laplacian corresponding to α as the operator on 2 (V ) given by where If α = 0, then c α = 1 and L 0 = L, the standard Laplacian, see [31] If α = 1, then c α (v, w) = kw kv and we call L 1 the hubs-attracting Laplacian, see [18].If α = −1, then c α (v, w) = kv kw and we call L −1 the hubs-repelling Laplacian, see [14].Set {e v } v∈V the standard 2 (V ) basis and define δ α (v) as (2.2) Then, we can see how L α acts on the vectors e v as: Nevertheless, this degree-biased Laplacian inherits some properties from its non-biased predecessors, proven in several previous articles [14,18]: If G is a locally-finite graph, then L α is a bounded operator.The spectra of this type of operators is real and 0 is an eigenvalue with associated eigenvector 1.By defining W α = diag(K α AK −α • 1), we can decompose the operators as: (2.4) Remark 2.2.Equation (2.4) allows us to extend easily the concepts in [14,15,18] for any α ∈ R, modifying the weight of node differences in the network.However, we will restrict ourselves to the ternary case α ∈ {−1, 0, 1}, as it is the sign of α which dictates the qualitative behaviour of the operator.

The adjoint of L α is an advection operator
We start here by defining the adjoint of L α : Definition 3.1.For a graph G without isolated nodes, we define the Adjoint degree-biased Laplacian as: It can be easily verified that and that in the case of finite graphs (L α ) * = (L α ) .
We prove now that this new operator on graphs corresponds to a degree-biased advection operator on undirected graphs.Theorem 3.2.The Adjoint operator of the degree-biased Laplacian, (L α ) * , is an advection operator where the time-invariant velocity field is biased by the degree of the nodes which create a flow in the direction of high or low degree nodes according to the respective value of α.
Before proceeding with the proof we should remind the concept of advection in the continuous setting: where ∇• is the divergence operator, u (t) is a conserved scalar quantity moving through a time-invariant velocity field v on the domain of u (t).We now will introduce a few concepts which will make possible to identify the operator (L α ) * with a degreebiased advection operator on undirected graphs.The necessity for using the double orientation of the graph emerges from the following.Here, we are considering undirected graphs which are then transformed such that every existing edge is replaced by weighted arcs pointing in both directions.Those weights are not necessarily the same in both direction, which oblige us to consider the double-orientation explicitly.In this context, the incidence matrix indicates how much of a quantity departs from a node to its nearest neighbor, as well as how much arrives at it.This obviously creates a gradient between every pair of connected nodes of the original graph.Note that this notation maintains the main property of the classical Laplacian with respect to the incidence matrix, that is L 0 = D * • D. Definition 3.6.Let G ↔ be a double-orientation of G.Then, the time-invariant velocity field V ∈ R m×n on the domain of u (t) is defined as We now proceed to prove Theorem 3.2: For that, recall that the product matrix can be written entry by entry as: Suppose i = j, then the product (D * ) il V lj = 0 if and only if the edge l goes from i to j, as the first term is non-zero if and only if l is connected to i and the second one is different from zero if l goes towards j.
Hence, there is only one term in (D * • V ) ij that is different from 0 on the edge l described before, leading to the calculation: On the other hand, for i = j, D * will be non-zero for each edge going towards i.This means that this term will only be equal to 1, so the calculation will be: Comparing this term with those produced in (3.1) for g = e i and v = j leads to the desired equality.
As the incidence matrix is the equivalent to the gradient operator in the graph setting, the previous decomposition shows that (L α ) * is the operator representing the flux of the "vector field" V , so its evolution equation is equivalent to that in equation (3.3).Remark 3.7.It should be remarked here that in previous works the only attempts to study advection processes on graphs have been on directed graphs.This was an obvious necessity due to the fact that a gradient was needed over the edges of the graph.However, here this gradient is created by the heterogeneities in the degrees of the nodes.
Before continue with the study of dynamical properties on graphs based on this new operator let us consider some of its more general spectral properties.
Proof.As it is the transpose of the degree-biased Laplacian, we can use its properties to write it as: Now, using the definition of W α and A, we can say that this Adjoint operator is similar to a symmetric matrix, thus, have real spectrum.
Proposition 3.9.The 0-eigenvector of the operator where k is the vector of degrees of each node.
Proof.Without any loss of generality we consider connected graphs.In case of disconnected ones we can consider the following results on a component-wise way, which is standard in network theory.
Let v be the zero-eigenvector.To find its expression, we must solve the system: Using the previous decomposition, we are going to find u := K α v so the system is now related to the matrix W α − A. Hence, we have to solve the system of equations: This can be written as: which has u j = (k j ) −α as solution.

Degree-biased advection equation on graphs
We start here by considering the degree-biased advection equation of undirected graphs: with initial conditions u(0) = u 0 .
To analyze the convergence of this system we first introduce the following.
We say that the process (3.7) converges towards the v set if: We now prove such convergence.
Theorem 3.11.Let G be any connected graph with n nodes, none of them isolated.Let k = (k i ) n i=1 be the vector of degrees of each node and define v = k −2α .Then, the dynamical system (3.7)converges towards the v set for any initial condition u 0 and the rate of convergence is given by the smallest nontrivial eigenvalue of (L α ) * , µ 2 ((L α ) * ).
Proof.The solution of the dynamical system (3.7) is given by: Using the previous decomposition of the adjoint operator and the fact that the matrix W α − A is symmetric, we can find a matrix U of eigenvectors such that where Λ is the diagonal matrix of eigenvalues of the adjoint operator.
If we take ψ j and φ j the j-th column of K −α U and (K −α U ) −1 , respectively, we can rewrite the solution as where the eigenvalues are ordered in a non-decreasing way: On the one hand, as µ 1 is the unique 0 eigenvalue, we can take the limit t → ∞ to obtain As we know that ψ 1 is an eigenvector associated to µ 1 , we can take it equal to v as shown in Proposition 3.9.Thus, we have: Moreover, the solution yields: On the other hand, we have that the second smallest eigenvalue is the slowest converging to 0 in the series of the exponential, hence dictating the rate of convergence.This result gives us the stationary state of a system governed by the evolution equation (3.7), but transient dynamics can be specific to particular graphs due to the fact that all nonzero eigenvalues of the operator would affect the dynamics.Consequently, the dynamics strongly depend on the structure of the network.

Advection-diffusion equation
For application to real-world problems, advection appears as a current leading a particle which is diffusing at the same time.For this reason we couple both phenomena into a degree-biased advection-diffusion equation on undirected graphs.Definition 4.1.We define the advection-diffusion equation as the differential equation: where u 0 is an initial condition, L a diffusion operator, which we will take as the usual standard Laplacian, and γ adv , γ dif are two constants modeling which action (diffusion or advection) acts stronger.We now show some properties of the operator (4.1).
Lemma 4.2.0 is the lowest eigenvalue of the operator (4.1).
Proof.First we have to see that the operator is positive semidefinite, which follows from the positive semidefiniteness of the Laplacians used and that both constants are non-negative: On the other hand, using that the spectrum of any matrix is equal to the spectrum of its transpose σ(A) = σ(A T ) , and using that the transpose of the Adjoint degree-biased Laplacian is the standard degree-biased Laplacian itself, we can calculate: T ) which shows that 0 is the lowest eigenvalue of the advection-diffusion operator.
Remark 4.3.Note that, as 1 is an eigenvector of the transpose, it directly becomes a left eigenvector of (γ dif L + γ adv L * α ).Following Theorem 3.11, we can see that the dynamical system controlled by the degree-biased advectiondiffusion operator converges for any γ adv , γ dif > 0 to a stationary state that depends on the eigenvector associated to the zero eigenvalue.We will prove it explicitly later, but first we are going to show the equation that this eigenvector should solve.
for any i ∈ V .
Proof.The vector v should solve the equation ( As we know the two terms forming the operators, this leads, for each i ∈ V , to the equation: Joining the terms of the equation for v i on the one hand and for all the v j 's on the other, leads to (4.2).Note that, using k i = j∈Ni 1, we can also reorder the equation to obtain: This shows that solutions behave continuously to the previously known cases, that are γ adv = 0, where the operator is the standard Laplacian and v = 1, and γ dif = 0, where v = k −2α , as shown in Proposition 3.9.This zero-eigenvector is far from trivial compared to the same eigenvector for both the classical Laplacian and the Adjoint hub-biased Laplacian.Before proving the general result we will need the following Lemma which gives an expression for this eigenvector for trees.
i=1 is the eigenvector of the operator γ adv L * α + γ dif L associated to the zero eigenvalue.Then, for any node i ∈ V , the coordinates of the vector are: where S(i, o) is the sequence of edges forming the shortest path between node i and node o, which is an objective node that can be chosen arbitrarily for each tree G.
For the sake of clarity, call For our objective, we have to prove the following equality: It suffices to show that they are equal addend by addend.For that purpose, consider that either S(i, o) passes through j or S(j, o) passes through i, which is provided by G being a tree.Consider the first case, so we can decompose S(i, o) = S(i, j) ∪ S(j, o), and note that P (S(i, j)) has only one term involving the two adjacent nodes.
Then, we do the calculation: Using this calculation on the left term, we arrive to the desired equality: If, on the other hand, it is S(j, o) the one that goes through i, we have to do the same calculation on the other term to achieve the equality.
As we are aiming for a more general result, we need to further understand the expression in equation (4.4), so we will manipulate it in the next Lemma to have a deeper understanding about what do the terms in the multiplication represent.
Lemma 4.6.The expression (4.4) for v i can be rewritten as where T i are all the paths in the tree departing from the node i to the leaves without repeating edges.
Proof.To prove it, it suffices to choose the objective node of (4.4) and divide each term in (4.5) by the term v o .Then, for any node i, the only different terms between T i and T o are those in the path S(i, o), so the operation will lead to a multiplication of the terms: and the terms k α l are canceled with the terms k −α j of the following term.This leaves us with a term k α o corresponding to the objective node o, which will be common to all the entries of v and thus, can be omitted and still represent the same eigenvector, and a term k −α i that already appeared in formula (4.4).Thus, equations (4.4) and (4.5) are equal.
The main goal of obtaining this result is not for the calculation of the corresponding eigenvector v i , for which many excellent numerical methods exist.Our main goal is to use it in terms its structural interpretation in order to establish structure-dynamics relations.
For that purpose, let us further explain what is the meaning of this result.Let us consider a node i in the bidirected graph G ↔ created from a tree graph G as described before.A directed walk on G ↔ is a sequence of edges (e 1 , e 2 , . . ., e h−1 ) for which there is a sequence of nodes (v 1 , v 2 , . . ., v h ) such that if e k and e k+1 are consecutive in the sequence, then e k = (v r , v s ) and e k+1 = (v s , v l ).A directed path is such a directed walk where all edges and vertices are distinct.A directed path departing from the node i is a directed path in which v 1 = v i .Let e l be an edge in a directed path departing from the node i, and let w l be the weight of this edge created as described before.Hence, equation (4.5) indicates that the ith component of the eigenvector corresponding to µ 1 = 0 is the product of the weights of all edges in the tree that are at least in one directed path departing from the node i.Let us consider V SPi to be the sum of the weights of all edges in at least on directed path departing from the node i, i.e., the volume of the corresponding directed tree.Let p l = w l /V SPi be the probability that the quantity moving across the network via a diffusive-advective process is located at the edge e l .Let us consider that the probabilities p i and p j are independent.Then, the total probability P i of finding a quantity that departed from the node i to any other node via advection-diffusion is given by the product of all the edge probabilities: P i = p l , such that v i ∝ P i .In plain words, the ith component of the eigenvector corresponding to µ 1 = 0 is proportional to the probability that a quantity leaves the node i.
We are now in position to prove the general expression for this eigenvector of the degree-biased advectiondiffusion operator based on the formula stated in [7] which is given in the following result.
Theorem 4.7.Let G be any graph.Suppose v = (v i ) n i=1 is the eigenvector of the operator γ adv L * α + γ dif L associated to the zero eigenvalue.Then, for any node i ∈ V , the coordinates of the vector are: where T (G) is the set of all spanning trees of the graph G and T i are the paths of each of those spanning trees departing from the node i and ending on the leaves without repeating edges.
Proof.As in the previous proof, we have to prove that j∈Ni Ti∈T (G) is the same if we change j and i.
Note that each tree in T (G) is the same in each term but for its root node, so it is enough to see that the terms are equal in each tree.Moreover, we can apply the previous theorem as we fulfil the assumption that the graph we are working on is a tree now.But, as we have changed the way of writing the eigenvector, we can offer a different proof to why both terms are equal for each tree.Consider one leaf of T i .Then, either the path passes through node j or it doesn't.If j is in the path, then the same path is in T j but for the term , which is exactly the term multiplying v j in (4.2), so it would appear in both cases.On the other hand, if the path does not pass through j, it happens exactly the opposite to the previous case, because then that path will be in T j but for a term γ dif + γ adv ki kj α , which appears multiplying v i in (4.2).
We can interpret this result for general graphs in the same way as we did for trees, using probability.As before, consider a quantity departing from the node i.Now, it can reach the other nodes following many different paths, one for each spanning tree of the graph.Hence, we need to use the total probability theorem to calculate P i .
A spanning tree departing from node i is a spanning tree in which, for every node l = i, there is a unique directed path departing from node i and ending in l.Let {SP i (j)} be the set of all spanning trees departing from i of G, which has the same size as T (G).For each j, we already know that the probability that a quantity departs from node i in the tree SP i (j) is the product l p l (j).The total probability theorem reads then as P i = j p(SP i (j)) l p l (j), where p(SP i (j)) is the probability of choosing SP i (j) as the spanning tree to follow.There are no reasons to choose one spanning tree over another so they are equiprobable.Thus, p l (j).This means that the ith component of the eigenvector corresponding to the 0 eigenvalue keeps being proportional to the probability that a quantity moving through the network via advection-diffusion departs from node i. Knowing this expression for the first eigenvector of the degree-biased advection-diffusion equation, we can proceed as in Theorem 3.11 to get the dynamics of the advection-diffusion equation for long times.
Theorem 4.8.Let G be any connected graph with n nodes, none of them isolated.Define v as in (4.6).Then, the dynamical system ruled by the advection-diffusion equation (4.1) converges towards the v set for any initial condition u 0 .
The convergence of the solution towards the stationary state depends directly on its second eigenvalue, and we can use simple estimates to show that this advection-diffusion case converges faster than the pure diffusion or the pure advection cases separately.Proposition 4.9.Let (L α ) * be the Adjoint biased Laplacian and L be a diffusion operator.Then, we have the bounds: Moreover, if L is the standard graph Laplacian, we have the bound: where δ and ∆ are the minimal and maximum degree of the graph.
Proof.Using Weyl's dual inequality, we get: Now, using that the first eigenvalue of the Laplacian is 0, the left term is simplified.Also, as γ adv is a fixed constant, we can just extract it multiplying.
The second equation is achieved applying the same steps to the other operator in the sum, as 0 is also an eigenvalue of the Adjoint biased Laplacian.
For the case where L is the Laplacian, proceed as before and then use the inequality 3.35 in [15]: to obtain the desired inequality.
As the second eigenvalue rules the convergence of solutions to the stationary state for the advection-diffusion, pure advection and pure diffusion operators, this result shows that solutions of an advection-diffusion system will converge towards its stationary state faster than in the other two cases.

Computational results
In this section we perform computational studies that put in practice the concepts developed in the previous part of this work.For that purpose, we start by studying some artificial graphs which allow us to visualize the main features of the flows of items under the new advection-diffusion dynamical system developed here.After this warming up, we develop a "proof-of-concept" example of this model in a real-life scenario.For this example we selected a landscape network, consisting of patches and corridors which allow the movement of a species through a vast geographic area.The selection of this example is twofold.On one hand, we have empirical data that allow us a "reverse engineering" exercise to estimate the parameters of the model.On the other hand, the foraging movement of species in a landscape is a classical example in which diffusive and advective motions are combined, giving us a realistic scenario to test our model.

System evolution on artificial networks
We can numerically show that Proposition 3.9 and Theorem 4.7 correctly represent the steady state of advective and diffusive-advective systems, respectively.For that purpose, we generated an small Erdős-Rényi graph with 8 nodes and 15 vertices [11].
In Figure 1 we illustrate the structure of this network, where it can be seen that although it is very small it contains enough structure-it is far from tree-like-as to illustrate the effects of advection and diffusion on a realistic networked scenario.On this graph, we solved numerically the dynamical system u = − Lu with random initial conditions for L = 0.5L + 0.5(L α ) * and L = (L α ) * .In Figure 2 we plotted the evolution of the solution along time (dashed lines) where each color represents a group of nodes with the same degree.At the same time, we computed the expressions of the 0 eigenvector v with the formulas given in Proposition 3.9 and Theorem 4.7.We plotted the values of v in the same figure as horizontal solid lines, whose color represents the degree of the nodes that should approach this steady state.As we can see in Figure 2, the analytical results perfectly fit the steady state of the system in all the cases.Moreover, it shows that the convergence of the advection-diffusion system (panels (B) and (D)) is faster than the convergence of the fully advective system (panels (A) and (C)), which is a consequence of Proposition 4.9.

In-and out-flows on networks
An important problem for the analysis of real-world networks is that, in general, we are able to "observe" them at a given fixed time t, or at best at a few fixed times t 1 , . . ., t h .In the case of a dynamics taken place on such networks we are then capable to "observe" the internode flow which is taking place only between nearest neighbors.Let us think for instance in the case of animals of a given species moving in a landscape network, which will be analyzed in the next subsection.If we take a photograph of the landscape at a give time t we will observe the animals located at every node and those moving between pairs of nearest neighbor patches.However, for those already in a given patch (node) we can say nothing about their previous states.That is, we cannot say whether they have came from a nearest neighbor, from a distant node, or even if they have remained at that patch for their entire life.On the contrary, for those animals at the edges of the graph we can infer that at a time τ > t they will almost surely arrive at the nearest neighbor nodes toward which they were moving at time t.
When we apply our mathematical models to a dynamical process taken place on a graph, supposing that we have a good estimation of the initial conditions, we are able to obtain a global picture of the movement of "items" across the whole network.For instance, if our network is p−q −r we can see that u p,r (t > 0) > 0 for certain time t.This indicates that some "items" initially located at node p have moved to node r at time t.But in a network such movement can only be realized through the edges of the graph (notice there is no edge between p and r), such that this "item" has had first to go from p to q and then from q to r.In our static photograph of the system at time t we possibly will observe such "item" in its transfer from q to r, together with those that has started their travel at the node q.Thus, what we really observe is the in-and out-flows of items to/from the nodes coming/going from/to its nearest neighbors.
The previously described situation prompts us to consider how can we compare the results of our computational experiments with those that can be observed at a given time in a network.Here, of course we consider a process controlled by the advection-diffusion Laplacian L (γ adv ) := γ adv L * α + γ dif L when γ adv = 1 − γ dif .At a given time t the flow between a pair of nodes v and w is given by (Γ (t, γ adv ) Therefore, the snapshot of the network flow at the time t is obtained by the matrix exp −t L A, where denotes the entrywise (aka Schur or Hadamard) product of the two matrices.In the previous example of the network p−q −r, this will always give values only for the pairs p − q and q − r, and not for p − r, where an edge does not exist.Consequently, we will focus here on the in-and out-flows of "items" to/from the nodes of the network.It is important to note here that we are using the standard graph terminology for the adjacency matrix of directed graphs [40], such that the (i, j) entry of A indicates that there is an arrow from i to j.Then, the in-and out flows of our operators are defined, respectively, as fin (t, γ adv ) := 1 T (Γ (t, γ adv ) A) , (5.1) where 1 is the all-ones column vector.These parameters can be compared with the "experimentally" observed in-and out-flows of items in the network at a given time as we will do in the next section.Here, we first will analyze these parameters in some artificially created networks.

The case of a superconcentrator (star graph)
We start by considering a star graph-a graph of n nodes with a central node connected to n − 1 nodes of degree one.The central node acts as a superconcentrator in term of the in-flow arriving at it from the rest of the nodes.A nice property of this graph is that fin and fout can be considered as independent of time.That is, let L = L (γ adv ) be the advection-diffusion Laplacian of a star graph with n nodes labeled such that the central node is the number 1. Let L = V ΛV −1 such that ψ j corresponds to the jth column of V and φ j to the jth column of V −1 .We can start by writing exp −t L where For instance, for n = 10 this terms is 4.54 • 10 −5 even for t = 1.Therefore, for sufficiently large n the term exp −t L vw is independent of t for all (v, w) ∈ E.
Let us focus on fin (i) where i is the central hub or the leaf of degree one and set α = −1.In this case: (5.5) In the purely diffusive regime, i.e., γ adv = 0, we have a vector fin (0) in which the entry corresponding to the hub, lets say node 1, is (n − 1) /n and the n − 1 resting entries are n −1 .Then, for γ adv > 0 we compare the fin,i (γ adv ) vector with the purely diffusive one.In order to compare both vectors we use the Pearson correlation coefficient.In Figure 3 we illustrate the results.The Pearson correlation coefficient displays an inversion of its sign from positive correlation in the diffusive case to anticorrelation in the advective one.We explain now the interpretation of these results.
In the purely diffusive regime we have: exp −t L .This term is also added to 1 n in the leaves (see Eq. (5.5)).Therefore fin,i (γ adv < 0.5) is similar to fin,i (γ adv = 0) as can be seen by the values of the cosine similarity, and more importantly they are linearly correlated because the hub still concentrates most of the in-flow, while the leaves have very little in-flow.In this regime, although γ adv > 0, we are still in a scenario mainly dominated by diffusion.However, when γ adv = 0.5 we have: fin,hub (γ adv = 1/2) = fin,leaf (γ adv = 1/2) = 1/2, which means that at this point the in-flow to the hub and to the leaves is exactly the same.From this point on we have an inversion in the in-flow arriving at the leaves in comparison with the hub.That is, when γ adv > 0.5 the in-flow arriving at the leaves in bigger than that arriving at the hub.Therefore, the similarity between the vectors fin (γ adv = 0) and fin (γ adv > 0.5) drops significantly as observed in Figure 3, and more importantly there is an anticorrelation between the the two in-flow vectors.This means that we are in a new regime where advection dominates, pushing the "items" toward the leaves of the star in spite of the counter-effect of diffusion which tries to push them to the hub.The same analysis can be made for the choice α = 1, but the results are less interesting to study due to the fact that advection pushes the "items" towards the hub, thus increasing the effect of diffusion in the in-flow and resulting in a high correlation between the fin (γ adv ) for all values of γ adv .

Influence of assortativity
The example for the star graph that we analyzed in the previous subsection asks the question about which structural characteristic of a network is responsible for the inversion of the Pearson correlation coefficient of fin in mainly diffusive to mainly advective processes.The star graph is characterized by a high degree heterogeneityone node of degree n − 1 and n − 1 nodes of degree one-as well as by a large degree disassortativity-the degree-degree correlation coefficient is negative.Intuitively, the degree assortativity seems more relevant for the advection-diffusion inversion.The reason is that fin takes into account the in-flow to a given node from its nearest neighbors.Thus, the kind of degree-degree connectivity pattern should make an important impact on the values of fin .
Here we consider Erdős-Rényi random graphs which have low degree heterogeneity, i.e., they have a normallike Poisson degree distribution.We then rewire the edges of these graphs maintaining the node degrees, such that we can tune their degree assortativity.The degree assortativity is measured by the assortativity coefficient, i.e., the Pearson coefficient of the degree-degree correlation, see [13,32].For every type of graph with a given assortativity we calculate the correlation coefficient between the vectors fin (γ adv = 0) and fin (γ adv > 0).The results, after the average of 10 random realizations per assortativity rewiring, are illustrated in Figure 4.As can be seen when γ adv is zero or close to zero, i.e., in the mainly-diffusive regime, the nodes with the largest degree concentrate most of the in-flow, while the low-degree nodes have low values of it.This produces high positive correlation coefficients between fin (γ adv = 0) and fin (γ adv > 0) in this region as observed in Figure 4.Note that, as advection when α = −1 tends to concentrate more flow on nodes with lower degree and diffusion tends to concentrate more flow on high degree nodes (as they are the more connected ones), the correlation between this in-flows monotonically decreases as γ adv .On the other hand, the most interesting thing is that this high correlation is irrespective of the degree assortativity of the network.By contrast, when γ adv > 0.5 we observe a large range of correlations between fin (γ adv = 0) and fin (γ adv > 0), which range from values close to 1 for highly assortative networks to values close to −1 for highly disassortative networks.That is, when there is a large contribution of advection in the diffusive-advective process, the degree-degree correlation of the network can change the in-flow of items from hubs (in the assortative network) to low-degree nodes (in the disassortative network).This transition does not behave in the same way as when we are changing the advection parameter.For instance, if we fix a large value γ adv , we can observe that the correlation is not monotone with respect to the assortativity index, which may indicate that other features of the network also affect the in-flows.Indeed, what we are reporting here is the existence of two completely different kinds of processes represented by exactly the same dynamical equation (including its parameters) for both classes of networks, i.e., assortative and disassortative ones.

Foraging in a patched landscape
As described in the Introduction, advection-diffusion processes are ubiquitous in biological, geophysical and technological systems.Here we perform the analysis of diffusive-advective dynamics on a landscape network.As described below this landscape network represents a good scenario for testing advective-diffusive movement of a species through the patches and corridors of a geographic region.Therefore, our goal here is to use this landscape network as a proof-of-concept for the new operators and dynamical systems developed in the current work on a realistic scenario.Consequently, we are not interested here in a microscopic analysis of the movement of a species across the landscape.Instead, we are interested in showing quantitatively how the newly developed operators capture some important ecological aspects of species movement which escape to the use of the standard Laplacian operator to model this kind of scenarios.This landscape network is a graph G = (V, E) where the nodes represent 161 small forest patches (with areas between l and 95 ha) that are scattered in an agricultural landscape in the Androy region, southern Madagascar [3] (see Fig. 5).This landscape network is used by the ring-tailed lemur (Lemur catta) to move across the environment.We consider here that this species forages at the patches of the landscape which contain fruits forming the basis of L. catta feeding.In studying foraging movement of species it is common to consider the existence of intensive and extensive search modes.Intensive search mode is typically modeled by diffusion, while movement in the extensive search mode is modeled by advection (see for instance [17,24,37,38]).Here we assume that the degree of a patch, which is the number of corridors that connects it to the rest of patches, is related to the amount of resources that the species needs for feeding.In general, it is expected that a patch with large area is connected to many other patches, while a small one does not.It is also plausible that a large patch is more abundant in food than a small one.A similar strategy was used in modeling advection-diffusion movement of species by Grünbaum [23] who considered that the behavior of predators was mediated by its gut fullness, a.k.a.satiety, which was a function of prey availability.Similarly, Fagan et al. [16] constructed a continuum diffusion-advection model in which consumers movement in a landscape is driven by local advection governed by intermediate-scale resource detection functions.Therefore, we consider the study of the movement of L. catta between patches of this fragmented landscape as an appropriate platform for applying the advection-diffusion model on networks developed here.
In the current work, the edges of the landscape network represent the corridors existing between pairs of patches.These edges were determined in the following way.First, it was determined that the dispersal flux rate (i.e., rate of movements) between any two habitat patches depends on: (1) the geographical distance between the habitat patches and (2) the area of the habitat patches.It was also determined that the geographical distance has a much stronger impact than the area on the dispersal flux rate.Therefore, a measure S ij , which is proportional to the dispersal flux rate, was calculated using a negative exponential dispersal kernel [25] where a j is the area of the patch j (in m 2 ), D ij is the distance between patch i and patch j (in m), and α and b are constants ( in 1/m) determined by [3], such that the model partly captured the underlying assumptions behind the Incident Function Model [25].The distances between all pair of patches were measured using the Landsat image.The minimum hospitable patch area for Lemur catta was estimated to be 1 ha [19].The parameters were set such that only 5% of the maximum S ij remained between patches of 1 ha, if separated by a distance D ij equal to the chosen level of vagility.This level of S ij was also set as the cut-off, i.e., an infinitesimal increase in distance or an infinitesimal decrease in patch j's area would disconnect these patches.Thus, pairs of patches with a S ij below the cut-off were not considered as connected, i.e. no edge connected such pairs of nodes.We can then consider that the weighted adjacency matrix W of the landscape network is the Hadamard (entrywise) product of the matrix S where the values of S ij are given for every pair of patches in the landscape and the adjacency matrix A of a directed weighted graph: W := S A. The adjacency matrix A contains information only about the topology of the system, that is about the corridors which are used by the ring-tailed lemur (Lemur catta) to move across the environment.This species is very important for the ecosystem as a local seed disperser since it forages on fruits of many different plant species in the area.
In the previous section we have shown the relevance of the degree-degree assortativity on the evolution of the diffusion-advection processes taking place on a network.Because the landscape network analyzed here is a weighted directed graph, we consider the assortativity coefficient defined specifically for this class of networks according to Yuan et al. [41]: where X and Y are two quantitative features for all the vertices, W := i,j∈V w ij is the sum of edge weights, Xsou = i,k∈V w ik X i /W is the average of feature X based on source vertices, Ȳtar = k,j∈V w kj Y j /W is the average of feature Y based on target vertices, and σ X = i,k∈V w ik X i − Xsou 2 /W and σ Y = k,j∈V w kj Y j − Ȳtar 2 /W are respectively the corresponding standard deviations.The features X and Y can be the in or out degree of the vertices.
Here we are interested in how the in-flow of species at a given node is correlated with the in-flow at another.Therefore, from the several combinations of in-and out-assortativity coefficients we are interested in the in-inassortativity one.In this case the landscape network is in-in-disassortativity with coefficient ρ ≈ −0.118 (the rest of assortativity coefficients are all positive ones).This value shows that in the network there is a certain tendency to nodes with high in-degree to be connected to nodes with low in-degree.This may give rise to a high concentration of the flow in low-degree nodes in the case that the hubs-biased advection with α = −1 is dominant, as shown in Section 5.4.
The in-flow of L. catta into each patch of the landscape can be measured by f in := 1 T W = 1 T (S A), which is equivalent to the in-strength-the term strength is used instead of degree in the case of weighted networks-of the weighted directed landscape network.In Figure 6 we present a scatterplot of the values of the f in vs. those of k in as well as its best fit.It can be seen that there is an exponential relation between the two variables: f in ≈ 0.6603e 0.2224kin .However, it is evident that k in cannot be used directly to estimate the values of the in-flow.The in-degree is able to explain only 45% of the variance of f in using the previous exponential model, i.e., r 2 ≈ 0.45.Consequently, we can observe that for single values of the in-degree there are wide ranges of values of the in-flow.For instance, for k in = 9 there are values of f in ranging from 2 to 25. Thence, the values of the in-flow cannot be derived directly from information coming purely from the adjacency matrix of the graph, i.e., by the in-degree.
Our goal here is to capture the information contained in the dispersal flux rate matrix S using the information provided by the time evolution operator based on the advection-diffusion Laplacian of the system.Here we proceed as follows.First, because the landscape network is directed, we remove all the directionality of the edges and create a graph G in which the nodes are the same as in the landscape network G, and two nodes in G are connected if and only if there is a directed edge in any direction between the corresponding nodes.Using this graph we build the matrix L := γ adv L * α + γ dif L for different values between 0 and 1 for the parameter γ adv and taking γ dif = 1 − γ adv .Then, we calculate the time evolution operator for different values of t, exp −t L and estimate the values: fin (t) := 1 T exp −t L A . (5.7) We do not consider here the steady state of the system as the one for which we determine fin (t).Instead, we investigate different values of t which better fit the values of fin (t) in relation to those of f in .The reason for doing that is that the foraging process of a species in a landscape is a dynamic process far from equilibrium.Therefore, our approach is a reverse engineering one in which we try to recover the value of t as well as of γ adv and γ dif from the existing data.For doing that we compare the values of fin (t) with those of f in by means of regression analysis (see further).In practice, what we are doing is estimating the values of the time t and the parameters γ adv and γ dif that better reproduce the dispersal flux rate matrix of the landscape system.Our goal here is not to reproduce exactly the values of f in observed in the real-world landscape network.Such "exact" reproduction would need the combination of network-theoretic parameters with other parameters that depend on vegetation patterns in the patches, geometric characteristics of the corridors, etc., all of which escape to the goal of this work.Therefore our aim is to use this experiment as a proof of concept showing that the diffusion-advection equation on networks proposed here is able to capture some of the dynamical characteristics of this foraging process in a real-world landscape.
In order to determine what is the set of parameters t, γ adv and γ dif that reproduce in a better way the values of f in we proceed as follow.First we observe that the kind of relation between fin (t) and f in is the kind of power-law.That is, a simple plot of the two parameters results in a straight line in a log-log plot.Then, we obtained the least square fit of fin (t) vs. f in using a time step of 100 and step 0.1 in γ adv and γ dif .With these values we determined the value of t which produces the largest Pearson correlation coefficient and smaller root of mean square error (RMSE) of the regression.We then fix t at 700, which produces the best results.At this point we drop the step in γ adv and γ dif to 0.01 and plot the results of the correlation coefficient and RMSE as illustrated in Figure 7 for the value of α = 1.A similar plot was obtained for α = −1, but in this case the correlation coefficient was smaller and the RMSE higher than for the value of α = 1.Note that, from the analysis given in Section 5.3, the fact that α = 1 matches better the dispersal flux rather than α = −1 indicates that the advection process on the landscape is increasing the amount of in-flow that receive the hubs of the network, compared to the in-flow that they would receive if only diffusion takes place.The largest correlation coefficients for α = 1 are obtained for values of γ dif > 0.7 (γ adv < 0.3) and in particular the minimum error in the prediction of the values of f in is obtained for α = 1 with γ dif = 0.8 and γ adv = 0.2 (see Fig. 7).The best fit between both parameters is f in = 1.993 • 10 4 f 2.59 in where r = 0.769 and RM SE = 2.524.In Figure 8 we show  the landscape where the size of each node is proportional either to our calculated in-flow (Fig. 8(a)) or to the empirical in-degree (Fig. 8(b)).Therefore, we have shown here under the scope and limitations of the current approach that advective motion of species in a landscape, together with the diffusive one, is a plausible mechanism for their foraging strategy.This agrees with the empirical evidence existing in the literature, and confirms that the current theoretical framework is valuable for modeling such realistic scenarios.

Conclusions
We have defined a new kind of advection operator for graphs.Unlike previous advection operators on graphs, this one is defined for undirected graphs, and the preferred direction of the motion in the network is defined by the degree of each node.This means that it only depends on the topological structure of the given graph.We have then introduced an advection-diffusion equation for the case of undirected graphs by combining the new advection operator and the standard graph Laplacian.This degree-biased advection-diffusion equation is shown to be very flexible for modelling all kinds of advection-diffusion phenomena, given that its parameters can be widely modified.Additionally, we tested this model in a real-life landscape network representing the habitat of Lemur catta in southern Madagascar.Here we have provided evidence that the combined contributions of advection and diffusion represents better the foraging movement of L. catta in a large geographic region represented by a patched landscape than the consideration of diffusion alone.This points out to the importance and relevance of the newly defined advection operator for undirected networks as well as of the generalized advection-diffusion equation proposed here for modeling such realistic scenarios.The new degree-biased advection operator as well as the corresponding advection-diffusion equation can be easily generalized to describe non-local temporal and spatial effects which are important in several real-life problems.

Definition 3 . 4 .Definition 3 . 5 .
Let G ↔ be a double-orientation of graph G = (V, E) with |V | = n nodes and |E| = m edges.Then, the gradient operator or oriented incidence matrix of G is the operator D ∈ R 2m×n whose entries are defined byD ve :=    1 if vertex v is the head of oriented edge e = (v, j), −1if vertex v is the tail of edge e = (j, v), Let G ↔ be a double-orientation of G = (V, E) with |V | = n nodes and |E| = m edges.Then, the divergence operator of G is the operator D * ∈ R n×2m mapping from the edges to the nodes.

Figure 1 .
Figure 1.Representation of the structure of the Erdős-Rényi graph that we used with 8 nodes and 15 edges, where the size of each node is proportional to its degree.

Figure 2 .
Figure 2. Illustration of of the results of simulations in an Erdős-Rényi random graph G (8, 15) (dashed lines) and the analytically calculated steady states (solid lines).For the choice α = −1, panel (A) shows the evolution and steady state of the advection equation while panel (B) represents the same for the advection-diffusion equation.Panel (C) represents the advection equation when α = 1 and panel (D) shows the results of the advection-diffusion equation for the same value of α.

Figure 3 .
Figure 3. Change of the Pearson correlation coefficient between the vectors fin (γ adv = 0) and fin (γ adv > 0) for a star graph with 10 nodes.
κ and µ n = n are the eigenvalues of L. Therefore because ψ j1 = 0 for all 2 ≤ j ≤ n − 1 we have exp −t L vw = ψ 1v φ w1 + ψ nv φ wn e −tn .(5.4)Consequently, for n 1 the last term in the expression for exp −t L vw in the star graph is nearly zero.
n for (v, w) ∈ E. Therefore, fin,hub (γ adv = 0) = n − 1 n and fin,leaf (γ adv = 0) = 1 n .This means that most of in-flow is concentrated on the central hub due to the diffusion from the leaves.For relatively small values of γ adv the value of fin,hub (γ adv ) is n − 1 n minus the small term γ adv n − 2 n

Figure 4 .
Figure 4. Contour plot of the Pearson correlation coefficient between the vectors fin (γ adv = 0) and fin (γ adv > 0) as a function of the assortativity coefficient and the advection parameter γ adv for random graphs created with the Erdős-Rényi mechanism for which the edges have been rewired assortatively or disassortatively maintaining the same degree of the nodes.

Figure 5 .
Figure 5. Representation of the foraging landscape studied in the article [3].

Figure 6 .
Figure 6.Illustration of the relation between the in-flow f in at every node of the real-world landscape network and the purely topological parameter of the in-degree k in .The solid line shown corresponds to the best fit, which is of exponential shape (see text for explanation).Such fit explains only 45% of the variance in the values of f in as given by the squared Pearson correlation coefficient.

Figure 7 .
Figure 7. Plot of the Pearson correlation coefficient r (left scale) and the root of the mean square error (RMSE) of the regression (right scale) for the fit of the type f in = a f b in , where a and b are parameters.The left plot is for the value α = 1 and the right one is for α = −1.The values of f in are obtained from the experimental data as described in the main text of the paper and fin (t) were calculated for t = 700 (see main text) and changing γ dif with step 0.01 (notice that γ dif = 1 − γ adv ).

Figure 8 .
Figure 8. Illustration of the in-flow of L. catta estimated by the advection-diffusion model presented here (a) as well as obtained according to the geographic characteristics of the landscape according to Bodin et al. [3] (b).