Are RNA networks scale-free?

A network is scale-free if its connectivity density function is proportional to a power-law distribution. It has been suggested that scale-free networks may provide an explanation for the robustness observed in certain physical and biological phenomena, since the presence of a few highly connected hub nodes and a large number of small-degree nodes may provide alternate paths between any two nodes on average-such robustness has been suggested in studies of metabolic networks, gene interaction networks and protein folding. A theoretical justification for why many networks appear to be scale-free has been provided by Barabási and Albert, who argue that expanding networks, in which new nodes are preferentially attached to highly connected nodes, tend to be scale-free. In this paper, we provide the first efficient algorithm to compute the connectivity density function for the ensemble of all homopolymer secondary structures of a user-specified length-a highly nontrivial result, since the exponential size of such networks precludes their enumeration. Since existent power-law fitting software, such as powerlaw, cannot be used to determine a power-law fit for our exponentially large RNA connectivity data, we also implement efficient code to compute the maximum likelihood estimate for the power-law scaling factor and associated Kolmogorov-Smirnov p value. Hypothesis tests strongly indicate that homopolymer RNA secondary structure networks are not scale-free; moreover, this appears to be the case for real (non-homopolymer) RNA networks.

Introduction

The connectivity (or degree) of a node v in a network (or undirected graph) is the number of nodes (or neighbors) of s, connected to v by an edge. A network is said to be scale-free if its connectivity function N(k), which represents the number of nodes having degree k, satisfies the property that , the unique solution of which is a power-law distribution, which by definition satisfies for some scaling factor (Newman 2006). Scale-free networks contain a few nodes of high degree and a large number of nodes of small degree, hence may provide a reasonable model to explain the robustness1 often manifested in biological networks—such robustness or resilience must, of course, be present for life to exist.

Barabasi and Albert (1999) analyzed the emergence of scaling in random networks, and showed that two properties, previously not considered in graph theory, were responsible for the power-law scaling observed in real networks: (1) networks are not static, but grow over time, (2) during network growth, a highly connected node tends to acquire even more connections—the latter concept is known as preferential attachment. In Barabasi and Albert (1999), it was argued that preferential attachment of new nodes implies that the degree N(k) with which a node in the network interacts with k other nodes decays as a power-law, following , for . This argument provides a plausible explanation for why diverse biological and physical networks appear to be scale-free. Indeed, various publications have suggested that the the following biological networks are scale-free: protein–protein interaction networks (Ito et al. 2000; Schwikowski et al. 2000), metabolic networks (Ma and Zeng 2003), gene interaction networks (Tong et al. 2004), yeast co-expression networks (Van Noort et al. 2004), and protein folding networks (Bowman and Pande 2010).

How scale-free are biological networks?

The validity of a power-law fit for previously studied biological networks was first called into question in Khanin and Wit (2006), where 10 published data sets of biological interaction networks were shown not to be fit by a power-law distribution, despite published claims to the contrary. Estimating an optimal power-law scaling factor by maximum likelihood and using goodness-of-fit tests, it was shown in Khanin and Wit (2006) that not a single one of the 10 interaction networks had a nonzero probability of being drawn from a power-law distribution; nevertheless, some of the interaction networks could be fit by a truncated power-law distribution. The data analyzed by the authors included data from protein–protein interaction networks (Ito et al. 2000; Schwikowski et al. 2000), gene interaction networks determined by synthetic lethal interactions (Tong et al. 2004), metabolic interaction networks (Ma and Zeng 2003), etc.

In Clauset et al. (2009), 24 real-world data sets were analyzed from a variety of disciplines, each of which had been conjectured to follow a power-law distribution. Estimating an optimal power-law scaling factor by maximum likelihood and using goodness-of-fit tests based on likelihood ratios and on the Kolmogorov–Smirnov statistic for non-normal data, it was shown in Clauset et al. (2009) that some of the conjectured power-law distributions were consistent with claims in the literature, while others were not. For instance, Clauset et al. (2009) found sufficient statistical evidence to reject claims of scale-free behavior for earthquake intensity and metabolic degree networks, while there was insufficient evidence to reject such claims for networks of protein interaction, Internet, and species per genus.

It is possible to come to opposite conclusions, depending on whether or Kolmogorov–Smirnov (KS) statistics are used to test the hypothesis whether a network is scale-free, i.e. follows a (possibly truncated) power-law distribution. Indeed, Khanin and Wit (2006) obtained a p value of for goodness-of-fit for a truncated power-law distribution for the protein–protein interaction data from Ito et al. (2000), while Clauset et al. (2009) obtained a p value of 0.31 for KS goodness-of-fit for a truncated power-law for the same data.

In this paper, we introduce the first efficient algorithm to compute the exact number of homopolymer RNA secondary structures having k neighboring structures, for each value of k, that can be reached by adding or deleting one base pair. Since there are exponentially many secondary structures, our time and space algorithm uses dynamic programming. By applying the Kolmogorov–Smirnov test, we then show that homopolymer RNA secondary structure networks are not scale-free. We also provide evidence that the same is true for real RNA networks. Prior to this paper, only fragmentary results were possible by exhaustively enumerating all secondary structures having free energy within a certain range obove the minimum free energy (Wuchty 2003).

Our work investigates properties of the ensemble of RNA secondary structures, considered as a network, and so extends results of Clote (2015), which described a cubic time dynamic programming algorithm to compute the expected network degree. The RNA connectivity algorithm described in Sect. 2.3 is completely unrelated to that of Clote (2015), and allows one to compute all finite moments, including mean, variance, skew, etc.

The plan of the remaining paper is as follows. Section 2 presents a brief summary of basic definitions, followed by a description of an efficient dynamic programming algorithm to determine the absolute [resp. relative] frequencies N(k) [resp. p(k)] for secondary structure connectivity of a given homopolymer, which allows non-canonical base pairs. Section 3 presents the statistical methods used to both fit RNA connnectivity data to a power-law distribution and to perform a goodness-of-fit test using Kolmogorov–Smirnov distance. Section 4 shows that RNA networks are not scale-free, by performing (computationally efficient) Kolmogorov–Smirnov bootstrapping tests. Section 5 presents concluding remarks, while the Appendix presents data that suggests that RNA networks satisfy a type of preferential attachment. The rigorous proof that RNA networks satisfy modified form of preferential attachment is suppressed for reasons of space, but is available in the preprint (Clote 2018).

