Longest Path in the Price Model

TL;DR

This paper analyzes the longest paths in directed acyclic graphs generated by the Price model, showing they scale logarithmically with network size and depend on attachment methods.

Contribution

It introduces a reverse greedy path concept and provides analytical and numerical evidence for the logarithmic scaling of longest paths in Price model variants.

Findings

Reverse greedy path scales with log of network size

Longest path length depends on attachment method

Numerical results confirm analytical predictions

Abstract

The Price model, the directed version of the Barab\'{a}si-Albert model, produces a growing directed acyclic graph. We look at variants of the model in which directed edges are added to the new vertex in one of two ways: using cumulative advantage (preferential attachment) choosing vertices in proportion to their degree, or with random attachment in which vertices are chosen uniformly at random. In such networks, the longest path is well defined and in some cases is known to be a better approximation to geodesics than the shortest path. We define a reverse greedy path and show both analytically and numerically that this scales with the logarithm of the size of the network with a coefficient given by the number of edges added using random attachment. This is a lower bound on the length of the longest path to any given vertex and we show numerically that the longest path also scales with the logarithm of the size of the network but with a larger coefficient that has some weak dependence on the parameters of the model.

Tables

Table 1. Table B1: Absolute and relative difference between chi-squared values obtained by fitting path data to a linear function, (3.1), χ lin 2 subscript superscript 𝜒 2 lin \chi^{2}_{\mathrm{lin}} and to a non-linear function f ​ ( t ) = a ​ ln ⁡ ( t ) + b + ( c / t ) 𝑓 𝑡 𝑎 𝑡 𝑏 𝑐 𝑡 f(t)=a\ln(t)+b+(c/t) , χ nonlin 2 subscript superscript 𝜒 2 nonlin \chi^{2}_{\mathrm{nonlin}} . The non-linear fit results in smaller χ 2 superscript 𝜒 2 \chi^{2} values, but the improvement is marginal.

m	$\alpha$	${\chi }_{\mathrm{gr},\mathrm{lin}}^2-{\chi }_{\mathrm{gr},\mathrm{nonlin}}^2$	$\frac{{\chi }_{\mathrm{gr},\mathrm{lin}}^2-{\chi }_{\mathrm{gr},\mathrm{nonlin}}^2}{{\chi }_{\mathrm{gr},\mathrm{lin}}^2}\times {10}^7$	${\chi }_{\max ,\mathrm{lin}}^2-{\chi }_{\max ,\mathrm{nonlin}}^2$	$\frac{{\chi }_{\max ,\mathrm{lin}}^2-{\chi }_{\max ,\mathrm{nonlin}}^2}{{\chi }_{\max ,\mathrm{lin}}^2}\times {10}^7$
2	1	1914.679	$1.8958$	52115.2	51.6
2	2	5500.539	$5.4468$	47292.36	46.9
2	3	10679.31	$10.574$	32425.11	32.1
2	4	15793.18	$15.638$	11345.38	11.2
3	1	4128.344	$4.0881$	53148.67	52.7
3	2	10891.1	$10.784$	30573.03	30.3
3	3	18021.77	$17.849$	14582.06	14.5
3	4	29150.27	$28.865$	535.4066	$0.531$
4	1	6381.005	$6.3185$	35694.73	35.4
4	2	15188.28	$15.041$	13707.66	13.6
4	3	30406.43	$30.109$	10.1752	$0.0101$
4	4	41108.91	40.72	3412.088	$3.39$
5	1	7572.032	$7.4983$	31364.85	31.1
5	2	23749.52	23.52	3088.395	$3.06$
5	3	36774.17	36.424	1128.841	$1.12$
5	4	66946.58	66.314	47772.77	47.3
6	1	11489.57	11.378	22020.32	21.8
6	2	26992.08	26.731	149.04	$0.148$
6	3	61386.35	60.782	45868.17	45.4
6	4	89429.041	88.563	93100.79	92.2

Table 2. Table B2: Variation in obtained fitting parameters for the greedy path scaling and the longest path scaling when using various initial graphs in the Price model with m = 6 𝑚 6 m=6 , α = 2 𝛼 2 \alpha=2 . The data from t = 1 , 000 𝑡 1 000 t=1,000 to t = 10 5 𝑡 superscript 10 5 t=10^{5} was fitted to the form a ​ ln ⁡ ( t ) + b 𝑎 𝑡 𝑏 a\ln(t)+b . The initial graph seems to have a small effect on the slop, parameter a 𝑎 a , but causes significant changes in the intercept, parameter b 𝑏 b .

Path type	Initial graph	$a$	$b$
Greedy	Star	$1.4973\pm 0.0052$	$-0.688\pm 0.04417$
Greedy	Complete	$1.5001\pm 0.0038$	$5.9741\pm 0.0309$
Greedy	Chain	$1.5213\pm 0.0035$	$8.5548\pm 0.0290$
Longest	Star	$3.0701\pm 0.0056$	$-7.9019\pm 0.0476$
Longest	Complete	$2.9970\pm 0.0034$	$0.5458\pm 0.0279$
Longest	Chain	$3.0303\pm 0.0034$	$2.9345\pm 0.02789$

Equations

$$\Pi (t,s)=p\frac{k^{(\mathrm{out})}(t,s)}{E(t)}+\bar{p}\frac{1}{N(t)}\text{mbox}ift\ge s\ge 1,$$

$$\ell (t)=\frac{m\bar{p}}{\ln (2)}\ln (t).$$

$$P(\ell ,t+1)$$

$$G(z,t)$$

$$\underset{t\to \infty }{\lim }\ell (t)$$

$$f(t)=a\ln (t)+b.$$

$${\mathcal{N}}^{(-)}(t)=\{s|(s,t)\in \mathcal{E}\}$$

$$P^{(-)}(t)$$

$$\Pi (t,s)$$

$${\Pi }_{\le }(t,s)$$

$$=$$

$${\Pi }_{\le }(t,s=t)=1,t\ge 1$$

$$k_{\mathrm{eff}}^{(\mathrm{out})}(t,s)$$

$$k_{\mathrm{eff}}^{(\mathrm{out})}(s,s)$$

$${\Pi }_{\le }(t,s)$$

$$=$$

$$k_{\mathrm{eff}}^{(\mathrm{out})}(t+1,s)-k_{\mathrm{eff}}^{(\mathrm{out})}(t,s)$$

$$y(t,s):=k_{\mathrm{eff}}^{(\mathrm{out})}(t,s)+\frac{m\bar{p}s}{p},\text{mbox}withy(s,s)=\frac{ms}{p}.$$

$$y(t+1,s)-y(t,s)$$

$$\Rightarrow y(t,s)$$

$${\Pi }_{\le }(t,s)$$

$${\Pi }_{\le }(t,s)=\frac{j(s)}{j(t)},j(s)=\frac{\Gamma (s+1)}{\Gamma (p+s)}.$$

$${\Pi }_{\le }(t,s)$$

$${\Pi }_{\max }(t,s)$$

$$\sum _{s=1}^t{\Pi }_{\max }(t,s)$$

$$P(\ell ,t+1)$$

$$G(z,t)=\sum _{\ell =0}^{\infty }{z}^{\ell }P(\ell ,t),$$

$$G(z=1,t)$$

$$G(z,t=1)$$

$${\frac{\partial G(z=1,t)}{\partial z}|}_{z=1}$$

$$G(z,t+1)$$

$${\Pi }_{\max }(t,s)$$

Full text

Introduction

