Distributed algorithms from arboreal ants for the shortest path problem

Significance Evolution has led to natural algorithms that regulate collective behavior in many biological systems. Here, we investigate natural algorithms that solve the shortest path problem, a basic optimization problem on graphs. A recent field study showed that arboreal ants solve a variant of the shortest path problem while creating and maintaining trail networks. They do this without any central control and with minimal computational resources. We propose a model for how trail networks change over time that explains this phenomenon and leads to surprising algorithms for the shortest path problem and its variants.


Introduction
Biological systems, such as ant trail networks, are fascinating examples of distributed algorithms in nature [Lyn96; NB14; FK17; Cou+05] , often finding globally optimum solutions using simple local interactions among individuals, devoid of central control. The study of natural algorithms has led to synergistic exchange between biology and computer science [NB11; Cha12a]. The algorithmic lens has enhanced our understanding of biological phenomena such as how birds flock [Cha09], how slime molds solve the shortest path problem [NYT00; BMV12; SV21] and how computation takes place in the brain [Pap+20; Val00]. Moreover, the process of evolution itself has been studied using an algorithmic lens [Cha+14; Val09; Liv+08; Vis14]. Also, inspiration from nature has led to new algorithmic ideas such as ant-inspired algorithms for distributed density estimation [MSL16], artificial neural networks in machine learning [LBH15], and algorithms for similarity search inspired by the fruit fly brain [DSN17], among others [Das+18; SDN20; Afe+11; CC12].
Here we investigate how the trail networks of the arboreal turtle ant (Cephalotes goniodontus) can solve variants of the shortest path problem, a basic optimization problem on graphs [For56;Bel58;Dij+59]. Textbook algorithms for this problem find optimum solutions using knowledge of the entire network [KT06; DPV08;Cor+09]. Turtle ants nest and forage in the tree canopy of the tropical forest; their trail network is constrained to lie on a natural graph formed by tangled branches and vines, and no ant has any global information about the network. Observations of turtle ants in the field show that a colony's trail network approximately minimizes the number of vertices [Cha+21]. We develop a model that gives a biologically plausible explanation for this observation, and outlines other intriguing phenomena as described in the next section.

Summary of model and results
Turtle ant colonies form trails on a graph whose vertices correspond to junctions in the vegetation, and edges correspond to branches connecting these junctions. A colony's network of trails connects many nests and food sources. The trail network minimizes the number of vertices [Cha+21] (compared to simulated random networks), approximately solving a variant of the Steiner tree problem [Win87; Cha+99; Lat+11], which is a generalization of the shortest path problem for multiple terminal vertices. Here we focus on a section of this network, considering two terminal vertices, such as a nest and a food source, and we seek to explain how a colony can find the path with minimum number of vertices between the two terminals. We model trails as a bidirectional flow of ants between the two terminal vertices; bidirectional flow is characteristic of the trails of this species [Gor17]. The flow in our model is similar to models of flow in traffic networks [War52; BMW55; RT02]. Ants lay trail pheromone on edges as they traverse them, and the next edge to traverse is chosen based on the level of pheromone. The pheromone decays with time. Some fraction of flow leaks as it passes through each vertex, modeling the loss of ants due to exploration. Chandrasekhar et al. [Cha+21] hypothesized loss of ants at the vertices to be the reason why ants prefer paths with fewer vertices. Our model leads to 4 main results: 1. We first consider the linear decision rule, which at each vertex, divides the flow among the next set of edges in proportion to their pheromone level. We show through simulations and analysis that when the incoming rate of flow remains unchanged, the dynamics converges to the path with the minimum leakage. This is also the path with the minimum number of vertices when all vertices have equal leakage. This result describes a biologically plausible process that explains how colonies can find paths with the minimum number of vertices.
2. We show that when the rate of flow increases with time, in the absence of leakage, the dynamics converges to the shortest path. Flow rate on ant trails can change over time [Den+86; Gor12; Bou+15; Rei+15; BCB17], for example, in turtle ants the flow rate can increase in response to new food sources [Gor12]. In other ant species, it has been shown that ant trails converge to the shortest path in certain simple graphs [Gos+89]. Our result shows a surprising link between these two phenomena: ant colonies can use their ability to increase the flow rate to find the shortest path.
3. We establish the utility of bidirectional flow by showing that it is necessary for convergence to the shortest or the minimum leakage path. In contrast, most flow-based problems considered in computer science and operations research have unidirectional flow [HR55; AMO88; Sch02; RT02].
4. We investigate the effect of increasing flow and leakage with decision rules other than the linear rule. For a general family of decision rules, we show that the linear rule is its unique member with guaranteed convergence to the shortest and the minimum leakage path. However, for various non-linear rules, we show that the dynamics still often converges to a path with smaller length and less leakage, compared to the path found in the absence of increasing flow and leakage respectively. Thus the utility of increasing flow and leakage is not limited to the linear decision rule.
Our model builds on a previous model by Chandrasekhar et al. [CGN18], adding components such as leakage and variation in flow rate. These components were not present in the model of Chandrasekhar et al. [CGN18] and are crucial for the phenomena we discuss above. Chandrasekhar et al. [CGN18] investigated how ants find alternative paths, not necessarily the ones with the minimum number of vertices, to route around ruptured links in a network. Here we ask how ants can find the path with the minimum number of vertices, and how the flow rate impacts the path found. We demonstrate that leakage at vertices can lead to convergence to the path with the minimum number of vertices, and increase in flow rate over time can lead to convergence to the shortest path.
Our work is different from traditional ant-colony optimization [DCM91; DB05; LSD18], in which the algorithms considered are not required to be biologically plausible. Models of ant colony optimization (ACO), inspired by ant behavior, solve combinatorial optimization problems, such as the traveling salesman problem [Yan+08] and the shortest path problem [ST12; NW06] . In ACO, individual agents, simulating ants, construct candidate solutions using heuristics, and then use limited communication, simulating trail pheromone, to lead other agents towards better solutions. The simulated ants have significantly more computational power than is biologically plausible for real ants. Unlike real ants, the simulated ants have the ability to remember, retrace and reinforce entire paths, and can use the quality of the global solution to determine the amount of "pheromone" to be laid.
Our model resembles the reinforced random walks introduced by Diaconis and others [DF80; Dav90; Pem+07], which have found applications in biology[SO97; CPB08;Smo+10]. Here, a single agent traverses a graph by choosing edges with probability proportional to their edge weight, with edge weight here analogous to the level of trail pheromone. This edge weight increases additively each time the edge is traversed. However, there are a few key differences between the model we study and the setup for reinforced random walks: 1) Our model involves many agents, modeled by a flow, and their behavior is affected by their collective action in putting down pheromone, while the model for reinforced random walks considers a single agent whose behavior is influenced by its past random choices. 2) Pheromone decays geometrically over time in our model, but edge weights do not decrease in reinforced random walks. 3) Our model has leakage at each vertex which is not present in the setup for reinforced random walks. Nevertheless, in a similar spirit to random walk-like processes studied before, the model we investigate is a Markov process of particular relevance to biology.

The Model
We consider bidirectional flow on a directed graph G = (V, E) with source vertex and destination vertex s and d respectively. For each vertex v, let f→ v (t) and f← v (t) denote the forward and backward flow on v present at time t. The forward flow moves along the direction of the edges, and the backward flow moves in the opposite direction. Let p uv (t) denote the pheromone level on edge (u, v) at time t. The pheromone level and flow together constitute the state of the system at any time t, which is updated as follows: 1. Flow movement: At each time step, the forward flow on a vertex moves along its outgoing edges, dividing itself based on a decision rule that depends on the pheromone levels on these edges. The simplest such rule divides the flow proportional to the pheromone level on the outgoing edges.
Formally, let f → uv (t) denote the forward flow moving along edge (u, v) at time t. Then according to this rule, which we call the linear decision rule, f → uv (t) = f→ u (t) p uv (t) ∑ z:(u,z)∈E p uz (t) . (1) The total forward flow on vertex v is the sum of flow along its incoming edges, multiplied by a leakage factor: The leakage parameter l v ∈ [0, 1] models the loss of ants due to exploration at each vertex.
Movement for backward flow takes place in exactly the same manner, with the direction of the flow reversed. Formally, let f ← uv (t) denote the backward flow moving along the edge (u, v) at time t. Then The total backward flow at u is the sum of backward flow along its outgoing edges, multiplied by the leakage parameter.
At each time t, new forward and backward flow f→ s (t) and f← d (t) appears on the source vertex and destination vertex s and d respectively.
2. Pheromone update: At each time step, the pheromone level on an edge increases by the amount of flow on it, and decays by a multiplicative factor of δ (similar to Chandrasekhar et al. [CGN18]): Note that unlike flow, the pheromone level present on an edge does not have any direction, and is influenced by flow from both the directions.