Computing degree frequency

Section 2.1 presents basic definitions and notation used; Sect. 2.2 presents an algorithm to compute the frequency of each degree less than K in the ensemble of all secondary structures with run time and memory requirements . Section 2.3 presents a more efficient algorithm, with run time and memory requirements , for the special case of a homopolymer, in which all possible non-canonical base pairs are permitted. We implemented both algorithms in Python, cross-checked for identical results, and call the resulting code RNAdensity. Since this paper is a theoretical contribution on network properties, we focus only on homopolymers and do not present the details necessary to extend the algorithm of Sect. 2.2 to non-homopolymer RNA, where base pairs are required to be Watson–Crick or GU wobble pairs—such an algorithm is possible to develop, using ideas of Sect. 2.2; however, since the resulting complexity is formidible, with time and space requirements, and since there are no obvious applications, we do not pursue such an extension.

Preliminaries

A secondary structure for a length nhomopolymer is a set s of base pairs (i, j), such that (1) there exist at least unpaired bases in every hairpin, where is usually taken to be 3, though sometimes 1 in the literature, (2) there are no basd triples, so for , if , then and , (3) there do not exist base pairs , such that ; i.e. a secondary structure is a type of outerplanar graph, where each base pair satisfies . The free energy of a homopolymer secondary structure s is defined to be times the number |s| of base pairs in s [Nussinov–Jacobson energy model (Nussinov and Jacobson 1980)]. Since entropic effects are ignored, this is not a real free energy; however it allows us to use the standard notation “MFE” for ‘minimum free energy’. Note that the MFE structure for a length n homopolymer has many base pairs.

For a given RNA sequence, consider the exponentially large network of all its secondary structures, where an undirected edge exists between any two structures s and t, whose base-pair distance equals one—in other words, for which t is obtained from s by either removing or adding one base pair. The connectivity (or degree) of a node, or structure, s is defined to be the number of secondary structures obtained by deleting or adding one base pair to s—this corresponds to the so-called move set (Flamm et al. 2000). At the end of the paper, we briefly consider the move set, where the degree of a structure s is defined to be the number of secondary structures obtained by adding, deleting or shifting one base pair (Bayegan and Clote 2015). The [resp. ] connectivity of the MFE structure for a homopolymer of length n is [resp. ]. ConnectivityN(k) is defined to be the absolute frequency of degree k, i.e. the number of secondary structures having exactly k neighbors, that can be obtained by either adding or removing a single base pair. The degree densityp(k) is defined to be the probability density function (PDF) or relative frequency of k, i.e. the proportion of all secondary structures having k neighbors, where Z denotes the total number of secondary structures for a given homopolymer. A network is defined to be scale-free, provided its degree frequency N(k) is proportional to a power-law, i.e. where is the scaling factor.

Computing the degree density

In this section, we describe a novel dynamic programming algorithm to compute the degree densityp(k) for the network of secondary structures for a homopolymer of length n. Note first that the empty structure of length n hasmany neighbors, each obtained by adding a base pair. Indeed,Using a simple induction argument, Eq. (1) implies that for all values of n, the maximum possible degree, , of a secondary structure for the length n homopolymer is

Let N(i, j) denote the number of secondary structures on interval [i, j], computed the following simple recurrence relation from Stein and Waterman (1978): for , set , and for setor more simplyAlthough Eq. (2) requires time and space, it can trivially be extended to compute the number of secondary structures for an arbitary RNA sequence , where base pairs are either Watson–Crick or wobble pairs. If no such extension is necessary, then the recurrence relation Eq. (3), first given in Stein and Waterman (1978), requires time and O(n) space, hence is more efficient by a factor of n. In a similar fashion, the recurrence relations (5–12) and pseudocode in Sect. 2.2 are given in a form that allows an extension (not given here) to the general case of computing the degree density for the ensemble of secondary structures of a given RNA sequence . The top-level pseudocode given in Algorithm 1 requires time and storage; however, in the next section, we improve this by a factor of n, both in time and space requirements.

Suppose that every hairpin loop is required to have at least unpaired positions; i.e. if (i, j) is a base pair, then . As described in the Eqs. (5–12) for Base Cases A–D, let Z(i, j, k, h, v) denote the number of secondary structures on the interval [i, j], for for the homopolymer model, that have exactly k neighbors, and for which there are exactly h unpaired positions (or holes) in , and for which there are exactly visible positions among . Concretely, this means that either(i) and the rightmost positions in the interval [i, j] are all unpaired, or(ii) that , and that position is paired to some . In base case D and inductive case D below, we will treat the two possible subcases of (i), in which the rightmost positions are unpaired – namely, the subcase (i) in which is unpaired (hence the rightmost positions are unpaired), and the subcase (i) in which position is paired with some .

Let denote the number of secondary structures on the interval [i, j] that have exactly k neighbors with respect to the move set (i.e. have degree k), so thatRecalling from Eq. (1) that , for any , we clearly have thatIn the sequel, we describe Base Cases A–D, which initialize the arrays Z(i, j, k, h, v) and , followed by Inductive Cases A–D, which treat the corresponding updates within the for-loops of the following pseudocode. Since arrays are initially set to zero, all updates to the arrays will be performed by adding a value val to the current value held in the array, so we will write and , which is a standard abbreviation for the assignments and . Explatory comments in the pseudocode are indicated by a double-slash. In Algorithm 1, assume a positive integer input of n to indicate a length n homopolymer.

In line 1, arrays are initialized to zero, then Base Cases A–D are applied; lines 2–9 then treat the Inductive Cases A–D. In this dynamic programming (DP) algorithm, the idea is to define for all intervals [i, j] where , after having computed and stored the values for for all intervals [i, j] where . All secondary structures of the interval [i, j] can be partitioned into structures having exactly degree k (i.e. k neighbors, in which k structures that can be obtained by either adding or removing a single base pair). To support an inductive argument, in proceeding from interval [i, j] to , we need additionally to determine the number of structures having degree k, which have a certain number h of positions that are visible (external to every base pair), which can be paired with the last position . Note that the position can not be base-paired with j in [i, j]; however, can be base-paired with j in . Thus in addition to keeping track of the number h of holes (positions in that are external to all base pairs, hence can be paired with j), we introduce the variable v to keep track of the number of visible positions in . This explains our need for the function Z(i, j, k, h, v) as defined in the Eqs. (5–12) for Base Cases A–D. We now proceed to the details, where for ease of the reader, some definitions are repeated.

