GENERALIZED NETWORK TRANSPORT AND EULER-HILLE FORMULA

. In this article we consider asymptotic properties of network ﬂow models with fast transport along the edges and explore their connection with an operator version of the Euler formula for the exponential function. This connection, combined with the theory of the regular convergence of semigroups, allows for proving that for fast transport along the edges and slow rate of redistribution of the ﬂow at the nodes, the network ﬂow semigroup (or its suitable projection) can be approximated by a ﬁnite dimensional dynamical system related to the boundary conditions at the nodes of the network. The novelty of our results lies in considering more general boundary operators than that allowed for in previous papers.

1. Introduction. One of the major challenges in mathematical modelling of real life phenomena is to strike a balance between the amount of data that is taken into considerations and the simplicity of the model so that its reliable analysis is feasible. In many cases we model large networks of interrelated dynamical systems that can be described at the microscale, by taking into account the detailed dynamics at each node and each edge of the network, or at the macroscale, by only looking at the overall, aggregated, interactions between the nodes and the edges. While models built from first principles largely depend on the intuition and general understanding of the field by modelers and could differ from one another, it is important that various models describing the same phenomenon are consistent; that is, that they could be linked by a rigorous mathematical procedure that informs under what circumstances they give the same results. In particular, since models at the macroscale are less complex than those at the microscale, it is of practical importance to develop a systematic approach to deriving the former from the latter in such a way that the essential features of the microdynamics are retained. Such a procedure has been known in applied mathematics under different names, such as aggregation, [3,17,18], or state lumping [2,14]. The process of aggregation is schematically presented on Fig. 1. It is important to note that in most cases such an aggregation, or lumping, can only be achieved is an approximate sense. In such a case the diagram on Fig. 1 will only be approximately commutative; that is, macrosolutions obtained in different ways will be close but not necessarily equal.
One of such systematic approximate aggregation methods is based on the existence of different time scales and is briefly explained below.
Many models describe intertwined processes occurring at different time scales. Such models can be reduced by nondimensionalization to systems of equations containing a small parameter, say , that encodes the ratio between the scales. We adopt the convention that macro effects begin to dominate when is small. Thus, to describe the processes at the macroscale it would be desirable to be able to set = 0 and hence simplify the equations. However, in many cases simply setting = 0 dramatically changes the properties of the model and a sophisticated limit procedure must be employed to correctly determine the appropriate macrosystem.
On the other hand, often we are faced with a reverse problem. We have a macromodel that is too crude for some purpose and we want to build a more detailed micromodel describing the same phenomenon which is consistent with the original macromodel in the sense that for some ranges of the parameters of the micromodel, its solutions can be approximated by the solution to the macromodel.
Problems of this type usually belong to the class of singularly perturbed problems and can be approached via a number of techniques such as asymptotic expansions (described in e.g. [9]), or methods based on the Trotter-Kato theorem for semigroups of operators, see [13], that are particularly suitable for infinite dimensional linear models.
In this paper we shall focus on models that originate from transport processes on graphs, developed in [15,22,10]. Example 1. Transport on networks. Let us consider a substance flowing through a network of channels. Due to the unidirectional character of the flow, following e.g. [15,22], the network can be treated as a directed metric graph (digraph) G consisting, say, of n nodes {w 1 , w 2 , ..., w n } and m edges {e 1 , e 2 , ..., e m } . For simplicity we assume that G has no multiple edges, no loops (edges with head and tail at the same node) and no sinks (nodes with no outgoing edge), see product of the unit intervals. Then the distribution of the substance on the graph can be described by a vector function u = (u 1 , u 2 , ..., u n ), where each u i , defined on [0, 1], describes the density of the material flowing along the edge e i . The wave that propagates along each edge is described by a linear scalar transport equation and the whole system takes the form . . , m. The second equation gives the amount of the material flowing into e i from all edges e j incoming to the node that is the tail of e i . To make it more clear, let us recall that the line digraph Q of G is defined as the graph whose nodes are the edges of G and the edges of Q are identified with the pairs (e i , e j ) of the edges of G whenever the head e i coincides with the tail of e j . In other words, each node of G generates an edge (or several edges) in Q whenever it joins an incoming and outgoing edge of G, see must have a special form -it must be a weighted adjacency matrix of the line graph of G. In particular, if the redistribution of the material at each node is governed by the Kirchhoff law, see e.g. [22,25], then B is column stochastic and the flow is conservative. The solution to (1) may be too detailed or too costly to obtain for particular applications. If we only are interested in a macro description of the process, which in this case may amount to giving the quantity of the material on each edge at a given time t; that is, by providing the vector v(t) = (v 1 (t), . . . , v n (t)) of masses at the edges at time t, then a standard mass conservation argument yields a system of ordinary differential equation where B is a n × n matrix describing the redistribution of the masses in the nodes of the line graph of G. We observe that if we have the micro solution u(x, t) = (u 1 (x, t), . . . , u n (x, t)) to (1), then the mass on each edge is given bȳ Since certainly solving (2) is easier than first solving (1) and then integrating the solution, a natural question arises whether, for a given problem (1), it is possible to find a matrix B such that the solution to (2) equals (or at least approximately equals) the integrated solution (3) to (1). In other words, we are interested whether one can aggregate (1) to obtain such a macrosystem (2) that the diagram on Fig.  1 is (approximately) commutative.
It turns out that mathematically similar models appear in other fields, not necessarily related to graph theory.