Definition 1.
For any path P from source s to destination d, we define its leakage l P as the fraction of flow that leaks out while moving through P, that is,

Discussion of Modeling Assumptions.
Our goal is to provide a simple and minimal biologically plausible model that explains how ants can find the path with the minimum number of vertices. We discuss our modeling assumptions below.
While we describe our model for unweighted graphs, it is general enough to capture the case with integral/rational edge lengths. For instance, an edge with length k can be represented using k unit length edges connected in series, with leakage at the vertices connecting these edges set to zero.
Our model uses a directed graph. Previous work shows that individual ants are unlikely to turn around on the trail, so that when an ant leaves a terminal, such as a nest or a food source, its distance from that terminal increases over time [Gor12]. To account for this, Chandrasekhar et al. [Cha+21] assign direction to each edge relative to a terminal vertex, where the outbound direction goes away from the terminal vertex and the inbound direction goes towards it. Similarly, we use directed edges in our model.
In our model, the backward flow at the destination vertex is not dependent on the forward flow reaching it at the previous time step. This is because an ant does not necessarily make a round trip from one terminal to the other and back. The flow between the two terminal vertices in our model represents only a section of the larger network, which includes many nests and food sources. Ants reaching a terminal vertex can go on to other parts of the network instead of turning back.

Results
In this section, we discuss the convergence properties for the model defined above. We defer all the proofs and simulation details to the supplementary material.

Convergence to the Minimum Leakage Path and the Shortest Path
Constant flow with time. Through simulations and analysis, we show that when the incoming forward and backward flow does not change with time, with the linear decision rule, the dynamics converges to the path with the minimum leakage (see Figure 3).
We run simulations for three families of directed graphs: 1. G(n, p): The G(n, p) model is a widely used random graph model which consists of graphs on n vertices where each pair of vertices has an edge with probability p.
2. G(n, p) with a locality constraint: The standard G(n, p) model allows edges between any two vertices (with probability p). However, the graphs formed by branches and vines in the natural vegetation have a local physical structure in which the edges are more likely between nearby vertices. To capture this, we consider the G(n, p) model with an additional locality constraint that an edge exists between vertex i and j only if |i − j| ≤ k for some parameter k. Here, the vertices are labelled from 1 to n with the source and the destination vertex labelled 1 and n respectively.
We generate a large number of instances with different values for n, p, k and other parameters, and observe convergence to the minimum leakage path in all the simulated instances. More details about the simulations can be found in supplementary material Appendix C. We complement our simulations on general graph models with a provable convergence result for graphs with two parallel paths, a case that has been experimentally investigated in the past [Gos+89; DS01].
Theorem 1. Consider a graph G consisting of two parallel paths P 1 and P 2 from s to d. Let the flow and the pheromone levels be updated according to the model in Section 2, and let P 1 be the path with the minimum leakage. If (i) the incoming flow values f→ s (t) and f← d (t) are non-zero and unchanging with time, and (ii) the initial pheromone level p uv (0) is positive for all edges (u, v) ∈ P 1 , then the flow dynamics governed by the linear decision rule converges to a state where all the flow goes through P 1 .
At a high level, the proof involves showing that relatively more pheromone accumulates on the path with less leakage as time progresses. Although the update rules governing our model are simple, mathematically understanding its dynamics is surprisingly non-trivial. Even for the seemingly simple case of two parallel paths, the progression of pheromone levels can be highly non-monotone, and the proof needs a careful construction of an appropriate potential function. We give a sketch of the proof in Section 4.
To connect this result to the observation of Chandrasekhar et al. [Cha+21] that ants form trails with approximately the minimum number of vertices, we need to connect leakage to the number of vertices. Note that the path with minimum leakage is also the path with the minimum number of vertices when all the vertices have equal leakage. Moreover, we can show that this connection between leakage and number of vertices degrades gracefully, and as long as the variation in leakage between different vertices is not too large, the path with the minimum leakage has approximately the minimum number of vertices. One way to formalize this is to assume that for any pair of vertices u and v, log(1 − l u ) and log(1 − l v ) are within a (1 + ϵ) factor of each other. We run the simulations for the same families of graphs considered for the last result. The differences in these set of simulations are the absence of leakage, and the incoming flow increases by a fixed factor in each step. We generate a large number of graph instances with different values of the parameters, and observe convergence to the shortest path in all the simulated instances.
We also show provable convergence to the shortest path for graphs consisting of two parallel paths.
Theorem 2. Consider a graph G consisting of two parallel paths P 1 and P 2 from s to d. Let the flow and the pheromone levels be updated according to the model in Section 2, and let P 1 be the shorter path. If (i) the initial pheromone level p uv (0) is positive for all edges (u, v) ∈ P 1 , (ii) leakage l P 1 = l P 2 = 0, and (iii) the incoming flow increases as follows: 1. Multiplicative increase: f→ s (t) = α t f→ s (0) and f→ then the flow dynamics governed by the linear decision rule converges to a state where all the flow goes through P 1 .
For ease of analysis, we consider only the cases when the flow increases by a fixed additive or multiplicative factor. Our analysis suggests that the outcome may be similar when the rate of increase is not fixed; further work is needed to demonstrate this.
For an intuitive explanation of this result, consider a graph with two parallel paths as shown in Figure 2 such that path P 1 is shorter than P 2 , and there is no leakage. Consider the forward flow on edges (d 1 , d) and (d 2 , d). Since P 1 is shorter, the forward flow on (d 1 , d) corresponds to the more recent forward flow that entered from s compared to the forward flow on (d 2 , d). Since the flow is increasing with time, more recent flow is larger ensuring that relatively more pheromone accumulates on (d 1 , d) than (d 2 , d) as time progresses. Similarly, relatively more pheromone accumulates on (s, s 1 ) than (s, s 2 ) as time progresses due to the increasing backward flow. Thus, as time progresses, relatively more pheromone accumulates on P 1 . However, as in the case with leakage, increase in relative pheromone levels on P 1 is not monotone and we require a more careful proof. We give a sketch of the proof in Section 4.
Previous studies have shown that ants are capable of finding the shortest path in certain simple graphs [Gos+89], and that factors such as detection of new food sources can positively reinforce the rate of flow of ants [Gor12; Den+86]. Our result shows a surprising connection between these two phenomena.
The main goal of our work is to demonstrate the intriguing connection between leakage and flow rate and the shortest path problem. We do this using simulations on general graph models, and analysis on graphs with parallel paths. Based on our simulations, we conjecture that the above results showing provable convergence hold for general graphs.

Conjecture 1.
The provable convergence results in Theorem 1 and 2 hold for general graphs.
We can view leakage and increasing flow as two possibly conflicting forces, leading to paths with minimum leakage and length respectively. An interesting direction for future work would be to investigate the dynamics when both these forces are active simultaneously.