Let denote the minimum number of unpaired positions required to be present in a hairpin loop. For a length n homopolymer, let , , , . Recall that Z(i, j, k, h, v) denotes the number of secondary structures on the interval [i, j], for for the homopolymer model, that have exactly k neighbors, and for which there are exactly h unpaired positions (or holes) in , and for which there are exactly visible positions among . For , this means that position is base-paired to some while positions are not base-paired to any position in [i, j]. When , this means simply that the rightmost positions in the interval [i, j] are all unpaired.

Recall as well our definition that . We begin by initializing for all values in corresponding ranges. Letting N(i, j) denote the number of secondary structures on [i, j] for the homopolymer model, as computed by Eq. (2), the following recurrences describe an algorithm that requires storage and time to compute the probability that a (uniformly chosen) random secondary structure has degree k for , where K is a user-defined constant bounded above by .

Base Case A considers all structures on [i, j], as depicted in Fig. 1, that are too small to have any base pairs, hence which have degree zero.

Fig. 1

Structures considered in base case A

Base Case A: For , define

Base Case B considers all structures on [i, j], as depicted in Fig. 2, that have only base pair (i, j), since other potential base pairs would contain fewer than unpaired bases. The degree of such structures is 1, since only one base pair can be removed, and no base pairs can be added. Moreover, no position in [i, j] is external to the base pair (i, j), so visibility parameters . The arrow in Fig. 2 indicates that the sole neighbor is the empty structure, obtained by removing the base pair (i, j).

Fig. 2

Structures considered in base case B

Base Case B: For and (i, j) is a base pair, define

Base Case C considers the converse situation, consisting of the empty structure on [i, j] where , whose sole neighbor is the structure consisting of base pair (i, j). The arrow is meant to indicate that the structure on the right is the only neighbor of that on the left, as depicted in Fig. 3. Since the size of the empty structure on [i, j] is and every position in [i, j] is visible (external to every base pair), and . the dotted rectangle in Fig. 3 indicates the unpaired positions at the right extremity as counted by .

Fig. 3

Structures considered in base case C

Base Case C: For and (i, j) not base-paired, define

Base Case D considers the empty structure on [i, j] where . The empty structure is the only structure having degree maxDegree, since maxDegree(i, j) many base pairs can be added to the empty structure. In Fig. 4, the dotted rectangle indicates the rightmost unpaired positions, corresponding to visibility parameter , while dotted circles indicate the holes, i.e. unpaired positions that could be paired with the rightmost position j.

Fig. 4

Structures considered in base case D

Base Case D: For all , the empty structure, as indicated by (so and ), has degree maxDegree(i, j) as defined by Eq. 1, where

Inductive Case A considers the case where left and right extremities i, j form the base pair (i, j), where . No position in [i, j] is visible (external to all base pairs), so visibility parameters . Recalling the definition of from Eq. 4, we have the following.

Inductive Case A: For and (i, j) base-paired in [i, j],From this point on, we use the operator , so that the previous equation would be written as (see Fig. 5).

Fig. 5

Structures considered in inductive case A

Inductive Case B considers the case where last position j base-pairs with the r, where . The value has already been considered in Inductive Case A, and values cannot base-pair to j, since the corresponding hairpin loop would constain less than unpaired positions. This situation is depicted in Fig. 6, where there are h holes (positions in that are external to all base pairs) and no visible positions in .

Fig. 6

Structures considered in inductive case B

Inductive Case B: For and (r, j) base-paired in [i, j] for some ,When implemented, this requires a check that .

For each value , inductive Case C(v) considers the case where position forms a base pair with position . The value is not considered here, since it was already considered in Inductive Cases A,B. Note that a structure s of the format has k neighbors, provided the restriction of s to has neighbors, and the restriction of s to has neighbors, where . The term vh is due to the fact that since base pair ensures that all holes are located in , hence located at more than distance from all visible positions in , a neighbor of s can be obtained by adding a base pair from any hole to any visible suffix position—there are vh many such possible base pairs that can be added. Finally, the last term is present, since one neighbor of s can obtained by removing base pair . This explains the summation indices and summation terms in Eq. (11). Figure 7 depicts a typical structure considered in case C(v).

Fig. 7

Structures considered in inductive case C(v)

Inductive Case C(v), for: For and base-paired in [i, j], for some , where are unpaired in [i, j],The first term handles the subcase where , so that is a base pair, while the second term handles the subcase where . Note that when implemented, this requires a test that .

Case D considers the case where there are h holes, and positions are unpaired, so that . Note that implies only that are unpaired, so Case D includes the addition requirement that position is unpaired. Structures s satisfying Case D can be partitioned into subcases where the restriction of s to has holes in , and visible positions in . Note that , accounting for the h holes in structure s in , and that it is essential that , since the case was considered in Case .

The term is due to the fact that the rightmost position in the restriction of s to can base-pair with position j, but not with , etc. since this would violate the requirement of at least unpaired bases in a hairpin loop. Similarly, the second rightmost position in the restriction of s to can base-pair with positions j and , but not with , etc.; as well, the third rightmost position can base-pair with positions j, and , but not with , etc. The number of neighbors of s produced in this fashion is thus . Finally, the term is due to the fact that each of the holes in the restriction of s to can base-pair to each of the positions in .

The argument just given shows the following. Let s be a structure that satisfies conditions of Case D with h holes and visible positions, and suppose that the restriction of s to has holes and w visible positions. Then s has k neighbors provided that the restriction of s to has neighbors on interval . The Eq. (12) now follows.

Inductive Case D: For and unpaired in [i, j], and ,

As in Case C(v), when implemented, this requires a test that (see Fig. 8).

Fig. 8

Structures considered in inductive case D

Our implementation of Eqs. (5–12) has been cross-checked with exhaustive enumeration; moreover, we always have that , so the degree density is correctly computed.

Faster algorithm in the homopolymer case

The algorithm described in Sect. 2.2 requires time and space, where K is a user-specified degree bound . By minor changes, that algorithm can be modified to compute the degree density function for any given RNA sequence . In the case of a homopolymer, any two positions are allowed to base-pair (regardless of whether the base pair is a Watson–Crick or wobble pair), provided only that every hairpin loop contains at least unpaired positions. For homopolymers, we have a faster algorithm that requires time and space. Since nucleotide identity is unimportant, instead of Z(i, j, k, h, v), we describe the function , where m corresponds to the length of interval [i, j].We begin by initializing for all , , , and . If , we assume that .