Example 2. Kimmel-Stivers and Lebowitz-Rubinow-Rotenberg models.
Consider a gene amplification process; that is, an increase or decrease of numbers of drug resistant gene copies observed in tumor cells and associated with their resistance to certain drugs. In the models considered in [20,21], the population is divided into a denumerable quantity of types according to the number of drug resistant genes, (V n (t)) n∈N0 (where N 0 = N ∪ {0}) and its evolution is modelled by a branching random walk. We assume that at death any cell in class n (apart from class 0) produces two offspring that only can belong to the classes n + 1, n, or n − 1. Then the expectations of the random variables (V n (t)) n∈N0 , denoted by (v n (t)) n∈N0 , satisfy the birth-and-death system with proliferation, where d n+1 , c n−1 are the rates of recruitment from the populations n + 1 and n − 1 into the population n, and a n is the net growth rate of the population n that incorporates birth, death and loss to other populations of cells of type n. We assume that the sequences (a n ) n∈N0 , (d n ) n∈N , (c n ) n∈N are bounded. However, the cells have their own vital dynamics and hence it is more accurate to describe them by the age specific density (u n (x, t)) n∈N0 . Also, the mutations can change the number of the resistant genes in an arbitrary way and cells can die at any time (at the rate µ n , n ∈ N 0 ). Further, the mutations can occur: • due to the replication errors occurring when the cell divides (at the rates r np µ n 0 1 Figure 6. Discrete Lebowitz-Rubinow-Rotenberg model. • due to external factors (mutagenes) that may happen at any moment of the cell's life cycle (at the rates (r ij ) i,j∈N0 ).
Taking these more detailed aspects of the process into account, we arrive at the wellknown discrete Lebowitz-Rubinow-Rotenberg model of cell maturation, [23,26], u n (x, 0) =ů n (x), n ≥ 0, where the summation is extended over all states j that, by mutation, can change the number of genes to n.
As before, a natural question is to which extent models (4) and (5) are related or, more precisely, under what conditions the integrals of solutions to (5) can be approximated by the solutions of (4).
Observe that models (1) and (5) are structurally the same and can be written in a unified way as   The problem is posed in the space X = L 1 ([0, 1], X), where X is an l 1 type space, precisely defined later (though for most results it can be an arbitrary Banach space). Further, u = (u 1 , . . . , u n , . . .) is a, possibly infinite, vector whose components u i , i = 1, . . . , give the density of some quantity at the state i (e.g. the density of a fluid on the i-th edge of a graph in the transport model, or the age specific density of cells with i copies of a particular gene in the mutation model). The matrix B : X → X gives the rules for the redistribution of the flow at the nodes (e.g. the Kirchhoff law in the graph flow, or the redistribution of daughter cells among different classes due to mutations occurring during mitosis) and M : X → X describes possible changes that occur outside the nodes (e.g. some loss in transport inside the channels, or mutations caused by external mutagenes). We note that not every system (6) can be interpreted as a graph transport model, see [6]. As mentioned above our interest is to determine, when the process described by (6) can be approximately aggregated to a macromodel; that is, to a system of ordinary differential equations modelling the transfer of mass from one edge to another. In recent papers, [7,8] (see also [8,14], where transport was replaced by diffusion), the authors were able to show the convergence of solutions to (6) with fast transport (that is, with ∂ x u(x, t) replaced by −1 ∂ x u(x, t)) balanced by slow exchange between the states (modelled by B = I + C, where C is an arbitrary matrix and I is the identity) to solutions of an appropriate system of differential equations.
These results also revealed an interesting interplay between the regular convergence of semigroups and the Euler-Hille formula where C is a bounded operator on a Banach space X. Note that (7), with suitable changes, is valid even if −C is an unbounded generator of a C 0 -semigroup, [16, Section III.5], see also [12], but we do not pursue this direction here.
To explain the relation between the regular convergence of semigroups and (7) in more detail, we recall that the regular convergence as → 0 of the semigroups solving (6) with the speed of transport −1 and B = I + C, that is equivalent to the convergence of their resolvents and thus relatively easy to prove, occurs on the space X (that is thus the space of the regular convergence of the semigroups, see [13]). Since the semigroup solving (6) is given as the composition of the iterates of B and the translation semigroup (see e.g. (12)), it follows that the regular convergence of the semigroups solving (6) is equivalent to (7). We emphasize the role played by the regular convergence of the semigroup solving (6) in the proof of (7). As illustrated in Example 4, a direct proof would require a detailed knowledge of the fine structure of eigenvalues of the perturbed operator that seems to be difficult to obtain in a general case.
The analysis of (6) with B = I + C and of (7) has been carried out in [8]. The ideas of that paper were extended in [7,12] to prove a similar result with I replaced by a cyclic operator T but the price to pay was that only the projections of the semigroups onto the eigenspace of T belonging to the eigenvalue 1 converge.
The main aim of this paper is to generalize the above results, both on the asymptotic behaviour of solutions to (6) and on the Euler-Hille type formula (7), to operators of the form B = K + C, where K is a contraction having 1 as an isolated eigenvalue, as well as to allow for the perturbation M in the transport terms. In particular, this generalization enabled an asymptotic analysis of a boundary perturbation of any network flow (1) with the node exchange rule given by the Kirchhoff rule. Such an analysis was beyond the theory developed earlier due to the presence of the identity matrix I in the boundary conditions there, that did not allow a network flow representation of (6), [6]. At the same time, as explained in Example 6, we only are able to approximate the long term behaviour of the flow.

Preliminaries. Let us consider the Banach space
j∈I be linear and bounded operators on X, while C( ) = (c ij ( )) i,j∈I be a family of bounded operators continuously depending on ∈Ī := (0, 0 ] for some 0 > 0; we denote In what follows we use the matrix norms induced from X.
Considered pointwise, these operators induce bounded linear operators on X that shall be denoted by the same symbols. We introduce the projection P : X → X by As discussed in Introduction, we consider the family of problems where we can reduce (9) to the equivalent problem, see e.g.
where C M ( ) = (c m,ij ( )) i,j∈I is given by Since the matrices K, C( ) and M are fixed throughout the paper, we introduce operators A 0 and A defined, respectively, by the right hand sides of the differential equations in (11) and (9) (11) and (9). We note that the domains are well defined as their elements are absolutely continuous functions. Then it is known, [5,Theorem 3.1], that (11), as well as (9) that is a bounded perturbation of (11), are well-posed in the sense of semigroups; that is, for each ∈ I there exists a unique semigroup solution v (x, t) = e tA0 v (x) of (11) (respectively, u (x, t) = e tA ů(x) of (9)). A crucial role in the considerations will be played by the representation formula, [15,Proposition 3.3], that takes the form and hence the properties of (e tA0 ) t≥0 are closely related to the behaviour of the powers of K + C M ( ). Then, the results for the solutions to (9) will be recovered from that for v (x, t) using where then, by e.g. [16, Proposition IV.1.13], there exists a spectral decomposition X = and Π i are the spectral projections given by Now we are ready to introduce the main assumption on K.
K is contractive and 1 is an isolated and semisimple eigenvalue of K.
Remark 1. By saying that 1 is a semisimple eigenvalue we mean that We note that if the multiplicity of 1 is finite, then X 1 is finite dimensional and then contractiveness of K ensures that 1 is semisimple, see e.g. [24,Section 7.10].
In further considerations we will use the 'diagonal-block' form of the operator K. Since 1 is an isolated element of σ(K), and, accordingly, we have the spectral projections Π i , i = 1, 2, and the decomposition X = X 1 ⊕ X 2 , as in (15), specified for K. Hence we can write where K 22 is a bounded operator with σ(K 22 ) = σ 2 (K). The plan of the paper is as follows. In Section 3 we will consider the convergence of the resolvents of A and apply the results of [13] on the regular convergence of semigroups to show that the solutions e tA ů converge as → 0 providedů belongs to the regularity space identified as Π 1 PX. In particular, in Corollary 3 we shall prove a version of the Euler-Hille formula (7) for the operator K + C M ( ). In Section 4 we show that the projections Π 1 Pe tA ů converge for any initial value u ∈ X, extending thus the results of [7,12].
It is worthwhile to emphasize that the results for (e tA0 ) t≥0 are obtained by a combination of classical spectral results of [19] with the results on the regular convergence from [13] without which the convergence on the subspace X 1 would be difficult to prove. This can be seen in the calculations in Example 4 that require a detailed knowledge of the form of the perturbed eigenvalue, which is difficult to achieve in a general case. On the other hand, the transfer of the results to (e tA ) t≥0 was relatively straightforward thanks to (13) -in particular, the lack of the commutativity of Π 1 and M would make it difficult to apply the Phillips-Dyson expansion as in [13,Chapter 29] to prove the second part of Theorem 4.1.
3. Regular convergence. This part is based on the Trotter-Kato theorem, see [13,Chapter 8], which gives a relation between the convergence of a sequence of semigroups and the convergence of the sequence of resolvents of their generators. More precisely, we consider a family of equibounded semigroups such that the resolvents of their generators converge to some operator R λ on X. Then it follows that R λ is the resolvent of the generator of a semigroup on X 0 = rangR λ ⊂ X, that is also the limit of the original sequence of semigroups if restricted to X 0 , see [13,Theorem 8.2]. Let us present an application of this result to (9). (17) is satisfied, then the family of semigroups (e tA ) t≥0 , ∈ I , is equibounded.
Proof. As noted above, it is sufficient to show the equiboundedness for (A 0 , D(A 0 )). The proof is similar to that of [7, Lemma 2.1] so that we only indicate major different points. The general solution of the resolvent equation is given by where f ∈ X, E λ (s) = e − λs I, 0 ≤ s ≤ 1 and r = (r ,i ) i∈I is an arbitrary vector. Using the boundary condition v(0) = (K + C M ( ))v(1), we obtain Let us fix > 0. Since we have where r is uniquely defined by the Neumann series and hence the resolvent of A 0 exists for λ > ω := C M , see the proof of [7, Lemma 2.1]. Without losing generality, we can assume K + C M ( ) ≥ 0 (otherwise we can define an analogous dominating problem with |K| + |C M ( )|, where | · | denotes the absolute value of the matrix, see [5, Theorem 2.1] for details). Then we can restrict our considerations to v ≥ 0. Adding together the rows in (21) and, using the contractiveness of K (note that |K| is contractive if and only if K is), we obtain j∈I r ,j ≤ j∈I (1 + c m,j ( ))e − λ r ,j + j∈I (1 + c m,j ( )) 1 0 e λ(s−1) f j (s)ds, (24) where c m,j ( ) = i∈I c m,ij ( ) for any j ∈ I. Integrating (20), we obtain Hence the estimate of R(λ, A 0 )f has the same form as if the operator K in (11) was replaced by I (compare [7, Eq. (11)]). Thus the equiboundedness of (e tA0 ) t≥0 follows as in [7, Lemma 2.1] and by, say, (13) and the equiboundedness of (e − xM ) ∈I , also the equiboundedness of (e tA ) t≥0 follows.

Remark 2.
We observe that the equiboundedness allows us to work with initial conditions independent of . In particular, it allows us to disregard the -dependence inv , introduced by transformation (10) in (11).
Although the equiboundedness is satisfied for any contractive K, it is not sufficient to obtain the convergence of the family of resolvents {R(λ, A )} >0 to the resolvent of some operator on the whole space. We have, however, in the uniform operator topology, where P was defined in (8).
The proof is based on the following lemma.
Proof. Let X 0 and X 2 be the spectral eigenspaces of S associated, respectively, with the eigenvalue 0 and with σ(S) \ {0}, X = X 0 ⊕ X 2 . Using the notation in (16), we rewrite the equation 1 where w i = Π i w, i = 0, 2 for any w ∈ X. Since B 00 is invertible, Then, by the invertibility of S 22 , the operator S 22 + B 20 − B 20 B −1 00 B 02 is also invertible for sufficiently small > 0 and Using again [1, Proposition 7.2] and the boundedness of the projection, we conclude that lim Finally which ends the proof.
Note that though both R 0 and R 1 depend on , they are equibounded operators as Recall that, by (20)- (21), the resolvent of A 0 takes the form Now, using (30), we obtain the first operator in the definition of r in the form appearing in Lemma 3.3 where I − K has 0 as an isolated semi-simple eigenvalue and Ω( ) contains higher order terms coming from the expansions of C M ( ) and E λ (1) so that it satisfies Ω( ) → 0 as → 0. Splitting, as before, X into the eigenspace X 1 related to the eigenvalue 1 of K and its complement X 2 , we have and we see that Π 1 (λK − C − KM)Π 1 is invertible for λ > ω. Hence, by Lemma 3.3 and the continuity of the inverse, [1, Proposition 7.2], we obtain lim →0 + hence Finally, by we conclude, using (33), (34), (36) and the equiboundedness of resolvent operator, that We can now state the result on the convergence of solutions to the problem (9), when tends to zero. for anyů ∈ Π 1 PX and the convergence is almost uniform on [0, ∞).
Proof. By (12), for u ∈ X we have for t − t < x ≤ 1.
we have In the following example we show that the assumption that 1 is a semi-simple eigenvalue (or, here, equivalently, that K is a contraction) is essential for the existence of a nontrivial limit semigroup.

Example 3. Let us consider (9) with
arbitrary C = (c ij ) i,j=1,2 and M = 0. Here 1 is an eigenvalue of K with algebraic multiplicity 2 and geometric multiplicity equal 1 and, clearly, K is not a contraction.
As in the proof of Theorem 3.2, first we consider r := I − e − λ (K + C) Hence we see that in general R(λ, A 0 ) does not converge as → 0. However, even if the limit exists, its only invariant subspace is {0} and hence the limit cannot be the resolvent of any nontrivial operator.
Corollary 1 gives us the convergence of semigroups (e tA ) t≥0 on a very restricted subspace Π 1 PX of functions that are not only x independent but also belong to the kernel of I − K. The example below suggests a way to overcome the latter restriction. The projection Π 1 onto the eigenspace of K associated with the eigenvalue 1 is given by We note that Π 1 and C do not commute, Π 1 C = CΠ 1 and Following (12), we examine the convergence (K + C) t when goes to 0. The eigenvalues of K + C are λ ±, = 3 + 2 ± √ 20 2 + 4 + 1 4 and λ +, → 1, λ −, → 1/2 as → 0. To find the spectral decomposition of K + C, first we find the corresponding eigenvectors. Noting that K + C is symmetric, the left and right eigenvectors are equal and are given by where by ·, · we denote the scalar product in R 2 . Hence To simplify calculations, we note that Since t − t is bounded and λ +, converges to 1, we have Hence, for any u ∈ R 2 , The example shows an important role played by the fact that the complement of {1} in σ(K), σ 2 (K), is contained in a disc of radius smaller than 1. Indeed, we have Corollary 3. Let (17) be satisfied and let sup{|λ| : λ ∈ σ 2 (K)} < 1 in (19). Then for any u ∈ X = PX and the convergence is almost uniform on (0, ∞).
Proof. Again, we begin with (e tA0 ) t≥0 . Let u ∈ PX and Π i , X i , i = 1, 2, be defined as in (15). According to [19,Theorem IV.3.16], the spectrum of K + C M ( ) can be decomposed in a similar way; that is, for sufficiently small > 0 the contours γ i enclose the disjoint components of σ(K + C M ( )) and there exist subspaces X 1 and X 2 such that X = X 1 ⊕ X 2 and dimX i = dimX i i = 1, 2.
Now the convergence of (e tA ) t≥0 follows as in the proof of Corollary 1.
As seen in the following example, one cannot expect convergence if there is an eigenvalue on the unit circle that is different from 1.
with some for 0 < α < π/2. Then 1 and e ±iα are the eigenvalues of K and the projection onto the eigenspace X 1 is given by Furthermore, the eigenvalues coincide with the one associated with operator K + C.
As in Example 4, we find the spectral representation of (K + C) We see that Π 1 R 3 is spanned by (1, 0, 0) and clearly (K + C) t u converges for u ∈ Π 1 R 3 , as in Corollary 2. On the other hand, (K + C) t u does not converge for any other u due to the properties of e ±iα t . However, The last result of the example is in agreement with [7,12]. Actually, it is even stronger since the spectrum of the operator K above may be not cyclic. In the next section we shall show how to combine [7,Theorem 4.1] with the earlier results of this paper to prove that (51) holds in a more general case. 4. The convergence of the projection. for anyů ∈ X, almost uniformly on [0, ∞).
Proof. As before, we begin with (e tA0 ) t≥0 . For t ≥ 0 and ∈ I , we denote n = t + 1. Then, using (12), for anyů ∈ X we have We see that the last term converges to zero due to (42) and thus we can restrict ourselves to the first two terms. Since Π 1 K = Π 1 , and the convergence of Π 1 Pe tA0 ů as → 0 reduces to that of Π 1 (K+ C M ( )) t Pů; that is, on X. Since we can write provided Π 1 w = 0. Since 1 is an isolated eigenvalue of K, we can consider the projection Π 1 as in (47). This projection commutes with K + C M ( ) and thus The last term equals 0 by (54) and the other two converge to 0 on account of (47) and the equiboundedness of (K + C M ( )) t .
Remark 3. Theorem 4.1 can be seen as a generalization of the concept of lumpability of linear evolution equations, [2], mentioned in the introduction. Roughly speaking, a semigroup (e tA ) t≥0 on a Banach space X is said to be lumpable by a bounded surjection M : X → Y if there exists a semigroup (e tÂ ) t≥0 on Y such that The operator M is viewed as a reduction of the state space and the process itself, under various names, has been discussed in many contexts, see e.g [3], where a comprehensive literature of the subject is provided. For instance, in the ecological context, see [17,18], lumpability in the sense of [2] is called a perfect aggregation. In op.cit. it is observed that the perfect aggregation very rarely occurs in applications and thus there is an interest in (56) occurring in some approximate sense. Thus, we can say that a family of semigroups (e tA ) t≥0 , ∈ I , on X is asymptotically lumpable by an operator M if there exists a semigroup (e tÂ ) t≥0 on rangM such that for all u ∈ X, t ≥ 0, lim →0 + M e tA u = e tÂ M u.
Example 6. Let us consider a transport problem on the graph described in Example 1, If we write (u 1, (x, t), . . . , u 4, (x, t)) = e tA ů for an arbitraryů ∈ X, then Formulae (57) and (58) should be looked at in the context of [10, Theorem 5.2] and [22,Theorem 4.5]. Namely, if we consider the flow described by (57) with = 1 and B = 0, then it will converge to the stationary distribution given by as t → ∞. Note that K is column stochastic and hence there is conservation of mass, so that the whole initial mass is eventually distributed according to the Perron eigenvector of K.
We observe that this result falls short of the expectation expressed in Example 1, as the dynamics of the macromodel is given by the single equation instead of the postulated system (2) of the same dimension as the adjacency matrix of the line graph of G. Moreover, the solution of the above equation approximates the total mass of the system instead of the masses on each edge of G, as suggested in (2). However, one should recognize that our macromodel consists of two elements: the Perron eigenvector e r giving the long term profile of the flow on the edges (the fraction of the total mass on each edge stabilizing after long time) and the time evolution of the total mass gathered on the edges corresponding to the nonzero entries of e r . Since in our example the long term behaviour of the flow is onedimensional, it is natural that the macroscopic model is one-dimensional as well. In other words, the aggregation method presented in this paper yields a macromodel approximating the long term dynamics of the given micromodel. Note also that the macro models retains some information of the micromodel -the structure of the original graph G partly determines the Perron eigenvector of K, the weighted adjacency matrix of the line graph of G.
We also note that the dimension of the eigenspace belonging to the eigenvalue 1 of K equals the number of terminal strong components of G, [11, p. 17]. If there are, say, k such strong components of G, then the long term flow consists of k independent flows, each in its own component, [10]. In such a case the limiting system of ordinary differential equations (2) consists of k differential equations describing the evolution of the material trapped in each terminal component.