General Rules and Fundamental Limits
In the previous subsection, we show that bidirectional flow with linear decision rule converges to the path with the minimum leakage when the flow is fixed, and to the shortest path when the flow is increasing and there is no leakage. How crucial is the bidirectional nature of the flow? How does the dynamics behave when the decision rule is non-linear? Next, we study these questions.
Necessity of bi-directional flow. We show that bi-directional flow is necessary to find the shortest or the minimum leakage path. This result holds independent of the decision rule, leakage levels and the change in flow rate; we formally state this result next. Let G be the set of all decision rules that distribute the forward (backward) flow at any vertex, onto the outgoing (incoming) edges only based on their pheromone levels.
Theorem 3. Consider any graph G with two parallel paths between s and d. Let there be unidirectional flow from s to d. For any decision rule in G, any setting of leakage parameters and with arbitrary incoming flow levels, there exists a setting of initial pheromone levels, such that the dynamics does not converge to the shortest or the minimum leakage path.
Thus, bi-directional flow is necessary for all pheromone based rules to guarantee convergence to the shortest or the minimum leakage path. The main idea behind this result is that with only unidirectional flow from s to d, pheromone on edges incident on s can not encode information about rest of the graph.
The conflict between leakage, change in flow and non-linearity. To understand the dynamics for other decision rules beyond the linear rule, we consider a family of decision rules F satisfying certain assumptions. This family is a subset of the family G discussed above. These rules distribute the forward flow among the outgoing edges of any vertex, and backward flow among the the incoming edges based on the normalized pheromone levels at these edges (as defined below). That is, allocation of flow does not depend on the absolute pheromone levels. We also assume that these rules are monotone in the sense that increasing the normalized pheromone level at any edge does not decrease the proportion of flow entering that edge. Most decision rules considered in the past [CGN18] satisfy these conditions. Here, the normalized pheromone level is defined as follows: Note that while the pheromone level on an edge has no associated direction, there is a forward and backward normalized pheromone level for each edge depending on whether the normalization is done with respect to the incoming edges or the outgoing edges.
More formally, we consider the family of rules given by functions g : [0, 1/2] → [0, 1] ∈ F . For any vertex with out-degree (in-degree) 2, a rule in this family uses a function g, which takes as input the minimum of the normalized pheromone levels on the two outgoing (incoming) edges, and returns the fraction of the forward (backward) flow on this edge. In other words, g(x) is the fraction of flow sent on edge with normalized pheromone level x, for x < 0.5. We assume that g is monotonically increasing, that is, g(x) ≤ g(x ′ ) for all x ≤ x ′ . Further, g satisfies g(0) = 0 and g(1/2) = 1/2. This condition says that if one of the edges has 0 normalized pheromone, there is no flow on it, and if both the edges have equal pheromone level, there is equal flow on them. For any vertex with out-degree (in-degree) 1, the forward (backward) flow goes to the next (previous) vertex, just as in the case of the linear rule. For our results below, we need to define these rules for only degree 1 and 2 vertices.
Note that the linear decision rule studied in previous subsection belongs to F , and corresponds to g(x) = x.
We study the dynamics when the decision rule belongs to the family of decision rules defined above and is non-linear. We show that for every non-linear rule g ∈ F , there exists a setting of leakage parameters in which the dynamics fails to converge to the path with the minimum leakage.
Theorem 4. Consider any graph G consisting of two parallel paths from s to d. When the incoming flow is fixed, for every non-linear decision rule g ∈ F , there exists a setting of leakage parameters and initial pheromone and flow levels dependent on g, such that the dynamics does not converge to the path with the minimum leakage.
We show an analogous result for the increasing flow case with no leakage.
Theorem 5. Consider any graph G consisting of two parallel paths from s to d with a unique shortest path. When the leakage is zero for all the vertices, for every non-linear decision rule g ∈ F , there exists a setting of initial pheromone and flow levels, with incoming flow increasing by a fixed multiplicative factor at each time step, such that the dynamics does not converge to the shortest path. The multiplicative factor and initial pheromone and flow levels are chosen as a function of g.
For an intuitive explanation for these results, consider the quadratic decision rule that distributes the flow in proportion to the square of the pheromone levels. Due to the square in the quadratic decision rule, given two edges incident on a vertex, this rule sends more than linearly proportional flow on the edge with the higher pheromone. Thus, if a path-that is not necessarily the shortest or the minimum leakage path-has relatively high pheromone initially, even more pheromone accumulates on it as time progresses, and the dynamics may not converge to the shortest or the minimum leakage path. Theorem 4 and 5 formalize this intuition for any non-linear decision rule belonging to the family F .
From the last subsection, we can view increasing flow and leakage as two conflicting forces, preferring the shortest and the minimum leakage path respectively. The above results suggest that non-linearity in the decision rule can be viewed as another force, in conflict with the forces of leakage and increasing flow, preferring certain states that may not correspond to the shortest or the minimum leakage path.  Usefulness of increasing flow and leakage not limited to the linear decision rule. Note that the above results only suggest that linear decision rule is necessary for guaranteed convergence to the shortest or the minimum leakage path. Can it be the case that even with some non-linear decision rules, the forces of increasing flow and leakage still help in finding shorter or smaller leakage paths respectively, compared to the paths found in the absence of these forces? To understand this question, we ran simulations for various non-linear decision rules, some of which have been previously used to model ant behaviour [CGN18;Den+90]. We observe that within each graph family, for a large fraction of the graph instances, the path found in the presence of these forces has length (respectively leakage) smaller than or equal to the length (respectively leakage) of the path found in their absence. Figure 4 demonstrates this for the quadratic decision rule for G(n, p) graphs with the locality constraint.
This rule divides the flow in proportion to the square of the pheromone levels.
In Figure 4a, we show the path length obtained by the quadratic decision rule with and without increase in flow, and by the linear rule with increase in flow, for 30 random graph instances. As discussed before, the linear rule finds the shortest path. But even with the quadratic rule, the path length obtained in the presence of increasing flow is smaller than or equal to the length obtained in its absence.
Similarly, Figure 4b shows the path leakage obtained by the quadratic decision rule with and without leakage present at the vertices, and by the linear decision rule with leakage present. There is a subtle distinction here between path leakage as an objective function and leakage as a process affecting the dynamics. For each graph instance, we assign leakage values to vertices (see supplementary material Appendix C for details). This gives us a path leakage objective function which we measure in all the three cases. However, in the case of the quadratic rule without leakage, the leakage process is not applied at vertices during the dynamics. This gives a baseline to which we compare the path leakage objective when the leakage process is applied. We observe that the quadratic rule with leakage applied leads to path leakage objective smaller than or equal to the baseline, while the linear rule with leakage applied minimizes the objective as discussed in Section 3.1.
We observe similar results for other graph families and non-linear decision rules. However, the extent to which the forces of increasing flow and leakage are effective varies with the non-linear decision rule and graph family. For instance, we observe that compared to the quadratic rule, these forces are more effective for a non-linear rule closer to the linear rule, dividing the flow in proportion to the 1.1 th power of the pheromone levels. Also, for most graph families and non-linear decision rules considered, there is a small fraction of instances where these forces end up increasing the path length (respectively path leakage). Nonetheless, for all the graph families and non-linear decision rules considered, for most (> 80%) graph instances, the path found in the presence of these forces has length (respectively leakage) smaller than or equal to the length (respectively leakage) of the path found in their absence. Thus the usefulness of the forces of leakage and increasing flow is not limited to the linear decision rule. We include more details and discussion of these simulations in supplementary material Appendix C.

Linear Decision Rule with Increasing Flow
For a graph consisting of parallel paths P 1 and P 2 , with len P 1 < len P 2 (see Figure 2), Theorem 2 says that the dynamics converges to P 1 when the incoming flow increases multiplicatively or additively with time.
To prove this, we show that as time progresses, relatively more pheromone is accumulated on P 1 than P 2 .
Let s 1 and s 2 be the neighboring vertices of s, and d 1 and d 2 be the neighboring vertices of d, on path P 1 and P 2 respectively. Note that only the pheromone level on edges (s, s 1 ), (s, s 2 ), (d 1 , d), (d 2 , d) affects the dynamics for this graph. Consider the ratio of pheromone levels on (s, s 1 ) and (s, s 2 ) at time (t + 1): where in the last equation, we set leakage l P 1 = l P 2 = 0, and write the backward flow in terms of the normalized pheromone level. For simplicity, let us assume f← d (t) = α t , for some α > 1. This gives As len P 1 < len P 2 , we know α (t−len P 1 +1) > α (t−len P 2 +1) . These backward flow terms, α (t−len P 1 +1) and α (t−len P 2 +1) , are the main reason why relatively more pheromone accumulates on ss 1 compared to ss 2 as time progresses. However, the ratio p ss 1 (t) p ss 2 (t) may not increase monotonically at each time step. To circumvent this issue, we carefully construct a potential function which increases monotonically with time. Our potential function is given by the minimum of the ratio of the pheromone levels p ss 1 (t) p ss 2 (t) and , and L def = max(len P 1 , len P 2 ). Our potential function is given by Using the definition of the linear decision rule, we know that These inequalities give us where we used α (t−len P 1 +1) > α (t−len P 2 +1) for the last inequality. This gives us p ss 1 (t+1) p ss 2 (t+1) = r ss 1 (t + 1) > r min (t). Thus the backward flow terms α (t−len P 1 +1) and α (t−len P 2 +1) ensure that the pheromone ratio at the edges incident on s at time t + 1 is greater than r min (t), the minimum of the pheromone ratios at the edges incident on s and d across last L time steps. Similarly, the forward flow ensures r d 1 d (t + 1) > r min (t). This implies that r min (t) never decreases, and strictly increases every L time steps. We use this to show convergence to the shortest path.
The proof for the case involving leakage with fixed flow (Theorem 1) is similar and uses the same potential function. The only difference is that in this case, the potential function goes up due to the leakage terms 1 − l P 1 and 1 − l P 2 (Equation 8), instead of the α (t−len P 1 +1) and α (t−len P 2 +1) terms.

