Infrastructure adaptation and emergence of loops in network routing with time-dependent loads

Network routing approaches are widely used to study the evolution in time of self-adapting systems. However, few advances have been made for problems where adaptation is governed by time-dependent inputs. In this work we study a dynamical system where the edge conductivities of a network are regulated by time-varying mass loads injected on nodes. Motivated by empirical observations, we assume that conductivities adapt slowly with respect to the characteristic time of the loads. Furthermore, assuming the loads to be periodic, we derive a dynamics where the evolution of the system is controlled by a matrix obtained with the Fourier coefficients of the input loads. Remarkably, we find a sufficient condition on these coefficients that determines when the resulting network topologies are trees. We show an example of this on the Bordeaux bus network where we tune the input loads to interpolate between loopy and tree topologies. We validate our model on several synthetic networks and provide an expression for long-time solutions of the original conductivities.

Several efficient methods have been proposed to solve this problem.A popular approach is that of messagepassing algorithms [31], where sources of mass are matched in sender-receiver pairs, and messages encode mass transfer between them [22,[32][33][34][35][36].Promising results have also been obtained with optimal transport theory [2,15,20,27,[37][38][39][40][41][42], the approach we consider in this work.The general idea behind this method is to describe the transport of mass as a process being regulated by edge capacities, quantities evolving with a dynamical system to allocate mass fluxes.
Despite their usage in modeling transportation problems across domains, a common drawback of all these methods is to consider only stationary loads, i.e. resources that are injected and travel through the network do not change with time.This assumption may not be valid in certain scenarios.For instance, blood vessels are known for adapting their structure continuously to meet changing metabolic demands [25,[43][44][45].Similarly, passengers in transportation systems enter stations with hourly, weekly, and seasonal time-varying rates [46].
A viable approach to model these systems is to control the network evolution considering an ensemble average of the stress generated by the loads [24,25,37].This relies on assuming stationary loads on nodes but with their positions varied stochastically.The ensemble average over the loads' locations is then computed as a proxy of a system with loads of fixed locations but time-varying amounts.This technique has also been employed to study network resilience to edge cutting [3] or for routing problems with spatially correlated loads [2].
Remarkably, adding stochasticity in the loads may lead to the emergence of loops in the resulting optimal networks topologies [2,3,24,25,37].This result is complementary to the hierarchical formation of trees since loops provide alternative routes to accommodate fluctuations, or guarantee robustness against broken links.Recently, loops formation has also been observed in multicommodity setups [27,28], where the loads are deterministic inputs of the problem.In this case loop generation is a consequence of having different commodities interacting in a unique shared infrastructure.
In all these works, the time-varying character of the transport network loads is neglected because the main problem variables are taken on average.Here we develop a model that considers the explicit time dependence of the mass inflows and investigate, both analytically and numerically, the long-term behavior of time-varying transport networks.This allows us to show that it is not the process uncertainty, inherent in any stochastic framework, but the non stationarity of the loads that promotes loops, which is fundamentally different from what can be concluded using stochastic formulations.
In particular, we generalize the routing problem in the work of Facca et al. [41] by considering periodic mass loads on nodes.We postulate an analytical relationship connecting physical quantities as the edge conductivities and the coefficients of the Fourier series expansion of the loads.We then define a dynamics that rapidly converges to the long-run average solutions of the original dynamics.
Our model relies on the idea of distinguishing slowvarying variables from fast ones.The first are capacities of edges that are regulated by the Fourier coefficients of the forcing; the second are, for example, loads of passengers entering and exiting network nodes.The physical intuition is that, while fluxes of passengers in a transportation system have the same rate of change of the network loads, roads do not.In fact, it is reasonable to assume that a network manager has a coarser observation scale of a transportation system than the users, whose paths rapidly fluctuate.In practice, this means that modifications in the network infrastructure occur on a much larger time scale than that of daily passengers' fluctuations.
Remarkably, we find that the Fourier decomposition of the loads yields a sufficient condition to determine whether the resulting optimal networks will contain loops or be a tree.Performing a numerical validation of our dynamics on synthetic networks we are also able to provide an analytical expression for the long-run conductivities.Precisely, we find that the conductivities start oscillating around constant values at large time scales, and at certain frequencies that can be expressed in terms of those of the input loads.Furthermore, we define a Lyapunov functional for our dynamical formulation, allowing us to interpret stationary topologies as optimal networks, i.e., structures minimizing the global cost to build the graph.Finally, we examine a case study with loads that are the sum of decoupled harmonic oscillators, finding that the condition on the Fourier coefficients can be equivalently reformulated in terms of the loads' amplitudes and phases.We numerically investigate this last setup on the Bordeaux bus network.

II. TIME-VARYING LOADS IN ROUTING OPTIMIZATION ON NETWORKS
Consider a network G with nodes v ∈ V and edges e ∈ E, each of length e > 0. The orientation of the edges is conventionally assigned by the signed incidence matrix of the graph, with entries B ve = ±1 if node v is the tail or the head of edge e, and B ve = 0 otherwise.We consider a routing optimization problem on G setting time-varying mass loads S(t) = {S v (t)} on nodes being the amount of mass either injected in (S v (t) > 0) or extracted from (S v (t) < 0) node v. Concretely, one could think of S(t) as a time-dependent origin-destination vector of passengers moving in a transportation network, where mass entries correspond to the fraction of passenger flowing though 1. Schematic visualization of the problem.We highlight conductivities (brown), the length of an edge (green), and the difference of pressure (purple) along an edge that triggers their fluxes.In the rightmost blue and yellow panels we depict two scenarios where the time-dependent loads S(t) generate fluxes that move from the central part of the network to its periphery (t = t1) and vice versa (t = t2 = t1).Node and edge widths are proportional to S(t) and F (t), respectively.
stations.This allows us to write Kirchhoff's conservation law as where µ = {µ e } are the non-negative edge conductivities, p(t) = {p v (t)} are pressure potentials on nodes and L vu (µ) := e B ve (µ e / e )B ue are the entries of the weighted Laplacian of the network [47].The conductivities can be interpreted as the capacities that the edges must have to allocate the mass loads acting on the nodes; thus we can consider them proportional to the edges' sizes.When considering passengers moving along a transportation network, µ can be seen as the width of a road or more generally a measure of the infrastructure's resources used to carry traffic flows.We propose a model in which the forcings S(t) dictate the time evolution of the conductivities by means of a feedback dynamics.In particular, we couple Eq. (1) with the system of with the system of ordinary differential equations (ODEs) − µ e (t) ∀e ∈ E (2) with m e > 0 initial values.For a solution trajectory µ(t), we define the fluxes F e (t) ≡ F e (µ(t), S(t)) := µ e (t)[p u (t)− p v (t)]/ e for e = (u, v), with p v (t) ≡ p v (µ(t), S(t)) := u L † vu (µ(t))S u (t) the solution of Eq. ( 1), where L † denotes the Laplacian pseudoinverse.We assume that the system is isolated, namely, v S v (t) = 0, ∀t ≥ 0, so p(t) is a well-defined potential (see [47], Lemma 0).Specifically, Klein and Randić [47] showed that L generally is not invertible, but Eq. ( 1) can be solved by the pseudoinverse within the subspace orthogonal to the unitary vector, that is, when v S v (t) = 0.
In Eq. (2), the growth in time of the conductivities is proportional to the flux forcing term F 2 e (t) with µ decaying exponentially when no flux flows though an edge.In practice, this corresponds to enlarging a road when many passengers travel along it, and reducing it when there is no traffic.We illustrate this intuition with a schematic representation in Fig. 1.
Our dynamical formulation assumes continuous variables for fluxes and conductivities, but the mass S(t) could be arbitrarily continuous or discrete.While this is valid in many scenarios (e.g., when modeling a large number of individuals), it may be limiting in cases where a discrete (or atomic) representation is necessary to capture fine-grain differences in the number of passengers.For this, one should consider alternative formulations and approaches, for instance, using message passing or belief propagation as in [22,[32][33][34][35][36].
Finally, we remark that Eq. ( 2) can be made scale independent with respect to the model variables by an opportune nondimensionalization that we describe in detail in Appendix A.

A. Slow adaptation of conductivities
In several biological systems the adaptation time of organisms is much slower (weeks) than the characteristic time of the mass injected in the system (seconds) [25,[43][44][45].In order to describe these organisms, a common approach is that of approximating the fast time-varying input loads with combinations of open and closed switchlike nodes with constant inflows, and to assume that the conductivities are regulated by an ensemble average of the pressures over different states of the loads [3,24,25,37].
Instead, here we want to model the evolution of these slow adapting conductivities taking in account the time dependence of the loads.We formalize this hypothesis by assuming: (i) the existence of a slow time scale τ , with τ = Kt and K 1, and that (ii) in a fixed time window ∆, small with respect to the slow variable τ but large with respect to the time t, holds for some conductivities μ = {μ e } with the natural time of evolution being τ .Time scales are depicted in Fig. 2. We can interpret t as seconds, ∆ as days, and τ as months.Such distinction between different natural time scales is observed in the interplay of rivers and tide loads in coastal delta formation [9], where the assumption is that tides cycle much faster than the river channel adaptation, a distinction analogous to that between t and τ .
Finally, we assume that (iii) the evolution in τ of μ is determined by the time integral average of the product of the mass loads.Assumptions (i), (ii) and (iii) together lead to the definition of for all e ∈ E and τ ≥ 0, where we introduced The functional Φ is the natural approximation of the right-hand side of Eq. ( 2), as shown in Appendix B.
We define then a trajectory μ(τ ) as a solution of the dynamics with me > 0 initial conditions.In general, Φ is difficult to manipulate as the loads S(t) may assume any arbitrary expression, possibly preventing the exact computation of the time integrals.For this reason, we investigate its behavior for a particular class of functions S(t) that allows for analytical tractability.

B. Periodicity of the loads
We consider periodic loads S(t), with period T small with respect to the fixed integration window ∆ introduced in Section III A: This allows us to express each S v (t) using its Fourier series S v (t) = nv∈Z c nv v exp(iωn v t), with ω = 2π/T .Substituting this into Eq.( 5) yields the pivotal result: holding for all u, v ∈ V and τ ≥ 0. The matrix C has entries

C. Periodic-loads dynamics
Combining Eq. ( 5) and Eq. ( 9), we can build a dynamics for some new conductivities μ = {μ e }, which evolve in the slow time scale τ .Precisely, we ignore negligible contributions in Eq. ( 9) and define with me > 0 initial conditions.The right hand side of Eq. ( 10) is such that Φ Φ for ∆ 1 and reads An important point is that the problem in Eqs. ( 10)-( 11) is not equivalent to the dynamics Eqs. ( 1)-( 3) with each S v (t) integrated over T .The latter case would imply that C vu had the form C uv = Su Sv , where Sv is the integral of S v (t) over the period.This is only a particular case of the dynamics in Eqs. ( 10)- (11).
Noticeably, in this case the condition rank(C) = 1 holds, i.e., since C is symmetric, there exists a vector y ∈ R |V | such that C uv = y u y v , ∀u, v ∈ V .This is a sufficient condition for Eqs. ( 10)- (11) to return a loopless network at convergence (see Appendix D for a proof), and confirms previous results observed for constant loads [23,26,29].
However, this condition does not hold generally for any arbitrary choice of the loads, as C may have a more general expression, in particular rank(C) > 1.Moreover, the case of constant loads is not the only one where rank(C) = 1.We provide an example of this in Section V B, where we explore the case of each S v (t) being the sum of a finite number of harmonic oscillators.

IV. CHARACTERIZATION OF THE FAST DYNAMICS
Finding an analytical expression for the fast conductivities µ(t) solutions of Eqs. ( 1)-( 3) cannot be done by directly solving the dynamics, because of the non linear dependence on µ(t) in the Laplacian pseudoinverse.Nevertheless, here we propose an argument to characterize their long-time behavior.
We support our findings with an empirical validation on synthetic networks built taking the Delaunay triangulation of |V | nodes placed at random in the unit square.In our experiments, we set |V | = 2 i , with i = 3, . . ., 9. The vector of loads S(t) is S(t) := 20 S 1 (t) + 10 S 4 (t) + 5 S 8 (t), where each factor is defined as S n (t) := q n cos(ωnt), with amplitudes extracted at random from a |V |-dimensional Dirichlet distribution as The period has been conventionally set to have ω = 2π.
We observe that the evolution of the fast conductivities is typically divided in two phases, as shown in Fig. 3a.First, the conductivities undergo a stabilization transient for t < t STAB , where they strongly depend on their initial conditions m e and significantly change their mean values.Then, when t > t STAB , the conductivities reach a plateau and oscillate around fixed values.More precisely, either they move around mean values that are far from zero and preserve their oscillatory nature for all times, or they decay to zero with negligible oscillations that are progressively damped as t increases.These experimental observations suggest the following ansatz for the stabilized solutions, for all t > t STAB and e ∈ E: ( We experimentally notice that also the fluxes start to oscillate around a constant value after a first stabilization time interval [see Fig. 3a].We use this evidence to deduce (see Appendix E) that the main oscillatory modes of the conductivities are resonant with the squared fluxes, and have the form  Remarkably, these numerical experiments also serve as a validation for hypothesis (ii) in Section III A. In fact, for any sufficiently slow time τ , the conductivities fluctuate around a constant value, suggesting the possibility of neglecting their fast oscillatory nature when studying asymptotics of Eqs. ( 1

A. Candidate Lyapunov functional
We empirically observe [see Fig. 3b] that our new dynamics Eqs. ( 10 where for each edge e we define the squared slow fluxes F 2 e := (μ 2 e / 2 e ) uv A eu (μ)A ev (μ)C uv .Noticeably, if rank(C) = 1 holds, it is possible to formally prove that Lγ (μ) is a well-defined Lyapunov functional for Eqs. ( 10)- (11) (see Appendix D for detailed derivations).In addition, we can interpret the functional as in Lonardi et al. [27] for multicommodity optimal transport.Namely, the Lyapunov is the sum of a dissipation cost, the first addend in Eq. ( 16), with an infrastructural cost, the price needed to build the transport network.
Additionally, we observe that the candidate Lyapunov functional Lγ converges to a value that is the same achieved by the running average functional over the period T : with µ that is evaluated along solution trajectories of Eqs.
(1)-(3).The functional Eq. ( 17) reaches a plateau at the stabilization time t STAB , when the fast conductivities µ(t) start oscillating around constant values.Remarkably, in Fig. 3b we see that t STAB t STOP , which is due to the fact that the time step δt for the numerical discretization of Eqs. ( 1)-(3) has to be set much lower than δt in order to capture the oscillatory nature of the loads.In our experiments we set it to δt = δt/10.A practical consequence of this is that the discretization of Eqs. ( 10)-( 11) is a fast and scalable alternative to extract the conductivities around which long-run solutions of Eqs. ( 1)-( 3) stabilize.
Because of this analogy between an optimal transport (functional minimization) setup and the solutions of our dynamical system, we can interpret the networks determined from the dynamics in Eqs. ( 10)- (11) as optimal topologies minimizing the infrastructural and dissipation cost.These networks can also be obtained by averaging long-run solutions of the original dynamics in Eqs. ( 1)-(3).In fact, as discussed in Section IV, long-run trajectories of Eqs. ( 1)-( 3) oscillate around asymptotics of the newly defined dynamical system in Eqs. ( 10)- (11).

A. Conditions for the generation of loops in closed form
If C has rank(C) = 1, i.e., C vu = y u y v for some y ∈ R |V | , the dynamics Eqs. ( 10)-( 11) produces trees at convergence.One trivial case where this holds is when the loads S(t) are static, i.e., constant for all times.However, this is not the only setting where rank(C) = 1 is satisfied.In particular, there are cases where such a condition holds but S do change in time.
Here, we explore a case of study proposing an ansatz where the loads are the sum of decoupled harmonic oscillators: with ω = 2π/T , n i v , N v ∈ N, and A i v , d v ∈ R. By construction, these loads are periodic in T , hence we compare them with their Fourier series representation . Equating this expression with Eq. ( 18) yields where we conventionally set φ 0 v = 0, ∀v ∈ V and where only a finite number of Fourier coefficients are different from zero, given that the sum in Eq. ( 18) is finite.
The goal here is to express rank(C) = 1 in terms of {A i v }, {n i v }, and {φ i v }, amplitudes, modes and phases of the harmonic oscillators, respectively.To do that, we start by noting that rank(C) = 1 is satisfied if and only if C uv = y u y v , ∀u, v ∈ V , with y v = ± √ C vv , and where the plus or minus signs have to be determined among 2 |V | possible choices in such a way that v y v = 0 (see Appendix F).
Defining the complex vectors ν v = {c nv v } with entries of the Fourier coefficients in Eq. ( 19), we rewrite where the centered dot denotes the complex dot product and || • || is its correspondent norm.Thus, the rank condition on C can be reformulated in terms of an equivalent linear dependence condition of the form ν v = λν u between the vectors ν v , v ∈ V , and for λ = 0. Finally, substituting Eq. ( 19) in this linear dependence condition leads to the following main result.Proposition 1.Let the time-dependent loads S(t) injected in the network nodes be as in Eq. (18).If the following hold, then, any γ ≤ 1, a stationary solution of Eqs. ( 10)-( 11) is a tree: e. sources and sinks are in phase,

A
For a formal justification of this result see Appendix F.

B. Numerical tests on the Bordeaux bus
In order to test the rank condition on C we design two experiments on the real network of the buses of Bordeaux.The network topology has been constructed focusing on a central region of the city, and using data collected from [48].Here we assume that the loads, representing passengers entering or exiting the network, vary much faster than the conductivities.These latter quantities can be thought of as the size of the roads that a network manager needs to design; thus we can safely assume their evolution to happen on a larger time scale with respect to that of S(t).
All the sinks u = v 1 , v 2 have loads S u (t) = −[S v1 (t) + S v2 (t)]/5 in both cases, to ensure conservation of mass.We expect that in the first case the network extracted from Eqs. ( 10)- (11) with γ ≤ 1 is a tree.In the second case the network can possibly contain loops.We run the dynamics setting γ = 0.9 and we display our findings in Fig. 4a.The empirical results reflect our predictions: the blue network (the first case) is a tree.In contrast, in the orange network (the second case) loops emerge.
We further validate our results on the bus network of Bordeaux with a second experiment.We assign the loads S v (t) = Exploiting the exact relation that the matrix C has with the modes of the loads (see Section V A), it is possible to see that our particular construction of S(t) gives ranks ranging in 1 ≤ rank(C) ≤ 6.
We show our results in Fig. 4b, where we plot the fraction of basis loops of the network against the rank of C. The dynamics is executed for γ = 0.5 and the random extraction of the forcings has been varied over 100 runs.In the plot, it is clearly visible that for all values of Q, the fraction of basis loops is zero at rank(C) = 1.Moreover, we can see that when we increase the complexity of the problem, i.e., when rank(C) grows, the values attained on the y axis also increase.This suggests that the rank of the C can be used as a qualitative proxy to predict the number of loops in the optimal transport network.Finally, as one could intuitively expect, the basis loop fraction increases with Q, i.e., with the number of nodes where mass is injected or extracted.

VI. CONCLUSIONS AND OUTLOOKS
Routing models on networks are relevant to study many real-world problems.While most of the works in the current literature consider stationary setups, i.e., the inflows injected in the network do not change in time [2,15,20,27,[37][38][39][40][41][42], few recent works investigate timevarying loads and the majority of these models study solely the averaged evolution of the networks' variables [3,9,24,37].
In this work we analyzed a dynamical system where the conductivities are regulated by time-varying mass inflows.Motivated by empirical observations [25,[43][44][45], we assumed the existence of auxiliary conductivities that have response times which are much slower than those of the loads.Furthermore, in order to make the problem analytically tractable, we supposed that all the loads injected in nodes are periodic, in a period that is substantially smaller than the adaptation time of the new conductivities.These two hypothesis together allowed us to deduce a dynamics where the evolution of the systems is solely regulated by an input matrix constructed using the Fourier series expansion of the loads.
The resulting dynamics allowed us to derive the main findings of our work.In detail, combining theoretical arguments with empirical evidence on synthetic networks, we found an expression for the long-run solutions of the original dynamics, which cannot otherwise be obtained by simply solving the original dynamics.These long-run solutions are the sum of stationary components, equal to the asymptotics of the dynamics we constructed, and an oscillatory one.This second contribution can be expressed as the sum of periodic signals, with modes related to those of the loads.Moreover, we discussed a sufficient condition on the loads that determines when optimal transport networks can be loopless.Such a condition was numerically validated on the Bordeaux bus network.Finally, our dynamics can be connected to an optimization setup, as shown by the proposed candidate Lyapunov functional.As a result, asymptotic trajectories of our dynamics minimize the total cost needed to build the network infrastructure.
Importantly, the numerical discretization of the dynamics we proposed in this work can be used as an efficient method to rapidly converge to average long-run solutions of the original dynamics.
Our results can be extended in several ways.For instance, it would be interesting to investigate different types of input loads that relax the periodicity hypothesis and use this to analyze the behavior of the conductivities in different problems' settings.Similarly, it would be interesting to explore how this formalism adapts to multilayer networks, where passengers can enter different stations corresponding to different transportation modes [20].Another relevant application could be that of integrating our findings with the recent work of Baptista et al. [19], where the authors studied how topological properties of the transport network change in time, as we approach stationary configurations, and how these reflect on the shape of the conductivities.
While our work constitutes a step towards extending the formalism of capacitated networks to time-dependent loads, it is important to remark that our findings are valid in a particular time limit.Specifically, this is the scenario where conductivities slowly evolve with the integral average of periodic forcings, as introduced in Section III.It is not clear how the theoretical analysis presented in this work could be adapted to scenarios where loads and conductivities evolve with the same time scale.This could be an interesting avenue for future work.Another interesting direction could be that of considering additional constraints on the evolution of the conductivities, which are not currently included in our model.For instance, one could introduce a threshold capacity above which the edge traffic saturates, causing blockage of roads.
Altogether, we believe that our results enrich the cur-rent knowledge on network routing problems with timevarying input loads and have immediate practical implications.In fact our model is deterministic, since there is only one single realization of the inputs, and thus adequate to model real-world scenarios where time-dependent loads are measured quantities, e.g., the amount of passengers traveling in a metro (which can be easily tracked), without the need of stochastic formulations that require the introduction of probability distributions that are hard to characterize.
To facilitate practitioners in using our model, we have made the algorithmic implementation publicly available [49].We enforce the hypothesis of periodicity of the loads, i.e. S v (t) = S v (t + T ), with T /∆ 1, and we parametrize the integration window ∆ as ∆ = KT , K 1.This allows us to split the integral in Eq. ( 9) into two separate contributions.In detail, making the reasonable hypothesis that S u (t)S v (t) is bounded by M < +∞ for all t ≥ 0 and for all u, v ∈ V , we can write Hence, we separate the first K integrals over the period T from the last one in ( K T, KT ).Since K 1, the first K contributions can be evaluated as with δ ij being the Kronecker delta for two indices i and j.
As for the second term, in the limit K 1 we can write showing that integrals over the small interval ( K T, KT ) are negligible for a large integration window.

Appendix D: Sufficient condition on the rank for optimal trees
We discuss in detail the sufficient condition rank(C) = 1 to obtain loopless optimal networks running the dynamics Eqs. ( 10)- (11).Our argument proceed as follows.
The matrix C is symmetric by construction; thus if its rank is 1 its eigenvalue decomposition is of the form C = N i=1 λ i x i x i , with all the eigenvalues equal to zero except one.We conventionally choose it to be λ 1 = v C vv > 0, with a unit norm eigenvector x 1 .Defining y := √ λ 1 x 1 and substituting the eigendecomposition of C in Eq. ( 12), we get that Φe is proportional to Fe := (μ e / e ) v B ev pv , with pv := u L † vu (μ)y u .In order to conclude, we need to show that p is a well-defined solution of Kirchhoff's law, u L uv (μ)p u = y v , i.e., y is a zero-sum vector [47].This comes as a consequence of conservation of mass.Indeed, since for all times v S v (t) = 0 holds, we have v S u (t)S v (t) = 0, ∀u ∈ V .Using Eq. ( 9) and ignoring negligible terms of O(∆), this yields v C uv = 0, ∀u ∈ V .Finally, substituting the eigendecomposition of C in this last relation gives v y u y v = 0, ∀u ∈ V .This is satisfied only if v y v = 0, i.e., y is a zero-sum vector.In this case, Eqs. ( 10)- (11) correspond to the standard dynamics Eqs. ( 1)-(3) with constant loads, which are S(t) = y, ∀t ≥ 0, and we recover the well-known result that optimal networks are trees for γ ≤ 1 [23,26,29].

Appendix E: Derivation of Eq. (15)
In order to discern the nature of the fast oscillating component b e (t) of the stabilized solutions, we need to investigate further the original dynamics Eqs. ( 1 Finally, substituting Eq. (E1) in Eq. (E9), we get the desired results, i.e., the main oscillatory modes of the conductivities, hence of b e (t), are resonant with the squared fluxes.Thus we obtain Eq. ( 15).

Validation on synthetic networks
We test these expressions numerically on networks generated as described in Section IV.We compute P e := From Eq. ( 15) we expect to have most of the spectral density of b e (t) concentrated on the modes in K, i.e., the ratio P e /P N ,e should be close to 1 for each edge.In Fig. 5a we plot P = {P e } versus P N = {P N ,e } for the example network considered in Fig. 3.The plot supports Eq. ( 15); indeed, the element-wise ratio P/P N is close to 1 for all points (each correspondent to a different edge) with a slight deviation only for small (thus negligible) values of the conductivities.
We further validate this result on an additional synthetic example network.We construct the Delaunay networks described in Section IV considering 100 combinations of seeds for the nodes' positions and for the random input loads.Then, we compute the spectral densities P and plot them against δP , with entries δP e := (P N ,e − P e )/P e .We show in Fig. 5b results for γ = 0.5, 1, 1.5 on 100 random graphs of size |V | = 8.Here we clearly see that δP are negligible for any edge with

<< 1 FIG. 2 .
FIG. 2. Schematic representation of the different time variables.The two arrows denote time scales t and τ .The time windows ∆, large with respect to t, are denoted with blue curved brackets and the fast period T in orange.Each window ∆ along which we integrate the dynamics Eqs.(1)-(3) contains a large number of periods T .
with c * denoting the complex conjugation of c.The term O(∆) contains all negligible contributions ε, decaying as ε/∆ → 0 for ∆ → +∞.For a detailed derivation of this result one can refer to Appendix C.
b e (t) = n,m∈N b m e b n e exp(iω(n + m)t) ∀e ∈ E , (15)with N := {n v } set of the Fourier modes of the loads.Hence, the conductivities oscillate with modes determined by those of the loads.This result is supported by several numerical experiments (see Appendix E for details).

FIG. 4 .
FIG. 4. Bordeaux bus optimal transport network.(a) Network visualization.Input loads have been built as described in Section V B. The tree network originated by C with rank 1 is plotted in blue, the loopy topology in orange.The yellow and the magenta stars denote the geographical location of the two loads, the green squares those of the sinks.Here the width of edges corresponds to slow conductivities at convergence μ∞ e .Results are plotted for γ = 0.9 (b) Basis loop fraction against rank(C).Points correspond to averages over 100 runs of the experiments where positions of the sources and sinks are extracted at random.Shaded regions denote their standard deviations.Results are displayed for γ = 0.5.

R
|F[b e ](f )| 2 df , the total spectral density of the oscillatory components b e (t), after the conductivities µ(t) stabilize.Here F[•](f ) is the Fourier transform operator.Additionally, we calculate P N , obtained summing the atomic contributions of the spectral density on the modes k ∈ K := {k s.t.k = n + m, for n, m ∈ N }.Namely, P N ,e := k∈K R |F[b e ](f )| 2 δ(f − k) df, ∀e ∈ E.

FIG. 5 .
FIG. 5. Spectral density validation.(a) Plot of PN versus P for the example network of Fig. 3.Each point corresponds to an edge; the color scale is that of Fig. 3a.(b) Plot of P versus δP .Each point corresponds to an edge, marker color denotes γ = 0.5, 1, 1.5.The dashed line is the cut-off used to build δPα.(c) Compatibility of δPα with δPα = 0 for different networks' sizes.Markers and bars correspond to averages and standard deviations over 100 random configurations of the problem, respectively.The networks have been obtained pairing ten seeds for node coordinate generation and ten seeds for mass and conductivity initialization µe ∼ U (0, 1).
FIG. 3. Characterization of the fast conductivities µ(t).All results are computed on a synthetic network with |V | = 8 and setting γ = 1.0.(a) Fast conductivities µ(t) and fluxes F (t) are drawn with solid lines, stationary solutions μ∞ are dashed.Labels on the x axis correspond to the number of iterations of the numerical discretization of Eqs.(1)-(3).Conductivities are depicted in two time windows, before and after their stabilization time t STAB .Fluxes are drawn only for t t STAB .Colors denote different edges.(b) Evolution of Lγ and of Lγ T in time.The green and the red circles denote t STAB and t STOP , respectively.

Derivation of Eq. (9)
In detail, in Eq. (B2) we used the definition of the fluxes F e (t) := µ e (t)[p u (t) − p v (t)]/ e , ∀e ∈ E and evaluated the pressure solving Kirch-hoff's law, i.e. p v (t) := u L † vu (µ(t))S u (t), ∀v ∈ V .The second important step is in Eq. (B4), where we used hypothesis (ii) in Section III A, namely the approximation in Eq. (4), to carry the conductivities out of the time integral, and we introduced A ev (µ(t)) := u B eu L † vu (µ(t)), ∀e ∈ E, ∀v ∈ V .
(3)3).From our numerical validation we observe that the fluxes start to oscillate around a constant value after a first stabilization time interval [see Fig.3a], analogously to the conductivities.This suggests the ansatz F e (t) = STAB and with the terms F ne e amplitudes of the Fourier series decomposition.We argue that pairing this expression with Kirchhoff's law, i.e., e B ve F e (t) = S v (t), yields Note that these definitions lead to F e (t) being a potentialbased flux and yield p v (t) = α v (t) + β v (t), ∀v ∈ V .Substituting Eq. (E7) and Eq.(E8) in Eq. (E5) and Eq.(E6), respectively, implies that ψ e (t) = 0, ∀e ∈ E, and for sufficiently large times.Hence, the only non zero terms in the Fourier decomposition of F e (t) have modes in N .This result is particularly useful to describe the behavior of µ(t) at large times.First, we recall that µ e (t) = μ∞ e + b e (t), ∀e ∈ E, as discussed in Section IV.Moreover, we observe that in our numerical experiments [see Fig.3a] the size of the amplitude of the oscillatory term b e (t) is negligible in size with respect to μ∞ e , unless µ e (t) decays to zeros.This allows us to approximate Eq. (2) as ne∈Z F ne e exp(iωn e t), ∀e ∈ E, for all times t sufficiently larger than t valid for all v ∈ V .Now, in order to guarantee that the fluxes {ϕ e (t), ψ e (t)} are well-defined, we suppose the existence of two time-dependent potentials α(t) = {α v (t)} and β(t) = {β v (t)}.These are defined on the network nodes and such that for all e ∈ E we have