The Price model [2, 3] is one of the oldest network models and it was motivated by the pattern of citations in academic papers. In a citation network, each node represents a document while every entry in the bibliography of a document $t$ is represented by a directed edge from an older document, node $s$, to node $t$. One of the key features of a citation network, one inherent in the Price model, is that there is a fundamental arrow of time in the network; bibliographies can only refer to older documents. This means that there are no cycles in the network, you can never find a path from a node that returns to that node. Thus a citation network is an example of a Directed Acyclic Graph (DAG).

Mathematically, DAGs have some distinctive properties and one of them is that for any pair of connected nodes there is a well defined and meaningful longest path length, for example see Fig. 1. Contrast this with, for example, undirected networks, where you can often find many paths between two given vertices that visit most of the nodes in a component so longest paths are often as long as the component is big, if all nodes in the path must be distinct, and infinite, if multiple visits to the same node were allowed. In directed graphs with cycles, the longest path is infinite, if multiple visits to a node are allowed. Both of these definitions of the longest path coincide if the network is acyclic, as the absence of cycles ensures that in any path, a node can only occur once.

In a citation network, it is not clear how useful the shortest path is. For instance, in writing this paper, the oldest citation we have is to a paper by Price [2]. The shortest path to Price’s paper from this work has length one. On the other hand, most of the knowledge of that work contained in this paper did not come directly from that paper. We only reread Price’s paper to check one detail while working on this project. So the length of the shortest path to that paper seems largely irrelevant. Rather, the information in this early bibliometrics paper by Price has reached us through a sequence of other work, much of it not explicitly referenced in our paper. We drew on much more recent documents such as the reference book by Newman [5] which in turn cites papers which developed various aspects of the Price model. Indeed there is much evidence [6, 4, 7] that typically 70% or so of a bibliography may not have been used directly when producing the work in an academic paper.

So our thesis is that for DAGs the longest path plays a much more important role than the shortest path. In simple models the longest path has been shown to be the best approximation to the geodesic for models of DAGs embedded in Minkowski space [8] where there is a single time direction222This has been exploited in real data sets where dimension and curvature of a DAG can be measured [26, 33] enabling us to embed DAGs such as citation networks in Minkowski space [31].. A similar rigorous link for undirected networks has only been made for the shortest path in networks embedded in Euclidean space where there is no arrow of time [9, 10].

The properties of the longest path have been investigated in the context of simple models known as Cube Spaces [11] which include those built from Poisson Point Processes in Minkowski spaces where all causally connected points are connected to form a network. However these are examples of transitively complete DAGs, that is if there is a path between two points then there is always an edge connecting those two points directly. However, that is not true for a citation network where the limited size of a bibliography means no document ever cites every older paper to which it has some connection. What we seek to do in this paper is to look at the properties of the longest path in a simple model, the Price model [3] and its variations, where the network is neither transitively complete nor (except for one parameter value) transitively reduced, the situation for most DAGs in a social context such as citation networks. Can we calculate the length of the longest path in the Price model? How does this length depend on the parameters of the model?

We start by outlining our analytic results. We then compare these predictions to numerical simulations and summarise our findings.

Analytic Results

In the Price model [3] (for instance see sec.14.1 of Newman [5]) we start from a network $G(t)$ defined at an integer ‘time’ $t$. We create a new graph $G(t+1)$ by first adding one new vertex, which we label with the time $(t+1)$. This new node, $(t+1)$, is connected to $m$ existing vertices $s$ in the graph $G(t)$. These $m$ existing vertices $\{s\}$ are each chosen with probability $\Pi(t,s)$. We will use a convention that these edges point from older to newer vertices, from $s$ to $(t+1)$. Once these edges have been added we have our new graph $G(t+1)$. The process is then repeated. For an example of how a network grows according to the Price model, see Fig. 2.

The mathematical and numerical simplicity of this model comes from the simple definition of $\Pi(t,s)$. To define the probability $\Pi(t,s)$ we first define $N(t)=N_{0}+t$ be the number of nodes in the graph $G(t)$ for some constant $N_{0}$. The number of edges in the graph $G(t)$, after all $m$ edges have been added to node $t$, is $E(t)=E_{0}+mt$ where $E_{0}$ is some constant. Finally in the graph $G(t)$ let the node created at time $s$ have out-degree $k^{(\mathrm{out})}(t,s)$, the number of edges leaving $s$ and connecting it to later nodes.

In this model, the connection of edges to new node $(t+1)$ is made in one of two ways. With probability $p$ node $(t+1)$ is connected to an existing vertex $s$ chosen in proportion to the number of edges leaving $s$ for later nodes in $G(t)$, $k^{(\mathrm{out})}(t,s)$. Price called this cumulative advantage and, after normalisation, we have that the probability of choosing $s$ is $k^{(\mathrm{out})}(t,s)/E(t)$. The second process happens with probability $\bar{p}=(1-p)$ and in this case we choose the source vertex $s$ uniformly at random from the set of vertices in $G(t)$, i.e. with probability $1/N(t)$. If we start the process at time equal to $1$, then the probability of connecting the vertex $(t+1)$ to existing vertex $s$ is $\Pi(t,s)$ where

[TABLE]

unless $s=t=1$ when $\Pi(t,s)=1$, otherwise $\Pi(t,s)=0$. Note that in his original paper, Price considered $p=m/(m+1)$ where $\Pi\propto k^{(\mathrm{out})}+1$. This more general form for the attachment probability $\Pi(s,t)$ in (1) has been used in many related contexts since Price, see Newman [5] for a review.

There is an issue about the starting point for this process. The usual form for $\Pi$, the $t\geq s>1$ form in (1), leaves us with a problem for $\Pi(t=1,s=1)$ when looking at the attachment to the second vertex, $t=2$. Our solution is to demand that $\Pi(t=1,s=1)=1$. This fixes the cumulative probability $\Pi_{\leq}$ to have a consistent value which is in fact all we need for this calculation. However, we will also assume that the number of nodes and number of edges are given by $N(t)=t$ and $E(t)=mt$ respectively. This is only needed for $t\geq 2$ so in principle we must allow multiple edges between nodes starting with $m$ edges added between node $2$ and node $1$. Again this cannot be true for at least the first node at $t=1$.

We also note that our analytic calculations allow our networks to contain multiple edges (node pairs linked by more than one edge). Of course, a real citation network and many numerical calculations of this model (though not our numerical calculations) do not have multiple edges. However, in the long time limit the effect of such edges becomes negligible as they form a small fraction of the edge population, a fraction that dies off as a power law in time [5].

Now we would like to define the longest path algebraically. Unfortunately, finding the longest path requires global knowledge of all the paths. This is extremely hard to do algebraically (though is surprisingly straightforward numerically). So the first stage of our calculation is to decide to calculate a path defined with local knowledge only. That is we define what is called a reverse greedy path using an iterative process where at each stage we only need to know about the properties of the next vertex in the path. We will denote the length of the reverse greedy path from the source vertex $s=1$ to a target vertex $t$ as $\ell(t)$. The length of the longest path from the source vertex to a target vertex $t$ will be denoted as $L(t)$.

The reverse greedy path to a node $t$ is a path running from the source node at the initial time $t=1$ to node $t$. This always exists and it is unique. To define it suppose that we have found the reverse greedy path to all earlier nodes. The last step on the reverse greedy path to node $t$ is made along the edge arriving at $t$ from its most recent predecessor node, say $s$. The idea is that the most recent predecessor of node $t$, furthest from the source node in terms of the time, is also the most likely to be the predecessor node furthest from the sink node in terms of network path lengths. There is no guarantee that our reverse greedy path is identical to the longest path, so the reverse greedy path length is a lower bound on the longest path length. A more formal definition is given in Appendix A.1 in the Supplementary Information.