General Rules and Fundamental Limits
In Theorem 3, we claim that for any pheromone based rule, bi-directional flow is necessary for convergence to the shortest or the minimum leakage path. When there is unidirectional flow from s to d, the pheromone levels on the edges incident on s is only a function of their initial pheromone levels and forward flow at s. It does not depend on the flow and pheromone levels on the rest of the graph. Therefore, in the case of two parallel paths, for a given decision rule and initial setting of the pheromone levels, if the dynamics converges to a particular path, then they will converge to the other path if we swap the initial pheromone levels on the two edges incident on s. Hence, we can always set the initial pheromone levels on the edges incident on s, such that the dynamics does not converge to the shortest or the minimum leakage path.
Next, we provide a proof sketch for Theorem 4. The proof sketch for Theorem 5 is similar. Let s 1 and s 2 be the neighboring vertices of s, and d 1 and d 2 be the neighboring vertices of d, on path P 1 and P 2 respectively (see Figure 2). Let P 1 be the minimum leakage path. We would show that the dynamics is not guaranteed to converge to path P 1 for non-linear g ∈ F . For any non-linear g ∈ F , consider some r ∈ (0, 1/2) such that g(r) ̸ = r. Such an r exists because g is non-linear. Consider the two cases: 1) g(r) < r, 2) g(r) > r.
Suppose g is such that g(r) < r. Consider an instance where initial normalized pheromone levels p → ss 1 (0) and p ← d 1 d (0) are at most r. And the initial flow values on the edges of path P 1 and P 2 are at most r and at least 1 − r respectively. Let the incoming forward and backward flow be equal to 1 at all times. The decision rule g sends at most r − c g,r amount of forward and backward flow on path P 1 (because g(r) < r and g is monotone) and at least 1 − r + c g,r on path P 2 , where c g,r is a positive constant dependent on g and r. The positive constant c g,r ensures that we can set the leakage levels to be small enough, satisfying l P 1 < l P 2 , such that even after leakage, the flow levels on P 1 and P 2 remain at most r and at least 1 − r respectively, throughout the future. And normalized pheromone levels p → ss 1 Using a similar idea, in the case when g(r) > r, we can set the initial flow and pheromone levels and the leakage parameters, such that the normalized pheromone levels p → ss 2 , p ← d 2 d and the flow levels on P 2 never fall below r and the flow levels on P 1 remain at most 1 − r. Thus, the dynamics never converges to the minimum leakage path P 1 .

Discussion
Like ant colonies, engineered systems such as molecular robots and swarm computing [RAN12; Lun+10; Bra+13; WPN14; HM15] involve a large population of individuals lacking central control and equipped with minimal computational resources. Searching for a target [HWN13; OB07], and in particular, finding the shortest path [Szy+06] is a basic task for such systems. Our algorithms based on leakage and increasing flow can also be applied to such swarms of robots equipped with the ability to release and detect pheromone [Na+21; Arv+15; Kur+09; Fuj+14; Rus99], to solve the shortest path problem and its variants.
Our result on convergence to the shortest path suggests that an ant colony has the ability to discover the shortest path merely by increasing its flow rate. An interesting direction for future research would be to empirically investigate the relationship between flow rate of ants and path length, and understand whether ant colonies increase flow rates to find short paths.
Our results also open up avenues for further theoretical investigation. While the algorithms designed by humans are often set up so as to be amenable to analysis, nature is not constrained in this way. For seemingly simple models of biological systems, it has thus been notoriously difficult to devise mathematical guarantees on the quality of the solutions produced [TKN07; MO+07; MO08; BMV12; Cha12a; Cha17; Bha+13; Cha14; Cha12b] . In our model, analyzing the dynamics is challenging as it involves understanding the progression of pheromone level with time, which is affected by the actions of a large number of agents (modeled by flow), and can be highly non-monotone even for the simple case of graphs with parallel paths. The behavior of this model in extensive simulations suggests that it should be possible to significantly generalize the results we prove here. In particular, we conjecture that provable convergence to the shortest or the minimum leakage path holds for general graphs (Conjecture 1). Another direction for future research is to extend our analysis to the case when multiple terminal vertices are present in the graph, as trail networks in nature usually include many nests and food sources.
In summary, our model for how ant trails change over time contributes to the synergistic exchange between biology and computer science, providing a plausible explanation for how turtle ant colonies can find paths that minimize the number of vertices, and suggesting a surprising algorithm for the shortest path discovery, by increasing the flow rate, applicable to distributed engineering systems. [Fuj+14] Ryusuke Fujisawa, Shigeto Dobata, Ken Sugawara, and Fumitoshi Matsuno. "Designing pheromone communication in swarm robotics: Group foraging behavior mediated by chemical substance". In: Swarm Intelligence 8.3 (2014), pp. 227-246. [Gor12] Deborah M Gordon. "The dynamics of foraging trails in the tropical arboreal ant Cephalotes goniodontus". In: PLoS One 7.11 (2012), e50472.
[NB14] Saket Navlakha and Ziv Bar-Joseph. "Distributed information processing in biological and computational systems". In: Communications of the ACM 58.

Supplementary Material A Proof of Convergence
Here we provide proof for all the results stated in Section 3.1.

A.1 Fixed flow with different leakage on each path (Theorem 1)
Here, we provide a proof of Theorem 1. We restate it below. Let s 1 , s 2 be the neighboring vertices of s that belong to paths P 1 and P 2 respectively. Similarly let d 1 and d 2 be the corresponding neighbors for d. Let f→ s , f← d > 0 be some fixed forward and backward flow values, such that f→ s (t) = f→ s and f← to be the relative pheromone levels at (s, s 1 ) and (d 1 , d) respectively. For notational simplicity, we define m def = len P 1 and n def = len P 2 to be the lengths of path P 1 and P 2 , and L def = max(m, n). We also define α def = 1 − l P 1 and Our potential function at any time t ≥ L is given by the minimum of the relative pheromone levels r ss 1 (t) and r d 1 d (t) across the last L time steps: We divide our proof into 3 steps: • Step 1: In this step, we show that r min (t) is non-decreasing at every time step and increases by a factor of γ(t) every L time steps, for all t ≥ L. Here, γ(t) is some appropriately defined function which is greater than 1 for all t ≥ L.
• Step 2: In this step, we will give a lower bound γ l > 1, on γ(t), to show that r min (t) increases sufficiently every L time steps.
• Step 3: We will find the rate of convergence based on the rate of increase shown in step 2. Now, we give proofs for each of these three steps.
Step 1. We will use the following lemma for the proof.