Base Case A: For , defineBase Case B: For , defineBase Case C: For , defineBase Case D: For all , defineInductive Case A: For and , defineInductive Case B: For , , and , defineWhen implemented, this requires a check that .

Inductive Case C(v): For and , defineInductive Case D: For , , and ,Note that h is strictly less than , since the case occurs only when additionally , which only arises in the empty structure. The general case for the empty structure was handled in Base Case D. When implemented, this requires a check that .

Statistical methods

Current software for probability distribution fitting of connectivity data, such as Matlab™, Mathematica™, R and powerlaw (Alstott et al. 2014), appear to require an input file containing the connectivity of each node in the network. In the case of RNA secondary structures, this is only possible for very small sequence length. To analyze connectivity data computed by the algorithm of Sect. 2.3, we had to implement code to compute the maximum likelihood estimation for scaling factor in a power-law fit, the optimal degree beyond which connectivity data is fit by a power-law, and the associated p value for Kolmogorov–Smirnov goodness-of-fit, as described in Clauset et al. (2009). We call the resulting code RNApowerlaw. This section explains those details.

Recall the definition of the zeta functionWe use both the generalized zeta function (22), as well as the truncated generalized zeta function (23), defined respectively byGiven a data set of positive integers in the range , the likelihood that the data fits a truncated power-law with scaling factor and range is defined byRather than sampling individual RNA secondary structures to estimate the connectivity of the secondary structure network for a given homopolymer, the algorithms from Sects. 2.2 and 2.3 directly compute the exact number N(k) of secondary structures having degree k, for all k within a certain range. It follows that the likelihood that secondary structure connectivity fits a power-law with scaling factor is given byhence the log likelihood is is given byThe parameter which maximizes the log likelihood is determined by applying SciPy function minimize (with Nelder-Mead method) to the negative log likelihood, starting from initial estimate , taken from equation (3.7) of Clauset et al. (2009)which in our notation yieldsIn results and tables of this paper, we often write the maximum likelihood estimate (MLE) simply as .

We compute the Kolmogorov–Smirnov (KS) p value, following (Clauset et al. 2009), as follows. Given observed relative frequency distribution D and a power-law fit P with scaling factor , the KS distance is defined to be the maximum, taken over all of the absolute difference between the cumulative distribution function (CDF) for the data evaluated at k, and the CDF for the power-law, evaluated at kwhere and are the actual and fitted cumulative density functions, respectively. The KS p value for the fit of data D by power-law P with scaling factor , is determined by (1) sampling a large number () of synthetic data sets from a true power-law distribution with scaling factor , (2) computing the KS distance between each synthetic data set and its power law fit with MLE scaling factor , (3) reporting the proportion of KS distances that exceed the KS distance between the original observed data set and its power-law fit with scaling factor .

Following (Clauset et al. 2009), is chosen to be that degree , such that the KS distance for the optimal power-law fit is smallest. In contrast, is always taken to be the maximum degree in the input data. Our computation of p value for goodness-of-fit follows Sect. 4.1 of Clauset et al. (2009), with the exception that we not generate any synthetic data for values , since the MLE scaling factor is determined for the (normalized) distribution of data values in the interval , a convention followed in Alstott et al. (2014). We have implemented Python code to compute , , , KS distance, p value, etc. as described above. In Sect. 4, we compare results of our code with that from powerlaw (Alstott et al. 2014) for very small homopolymers. Though our code does not do lognormal fits, this is performed by powerlaw, where the density function for the lognormal distribution with parameters is defined byIn computing the p value for power-law goodness-of-fit using Kolmogorov–Smirnov statistics, it is necessary to sample synthetic data from a (discrete) power-law distribution with scaling factor , a particular type of multinomial distribution. Given an arbitrary multinomial distribution with probability for each , it is straightforward to create M synthetic data sets, each containing N sampled values, in time O(mNM); however, since and N is the (exponentially large) number of all secondary structures having degrees in , the usual sequential method would require prohibitive run time. Instead, we implemented the much faster conditional method (Malefaki and Iliopoulos 2007). Our goal is to sample from a multinomial distribution given bywhere is the number of degrees in the synthetic data, and in the sample set of size N there are many occurrences of degree . To do this, we sample from the binomial distribution of N coin tosses with heads probability , then from the binomial distribution of coin tosses with heads probability , then from the binomial distribution of coin tosses with heads probability , etc. where each is determined with the function binom from Python Scipy.stats.

Results

Below, we use the algorithms described in previous sections to compute RNA secondary structure connectivity, determine optimal scaling factor and minimum degree for a power-law fit, then perform Kolmogorov–Smirnov bootstrapping to determine the goodness-of-fit for parameters . In Appendix A, we show that preferential attachment appears to hold for the network of RNA structures, at least for our definition of preferential attachment.

Analysis of RNA networks using RNAdensity and RNApowerlaw

The algorithm RNAdensity described in Sect. 2.3 was used to compute absolute and relative degree frequencies for the following cases: (1) homopolymers of length with for maximum possible degree upper bound , (2) homopolymers of length with , where degree upper bound for and for , (3) homopolymers of length with , where degree upper bound for and for . For small homopolymers of length at most 30, optima values for , power-law scaling factor , Kolmogorov–Smirnov distance were determined using software powerlawpowerlaw (Alstott et al. 2014) as well as RNApowerlaw from Sect. 3. Table 1 summarizes these results, which show the agreement between powerlaw and RNApowerlaw. Moreover, both both programs indicate that formal hypothesis testing rejects the null hypothesis that a power-law distribution fits connectivity data; indeed, powerlaw determines a negative log odds ratio R for the logarithm of power-law likelihood over lognormal likelihood, indicating a better fit for the lognormal distribution, and RNApowerlaw determines small p values for Kolmogorov–Smirnov goodness-of-fit of a power-law distribution. Figure 9a shows connectivity density function for a 100-mer, with overlaid Poisson and lognormal distributions—since Erdös-Rényi random graphs have a Poisson degree distribution (Albert and Barabási 2002), it follows that RNA secondary structure networks are strikingly different than random graphs. Figure 9b shows a portion of the power-law fit for degrees in , where scaling factor and . Although maximum degree probability at is less than 0.05 for the raw data, the connectivity density for is normalized, which explains why the degree probability for is . Visual inspection might suggest an excellent fit for the power-law distribution; however, a Kolmogorov–Smirnov p value of indicates that the distribution is not power-law. The seemingly good power-law fit for RNA connectivity data, suggested by visual inspection, however motivated our initial investigation of preferential attachment.