Of course in any one instance of the Price model, this reverse greedy path length will fluctuate if we look at nodes of similar ages, not least because $\ell$ is integer valued. We will use a mean field approach so our $\ell(t)$ is an average over many realisations of the model though for simplicity we will not include the expectation value notation $\langle\ldots\rangle$. For that reason our $\ell(t)$ will be a real valued monotonically increasing function of time $t$.

We can find the long-time behaviour using the following simple argument. On average, there are $pm$ edges added with cumulative advantage at each time step. Suppose we are adding a new node at time $(t+1)$ and we are looking for source nodes $s$ for the $m$ edges we are going to add. The source nodes chosen with cumulative advantage are those with largest degree and those are the oldest nodes created at small values of time $s$ (since $k^{(\mathrm{out})}(t,s)\propto(t-s)^{p}$, for instance see Newman [5]). So the youngest source node chosen, nodes created at the largest value of time $s$, is likely to be one of the $m_{\mathrm{rnd}}=m\bar{p}$ nodes chosen uniformly at random. The probability that all these $m_{\mathrm{rnd}}$ randomly chosen source nodes are chosen between time $1$ and time $\hat{s}$ inclusive is $(\hat{s}/t)^{m_{\mathrm{rnd}}}$. Suppose we consider the time $\hat{s}_{1/2}$ where with probability one half the time coordinate of the largest randomly chosen source node is $\hat{s}_{1/2}$ or less, then this sets the scale for the birth date of the youngest source node connected to node $(t=1)$, namely that $\hat{s}_{1/2}=\mu t$ where $\mu=2^{-m_{\mathrm{rnd}}}$. This is the previous node on the reverse greedy path from the initial node to node $t$. We can then estimate the numbers of steps it takes to get back to the source node at $t=1$ as $\mu^{\ell}t\approx 1$ which leads to

[TABLE]

The simplicity of the attachment probability in the Price model means we can also produce more detailed derivation using a mean-field approach. Let the probability that the length of the reverse greedy path, $\ell$, from new node $(t+1)$ to the initial node at $t=1$, be $P(\ell,t)$. Then the master equation is of the form

[TABLE]

Here $\Pi_{\mathrm{max}}(t,s)$ is the probability that of the $m$ predecessor nodes connected to a new node at $(t+1)$, the oldest of them is $s$. In terms of the generating function $G(z,t)=\sum_{\ell=0}^{\infty}z^{\ell}P(\ell,t)$ we find that the exact solution in the Price model is (see Appendix A.3 in the Supplementary Information for details)

[TABLE]

Exact forms for the expected reverse greedy path length can be found from this expression, especially for specific small values of $m$. However, the leading order contribution for large times is always of the form

[TABLE]

where it is implicit that there is no contribution from the term with the sum for the case of $m=1$. Here $\psi(z)$ is the digamma function and $\zeta(z)$ is the Riemann zeta-function. The details of the calculation are given in Appendix A.3 of the Supplementary Information.

Finally, the scaling properties of the longest path in the Price model suggests that the properties of height antichains are also very simple. The height of a node in a DAG is the length of the longest path to a node from a source node, any node with zero in-degree. Thus in the Price model, the height of a node is simply the length of the longest path length from the initial node to the given node, our $L$.

Nodes connected by a path cannot be of the same height. Thus the subset of all nodes at the same height form an antichain, a set of nodes in which no two are connected by a path [34]. The scaling properties of these height antichains are simple to estimate if we conjecture that the average longest path $L$ of a node $t$, its average height, scales as $\ln(t)$. This suggests that if the median index of a node in an antichain of integer valued height $h$ is $t_{\mathrm{mid}}=(\mu)^{h}$ then the mean index of nodes in the antichain will scale as $\cosh(\sqrt{\mu})(\mu)^{h}$, the variance in the index of nodes in the antichain will be roughly $(1/\sqrt{3})\sinh(\sqrt{\mu})(\mu)^{h}$, and the number of nodes in the antichain will vary as in the $2\sinh(\sqrt{\mu})(\mu)^{h}$.

Numerical Methods and Results

In the master equation (3), multiedges (attaching two edges from the new vertex to the the same vertex) are not excluded. In our numerical implementation code we also allowed multiedges to be created. However the probability of attaching one edge from new vertex $(t+1)$ to any existing vertex $s$ is decreasing as $(s/t)^{p}$ (for instance see p. 489-90 Newman [5]). So the creation of a multiedge becomes negligible at large times hence our networks are essentially the same as implementations of the Price model in which multiedges are excluded.

The first few steps of the numerical implementation of the Price model have some subtleties which are worth mentioning. The problems noted analytically with the initial node at $t=1$, which is the only node with zero-in-degree, exemplify the issue. The earliest nodes are those with the shortest values of our $\ell$ and $L$ path lengths to the first node. Since the path lengths of the first few nodes will be added to any other path routed via one of these early nodes, we expect the initial graph to give a constant contribution to the path lengths we measure but not to alter the growth in length scales over long-times.

We chose to start our simulations from a complete graph of $(2m+1)$ nodes, labelled $t=1$ to $t=(2m+1)$. All pairs of nodes are connected in this initial graph, with the edge direction from earlier to later node. This initial graph ensures that $E(t)=mt$ for all graphs generated numerically, $G(t)$ for $t\geq(2m+1)$. The out-degree distribution is fairly even but the in-degree is not fixed to be $m$ for nodes in this initial graph. Further notes on the effect of the initial graph are given in Appendix B.4 in the Supplementary Information.

In order to simplify and accelerate the numerical analysis, for each new node $(t+1)$ added we drew nodes uniformly at random from an “attachment list” which we maintain. After we have chosen the source edges for the $m$ edges attached to the new vertex, we update our attachment list by adding each source node once for every edge, and we add $(\bar{p}/pm)$ references to the new node. This means we restrict our results to cases where $(\bar{p}/pm)$ is an integer. Drawing nodes uniformly from our attachment list means that we are choosing vertices according the the probability (1). For the special case where there was no cumulative advantage, $p=0$, the attachment list was simply a list of the existing vertices where each is referenced once. A more detailed explanation is given in Appendix B.1 of the Supplementary Information.

For each node $s$, we also record values for the lengths of the reverse greedy path $\ell(s)$ and the longest path $L(s)$. When adding a new node $t$, it is simple to look at the values of the lengths of these paths to the $m$ nodes attached to the new node. From that information, it is simple to record the lengths of these paths to the new vertex, $\ell(t)$ and $L(t)$. Storing and manipulating these results proved to be more of a limiting factor than the speed to produce them. We produced results for networks of up to $10^{8}$ nodes.

The results for these path lengths are quite noisy for any one node as shown by an exemplary run in Fig. 3. Despite the relatively large fluctuations in results for any one node,there is a clear trend in the nodes created at later times. The fluctuations of the path length scaling are greatly reduced when averaged over multiple networks as shown in Fig. 3. So we use 100 runs for each set of parameter values in our work.

In order to compare our numerical data with the analytical results we fitted the path lengths found to the function $f(t)$ where

[TABLE]

The fit was made by using a non-linear fitting routine based on the optimisation of the chi-squared measure of goodness of fit (for instance see [12]) as described in more detail in Appendix B.3 of the Supplementary Information. Errors on parameters were estimated from the covariance matrix produced by such a method. Given that our analytical work only studies the long-time limit and that the early times in the numerical simulation do not satisfy all the conditions of the analytical work, it is not surprising that in Fig. 4 we still see significant difference between the numerical results and analytical predictions for the length of paths from $t=1$ to those nodes created at early times. So when fitting to our numerical data we only use data for nodes created from time $t_{0}=1,000$ up to the last node at $t=10^{8}$. The effect of this cutoff is discussed in Appendix B.3 but we found varying this lower cutoff had little effect on our results since we had so many data points from the region where the asymptotic growth dominates.