Lemma 1.
For any time t ≥ L, the following is true, Proof. Under the linear decision rule we know that, Let and from the definition of r min (t), we know that min(r d 1 d (t − m + 1), r d 1 d (t − n + 1)) ≥ r min (t).
Now we show that r min (t) is non-decreasing at every time step and increases by a factor of γ(t) every L time steps, for all t ≥ L. To show this, we would show that r ss 1 (t + 1) ≥ r min (t)γ s (t), and r d 1 d (t + 1) ≥ r min (t)γ d (t). Here, γ s (t) and γ d (t) are appropriately defined functions such that γ s (t)) > 1 and γ d (t)) > 1 for all t ≥ L.
This would give us r min (t For all time t ≥ L, r ss 1 (t + 1) ≥ r min (t)γ s (t), for some appropriately defined function γ s (t), such that γ s (t) > 1 for all t ≥ L.
Proof. The pheromone levels on edges (s, s 1 ) and (s, s 2 ) at time t + 1 is provided by the following expressions.
Therefore, for t ≥ L, we get By the definition of r ss 1 (t) and linear decision rule we know that, From Lemma 1, we know that For t ≥ L, define Note that, from the definition of r min (t) and using Lemma 1, we know that a(t) ≥ 1 and b(t) ≥ 1, for all t ≥ L. Now, we can write Substituting this in Equation (17), we get where we define Since a(t) ≥ 1, b(t) ≥ 1 and α > β, we get that γ s (t) > 1 for all t ≥ L.
This completes the proof of the Lemma.
Step 2. In this step, we will give a lower bound γ l > 1, on γ(t), to show that r min (t) increases sufficiently every L time steps. If the pheromone levels on the edges are too high as compared to the flow, it will take more time for the relative pheromone levels to change. We will first show that there exists a time T 1 , such that for t ≥ T 1 , the pheromone levels and flow are comparable. Our lower bound γ l ≤ γ(t) will hold for all t ≥ T 1 .

Lemma 3.
In the flow dynamics governed by the linear decision rule, the pheromone level on any edge e = (u, v) is always bounded as follows: Proof. Consider any edge e = (u, v) and time t > 0, Since , we get that for t ≥ T 1 , In Equation (19), we defined a(t) and b(t) such that Similarly, for t ≥ L, we define c(t) r min (t) . From the definition of r min (t), we know that c(t) ≥ 1 for all t ≥ L. We will show a upper bound on c(t) in terms of b(t) which will be useful later.
Proof. We know which gives which gives Finally, from the definition of r min (t), we know Since p ← Substituting the value p ← This finishes the proof of the Lemma. Now, we come to the main part of step 2 where we show γ(t) > γ l for all t ≥ T 1 + L, for some γ l > 1.
, · · · , γ d (t)}. Therefore if we show that γ s (t) ≥ γ s l and γ d (t) ≥ γ d l for all t ≥ T 1 + L, for some γ s l , γ d l > 1, we can set γ l = min(γ s l , γ d l ), and we will be done.
Below we will prove γ s (t) ≥ γ s l . The proof for γ d (t) ≥ γ d l is similar.
Proof. From Lemma 2, we know that where We want to show that γ s (t) ≥ γ s l for some fixed constant γ s l > 1. Using a(t) ≥ 1 and b(t) ≥ 1, we get Now, to lower bound γ s (t), we need to upper bound . Using the definition of normalized pheromone level, we can write where we used the upper bound for pheromone level shown in Lemma 3. Similarly, we can write Combining Equation (37) and Equation (38), we get where we define C f→ . Now, we upper bound . From the definition of normalized pheromone level, we get We know r d 1 d (t − n + 1) = c(t)r min (t) where c(t) ≥ 1, and we earlier defined r ss 1 (t) = a(t)r min where a(t) ≥ 1. This gives us In Lemma 4, we show that (1 + c(t)r min (t)) ≤ b(t)(1 + r min (t)), which gives us Substituting this in Equation (35), we get , f← d ,δ +β . Since γ s l > 1, this completes the proof of the lemma.
Using a very similar proof, we can bound γ d (t) by γ d l . Finally, setting γ l = min(γ s l , γ d l ), we will get a lower bound γ l > 1 on γ(t), for all t ≥ T + L. This finishes the proof for step 2.
Step 3. In step 1, we show that r min (t) is non-decreasing for t ≥ L, therefore r min (L + T 1 ) ≥ r min (L). And from step 2, we know that r min (t) increases at least by a factor of γ l every L time steps, for all t ≥ L + T 1 .

This gives us r min
. For t ≥ L + T 1 + T 2 time steps, we would get that r min ≥ 2 ϵ , which implies p → ss 1 Since the pheromone level on all edges on P 1 is always non-zero, we trivially know that p → uv (t) = 1 for (u, v) ∈ P 1 \ (s, s 1 ) and p ← uv (t) = 1 for (u, v) ∈ P 1 \ (d 1 , d). This gives us that for t ≥ L + T 1 + T 2 , normalized pheromone level on all edges on P 1 are at least 1 − ϵ, where This completes the proof of Theorem 1.

A.2 Increasing flow with no leakage on both paths (Theorem 2)
Here, we provide a proof of Theorem 2. We restate it below. Theorem 2. Consider a graph G consisting of two parallel paths P 1 and P 2 from s to d. Let the flow and the pheromone levels be updated according to the model in Section 2, and let P 1 be the shorter path. If (i) the initial pheromone level p uv (0) is positive for all edges (u, v) ∈ P 1 , (ii) leakage l P 1 = l P 2 = 0, and (iii) the incoming flow increases as follows: Let s 1 , s 2 be the neighboring vertices of s that belong to paths P 1 and P 2 respectively. Similarly let d 1 and d 2 be the corresponding neighbors for d. Define r ss 1 (t) def = p ss 1 (t) p ss 2 (t) and r d 1 d (t) def = p d 1 d (t) p d 2 d (t) to be the relative pheromone levels at (s, s 1 ) and (d 1 , d) respectively. For notational simplicity, we define m def = len P 1 and n def = len P 2 to be the lengths of path P 1 and P 2 , and L def = max(m, n).
Our potential function at any time t ≥ L is given by the minimum of the relative pheromone levels r ss 1 (t) and r d 1 d (t) across the last L time steps: The proof is similar to the proof for the case of leakage with constant flow. We divide our proof into 3 steps: • Step 1: In this step, we show that r min (t) is non-decreasing at every time step and increases by a factor of γ(t) every L time steps, for all t ≥ L. Here, γ(t) is some appropriately defined function which is greater than 1 for all t ≥ L.
• Step 2: In this step, we will lower bound γ(t) to show that r min (t) increases sufficiently every L time steps.
• Step 3: We will find the rate of convergence based on the rate of increase shown in step 2. Now, we give proofs for each of these three steps.
Step 1. We will use the following lemma for the proof.
Lemma 6. For any time t ≥ L, the following is true, Proof. Under the linear decision rule we know that, and from the definition of r min (t), we know that min(r d 1 d (t − m + 1), r d 1 d (t − n + 1)) ≥ r min (t).
Now we show that r min (t) is non-decreasing at every time step and increases by a factor of γ(t) every L time steps, for all t ≥ L. To show this, we would show that r ss 1 (t + 1) ≥ r min (t)γ s (t), and r d 1 d (t + 1) ≥ r min (t)γ d (t). Here, γ s (t) and γ d (t) are appropriately defined functions such that γ s (t)) > 1 and γ d (t)) > 1 for all t ≥ L.
Lemma 7. For all time t ≥ L, r ss 1 (t + 1) ≥ r min (t)γ s (t), for some appropriately defined function γ s (t), such that γ s (t) > 1 for all t ≥ L.
Proof. The pheromone levels on edges (s, s 1 ) and (s, s 2 ) at time t + 1 is provided by the following expressions.
Therefore, for t ≥ L, we get where we used l P 1 = l P 2 = 0 in Equation (50) and express the flow in terms of normalized pheromone level in Equation (51). By the definition of r ss 1 (t) and linear decision rule we know that, From Lemma 6, we know that For t ≥ L, define Note that, from the definition of r min (t) and using Lemma 6, we know that a(t) ≥ 1 and b(t) ≥ 1, for all t ≥ L. Now, we can write Substituting this in Equation (51), we get where we define Since the flow is flow is increasing with time, f← This completes the proof of the Lemma.
Step 2. In this step, we will lower bound γ(t) to show that r min (t) increases sufficiently every L time steps.
For the multiplicative increase case, we will show that γ(t) ≥ 1 + c 1 α len P 2 for some constant c 1 > 0. For the additive increase case, we will show that γ(t) ≥ 1 + c 2 t for some constant c 2 > 0. If the pheromone levels on the edges are too high as compared to the flow, it will take more time for the relative pheromone levels to change. So we will first show that there exists a time T 1 , such that for t ≥ T 1 , the pheromone levels and flow are comparable. Our lower bounds for γ(t) will hold for all t ≥ T 1 .