Table 1

Table comparing goodness-of-fit computations for software powerlaw (Alstott et al. 2014) and RNApowerlaw for homopolymer lengths less than 30 nt

n	S_n	k_min	α (PL)	α (RNAPL)	KSdist (PL)	KSdist (RNAPL)	⟨KSdist⟩	log odds ratio R (PL)	p-val for R (PL)	p-val (RNAPL)
10	65	3	3.13752	3.13753	0.05576	0.05576	0.02721	-0.15	0.765	0.813
12	274	4	3.23011	3.23011	0.03650	0.03650	0.01277	-0.81	0.482	0.746
14	1184	5	3.38933	3.38935	0.02021	0.02021	0.00669	-1.70	0.270	0.699
16	5223	6	3.51285	3.51289	0.02252	0.02253	0.00603	-6.78	0.029	0.051
18	23,434	9	3.79069	3.79073	0.02333	0.02333	0.00624	-16.00	0.001	0.001
20	106,633	10	3.87168	3.87165	0.02116	0.02116	0.00581	-82.12	0.000	0.000
22	490,999	10	3.85806	3.85809	0.02304	0.02304	0.00523	-670.64	0.000	0.000
24	2,283,701	14	4.16480	4.16477	0.02242	0.02242	0.00484	-1452.24	0.000	0.000
26	10,713,941	15	4.24485	4.24486	0.02298	0.02298	0.00417	-7129.42	0.000	0.000
28	50,642,017	16	4.33086	4.33089	0.02167	0.02168	0.00347	-33,020.89	0.000	0.000
30	240,944,076	–	–	4.33681	–	0.02393	0.00298	–	–	0.000

Given homopolymer length n, the connectivity density is computed over all secondary structures for (all possible) degrees using the algorithm described in Sect. 2.3. Program powerlaw requires an input file containing the degrees of all structures (i.e. containing values, where is the exponentially large number of all secondary structures), while our program RNApowerlaw requires as input a list of degrees and their (absolute) frequencies. Table headers as follows: n is homopolymer length, is the number of all secondary structures, is the maximum likelihood value for the scaling factor of the optimal power-law fit, as computed by powerlaw (PL) and RNApowerlaw (RNAPL), KSdist is the Kolmogorov–Smirnov (KS) distance using Eq. (29), is the mean KS-distance obtained by replacing ‘max’ by ‘mean’ in Eq. (29), R is the log-odds ratio with associated p value as computed by powerlaw, and the p value in the last column is computed by RNApowerlaw as described in Sect. 3. Since powerlaw required more than 24 h for the computation when , we did not attempt a computation for ; in contrast, RNApowerlaw requires a few seconds computation time. Since the log-odds ratio R is the logarithm of the power-law likelihood divided by lognormal likelihood, a negative value indicates that the lognormal distribution is a better fit for the tail of RNA secondary structure connectivity data. A small p value computed by RNApowerlaw indicates that RNA connectivity data is not well-approximated by a power-law distribution. While our code RNApowerlaw computes the p value for the power-law fit, Alstott’s code powerlaw only computes the p value for the log-odds ratio test

Fig. 9

a Connectivity degree distribution for homopolymer of length 100 where , computed with the algorithm described in Sect. 2.3 for all degrees bounded by . There are secondary structures for the 100-mer (exact number 6.31986335936396855341222902079183), and of the structures have degree bounded by K. Using the output degree densities, the degree mean [standard deviation] is [resp. ]; note that the mean computed from the algorithm in Sect. 2.3 is very close to the exact degree mean of , computed over all structures using the different dynamic programming algorithm in Clote (2015). The Poisson distribution (blue curve) with same mean is shown, as well as the lognormal distribution (red) with parameters and ; i.e. [resp. ] is the mean [resp. standard deviation] for logarithms of the connectivity degree—see Eq. (30). b Power-law fit of tail with scaling factor and , determined by maximum likelihood. Kolmogorov–Smirnov (KS) distance for the fit is 0.01213—see Eq. (29), while average KS distance for the alpha power-law fit 0.00400. Nevertheless, since the p value 0 (to 10 decimal places), hypothesis testing would reject the null hypothesis that the power-law distribution is a good fit for the tail (color figure online)

Since powerlaw requires input files of (individually observed) connectivity degrees, when creating Table 1, we could not run powerlaw for homopolymer length greater than 28, for which latter the input file contained 50, 642, 017 values. A potentially attractive alternative is to generate input files consisting of many occurrences of the value k, where denotes the total number of samples, and where relative frequency p(k) is the proportion of structures having degree k. However, Table 2 shows that neither scaling factor nor are correct with this alternative approach, even for small homopolymers of length 20, 30 and 40. This table justifies the need for our implementation of RNApowerlaw as described in Sect. 3. Table 3 shows maximum likelihood scaling factors and Kolmogorov–Smirnov p values for optimal power-law fis of connectivity data for homopolymers of lengths from 30 to 150.

Table 2

Table showing that approximate [resp. exact] scaling factor [resp. ] and minimum degree for optimal power-law fit of homopolymer connectivity data can not be reliably computed by using software powerlaw (Alstott et al. 2014) on data sampled from relative frequencies

N	10^2	10^3	10^4	10^5	10^6	10^7	RNAPL	S_n
α_0, n=20	6.58318	3.66505	3.93389	3.86017	3.84749	3.84657	3.84648	106633≈1.1·10^5
k_min	10	7	10	10	10	10	10	–
α_0, n=30	5.27581	4.42183	4.46307	4.35008	4.32651	4.32272	4.32213	240944076≈2.4·10^8
k_min	12	13	16	16	16	16	16	–
α_0, n=40	5.15978	5.09714	5.03719	5.24488	5.16985	5.70916	5.94561	633180247373≈6.3·10^11
k_min	15	19	23	29	29	42	49	–
α, n=20	6.76575	3.70988	3.96139	3.88570	3.87271	3.87180	3.87165	106633≈1.1·10^5
k_min	10	7	10	10	10	10	10	–
α, n=30	5.33162	4.44651	4.47963	4.36511	4.34122	4.33739	4.33681	240944076≈2.4·10^8
k_min	12	13	16	16	16	16	16	–
α, n=40	5.19197	5.11604	5.049419	5.25365	5.17824	5.65206	5.95033	633180247373≈6.3·10^11
k_min	15	19	23	29	29	41	49	–