The dependence of the coefficient of the $\ln(t)$ term found from the fit, $a$, on the model parameters is shown in Fig. 5 and Fig. 6.

The next-to-leading order coefficient, $b$ of (6), showed no clear trends. We also considered a non-linear fit with a term of $c/t$ added to the expression in the (6). We found that in practice, this term had little influence on the remaining parameters of interest, namely, $a$ and $b$. Furthermore, the errors in $c$ were found to be relatively large in comparison to the errors of the parameters $a$.

Finally, it is clear that the longest path length is scaling as $\ln(t)$ to a good approximation. As noted above this then implies that the properties of height antichains in the Price model should follow a regular pattern which depends on the height of nodes in each antichain. Numerical confirmation of these patterns are given in Appendix B.6 of the Supplementary Information.

Discussion

The numerical results for the leading behaviour of the path length scales are striking. Within the margin of numerical error, our results in Fig. 5 show that the length of the reverse greedy path scales asymptotically as $m\bar{p}\ln(t)$ for a wide range of parameter values. This is consistent with both the simple argument and the detailed analytical calculation presented in the Analytic Results (see also Appendix A of the Supplementary Information). The analytical approach also shows that for long times, the distribution of lengths of reverse greedy paths in the Price model is Poisson distributed with mean equal to $m\bar{p}\ln(t)$ (see Appendix A).

The reverse greedy path length is a lower bound on the length of the longest path so it is no surprise that the longest path length also scales as $\ln(t)$ with a coefficient, $a_{\mathrm{max}}$, which is larger than the corresponding scaling factor for the reverse greedy path length, $a_{\mathrm{gr}}$. Interestingly this coefficient of the $\ln(t)$ term, $a_{\mathrm{max}}$, for the longest path shows some additional weak dependence on the parameters beyond the $m\bar{p}$ found for the reverse greedy path, as both Fig. 5 and Fig. 6 clearly show.

The Price model is not in itself a very realistic model for any particular context. For instance, a true citation network often shows many other features such as a preference to cite recent papers, for example see [13, 14, 15, 16, 7, 17, 18]. The choice of a simple linear form for the attachment probability (1) appears to be part of this simplification, a form linear in degree motivated by the need for mathematical simplicity. At first sight, this form seems unrealistic since it requires authors of papers to have global information about the citation network because of the normalisation factors. No author can know exactly how many citations a paper has let alone the total number of citations in the network. However, this form emerges naturally in many situations as the result of doing local searches on the network, see [19, 20, 21, 22, 23, 24, 7, 18] and references therein. In more realistic models, the cumulative advantage, the $p$ term in $\Pi$ emerges from doing a local search back through the current citation network, while the random attachment, the $\bar{p}$ term in $\Pi$, represents a simple model of other possible processes. So like all good models, the emphasis in the Price model on the linear form for $\Pi$ in terms of degree does capture an important and realistic feature of many real situations. This linear form of the attachment probability $\Pi$ is also the critical feature in the analysis of undirected versions of this model, such as the Barabási-Albert model [25] where the cumulative advantage aspect is known as preferential attachment and the original example worked with $p=1/2$ in our notation.

However, the Price model, simple as it is, also emphasises another critical aspect of a citation network, and that is the inherent arrow of time in this context. Citations (almost) always point backwards in time333Typical data sets [4, 26] suggest that less than 1% of citations are to documents which are labelled as being published later than the citing document. and the networks created in the Price model are realistic in this way, they always produce directed acyclic graphs. This acyclic property is lost when the edge direction is ignored, as in the Barabási-Albert implementations of this model. Since many analyses work in the undirected version, they have missed this key feature of the Price model and of real-world citation networks.

For instance, the length of the shortest path between two nodes is a natural measure for undirected networks since in some circumstances it can be related to the geodesic of networks embedded in Euclidean space, for example see [9, 10]. For an undirected version of the Price model, the LCD model of Bollobás and Riordan [27, 28] (a more precisely defined version of the Barabási-Albert model [25]), it is known that the diameter, the largest length of any shortest path between two nodes, scales as $\ln(t)/(\ln\ln(t))$ if $m>1$ while the diameter scales as $\ln(t)$ for the special case of $m=1$ [27] (see also theorem 18 of [28]). For the case $m>1$, Bollobás points out that while in any random graph we expect to see the small-world effect [29] and a $\ln(N)$ scaling of lengths (for $N$ nodes in the network), for the undirected version of this model “one might expect the diameter to be even smaller: the unbalanced degree distribution pushes up the number of small paths, and thus, perhaps, pushes the diameter down” (see Bollobás [28] section 13, p25). That is, the unusual slow scaling of the shortest path distance scale in this undirected version of this model is due to the effect of the very high degree nodes created because of the cumulative advantage (preferential attachment) process.

However, the situation is completely different when we take account of the direction of edges in this model. First, the link between shortest path lengths and geodesics in Euclidean space used in [9, 10] is lost. The natural order of nodes in a DAG, the arrow of time, means we should compare the path lengths network against geodesic lengths for network models embedded in Minkowski space, and indeed there is a proven relationship between these two [8, 30, 31]. Following on from this, when using the longest path in the directed form of these models, our analysis has shown that the longest path is likely to be created by edges created from random attachment not those formed using the cumulative advantage mechanism, the opposite of what is suggested for the shortest path in the undirected form of these models. Thus the fat-tailed nature of the degree distribution in the Price model (or directed versions of the Barabási-Albert/LCD models) is not a factor for the longest path and so, using Bollobás’ insight [28], we should expect the longest path to scale simply as $\ln(t)$, and not something slower than that. That is, indeed what we have shown in our work here.

Looking more widely, we note that Bollobás [28] (p10) suggested that “For these models the orientation is not very interesting”. Our conclusion is the opposite. Namely that for any directed network in which vertices are added sequentially, the arrow-of-time inherent in these growing network models is both physically relevant and this vertex order produces new and distinctive features. Our analytical and numerical analysis of the longest path length is just one illustration of what is possible.

Acknowledgements

V.V. acknowledges support from EPSRC, grant number EP-R512540-1.

Author contributions statement

T.S.E. derived the analytical solution. L.C. and V.V. conducted the data analysis. T.S.E. and V.V. wrote the manuscript. All authors developed computational framework, analysed and interpreted the results and reviewed the manuscript.

Additional information

The authors declare no competing interests.

Supplementary Information.

A Full Analytic Calculation

In this appendix we will work with a DAG $\mathcal{G}$ with vertex set $\mathcal{V}$ and edge set $\mathcal{E}$. The vertices will be labelled by sequential integers $t$ which take values from $1$ to the number of vertices $N=|\mathcal{V}|$. The directed edges are defined to run from the lower value vertex to the higher value vertex, so if $(s,t)\in\mathcal{E}$ then $s<t$.

A.1 Definition of reverse greedy paths

Every vertex $t$ in a DAG $\mathcal{G}$ has two reverse greedy paths: the forward reverse greedy path to nodes of times later than $t$, the other the reverse greedy path coming from nodes of earlier times. Here we will define the reverse greedy path as that is what is used in this work.