Lemma 8.
In the flow dynamics governed by the linear decision rule, the pheromone level on any edge e = (u, v) is always bounded as follows: for all Proof. Consider any edge e = (u, v) and time t > 0, where we used that f→ s (t) and f← , we get where we used that f→ s (t) and f← d (t) are monotonically increasing. Since , we get that for t ≥ T 1 , In Equation (53), we defined a(t) and b(t) such that . From the definition of r min (t), we know that c(t) ≥ 1 for all t ≥ L. We will show a relationship between b(t) and c(t) which will be useful later.
Proof. We know which gives which gives Finally, from the definition of r min (t), we know Substituting the value p ← This finishes the proof of the Lemma. Now, we come to the main part of step 2 where we lower bound γ(t) for all t ≥ T 1 + L. From step 1, we know that Therefore, to lower bound γ(t), we need to lower bound γ s (t) and γ d (t).
Below we will prove a lower bound for γ s (t).
Lemma 10. For all time t ≥ L + T 1 , r ss Proof. From Lemma 7, we know that Using a(t) ≥ 1 and b(t) ≥ 1, we get Now, to lower bound γ s (t), we need to upper bound . Using the definition of normalized pheromone level, we can write where we used the upper bound for pheromone level shown in Lemma 8. Similarly, we can write Combining Equation (71) and Equation (72), we get . From the definition of normalized pheromone level, we get where c(t) ≥ 1, and we earlier defined r ss 1 (t) = a(t)r min where a(t) ≥ 1. This gives us In Lemma 9, we show that (1 + c(t)r min (t)) ≤ b(t)(1 + r min (t)), which gives Substituting this in Equation (74), we get Using a similar argument, we can bound γ d (t) as ) . Now, using these lower bounds on γ s (t) and γ l (t), we can lower bound γ(t) for the multiplicative increase and additive increase case.
Lemma 11. Consider the multiplicative increase case, that is, when f→ s (t) = α t f→ s (0) and f← Proof. From Lemma 10, we know that Using n > m, this gives us Lemma 12. In the additive increase case, that is, when f→ s (t) = f→ s (0) + αt and f← for all t ≥ L + T 1 .
Proof. From Lemma 10, we know that Step 3. In this step, we use the lower bounds on γ(t) shown in step 2, to find the rate of convergence.
Lemma 13. Consider the multiplicative increase case, that is, when f→ s (t) = α t f→ s (0) and f← Here, Proof. In step 1, we show that r min (t) is non-decreasing for t ≥ L, we know that r min (L + T 1 ) ≥ r min (L). And from step 2, we know that r min (t) increases at least by a factor of C f→ . For t ≥ L + T 1 + T 2 time steps, we would get that r min ≥ 2 ϵ , which implies p → ss 1 Since the pheromone level on all edges on P 1 is always non-zero, we trivially know that p → uv (t) = 1 for (u, v) ∈ P 1 \ (s, s 1 ) and p ← uv (t) = 1 for (u, v) ∈ P 1 \ (d 1 , d). This gives us that for t ≥ L + T 1 + T 2 , normalized pheromone level on all edges on P 1 are at least 1 − ϵ.
To get the desired upper bound on T 2 , we use x 1+x ≤ log(1 + x) for all x ≥ −1. Using this we can write Substituting this in the expression for T 2 , we get the desired bound.
Note that the expression for T 2 involves a α m α−1 term. So if α is a constant, this would give exponential dependence on the length of the shortest path m. But by setting α ≈ 1 + 1 m , we would get α m α−1 ≈ m (for m large enough), making the dependence on the length of the shortest path linear.

Lemma 14.
Consider the additive increase case, that is, when f→ s (t) = f→ s (0) + αt and f← Here, Proof. From step 2, we know that r min (t) increases at least by a factor of γ(t) every L time steps for t ≥ L + T 1 , where Note that for t ≥ T 2 , This gives us that r min (t) increases at least by a factor of exp C δ 2t every L time steps for t ≥ L + T 1 + T 2 . In step 1, we show that r min (t) is non-decreasing for t ≥ L, so we know that r min (L + T 1 + T 2 ) ≥ r min (L).
Define C δ,L = 2L C δ and T def = L + T 1 + T 2 . Also, suppose t is of the form T + kL for some positive integer k. Then this gives us 1 c+i ≥ log( x+c c ) for the last inequality. From here, we get that for t ≥ T 2 r min (L)ϵ C δ,L , r min (t) ≥ 2 ϵ . In the calculations above, we also assumed that t is of the form T + kL for some positive integer k. But as we know that r min (t) is monotonically increasing for all t ≥ L, r min (t) ≥ 2 ϵ holds for all Since the pheromone level on all edges on P 1 is always non-zero, we trivially know that p → uv (t) = 1 for (u, v) ∈ P 1 \ (s, s 1 ) and p ← uv (t) = 1 for (u, v) ∈ P 1 \ (d 1 , d). This gives us that for t ≥ T This completes the proof of Theorem 2.

A.3 Connecting leakage and number of vertices
Here, we show that for any graph (not necessarily the one with parallel paths), the path with minimum leakage also has approximately the minimum number of vertices as long as variation in leakage between different vertices is not too large. (The leakage we consider in the lemma below are only for non-terminal vertices. By convention, we do not have leakage at terminal vertices in our model.) Lemma 15. Suppose for all pairs of vertices u and v, log(1 − l u ) and log(1 − l v ) are within a (1 + ϵ) factor of each other, where l v ∈ (0, 1) for all v. Then the path with the minimum leakage has number of vertices at most (1 + ϵ) times the path with the minimum number of vertices.
Proof. Let max v (1 − l v ) = α for some α ∈ (0, 1). Then, since log(1 − l u ) and log(1 − l v ) are within a factor of (1 + ϵ) for all pairs of vertices u and v, we get Let P 1 be the path with the minimum number of vertices between s and d. And let the number of vertices on P 1 between equal to m (not counting s and d). Let P 2 be the path with the minimum leakage among all paths between s and d, and let the number of vertices on P 2 be equal to n (not counting s and d).
Then, by the bounds on leakage, we get bounds on path leakage for P 1 and P 2 , Since l P 2 ≤ l P 1 (P 2 is the path with the minimum leakage), we get Therefore, the path with the minimum leakage has number of vertices at most (1 + ϵ) times the path with the minimum number of vertices.
We do not consider the degenerate case of l v = 0 or l v = 1 in the above lemma as in that case, the condition on leakage would imply that either all Proof of Theorem 3. Let P 1 and P 2 denote be two parallel paths between s and d, and s 1 and s 2 be neighboring vertices of s on P 1 and P 2 respectively. Without loss of generality, let P 1 be the minimum leakage or the shortest path, and let the flow be unidirectional from s to d.
As the flow is unidirectional, we get that the pheromone levels on the edges incident on s is only a function of their initial pheromone levels and forward flow at s. That is, for some function F, p ss 1 (t) and p ss 2 (t) can be written as p ss 1 (t) = F(p ss 1 (0), p ss 2 (0), f→ s (0), f→ s (1), . . . f→ s (t − 1)) p ss 2 (t) = F(p ss 2 (0), p ss 1 (0), f→ s (0), f→ s (1), . . . f→ s (t − 1)) Given (p ss 1 (0), p ss 2 (0)), and f→ s (t) for all t ≥ 0, suppose the dynamics converges to P 1 . Now if we swap the initial pheromone level on the two edges incident to s, then the dynamics would converge to P 2 . Therefore, for one of these two initial pheromone settings, the dynamics does not converge to the minimum leakage or the shortest path.