Approximate value is computed from Eq. (27), while is the maximum likelihood estimate (MLE) of the optimal power-law scaling factor. Given homopolymer length , connectivity density is computed over all secondary structures for (all possible) degrees using the algorithm described in Sect. 2.3. Since powerlaw requires input files of (individually observed) connectivity degrees, rather than a histogram of (absolute) frequencies F(k) of connectivity degrees, we generated a file consisting of many occurrences of the value k, where denotes the total number of samples, and where relative frequency p(k) is defined by . In contrast to powerlaw, our program RNApowerlaw (RNAPL) computes exact values from connectivity degree (absolute) frequencies. When using powerlaw, it is clearly necessary to create input files of ever-increasing sizes N, in order to have more accurate values of , and . Since the number of RNA secondary structures is exponential in homopolymer length n, it rapidly becomes impossible to use powerlaw for large RNAs—for instance, table values for required an overnight run of powerlaw, while our software returned the exact value within a few seconds

Table 3

Table showing maximum likelihood scaling factors with associated p values for optimal power-law fits of RNA secondary structure connectivity data for homopolymers of length to 150

n	k_max	% of S_n	k_peak	k_mfe	k_min	α(k_min,k_max)	KS(k_min,k_max)	⟨KS⟩	p val
30	60	0.998861	10	13	16	4.223674	0.014393	0.004877	0.000000
35	70	0.999174	12	16	18	4.395936	0.015391	0.005284	0.000000
40	80	0.999404	14	18	23	4.736679	0.015298	0.004859	0.000000
45	90	0.999563	16	21	30	5.146670	0.012075	0.004291	0.000000
50	100	0.999681	18	23	32	5.310801	0.012421	0.004345	0.000000
55	110	0.999762	20	26	39	5.674231	0.011649	0.003979	0.000000
60	120	0.999823	22	28	41	5.829310	0.012328	0.003979	0.000000
65	130	0.999866	24	31	49	6.200720	0.010772	0.003572	0.000000
70	140	0.999899	26	33	52	6.386452	0.010836	0.003464	0.000000
75	150	0.999923	28	36	60	6.721588	0.009719	0.003151	0.000000
80	160	0.999941	31	38	63	6.897103	0.009818	0.003067	0.000000
85	170	0.999955	33	41	67	7.097544	0.009737	0.002940	0.000000
90	180	0.999965	35	43	74	7.373569	0.008916	0.002726	0.000000
95	190	0.999973	37	46	78	7.564208	0.008755	0.002615	0.000000
100	200	0.999979	40	48	83	7.775022	0.008444	0.002476	0.000000
105	135	0.999388	42	51	67	7.204937	0.010712	0.003853	0.000000
110	140	0.999432	44	53	70	7.360192	0.010810	0.003854	0.000000
115	145	0.999474	46	56	73	7.513728	0.010889	0.003852	0.000000
120	150	0.999512	49	58	77	7.706717	0.010405	0.003703	0.000000
125	155	0.999549	51	61	80	7.856962	0.010504	0.003696	0.000000
130	160	0.999582	53	63	84	8.045458	0.010102	0.003556	0.000000
135	165	0.999614	55	66	88	8.231267	0.009724	0.003425	0.000000
140	170	0.999643	58	70	91	8.377410	0.009809	0.003418	0.000000
145	175	0.999670	60	71	94	8.522170	0.009884	0.003413	0.000000
150	180	0.999695	62	75	98	8.703041	0.009515	0.003289	0.000000

Absolute and relative connectivity degree frequencies were computed by RNAdensity from Sect. 2.3, while the optimal parameters and p values were computed by RNApowerlaw from Sect. 3. Column headers are as follows: n is sequence length, is the degree upper bound K for RNAdensity, of indicates the proportion of all secondary structures having degree bounded by , is the location of the density maximum, is the degree of the minimum free energy structure (having largest number of base pairs), is the optimal lower bound for a power-law fit, is the maximum likelihood scaling factor for power-law fit, is the Kolmogorov–Smirnov (KS) distance between connectivity data and power-law fit, p val is goodness-of-fit p value for Kolmogorov–Smirnov statistics, and is the average KS distance, obtained by replacing ‘max’ by ‘mean’ in Eq. (29)

Figure 10a shows a scatter plot with regression line for the cut-off values , defined to be the least value such that the probability that a secondary structure for length n homopolymer has degree greater that is at most 0.01. From this figure, we determined that for homopolymer length , it more than suffices to take degree upper bound . Figure 10b shows the connectivity degree distribution for a homopolymer of length 20, where degree dg(s) is redefined to be the number of structures t that can be obtained from s by adding, removing, or shifting a base pair in s. The so-called move set, consisting of an addition, removal or shift of a base pair is the default move set used in RNA kinetics software kinfold (Lorenz et al. 2011). Although a dynamic programming algorithm was described in Clote and Bayegan (2015) to compute the average network degree, the methods of this paper do not easily generalize to connectivity densities. Figure 11 shows a least-squares regression line for the log-log density plot for connectivity (computed by brute-force) for a homopolymer of length 20, together with an optimal power-law fit computed by RNApowerlaw. Since there are only 106.633 secondary structures for the 20-mer with , we ran powerlaw on connectivity data, which determined , , and a log odds ratio with p value of 0.248. Since RNApowerlaw determined , , and a Kolmogorov–Smirnov p value of 0.219, we can not reject the null hypothesis that a power-law distribution fits the tail of connectivity data for a 20-mer.

Fig. 10

a Plot of the least cut-off value as a function of homopolymer length n, for . Here is defined as the least value such that the probability that a secondary structure for length n homopolymer has degree greater that is at most 0.01. For the least-squares fit, the regression equation is , with p value of for slope value, and p value of for the y-intercept. b connectivity for the 106,633 secondary structures for a 20-nt homopolymer with (green shaded curve), with Poisson distribution of the same mean. Connectivity values range from (with many intermediate gaps before the max degree). The distribution mean [resp. standard deviation] is [resp. ]; these values should be contrasted with the corresponding values of [resp. ] for connectivity for the same 20-nt homopolymer (data not shown) (color figure online)

Fig. 11