The reverse greedy path arriving at a vertex $t$ is defined by following the edge arriving at that vertex $t$ which links back to the most recent predecessor vertex $s$. More formally, suppose $\mathcal{N}^{(-)}(t)$ is the set of predecessors of a vertex $t$, the source vertices of edges whose target vertex is $t$:

[TABLE]

where $\mathcal{E}$ is the edge set of the DAG. Then the reverse greedy path to vertex $t$ is a sequence of vertices denoted as $P^{(-)}(t)$ which is defined to be

[TABLE]

Note that the length of the reverse greedy path to $t$ is denoted $\ell(t)$ so $\ell(t)=|P(t)|-1$. The reverse greedy path will always terminate at a global source node, the only node with zero in-degree.

A.2 Price Model Definition

The Price model [3] (see also sec.14.1 of [5]) is a growing network model giving us a sequence of graphs $\{\mathcal{G}(t)\}$ where these are indexed by an integer $t$, the ‘time’. Each graph has a vertex set $\mathcal{V}(t)$ and an edge set $\mathcal{E}(t)$. To create the next graph $\mathcal{G}(t+1)$ in the Price model sequence, we add one new vertex at time $(t+1)$ to $\mathcal{V}(t)$. We will overload our notation so that the positive integer $t$ is both the ‘birth’ time of a vertex and it is also used to represent that vertex. So formally $\mathcal{V}(t+1)=\mathcal{V}(t)\cup\{(t+1)\}$. We choose to index our nodes from one upwards so that at time $t$ the total number of nodes in the graph $\mathcal{G}(t)$, denoted as $N(t)$, is $N(t)=|\mathcal{V}(t)|=t$.

A new node $(t+1)$ is connected by $m$ new edges to existing vertices $s$ chosen with probability $\Pi(t,s)$. The number of edges in the edge set $\mathcal{E}(t)$ of the DAG $\mathcal{G}(t)$, edges between nodes created at time $t$ or earlier, is $E(t)=E_{0}+mt=|\mathcal{E}(t)|$ where $E_{0}$ is some constant. Note this is the total number of edges after we have finished adding all the incoming edges to vertex, $t$.

Finally the node created at time $s$ in the graph $\mathcal{G}(t)$ has out-degree $k^{(\mathrm{out})}(t,s)$.

The connection of edges to new node $(t+1)$ as encoded by $\Pi(t,s)$ is made in one of two ways. With probability $p$ the node $(t+1)$ is connected to an existing vertex $s$ chosen with cumulative advantage $k^{(\mathrm{out})}(t,s)/E(t)$. Otherwise, so with probability $\bar{p}=(1-p)$, we choose an existing vertex $s$ uniformly at random from the set of existing vertices, i.e. with probability $1/N(t)$. That is the probability of connecting new vertex $(t+1)$ to existing vertex $s$ is444In Newman’s notation (equation (14.1), sec.14.1 of [5]) $\Pi(t,s)=(k^{(\mathrm{out})}(t,s)+a)/(t(c+a))$ with $a=m\bar{p}/p$ and $c=m$ (with $c+a=m/p$).

[TABLE]

The usual form, the first equation, leaves us, however, with a problem for $\Pi(t=1,s=1)$ when looking at the attachment to the second vertex, $t=2$, if we think of the initial vertex at $t=1$ as having no incoming edges. The solution is to insist that $\Pi(t=1,s=1)=1$. We will see that this fixes the cumulative probability $\Pi_{\leq}$ to have a consistent value which is in fact all we need for this calculation.

For simplicity we will almost always assume that $E(t)=mt$ and $N(t)=t$. This is impossible to satisfy at early times as at the end of time step $1$ we have one node which cannot have the $m$ edges required since self-loops are not allowed for a DAG. In turn, this early time issue is a clear signal that there are some initial graph effects in the Price model. In practice, we can deal with this issue in various ways. If multiple edges between pairs of vertices are allowed then we can enforce this $E$ and $N$ constraint from step $t=2$ onwards (set $2m$ edges from node $1$ to node $2$) but at the cost of having a very skewed initial degree distribution. Without multiple edges, then the earliest point where we can have our $E(t)=mt$ and $N(t)=t$ is at $t=2m+1$ if at that time we have a complete directed graph in which the first $2m+1$ nodes have edges leaving them to connect to every later node, i.e. where the edge set is $\mathcal{E}(t=2m+1)=\{(i,j)|1\leq i<j\leq 2m+1\}$. Unless otherwise stated, we used this complete graph graph at time $t=2m+1$ as the starting point for our numerical simulations.

We will need the cumulative probability function which is the probability that you attach all $m$ edges from the node created at time $(t+1)$ to a node created at time $s$ or earlier

[TABLE]

It is clear that we must have

[TABLE]

as all existing vertices are then of time $t$ or less. This shows that if $\Pi_{\leq}(t=1,s=1)=1$ then we need $\Pi_{\leq}(t=1,s=1)=\Pi(t=1,s=1)=1$ for consistency.

To find this cumulative advantage probability $\Pi_{\leq}(t,s)$ we can work in terms of an effective node, a single ‘super’ node, which represents the $s$ nodes created from the initial time to time $s$. If we coarse grain these nodes from $t=1$ to $s$ into one super node, then we should attach all their edges to this one super node (we can picture this as having multiple self-loops) so that algebraically the degree of this super node, $k^{(\mathrm{out})}_{\mathrm{eff}}(t,s)$, is the sum of all the degrees of the individual nodes created at or before $s$,

[TABLE]

Now we can use this to write the cumulative probability for attachment as

[TABLE]

The master equation for this effective node is therefore

[TABLE]

We can rewrite this in terms of a function $y(t,s)$ where

[TABLE]

The master equation then becomes

[TABLE]

using $y(s,s)$ value in (A11). This gives

[TABLE]

We note that this is separable with

[TABLE]

For large $t$ and $s$ we have that

[TABLE]

A.3 Master Equation for Reverse Greedy Path Length $\ell$

We are interested in looking at the reverse greedy path in our growing network model. For this we need to know $\Pi_{\mathrm{max}}(t,s)$, the probability that of the $m$ predecessor nodes connected to a new node at $(t+1)$, the oldest of them is $s=\max(\mathcal{N}^{(-)}(t+1))$. The probability that we connect a new node $(t+1)$ to nodes of age $s$ or less is simply $(\Pi_{\leq}(t,s))^{m}$ where $\Pi_{\leq}(t,s)$ is the cumulative probability of attachment of (A5). This gives us that

[TABLE]

Note the case $s=1$ is covered as we defined $\Pi_{\leq}(t,0)=0$. We may check that

[TABLE]

where we use the definition $\Pi_{\leq}(t,0)=0$ and (A6).

Let the probability that the length of the reverse greedy path, $\ell$, from new node $(t+1)$ to the initial node at $t=1$, be $P(\ell,t)$. This will satisfy the equation

[TABLE]

We can rewrite this master equation for $P(\ell,t)$ in terms of the generating function $G(z,t)$ defined as

[TABLE]

where

[TABLE]

This last follows because the only source node is the initial node at $t=1$ with reverse greedy path length zero, $\ell(t=1)=0$. Ultimately, we want to look at the average length of the reverse greedy path, $\langle\ell(t)\rangle$, which is obtained as

[TABLE]

In terms of the generating function, the master equation (A19) becomes

[TABLE]

Peeling off the top term of the sum in (A24) gives us

[TABLE]