B.2 Necessity of the linear rule for convergence to the minimum leakage path (Theorem 4)
Theorem 4. Consider any graph G consisting of two parallel paths from s to d. When the incoming flow is fixed, for every non-linear decision rule g ∈ F , there exists a setting of leakage parameters and initial pheromone and flow levels dependent on g, such that the dynamics does not converge to the path with the minimum leakage.
Proof of Theorem 4. For notational convenience, we define α = 1 − l P 1 and β = 1 − l P 2 . Without loss of generality we assume P 1 to be the path of minimum leakage and therefore we let α ≥ β. Let n and m be the number of vertices between s and d on paths P 1 and P 2 respectively. We name the vertices from left to right on path P 1 by v 0 to v n+1 and P 2 by u 0 to u m+1 , with the convention v 0 = u 0 = s and v n+1 = u m+1 = d.
We also let s 1 = v 1 , s 2 = u 1 , d 1 = v n and d 2 = u m . Let f→ s and f← d denote the fixed incoming forward and backward flow. We divide the analysis into two cases. For any non-linear g ∈ F , there exists an r ∈ (0, 1/2) such that g(r) ̸ = r. For such an r, one of the following two conditions holds: • g(r) < r.
Case 1: Suppose g(r) < r, then pick the following initial configuration: • Assign pheromone value on the edges (s, s 1 ), (s, s 2 ), (d 1 , d), (d 2 , d) such that the normalized pheromone level on (s, Further assign the flow values on the edges such that they satisfy the following, (Since we do not assume any leakage at the terminal vertices while defining the path leakage (see definition 1), we set l v 0 = l u 0 = l v n+1 = l u m+1 = 0 in the expressions above. ) • Let the leakage value at each vertex be such that the parameters α and β satisfy the following inequality, α β ≤ 1 + c g,r 4·r · min , where c g,r def = r − g(r) > 0. As the forward and backward flow is fixed, the above constraint on the leakage parameters only depend on g and therefore satisfy the conditions of the theorem.
We show by induction that the above inequalities on flow and pheromone levels hold for all time t ≥ 0 and therefore the system does not converge to the minimum leakage path.

Hypothesis:
At time t ≥ 0, pheromone values on the edges (s, s 1 ), (s, s 2 ), (d 1 , d), (d 2 , d) are such that the normalized pheromone level on (s, Further the flow values on the edges satisfy the following, Base case: The conditions trivially hold at time 0 because of the initial setting of flow and pheromone levels described above.

Induction
Step: To prove the hypothesis for time t + 1, all we need to show is that, p → ss 1 all the remaining inequalities follow from these basic inequalities. Also note that from the definition of our case, that is g(r) < r, we get that the inequalities p → ss 1 and similarly, In the above we used monotonically non-decreasing property of decision rule g. Therefore it is enough to show that p → ss 1 We know by the induction step that, Also note that, In the above we used the monotonically non-decreasing property of x/(1 + x) and the conditions provided by the induction step at time t, that is f ← In the fourth inequality, we used α − β ≥ 0 and r ≥ 0.
Let c g,r def = r − g(r) > 0 and note that we have the following upper bound on the flow value, In the above we used p → ss 1 (t) ≤ r and monotonicity of decision rule g. Using these bounds, we provide an upper bound for the normalized pheromone level.
In the second inequality we used Equations (84) and (85). In the third equality, we rearranged the terms.
In the fourth inequality, we used f ←  (t + 1) ≤ r , and we conclude the first case.

Case 2:
Suppose g(r) > r, then pick the following initial configuration: • Assign pheromone value on the edges (s, s 1 ), (s, s 2 ), (d 1 , d), (d 2 , d) such that the normalized pheromone level on (s, Further assign the flow values on the edges such that they satisfy the following, (Since we do not assume any leakage at the terminal vertices while defining the path leakage (see definition 1), we set l v 0 = l u 0 = l v n+1 = l u m+1 = 0 in the expressions above. ) • Let the leakage value at each vertex be such that the parameters α and β satisfy the following , where c g,r def = g(r) − r > 0. As the forward and backward flow is fixed, the above constraint on the leakage parameters only depend on g and therefore satisfy the conditions of the theorem.
We show by induction that the above inequalities on flow and pheromone levels hold for all time t ≥ 0 and therefore the system does not converge to the minimum leakage path.
Base case: The conditions trivially hold at time 0 because of the initial setting of flow and pheromone levels described above.

Induction
Step: To prove the hypothesis for time t + 1, all we need to show is that, p → ss 2 all the remaining inequalities follows from these basic inequalities. Also, from the definition of our case, that is g(r) > r, we get that the inequalities p → ss 2 To see this, note that f → ss 2 when p → ss 2 (t + 1) ≥ 1/2, where we used monotonically non-decreasing property of g for the above inequalities. Also, we used g(1/2) = 1/2 and r ≤ 1/2 for the last inequality. Here, we had to consider two cases because function g takes as input the minimum of the two normalized pheromone levels. Similarly we can show f ← Therefore it is enough to show that p → ss 2 (t + 1), p ← d 2 d (t + 1) ≥ r. As the proof for the bound on p ← d 2 d (t + 1) is analogous to that of p → ss 2 (t + 1), in the remainder we focus our attention on the proof for p → ss 2 (t + 1). Recall the definition of p → ss 2 (t + 1), We know by the induction step that, Also note that, In the above we used the monotonically non-decreasing property of x/(1 + x) and the conditions provided by the induction step at time t, that is f ← In the fourth inequality, we used α − β ≥ 0 and 1 ≥ r ≥ 0.
Using these bounds, we provide an upper bound for the normalized pheromone level.
In the second inequality we used Equations (99) and (100). In the third equality, we rearranged terms. In the fourth inequality, we used

B.3 Necessity of the linear rule for convergence to the shortest path (Theorem 5)
Theorem 5. Consider any graph G consisting of two parallel paths from s to d with a unique shortest path. When the leakage is zero for all the vertices, for every non-linear decision rule g ∈ F , there exists a setting of initial pheromone and flow levels, with incoming flow increasing by a fixed multiplicative factor at each time step, such that the dynamics does not converge to the shortest path. The multiplicative factor and initial pheromone and flow levels are chosen as a function of g.
Proof of Theorem 5. Let n and m be the number of vertices between s and d on paths P 1 and P 2 respectively. Without loss of generality we let P 1 to be the shortest path (that is, n < m). We name the vertices from left to right on path P 1 by v 0 to v n+1 and P 2 by u 0 to u m+1 , with the convention v 0 = u 0 = s and v n+1 = u m+1 = d.
We also let s 1 = v 1 , s 2 = u 1 , d 1 = v n and d 2 = u m . We divide the analysis into two cases.
For any non-linear g ∈ F , there exists an r ∈ (0, 1/2) such that g(r) ̸ = r. For such an r, one of the following two conditions holds: g(r) < r or g(r) > r.

C Simulation Details
In Section 3, we discussed various simulation results for linear and non-linear decision rules. We give more details of these simulations here.
Graph families considered. We ran the simulations for three kinds of directed graphs: graphs sampled from the G(n, p) model, graphs sampled from the G(n, p) model with the additional locality constraint that an edge can exist between vertex i and j only if |i − j| ≤ k for some parameter k (with the source vertex and the destination the vertex numbered 1 and n respectively), and the grid graph. G(n, p) model with the additional constraint ensures that the graph has long paths between the source and the destination. We only considered instances where there was at least one path from the source to the destination.
For the G(n, p) model and its locally constrained version, we consider two kinds of graph families: one where an edge is allowed from i to j only if i < j (resulting in a DAG), and the other where the edges can go both ways. We will use G(n, p) and G(n, p) local to denote the G(n, p) graph and it's locally constrained version where edges can go both ways, and G(n, p) DAG and G(n, p) local DAG to denote their acyclic versions respectively.

C.1 Linear Decision Rule
In Section 3.1, we discussed that with linear decision rule, in the fixed incoming flow setting, the dynamics converges to the path with the minimum leakage. With increasing flow, and no leakage, the dynamics converges to the shortest path. We generated a large number of instances with different values for n, p, k and other parameters and observed the desired convergence in all the simulated instances. Below, we describe the details of the parameter settings we considered.
For all the simulations, the decay parameter δ was set to 0.9. The initial forward flow level at s and backward flow level at d was chosen uniformly at random from (0.5, 1), and the initial flow at all other vertices was set to 0. We ran each simulation instance with two settings for initial pheromone level: (i) setting initial pheromone levels at all edges to 1, (ii) choosing initial pheromone level on the edges uniformly at random from (0, 1).

C.1.1 Leakage with fixed flow
For all the randomly generated instance, we chose the leakage value of all the vertices uniformly at random from (0, 0.1). For all these instances, the dynamics converged to the path with the minimum leakage. Below, we describe the parameter values for different graph families.
3. Grid graph: We ran the simulations on a 10X10 grid graph with source and destination at the diagonally opposite extreme vertices. We ran 100 instances. The graph structure was same across all these instances, but other parameter values such as leakage levels, initial pheromone levels etc. were chosen randomly as described above.

C.1.2 Increasing flow with no leakage
In this case, the leakage was set to 0 for all the vertices. We increase the forward flow level at s and backward flow level at d by a factor of 1.1 every time step. To avoid floating point overflow, instead of multiplying the forward and backward flow values at s and d by 1.1, we keep these two values the same, and divide all other flow values and pheromone levels in the graph by a factor of 1.1. For any decision rule that only depends on the normalized pheromone levels (which holds for the linear rule), it is not difficult to see that only the relative values of flow and pheromone level matter, and this gives rise to exactly the same dynamics (up to scaling) as when we multiply the forward flow at s and backward flow at d by 1.1. When the value of pheromone level or flow at any edge became too small (smaller than the minimum allowed value for a 64 bit floating point number which is ≈ 10 −323 ), we rounded it to zero.
For all the graphs, we planted a short path in the graph so as to ensure that the shortest path is unique. In all the simulations, we observed that the dynamics converges to this shortest path. Below, we describe the details for different graph types.

C.2 Nonlinear Decision Rules
In Section 3.2, we showed that for a reasonably general family of decision rules, the linear decision rule is its unique member with guaranteed convergence to the path with minimum leakage or the shortest path. Note that these results only suggest that linear decision rule is necessary for guaranteed convergence to the shortest or the minimum leakage path. Can it be the case that even with non-linear decision rules, the forces of increasing flow and leakage still help in finding shorter or smaller leakage paths respectively, compared to the paths found in the absence of these forces? To understand this question, we ran simulations for various non-linear decision rules. We observe that within each graph family, for a large fraction of graph instances, the path found in the presence of these forces has length (respectively leakage) smaller than or equal to the length (respectively leakage) of the path found in their absence. These simulations suggest that the usefulness of the forces of leakage and increasing flow is not limited to the linear decision rule. Below, we discuss more details of these simulations.

C.2.1 Decision rules details
We consider the following non-linear decision rules: • Quadratic: This rule divides the flow in proportion to the square of the pheromone levels.
• 1.1 power: This rule divides the flow in proportion to the 1.1 th power of the pheromone levels. We study this rule to understand the effect of strength of non-linearity. Since this rule is closer to the linear rule compared to the quadratic rule, we would expect that the forces of leakage and increasing flow are more effective with this rule compared to the quadratic rule.
• Quadratic-with-offset: This rule is a slight variant of the quadratic rule and has been used previously [Den+90] to model ant behaviour. This rule adds a fixed positive constant c to the pheromone levels, and divides the flow in proportion to the square of the offsetted pheromone levels. The parameter c was added to encourage exploration.
• Rank-edge: This rule was introduced in [CGN18]. This rule ranks the edges from highest to lowest pheromone levels. If multiple edges have the same pheromone level, they get the same rank. Let q ∈ (0, 1) be some parameter. This rule sends (1 − q) fraction of flow to the first ranked edge, q(1 − q) fraction of flow to the second ranked edge, q 2 (1 − q) fraction of flow to the third ranked edge and so on. In general, the i th ranked edge gets q i−1 (1 − q) fraction of flow. If there are k edges at the i th rank, then each of them gets fraction of flow. One exception to this rule is the lowest ranked edge which gets all the remaining flow that has not been assigned to any other edge. If there are multiple edges at the lowest rank, then this remaining flow gets divided equally among them.

C.2.2 Results
As in the linear decision rule case, we considered the G(n, p), G(n, p) local, G(n, p) DAG, G(n, p) local DAG and the grid graph for our simulations with non-linear decision rules. In the increasing flow case, similar to the linear decision rule setting, we add a random shortest path to these graphs so that the shortest path is unique. For each non-linear decision rule and each graph family, we ran simulations with 100 random graph instances.
In Table 1, we show the effect of leakage with non-linear decision rules. There is a subtle distinction here between path leakage as an objective function and leakage as a process affecting the dynamics. For each graph instance, we assign leakage values to vertices. This gives us a path leakage objective function for each graph instance. For each graph instance, we compare two settings, the baseline setting in which the leakage process is not applied at vertices during the dynamics and the main setting in which the leakage process is applied. Now for each graph instance, we compare the path leakage objective function of the path found in the baseline setting and the main setting. We report the fraction of instances where the path leakage objective in the main setting is unchanged or smaller as compared to the baseline setting. We also report the average percentage difference in the path leakage objective between the baseline and the main setting, that is 1 n n ∑ i=1 100 * path_leakage(main i ) − path_leakage(baseline i ) path_leakage(baseline i ) .
Here, path_leakage(main i ) and path_leakage(baseline i ) denote the path leakage objective obtained in the main setting and the baseline setting respectively, for the i th graph instance. In the last column, we report the average percentage difference in path leakage for the baseline compared to the optimum (minimum) path leakage, that is 1 n n ∑ i=1 100 * path_leakage(opt i ) − path_leakage(baseline i ) path_leakage(baseline i ) .
Here, path_leakage(opt i ) denotes the optimum path leakage for the i th graph instance. Due to symmetry in the grid graph, the dynamics does not converge to a path in the baseline setting, therefore we do not consider grid graphs while studying the effect of leakage with non-linear decision rules.
In Table 2, we show the effect of increasing flow with non-linear decision rules. Again, we consider the baseline setting where there is no increasing flow, and the main setting with increasing flow. We report the fraction of instances where the path length in the main setting is unchanged or smaller as compared to the baseline setting. We also report the average percentage difference in the length of path found between the baseline and the main setting, that is 1 n n ∑ i=1 100 * path_length(main i ) − path_length(baseline i ) path_length(baseline i ) .
Here, path_length(main i ) and path_length(baseline i ) denote the path length obtained in the main setting and the baseline setting respectively, for the i th graph instance. In the last column, we report the average percentage difference in path length for the baseline compared to the shortest length path, that is Here, path_length(opt i ) denotes the shortest path length for the i th graph instance.
We discuss our major observations below: • Usefulness of forces of leakage and increasing flow not limited to the linear rule: For a large fraction of graph instances in each graph family, the path found in the presence of these forces has length (respectively leakage) smaller than or equal to the length (respectively leakage) of the path found in their absence. This holds for all non-linear decision rules considered.
• Forces of leakage and increasing flow more effective with weaker non-linearity: One would intuitively hope that the closer a non-linear decision rule is to the linear rule, the more effective the forces of leakage and increasing flow would be. To study this, we compare the quadratic decision rule with 1.1 power rule. For all graph families, we observe that the percentage of instances with strictly smaller path leakage (path length resp.) is higher for 1.1 power rule. We also observe that the average percentage change in path leakage (path length resp.) is more negative for 1.1 power rule. This shows that the forces of leakage and increasing flow are more effective for 1.1 power rule (which is closer to the linear rule) compared to the quadratic rule.
• Effectiveness of the forces of leakage and increasing flow varies across graph families: We observe that the effectiveness of leakage and increasing flow varies across graph types. For instance, across all non-linear rules considered, we observe that the percentage of instances with strictly smaller path leakage (path length resp.) is higher for G(n, p) local and G(n, p) local DAG graphs compared to G(n, p) and G(n, p) DAG graphs. We also observe that there are a few instances where these forces end up increasing the path leakage (path length resp.). For instance, with 1.1 power rule and for G(n, p) local graph family, 1% of instances end up with increased leakage compared to the baseline. It is an interesting direction for future research to understand what graph properties affect the effectiveness of leakage and increasing flow with non-linear decision rules.
level at s and backward flow level at d was chosen uniformly at random from (0.5, 1), and the initial flow at all other vertices was set to 0.
• Leakage: For the case when there is leakage at vertices, the leakage at each vertex is chosen uniformly at random from (0, 0.1).
• Increasing flow rate: In the case of increasing flow, we increase the incoming forward and backward flow by a factor of 1.1 at each time step. Similar to the linear decision rule case, to avoid floating point overflow, we implement this by decreasing all other flow and pheromone levels by a factor of 1.1 except the incoming forward and backward flow. For any decision rule that only depends on the normalized pheromone levels , it is not difficult to see that only the relative values of flow and pheromone level matter, and our implementation gives rise to exactly the same dynamics (up to scaling) as when we increased the incoming forward and backward flow by a factor of 1.1. All the non-linear decision rules we consider here only depend on the normalized pheromone levels except the quadratic-with-offset rule. For the quadratic-with-offset rule, we additionally also scale down parameter c by 1.1 in each step, which makes the dynamics equivalent (up to scaling) to the case where incoming flow is multiplied by 1.1. When the value of pheromone level or flow at any edge becomes too small (smaller than the minimum allowed value for a 64 bit floating point number ≈ 10 −323 ), we round it to zero.
• Decision rule parameters: The quadratic and 1.1 power rules do not have any parameters to be set. Quadratic-with-offset rule has a parameter c. Note that when c is large, the dynamics would not converge to a single path. We set c to be small enough such that the dynamics for most instances end up converging to a single path (see definition of convergence below). For the leakage simulations, we set the value of c = 0.25 for all graph instances . For the increasing flow simulations, we set the value of c = 1 for the grid graph instances and c = 0.5 for all other graph instances. Rank-edge decision rule also has a parameter q. We set q = 0.01 for all our simulations.
• Convergence criterion: For all simulations (including linear rule simulations) except the quadraticwith-offset rule simulations, we consider the dynamics to have converged to some path when for each vertex on the path, at least 99% of the forward (backward) flow entering it passes through its outgoing (incoming) edge on the path. For linear rule, this corresponds to each edge on the path having backward and forward normalized pheromone level at least 0.99. With quadratic-with-offset rule, we keep the threshold slightly lower to 95%. If the threshold is too high for this rule, most instances would only converge to a path when c is fairly small. However, with small c, this rule would be similar to the quadratic rule. To differentiate the dynamics with this rule from the quadratic rule, we keep the threshold slightly lower to 95% which allows us to choose a slightly larger value of c.
• Initial Pheromone levels: We set initial pheromone level on all edges equal to one. Recall that with the linear rule, our simulations worked as expected even when the pheromone level at each edge is chosen uniformly at random from (0, 1). However, for non-linear rules, we observed that for most instances, the simulations do not converge to a path when pheromone level is randomly initialized.
In that sense, the linear rule seems more robust to pheromone initialization. Therefore, we run our simulations with non-linear rules with pheromone levels uniformly set to one so that the dynamics converges to a path.
• Ignoring instances that do not converge to a path: We discussed that we set the initial pheromone level uniformly to one so that the dynamics converges to a path. But even in this case, with non-linear decision rules, we observe that for a few instances the dynamics does not converge to a path. In our simulations, we ignored such instances, and only considered the instances where the dynamics converges to a path in both the baseline setting and the setting with leakage or increasing flow.