a Plot of as a function of for the degree distribution for connectivity of the 20-nt homopolymer with , for degrees . The distribution tail appears to satisfy a power-law with exponent , i.e. , where x is degree and p(x) is the relative frequency of the number of nodes having degree x (regression equation log-log plot is ). b It is well-known that linear regression of the log-log plot is less reliable than using maximum likelihood when establishing whether the tail of empirical data is fit by a power-law distribution. For the connectivity data of a 20-nt homopolymer, the maximum likelihood estimation (MLE) of optimal power-law scaling factor is with p value is 0.219 when and . Since the p value is not less than 0.05, we can not reject the null hypothesis that connectivity is well-fit by a power-law distribution (color figure online)

Conclusion

Since the pioneering work of Zipf on the scale-free nature of natural languages (Zipf 1949), various groups have found scale-free networks in diverse domains ranging from communication patterns of dolphins (McCowan et al. 2002), metabolic networks (Jeong et al. 2000), protein–protein interaction networks (Ito et al. 2000; Schwikowski et al. 2000), protein folding networks (Bowman and Pande 2010), genetic interaction networks (Tong et al. 2004; Van Noort et al. 2004) to multifractal time series (Budroni et al. 2017). These discoveries have galvanized efforts to understand biological networks from a mathematical and topological standpoint. Using mathematical analysis, Barabasi and Albert (1999) established that scale-free networks naturally emerge when networks are dynamic, whereby newly accrued nodes are preferentially connected to nodes already having high degree. On such grounds, one might argue that protein folding networks and protein–protein interaction (PPI) networks should exhibit scale-free properties, since nature is likely to reuse and amplify fast-folding domains—cf. Gilbert’s exon shuffling hypothesis (Gilbert 1978). Indeed, Cancherini et al. (2010) have established that in 4 metazoan species analyzed (H. sapiens, M. musculus, D. melanogaster, C. elegans) those genes, which are enriched in exon shuffling events, displayed a higher connectivity degree on average in protein–protein interaction (PPI) networks; i.e. such genes had a larger number of interacting partners. On similar grounds that nature should reuse and amplify successful metabolic networks, one might argue that metabolic networks should exhibit scale-free properties. However, rigorous statistical analysis has shown that metabolic networks fail a goodness-of-fit test for scale-free distribution, while PPI satisfy a goodness-of-fit test for scale-free distributions over a certain range of connectivity (Khanin and Wit 2006; Clauset et al. 2009).

There appears to be a current polemic whether certain naturally occurring networks are scale-free. Broido and Clauset (2019) provide statistical arguments that less than 45 of the 927 real-world network data sets (i.e. ) found in the Index of Complex Networks exhibit the “strongest level of direct evidence for scale-free structure”. In a response to a preprint of Broido and Clauset (2019) dated March 6, 2018 and posted on the Barabási Lab web site https://www.barabasilab.com/post/love-is-all-you-need, A.L. Barabási argued against the conclusions of Broido and Clauset (2019). Here, it should be noted that this is not the first time a polemic has arisen about the power-law distribution—indeed, there was a heated exchange between Mandelbrot and Simon almost 60 years ago in the journal Information and Control. For details, references, and a history of the power-law distribution, see Mitzenmacher (2004).

Given the current interest in whether certain naturally occurring networks are scale-free, we have introduced a novel algorithm to compute the connectivity density function for a given RNA homopolymer. Our algorithm requires run time and storage, where K is a user-specified degree bound . Short of exhaustively listing secondary structures by brute-force, no such algorithm existed prior to our work. Since existent software appears unable to perform power-law fitting for exponentially large RNA connectivity data, we have also implemented code to compute and statistically evaluate the maximum likelihood power-law fit for an input histogram, using a very fast method to the density function of a sampled power-law distribution with given scaling parameter. Using the resulting code, called RNAdensity and RNApowerlaw, we have computed the connectivity density function for RNA secondary structure networks for homopolymers of length up to 150. Statistical analysis categorically shows that there is no statistically significant power-law fit for homopolymer RNA secondary structure network connectivity, despite the seemingly good visual fit shown in Fig. 9. Figure 12 shows that secondary structure network connectivity is not scale-free for the (real) 32 nt selenocysteine insertion sequence fruA. Figure 13 shows that the and degree distributions for other naturally occurring RNAs are not scale-free, in particular for the 65 nt RNA thermometer from Klebsiella pneumoniae subsp. pneumoniae with EMBL accession code CP000647.1/1773227-1773291 and the 76 nt alanine transfer RNA from Mycoplasma mycoides with accession code tdbR00000006 from tRNAdb Juhling et al. (2009) (accession code RA1180 from the Sprinzl tRNA database). While the density plot in Fig. 12 was produced by exhaustively enumerating all 971,299 secondary structures of the 32 nt fruA, Figure 13 was produced by enumerating all secondary structures having free energy within 13 kcal/mol of the minimum free energy, as computed by RNAsubopt from the Vienna RNA Package (Lorenz et al. 2011); this procedure generated 1,079,102 secondary structures (out of a total of structures) for the 65 nt fourU RNA, and 408,414 secondary structures (out of a total of structures) for the 76 nt tRNA.

Fig. 12

a- and -degree distribution for the 32 nt selenocystein insertion (SECIS) element fruA with sequence CCUCGAGGGG AACCCGAAAG GGACCCGAGA GG, obtained by brute-force computation from an enumeration of all secondary structures (exact number 971299), ranging in degree from 4 to 123. Average -degree 13.10; average -degree 33.25. Using notation from Table 9, the MLE power-law fit for -degree data has values of , , , , p value of 0.0000. In contrast, the MLE power-law fit for -degree data has values of , , , , p value of 0.729. Summarizing, Kolmogorov–Smirnov statistics indicate that the data is not scale-free, while it cannot be refuted that the data is scale-free. However, the range of degrees for which the data might be scale-free is from 93 to 123, which accounts for only of the density. As shown in (b), even log-log regression suggests that the data is not scale-free. b Log–log plot of -density of fruA with regression equation , determined from the relative frequency of structures having -degree in the range of 29 to 4123, corresponding to the portion of the density starting after the peak of 0.04987 in previous panel at degree (color figure online)

Fig. 13