For now, we will simply assume that cumulative probability $\Pi_{\leq}(t,s)$ is separable, i.e. $\Pi_{\leq}(t,s)=h(t)j(s)$ for some functions $h(t)$ and $j(s)$. We can use the fact that $\Pi_{\leq}(t,s=t)=1$ from (A6) to see that $h(t)=1/j(t)$ so that $\Pi_{\leq}(t,s)={j(s)}/{j(t)}$. We have the precise form for $j(s)$ in the Price model in (A15). That gives us that $\Pi_{\mathrm{max}}(t,s)$ is also separable and this can be written as

[TABLE]

Then

[TABLE]

Using (A24) leaves us with

[TABLE]

which has the formal solution

[TABLE]

given $G(z,1)=1$. Now we can check $z=1$ value which is clearly $G(z=1,t)=1$ as required.

We can obtain the precise form for the Price model by using the form for $H(s)$ in (A27) given $j(s)$ in the Price model from (A15). This gives us that

[TABLE]

The solution for the generating function then becomes

[TABLE]

In principle this is a ratio of $m$-th order polynomials in $s$ so you can write this product as the product of $m$ ratios of Gamma functions, one term per root of the polynomial.

Average Path Length

We can find $\ell(t)$, the average length of the reverse greedy path in the Price model, from this general solution (A30) using (A23) and $G(z=1,t)$ from (A21). This gives us that

[TABLE]

so then

[TABLE]

and finally

[TABLE]

For the Price model we have $H(s)$ from (A33) and so we find that

[TABLE]

Now we can expand the product using the binomial expansion to find that

[TABLE]

For each value of $n$ we can write the $s^{-n}$ series as the difference of two Hurwitz zeta functions $\zeta(m,t)$

[TABLE]

The Hurwitz zeta functions are finite for $t>0$ and $n>1$ so we see that the terms $n\geq 2$ only contribute a constant plus terms from $\zeta(n,t)$ which fall off as $t^{1-n}$ or faster.

The leading term in the large time limit of $\langle\ell(t)\rangle$ comes only from the $n=1$ term in (A42) and this is

[TABLE]

where $\psi(t)$ is the digamma function.

The constant in an asymptotic expansion in $t$ of $\langle\ell(t)\rangle$ picks up further contributions from the $n\geq 2$ terms so we have that

[TABLE]

where it is implicit that there is no contribution from the term with the sum for the case of $m=1$. Here $\zeta(n)$ is the Riemann zeta-function.

Variance

We can also use the generating function solution (LABEL:aGsol) to find the variance $\sigma^{2}(t)$ of the reverse greedy path length. We have that

[TABLE]

We have that

[TABLE]

so that

[TABLE]

With $\sigma^{2}(t)=\langle(\ell(t))^{2}\rangle-(\langle\ell(t)\rangle)^{2}$ we have that

[TABLE]

In the long-time limit we have that

[TABLE]

Note that these correction are all finite so we have basically shown that the variance grows with the mean i.e. Poisson-like behaviour.

Long-time Generating Function Result

From (LABEL:aGsol) the logarithm of the generating function $G(z,t)$

[TABLE]

Expanding around $z=1$, or with $\bar{p}$ small, we have that

[TABLE]

We recognise the first term as $\langle\ell(t)\rangle$ from (A41) so that

[TABLE]

That is the leading term comes from the only term with a sum over $1/s$ which diverges as $\ln(t)$ and is completely captured by the $\langle\ell(t)\rangle$ contribution. The higher order terms in an expansion around $z=1$ are all finite which may be expressed as polynomials in $\bar{p}$ with terms from $\bar{p}^{q}$ to $\bar{p}^{qm}$ for the coefficient of the $(1-z)^{q}$ term.

From this we wee that the leading order term in $G(z,t)$ in the long time limit

[TABLE]

which is the generating function for a Poisson distribution. So we deduce that in the long-time limit, the distribution of lengths of the reverse greedy path in the Price model is a Poisson distribution with mean equal to $m\bar{p}\ln(t)$.

[TABLE]

A.4 Reverse Greedy Path for $m=1$

The Price model for $m=1$ is a special case in this model as then the graphs produced are directed trees. Various aspects of this $m=1$ case of the Price model and its variants can be studied analytically, for example [32, 35]. Here we just note that many of the equations in our calculations become simple first order polynomials, such as for $\Pi_{\mathrm{max}}(t,s)$ in (A17), which allows for a simple direct solution in the $m=1$ that case. The solutions for the reverse greedy path length also takes on a much simpler form.

For $m=1$ we have from (LABEL:aGsol) that the generating function is just

[TABLE]

Similarly, the expression for the reverse greedy path length is also simple:

[TABLE]

where $\gamma\approx 0.577$ is the Euler-Mascheroni constant.

A.5 Vertex Partitions

A partition $\mathcal{P}$ of the vertex set $\mathcal{V}$ is a set of non-overlapping non-empty subsets which contain each and every vertex once and only once. That is $\mathcal{P}=\{\mathcal{P}(i)\}$ where the blocks of the partition, $\mathcal{P}(i)$, are such that $\mathcal{P}(i)\subset\mathcal{V}$, $\mathcal{P}(i)\neq\emptyset$, $\mathcal{P}(i)\cap\mathcal{P}(j)=\emptyset$ unless $i=j$, and $\cup_{i}\mathcal{P}(i)=\mathcal{V}$. The definition of a unique integer length scale associated with each node in any one instance of the model gives a natural partition of the set of vertices. reverse greedy path to each node $t$ from a source node, $\ell(t)$.

Of particular interest here is a partition in terms of the longest path $L(t)$ to a given node $t$. In any DAG, the height of a node is the length of the longest path to a node from any global source node, any node with zero in-degree. Thus in the Price model, the height of a node is simply the longest path length from the initial node to any node, that is our $L(t)$ is the height. A natural partition of the set of vertices in a DAG is where each block, $\mathcal{P}_{\mathrm{height}}(h)$, contains all the vertices of height $h$,

[TABLE]

This partition by height has the special property that no two vertices in any one block are connected by any path. In a DAG nodes connected by a path cannot be of the same height. A set of disconnected vertices is known as an antichain so the blocks of the height partition are all antichains. In some sense, the vertices in an antichain can be considered to be of equivalent ages. There is no relationship between them which says that any one vertex in the antichain need to come before or after another. The question here is can we estimate some of the properties of these height antichains? To do that we need to define a slightly different partition.

For the Price model, we can also define another type of partition in terms of the average path length $\langle\ell(t)\rangle$ associated with a node created at time $t$ averaged over all possible instances of the Price model. We will discuss this partition in terms of the reverse greedy path length $\ell$ where we have analytical results but the construction is identical for the longest path length $L$. Since our average path lengths are monotonically increasing functions of time $t$, we can define the blocks $\widehat{\mathcal{P}}(\hat{\ell})$ of one partition in terms of an integer value $\hat{\ell}\in\mathbb{Z}$ as follows. All nodes in one block of the partition have an expectation value for their reverse greedy path length which rounds to the integer $\hat{\ell}$,

[TABLE]

The idea here is that we can estimate the properties of the blocks in this partition, $\widehat{\mathcal{P}}(\hat{\ell})$, by using the leading behaviour for the lengths scales, namely they grow logarithmically with $t$. We would then expect that the partitions of a single instance based on the measured values of an integer valued path length would show similar features.

The leading behaviour for the mean reverse greedy path length of a node $t$ can be expressed as

[TABLE]

where $\bar{\ell}=\langle\ell(t)\rangle$ while $\mu_{\mathrm{gr}}$ and $\tilde{t}$ are some constants. Our leading order results give $\mu_{\mathrm{gr}}=\exp(1/(m\bar{p}))$ for the reverse greedy path, while we can estimate $\mu=\exp(1/a)$ and $\tilde{t}=\exp(-b/a)$ using the numerical results for the coefficients $a$ and $b$ obtained from the fit to $a\ln(t)+b$.