a- and -degree distribution for the 65 nt fourU RNA from Klebsiella pneumoniae subsp. pneumoniae with sequence GGACAAGCAA UGCUUGCCUU UAUGUUGAGC UUUUGAAUGA AUAUUCAGGA GGUUAAUUAU GGCAC and EMBL accession code CP000647.1/1773227-1773291. FourU RNA is a class of thermometers found in bacteria such as E. coli, Salmonella, V. cholerae, etc. that regulate protein expression by undergoing a conformation change triggered by temperature—for instance, the conformational change of the V. cholerae fourU thermometer at C permits the transcription of a virulence factor. All 1,079,102 secondary structures having free energy within 13 kcal/mol of the minimum free energy (MFE) of this RNA were generated using RNAsubopt from the Vienna RNA Package (Lorenz et al. 2011). The and degree of each secondary structure were determined in order to produce the degree relative frequency histogram. Although the collection of structures having free energy within C of the MFE contains over one million structures (computation required 1–2 days), there are 1995457849526533 () many secondary structures altogether. The average degree is 38.0, while the average degree is 64.2. FourUanalysis: Using RNApowerlaw, , , and p value is 0 (to 10 decimal places). Using powerlaw, , , and the log ratio of power-law fit to log-normal fit is with corresponding p value of —in other words, a log-normal distribution provides a significantly better fit than a power-law distribution for the degree data of this fourU RNA. FourUanalysis: Using RNApowerlaw, , , and p value is 0 (to 10 decimal places). Using powerlaw, , , and the log ratio of power-law fit to log-normal fit is with corresponding p value of —in other words, a log-normal distribution provides a significantly better fit than a power-law distribution for the degree data of this fourU RNA. b- and -degree distribution for the 76 nt alanine transfer RNA from Mycoplasma mycoides with accession code tdbR00000006 from tRNAdb (Juhling et al. 2009) (accession code RA1180 from the Sprinzl tRNA database) with sequence GGGCCCUUAG CUCAGCUGGG AGAGCACCUG CCUUGCACGC AGGGGGUCGA CGGUUCGAUC CCGUUAGGGU CCACCA. All 408414 secondary structures having free energy within 13 kcal/mol of the minimum free energy of this RNA were generated using RNAsubopt from the Vienna RNA Package (Lorenz et al. 2011). The MS1 and MS2 degree of each secondary structure were determined in order to produce the degree relative frequency histogram. Although the collection of secondary structures having free energy within C of the MFE contains about one-half million structures (computation required 1-2 days), there are 877346780605139050 () many secondary structures altogether. The average degree is 38.1, while the average degree is 68.3. tRNAanalysis: Using RNApowerlaw, , , and p value is 0 (to 10 decimal places). Using powerlaw, , , and the log ratio of power-law fit to log-normal fit is , with corresponding p value of —in other words, a log-normal distribution provides a significantly better fit than a power-law distribution for the degree data of this fourU RNA. tRNAanalysis: Using RNApowerlaw, , and p value is 0 (to 10 decimal places). Using powerlaw, , , and the log ratio of power-law fit to log-normal fit is with corresponding p value of —in other words, a log-normal distribution provides a better fit than power-law for the degree data of this tRNA (color figure online)

Since (Day et al. 2003; Kihara and Skolnick 2003) have presented data that suggests that existent protein structures can be explained using only a small number of protein folds, we presented data in Table 4 that suggests that RNA secondary structures may satisfy a type of preferential attachment—a rigorous combinatorial argument establishes this fact for a modified notion of preferential attachment [data not shown, but available in the Appendix of Clote (2018)]. Finally, Python implementations of the algorithms from this paper are publicly available at http://bioinformatics.bc.edu/clotelab/RNAnetworks.

Table 4

Table showing secondary structure preferential attachment probabilities

n	n+1	S_n	S_n+1	Succ	Fail	Succ/(Succ+Fail)	⟨p(s^′,t^′|s,t)⟩
5	6	2	4	5	1	83.33%	0.8333±0.1667
6	7	4	8	18	8	69.23%	0.7222±0.4157
7	8	8	16	90	37	70.87%	0.7748±0.3260
8	9	16	32	419	131	76.18%	0.8105±0.2941
9	10	32	65	1891	575	76.68%	0.8122±0.2887
10	11	65	133	7883	2498	75.94%	0.8125±0.2891
11	12	133	274	33,069	9763	77.21%	0.8300±0.2730
12	13	274	568	142,968	40,797	77.80%	0.8322±0.2709
13	14	568	1184	621,884	171,384	78.40%	0.8366±0.2646
14	15	1184	2481	2,723,993	723,887	79.00%	0.8428±0.2587
15	16	2481	5223	12,041,929	3,108,978	79.48%	0.8478±0.2556
16	17	5223	11,042	53,730,451	13,544,005	79.87%	0.8518±0.2523
17	18	11,042	23,434	241,738,083	59,258,399	80.31%	0.8561±0.2485
18	19	23,434	49,908	1,096,087,115	261,730,198	80.72%	0.8598±0.2455

The first two columns contain homopolymer length n and , followed by the number of secondary structures in and , then the total number of 4-tuples that succeed in demonstrating [resp. fail to demonstrate] preferential attachment, denoted by Succ [resp. Fail]. The next column contains the proportion Succ/(Succ+Fail) of 4-tuples that demonstrate preferential attachment, defined by Eq. (33), while the last column contains the expected preferential attachment , defined by Eq. (35). This expectation is obtained by computing the arithmetical average of the conditional probabilities , defined by

As an afternote, our personal opinion is that it doesn’t much matter whether a naturally occurring network arising from physical phenomena is precisely scale-free or not. If network connectivity appears to follow a power-law distribution, even approximately, then by results of Barabasi and Albert (1999), this suggests that preferential attachment could play a role in how the network may have been constructed by nature. Preferential attachment might well have been a factor in how protein and RNA structures have been formed by evolutionary forces—even in the emergence of stable folds in prebiotic times (Abkevich et al. 1997). It is noteworthy that only a small number of protein folds suffice to explain the diversity of all protein folds found in the Protein Data Bank (PDB) (Kihara and Skolnick 2003): “The number of proteins required to cover a target protein is very small, e.g. the top ten hit proteins can give 90% coverage below a RMSD of 3.5 Å for proteins up to 320 residues long.” As well, the 30 most populated metafolds represent “about half of a nonredundant subset of the PDB” (Day et al. 2003). However, other evolutionary factors seem to be present in the evolution of protein folds, such as kinetic accessibility (Cossio et al. 2010), as well as the ability to switch between alternate conformations (Porter and Looger 2018).