From this behaviour we can estimate various properties of the blocks $\widehat{\mathcal{P}}(\hat{\ell})$. It is useful to express these results in terms of the time scale $\hat{t}(\hat{\ell})$ (not necessarily an integer) for a node whose path length is expected to be the integer $\hat{\ell}$, namely $\hat{t}(\hat{\ell})=\tilde{t}(\mu_{\mathrm{gr}})^{\hat{\ell}}$ for $\hat{\ell}\in\mathbb{Z}^{+}$. Other properties of the path length partition $\widehat{\mathcal{P}}(\hat{\ell})$ can be estimated in terms of $\hat{t}$: the average time of nodes in the partition $\bar{t}(\hat{\ell})$, the standard deviation of the times of nodes in the partition $\sigma_{t}(\hat{\ell})$ and the number of nodes in the block $|\widehat{\mathcal{P}}|(\hat{\ell})$. We find that

[TABLE]

The results above are for the partition into blocks $\widehat{\mathcal{P}}(\hat{\ell})$ of (A64) which are defined in terms of averages over an ensemble. More interesting is to look at the partition into height antichains, $\mathcal{P}_{\mathrm{height}}(h)$ of (A63) since this based on the topology of a single network and so much more relevant to studies of real data sets. Our suggestion is that in the Price model all the length scales behave in the same way, following our description in (A65). So we should expect that all partitions in the Price model based on such length scales should show the scaling behaviour suggested in (A68). For example that suggests that for the height antichains, we should compare the ratios of a particular property for adjacent blocks, for instance $\bar{t}(h+1)/\bar{t}(h)$, since this analysis suggests such ratios will be constant.

B Numerical Implementation

B.1 Price Model Algorithm

To reduce the number of random numbers drawn, we use the following algorithm. First we observe that the probability $\Pi(t,s)$ of (A3) used to choose existing nodes $s$ to attach to a new node $(t+1)$ is proportional to $\Pi(t,s)\propto k^{(\mathrm{out})}(t,s)+(m\bar{p}/p)$. If we limit ourselves to parameter values where $\alpha=(m\bar{p}/p)$ is a non-negative integer, then we can choose existing vertices with probability $\Pi(t,s)$ by drawing uniformly at random from a list, attachment_list, if it is properly formed. We define our algorithm as follows.

Set initial graph G with $t$ nodes. 2. 2.

Initialise attachment_list to match initial graph.

• •

For every edge in the initial graph, we add the source end of the edge to attachment_list. That is is $(s,t)\in\mathcal{E}(t)$ where $s<t$ then we add $s$ only to attachment_list.

• •

Each node in the initial graph is added $\alpha=(m\bar{p}/p)$ times to attachment_list. 3. 3.

Increment $t$. 4. 4.

Draw $m$ times uniformly at random from attachment_list, to give the list (sequence) $S(t)=[s_{1},s_{2},\ldots,s_{m}]$.

• •

If drawn with replacement, this is simple and fast and it matches the assumptions implicit in the algebraic calculations. This allows multiple edges between nodes.

• •

If drawn without replacement we avoid multiple edges but this is liable to be slower. However this constraint is not enforced in the analytic equations. 5. 5.

Append the $m$ nodes in the list $S(t)$ to attachment_list.

This captures the cumulative advantage process. 6. 6.

Append $\alpha=(m\bar{p}/p)$ copies of $t$ to attachment_list.

This encodes the uniform random process. 7. 7.

If adding more nodes, return to step 3.

We allowed multiple edges in our simulations and in this approach one random number was needed for every edge added, a total of $mN$ random numbers if the final network had $N$ nodes. For the attachment list $(m+\alpha)N$ memory locations are needed. Note that in our work here we only created the attachment_list along with lists for the length of paths to each vertex. There was no need for use to create the full network structure. However, should the full network structure be required, it can be deduced later from attachment_list as defined here.

Note that the $p=0$ case is treated separately since the attachment list is not needed. In that case we draw uniformally at random from the set of nodes, with replacement in our case since we allow multiple edges.

In our work we choose to allow multiple edges in the networks created in our simulations. This will have relatively little effect on results as the model is inherently a sparse graph model. As the number of nodes increases, as $m$ is fixed, the chance of choosing the same existing vertex twice when connecting it to a new node decreases.

As noted elsewhere, our initial graph was a DAG of $2m+1$ nodes with an edge between every node, that is the initial edge set was $\mathcal{E}(2m+1)=\{(s,t)|1\leq s<t\leq 2m+1\}$. This has $E(t)=mt$ for $t\geq 2m+1$. However, we did investigate the effect of the initial graph and this is discussed in Section B.4.

An alternative algorithm for the Price model, which we did not use here, would be to carry two separate lists of existing vertices, one used for the cumulative advantage process, and one for the selecting uniformly at random from the set of vertices. This is equally easy to code. The disadvantage is that we would now need to draw an additional random number for every edge added, a total now of $2mN$, as this extra random number has to be compared against the parameter $p$ in order to decide which list to draw from. The extra random numbers would slow the process but the gain would be flexibility as there would be no restriction on the parameters to ensure that $\alpha=(m\bar{p}/p)$ was integer. This alternative method would also only require $(m+1)N$ memory locations for the lists.

B.2 Numerical Reverse Greedy Path Algorithm

Numerically, the reverse greedy path may be found as follows

Set initial graph G with $t$ nodes. 2. 2.

Set the reverse greedy path length for the initial $t$ nodes, greedy_length[s] for $s=1,2,\ldots,t$. 3. 3.

Set the nearest neighbour along this reverse greedy path length for the initial $t$ nodes, greedy_neighbour[s] for $s=1,2,\ldots,t$. 4. 4.

Increment $t$. 5. 5.

Add a new node index $t$ to G. 6. 6.

Update graph G by adding $m$ new edges, from vertices $\{s\}$ to new node $t$, which then define the past neighbour set of vertex $t$, i.e. $\mathcal{N}^{(-)}=\{s\}$. 7. 7.

Find the past neighbour with the largest time $s_{\mathrm{gr}}=\max(\mathcal{N}^{(-)}(t))$. 8. 8.

Record this closest past neighbour as the last step on the reverse greedy path to $t$,

i.e. set greedy_neighbour[t] $=s_{\mathrm{gr}}$. 9. 9.

Set the length of the reverse greedy path of new node $t$ to be one more than the length of the reverse greedy path to this nearest past neighbour $s_{\mathrm{gr}}$, so set $\ell(t)=\ell(s_{\mathrm{gr}})+1$. Numerically we just need to set

 greedy_length[t]=

greedy_length[greedy_neighbour[t]]+1. 10. 10.

If adding more nodes, return to step 4.

The length of the reverse greedy path from the first node $t=1$ to a given node $t$ is now stored for each node. The vertices on the unique reverse greedy path can be found by iterating back though the greedy_neighbour property of the nodes. If the specific path used is not required we need not record this information.

We also note that the length of the longest path, the longest path from the initial node $t=1$ to any node $t$, may be tracked numerically in a similar way. Such a longest path is always exists but is not unique so we could only ever record one example longest path using this approach.

B.3 Fitting

For each set of parameters $p,m,N$, we ran the model $R=100$ times. We collected the length of the reverse greedy path and of the longest path for each node created at time $t$. In this case, each data point $y_{tr}$ is the relevant path length, $\ell(t)$ or $L(t)$, measured for node $t$ on run $r$.

We make several assumptions about our data. First, we assume that the measurements we make on each individual run are independent. This assumption is not completely true in our data.

Secondly, we assume that the set of measurements at the same point $t$ across many runs have a distribution that is normally distributed. As Fig. B7 shows, this is approximately true.

To evaluate the quality of our fit, we looking for parameter values which minimised the value of the chi-squared function (for instance see [12])

[TABLE]

where $y_{tr}$ is a path length of the node created at time $t$ as measured in run $r$. The fitting function used in the main text is $f_{t}=a\ln(t)+b$ but we also tried $f_{t}=a\ln(t)+b+(c/t)$ where $a,b,c,$ are constants to be found by fitting.

For numerical convenience (minimising memory requirements), for each of the path length we kept track of two values for each node occurring at $t$, the total sum of values $T_{t}$ and the sum of squares $S_{t}$, as given by

[TABLE]

We do not expect the fit to be valid for small system sizes. For instance, in most of our work the initial graph is a complete DAG (see Appendix B.4 for a discussion of initial graph effects) in within this initial the path lengths scale as $y_{t}=L(t)=\ell(t)=t$. To deal with this we fit our data from time $t_{0}$, an additional parameter in our fit, to the largest time which is equal to the number of nodes $N$. So, in terms of the totals and the squares of totals, our chi-squared function becomes

[TABLE]

Under our assumptions of normal fluctuations and independence as noted above, the probability of obtaining a particular $\chi^{2}$ value is given by the integral of the complementary cumulative distribution function of the $\chi^{2}$ distribution for the given degrees of freedom. The number of degrees of freedom in the model is equal to $(N-t_{0}+1).R$. We varied the $t_{0}$ parameter between $100$ and $10,000$ but found that this cutoff parameter $t_{0}$ had no significant influence on the resulting fits. This is to be expected, as even the largest $t_{0}$ value considered, $10,000$, constitutes a mere $0.1\%$ of the data we were using. So, in our work we used a fixed value of $t_{0}=1,000$.

The errors in the parameters of the linear fit to $f_{t}=a\ln(t)+b$ come from the square root of the diagonal elements of the covariance matrix. We found that there was no consistent pattern in the next-to-leading order terms, that is in $b$, so these are not shown. As noted elsewhere, this constant factor $b$ in our fit is likely to be effected by the form of the initial graph since this can add a significant constant to the lengths scales and the initial graph will effect lengths scales later in different ways depending on the parameters chosen. We also tried looking at higher order terms, using non-linear fits to $f(t)=a\ln(t)+b+(c/t)$, but again saw no improvements in our results. For the results of the comparison between the linear and non-linear fits see Table B1.

B.4 Initial graph effect

Throughout this work, we started our simulations from an initial graph which was the directed complete graph composed of $(2m+1)$ nodes. That is $\mathcal{G}_{\mathrm{compl}}(t)$ where $\mathcal{V}_{\mathrm{compl}}(t)=\{1,2\ldots,t\}$ and $\mathcal{E}_{\mathrm{compl}}(t)=\{(s,r)|1\leq s<r\leq t\}$. This is a dense, transitively complete graph but at each step in our simulation we always have exactly $m$ more edges than nodes which is an assumption in the analytical work. A possible drawback is that this initial graph has long path lengths, $\ell(t)=L(t)=(t-1)$ for nodes in this graph. As both our longest and reverse greedy paths will start with paths in the initial graph, these long initial graph path lengths will produce a significant addition to those we measure. For instance, for $m=5$ this initial graph contributes up to 10 to any path we measure while the typical length scales we measure at late times are less than 100 in general.

We can, however, choose alternative initial graphs. To match the standard Price model and our analysis, we confine ourselves to initial graphs which are weakly connected and which have a single global source node555In fact, if we have several global source nodes we could add just one more “global parent node” with an edge from this extra node to all of the original global source nodes. This adds one to all all the lengths we measure but creates a DAG with a single node of $k^{(\mathrm{in})}=0$. where $k^{(\mathrm{in})}=0$. A transitively reduced [4] version of our chosen initial graph is a single chain of nodes, so $\mathcal{G}_{\mathrm{chain}}(t)$ where $\mathcal{V}_{\mathrm{chain}}(t)=\{1,2\ldots,t\}$ and $\mathcal{E}_{\mathrm{chain}}(t)=\{(r,r+1)|1\leq r<t\}$, would also be a valid choice. This type of initial graph is as sparse as possible, all nodes have $k^{(\mathrm{in})}$ and $k^{(\mathrm{out})}$ equal to zero or one, but this initial graph has the same long path lengths as the complete graph. The reason the chain and the complete graph lead to different behaviour in the shift in the average path length at long times is because of their different initial degree distributions which alters the pattern of attachment at early times. That in turn alters the likelihood that the paths we measure join the initial graph at a particular node. In particular, we expect that the paths we measure in models starting with the complete graph are more likely to contain a smaller fraction of the initial graph as the earliest nodes will have higher out degrees but add shorter path lengths to our measurements. By way of contrast, in the chain, most nodes have the same initial degree, they are likely to have similar degree over time so the paths we measure are more likely to leave the initial graph at a later node so giving a longer contribution to the path length coming from this initial chain graph.

Alternatively, we could also use an initial graph, in which the single global source node points to the remaining nodes, which in turn are not pairwise connected. This is a “directed star graph” $\mathcal{G}_{\mathrm{star}}(t)$ where $\mathcal{V}_{\mathrm{star}}(t)=\{1,2\ldots,t\}$ and $\mathcal{E}_{\mathrm{star}}(t)=\{(1,r)|1<r\leq t\}$. Now the degree distribution is highly skewed, with the central node going to dominate attachment through the the cumulative advantage mechanism. However, connecting to one of the other nodes only adds one to the paths we measure so we expect our models using this initial star graph will have the smallest paths lengths of the three cases considered.

As Fig. B8 shows, the results obtained using either of the three studied variants of the initial graph are distinguishable. The asymptotic scaling of the path lengths as $a\ln{t}$ remains the same with the same value of the parameter $a$. However the constant contribution, $b$ of (3.1), seems to depend on the choice of the initial graph. As Fig. B8 and Table B2 show, the $b$ coefficient behaves exactly as suggested above: the star graph gives the shortest path lengths while the chain gives the longest.

B.5 Next-to-leading Order Coefficient

The next-to-leading order coefficient, $b$ of (3.1) did not reveal as clear trends as the coefficient $a$ of the same equation. There does not seem to be an obvious relation between $b_{\textrm{max}}$ and $b_{\textrm{gr}}$, as Fig. B9 shows.

B.6 Height Antichain Properties

The height of a node in a DAG is the length of the longest path to a node from any global source node, any node with zero in-degree. Thus in the Price model, the height of a node is simply the longest path length from the initial node to any node, our $L$. We can also define antichains, sets of nodes which are not connected by any path, for example see [34]. A simple example in a DAG is a height antichains, a set of nodes which are all at the same height since nodes connected by a path cannot be of the same height. Based on the $\ln(t)$ scaling of the reverse greedy path length, we can speculate that the average height $L$ of a node $t$ will scale as $t=\tilde{t}_{\mathrm{max}}(\mu)^{L}$. Then, as discussed in Section A.5, we find that various properties of the antichains should scale in a simple way with height: the average time of nodes in the partition $\bar{t}(L)$ as in (A66), the standard deviation of the times of nodes in the partition $\sigma_{t}(L)$ as in (A67), and the number of nodes in the block $|\widehat{\mathcal{P}}|(L)$ as in (A68). Results for ratio of these quantities for adjacent blocks are shown in Fig. B10.