# Whip dynamics in a freely falling chain

Tianyi Guo1, Xiaoyu Zheng2, Peter Palffy-Muhoray1,2
_1 Advanced Materials and Liquid Crystal Institute, Kent State University, OH,
USA_
_2 Department of Mathematical Sciences, Kent State University, OH, USA _

###### Abstract

In this work, we examine the whip dynamics in a freely falling chain. We
consider an inextensible chain with free ends in the presence of gravity
hanging from a fixed pulley. If the configuration is unbalanced, the chain
begins to accelerate and eventually may lift off the pulley. We show that
after the liftoff, the dynamics of the chain resembles that of a whip whose
tip speed diverges. The curved region of the chain plays a key role in the
dynamics, lifting the chain. The chain is in free fall, yet it exhibits whip
dynamics, with one end snapping upward periodically as the chain falls. We
have conducted experiments demonstrating this phenomena.

Keywords: Chain dynamics, whip dynamics, free fall, pulley

##  1 Introduction

Chains and ropes are among our oldest tools, yet, remarkably, they are still
capable of exhibiting surprising behavior. One example is the Mould effect,
discovered by Stephen Mould [1], where a simple ball chain, pulled by gravity
from a cup, forms an astonishing chain fountain, where the chain first rises,
then turns and falls, forming essentially an inverted catenary [2]. Ropes can
show similar behavior [3]. Considerable effort has gone into unraveling the
detailed dynamics in this process [4, 5, 6, 9, 10]. One key ingredient in the
unexpected behavior of chains and ropes is the acquisition of linear momentum
by rigid objects from their environment via hindered rotation [7, 11].

In this paper, we report observations and analysis of a different kind of
unusual behavior: whip dynamics in a freely falling chain. In Section 2, we
introduce three facets of simple chain whip dynamics which lays the foundation
of the essence of whip dynamics. In Section 3, we study in detail the the
dynamics of the free falling chain, both when the chain is on the pulley, and
after liftoff and detachment from the pulley. We establish the analogy of the
free falling chain with simple whip dynamics. We conclude the paper in Sect.
4.

##  2 Prelude: three facets of simple chain whip dynamics

A traditional whip is a handle to which a flexible strip of leather, of a
length of cord had been attached. Chains instead of leather or cord have been
used as whips, frequently as weapons. Unlike leather or cord whips, which are
typically tapered, chain whips are distinguished by uniform mass density and
low dissipation. Below we discuss three aspects of whip dynamics which are
clearly illustrated by simple 1D models. For simplicity, we ignore gravity
here, and consider energy and momentum conservation.

###  2.1 Moving chain with one end arrested

We consider an inextensible uniform chain, such as a bicycle chain, of length
$L$, in a straight configuration, moving at constant velocity $\mathbf{v}_{0}$
in a direction tangent to the chain. Suppose that a hand appears at the front
of the chain, folds the first link clockwise through the angle $\pi$, and
holds it at rest. (This is essentially how the cracking of a whip is
initiated.) Now the chain is folded, as shown in Fig. 1, with a moving
straight portion on top, and a stationary straight portion below; the two
parts are joined with a (negligibly small) semicircular section, which we call
the ’virtual pulley’ (VP), which is also moving in the same direction.

Figure 1: Schematic of moving chain with one end arrested

As time goes on, the top part of the chain, of length $x(t),$ moving with
velocity $\mathbf{v}_{x}(t)$, becomes shorter, while the standing part of the
chain, of length $y(t)=L-x(t)$ at rest, becomes longer; the semicircular
section, which we term the ‘VP’, is moving forward with velocity
$\mathbf{v}_{p}(t)$. It follows that

$\dot{x}=v_{p}-v_{x},$ (1)

and

$\dot{y}=v_{p}.$ (2)

and since the chain is inextensible, $\dot{x}+\dot{y}=0$, we have

$-\dot{x}=v_{p}=\frac{1}{2}v_{x}.$ (3)

To obtain an explicit expression for the dynamics, conservation of momentum
needs to be considered.

There is tension in the chain which is responsible for changing the momentum
of the chain passing through the ‘VP’. The speed of the chain relative to the
pulley is $-\dot{x}$. The change of momentum for mass $\rho\Delta L$ in the VP
is $-2(\rho\Delta L)\dot{x}$ which takes place in time interval $\Delta
t=\Delta L/\dot{x}$ ; this is due to the tension $T$ in the chain at the two
sides of the pulley [4]. That is,

$T=\rho\dot{x}^{2}.$ (4)

The rate of change of momentum in the top part of the chain is111The tension
at one point in the chain which changes the momentum of the attached part does
not involve a change in the length of the attached part. $\rho
x\dot{v}_{x}=-2\rho x\ddot{x}$, and the rate of change of momentum is the
tension $T$, therefore

$\dot{x}^{2}+2x\ddot{x}=0.$ (5)

With initial conditions $x(0)=L$ and $\dot{x}(0)=v_{0}$, the solution is
explicitly given by

$x=L(1-\frac{t}{t_{0}})^{\frac{2}{3}},$ (6)

where

$t_{0}=\frac{4}{3}\frac{L}{v_{0}}$ (7)

is the time when the chain becomes straight again. We see that momentum
conservation is sufficient to arrive at the equation of motion. This result
can also be obtained by conservation of kinetic energy.

The speed $v_{p}$ of the bend, or the ‘VP’, is then

$v_{p}=-\dot{x}=\frac{1}{2}v_{0}(1-\frac{t}{t_{0}})^{-\frac{1}{3}},$ (8)

and both $v_{p}$ and the speed $v_{x}$ of the moving part of the chain diverge
as $t\rightarrow t_{0}$ or as $x\rightarrow 0$. _This is the essence of the
whip dynamics; the conserved kinetic energy becomes localized in the ever
shortening moving part of the whip_. The graphs $x$ and $v_{p}/(2L/(3t_{0}))$
as function of $t$ are shown in Fig. 2. It is interesting note that $v_{p}$
remains nearly constant until $t\simeq t_{0}$, and after this, increases
rapidly. We see that the velocity of the tip of an ideal whip diverges.

(a) (b)

Figure 2: (a) The length of the top portion of the moving chain $x$ vs. $t$.
(b) Speed of virtual pulley $v_{p}$ vs. $t.$

We now calculate the magnitude of the applied force needed to hold the bottom
chain fixed. Since the net force on the standing part of the chain is zero, so
the force and the tension in the chain at the two sides of the ‘VP’ are
balanced, which gives

$F=-T=-\rho\dot{x}^{2}.$ (9)

Substituting Eq. (8) for $\dot{x}$, the applied force is

$F=-\frac{1}{4}\rho v_{0}^{2}(1-\frac{t}{t_{0}})^{-\frac{2}{3}}.$ (10)

The consequence of energy and momentum conservation222We note that the
assumption that the chain is lossless has been implicitly included in the
assumption that the tension in the chain on both sides of the ‘virtual pulley’
is the same. is that the velocity of the tip an ideal lossless whip is
infinite. It appears in this case that the effect of the bend in the chain is
to constrain the kinetic energy, initially distributed in the entire chain, to
the moving part part on one side of the bend. As the bend reaches the end of
the chain, the energy and momentum densities diverge.

###  2.2 Stationary chain with one end pulled

Next, we consider the above scenario viewed from a moving inertial frame of
reference, moving with constant speed $\mathbf{v}_{0}$. In this moving frame,
the chain is initially at rest; then a moving hand appears, and pulls the
stationary chain with velocity $-\mathbf{v}_{0}$. The situation is illustrated
in Fig. 3.

Figure 3: Stationary chain with one end pulled

From this new perspective, the only change is that all the velocities have a
constant $-\mathbf{v}_{0}$ added. In this frame,

$v_{p}=v_{0}(\frac{1}{2}(1-\frac{t}{t_{0}})^{-\frac{1}{3}}-1)$ (11)

and the ‘VP’ starts to move with velocity $-\frac{1}{2}v_{0}$ at $t=0$, and
speeds up reaching $v_{0}$ at $t=t_{0}$. More interestingly, the top part of
the chain has velocity $v_{x}$ given by

$v_{x}=v_{0}((1-\frac{t}{t_{0}})^{-\frac{1}{3}}-1).$ (12)

As soon as the bottom part of the chain is pulled to the left with constant
velocity $v_{0}$, remarkably, the entire top part of the chain moves and
accelerates to the right! Although not easy to observe on a tabletop due to
friction, the effect can clearly be seen. Applying a force in one direction to
one part of the chain, gives rise to an equal and opposite force acting on the
other part. Forces are unchanged when viewed from inertial frames, and the
force accelerating the top part is just the applied force. The folded chain
appears to guide tensile stress; the stress, initiated by the pull, traverses
the bend and the sample, and ‘pulls’ it, causing the far end to move in the
opposite direction from the external pull.

###  2.3 Chain dynamics with no external force

In the third scenario, we consider a chain in space, at $t=0$, folded exactly
in half. The folded chain is horizontal, the fold is at the right. The top
half is moving to the right, with velocity $v_{x}(t)$, the bottom half is
moving to the left, with velocity $v_{y}(t)$, as shown in Fig. 4. The center
of mass of the chain is at rest, since there are no external forces acting on
the chain.

Figure 4: Moving chain without external force

As before, we have that

$\dot{x}=-v_{x}+v_{p},$ (13)

and

$\dot{y}=v_{y}+v_{p}.$ (14)

Since $\dot{x}+\dot{y}=0$, we have

$v_{p}=\frac{1}{2}(v_{x}-v_{y}).$ (15)

Since there are no external forces acting and the initial total momentum is
zero, we have

$\rho xv_{x}-\rho yv_{y}=0.$ (16)

Substituting Eqs. (13) and (14), and $x+y=L$ into Eq. (16), gives

$v_{p}=\dot{x}(\frac{2}{L}x-1).$ (17)

The tension in the chain at the two sides of the ‘VP’ is, as before

$T=\rho\dot{x}^{2}$ (18)

and force balance on the top part gives

$\rho x\dot{v}_{x}=\rho\dot{x}^{2}.$ (19)

From (13), we have $\dot{v}_{x}=\dot{v}_{p}-\ddot{x},$ and from (17)

$\dot{v}_{p}=\ddot{x}(\frac{2}{L}x-1)+\frac{2}{L}\dot{x}^{2}.$ (20)

Substitute those relations in Eq. (19), we arrive at

$\ddot{x}x(\frac{x}{L}-1)+\dot{x}^{2}(\frac{x}{L}-\frac{1}{2})=0.$ (21)

This is the equation of motion. After some manipulation, this can be solved to
give an implicit relation between $x$ and $t,$

$t=\frac{L}{2v_{0}}((1-2\frac{x}{L})\sqrt{\frac{x}{L}(1-\frac{x}{L})}+\frac{1}{2}\sin^{-1}(1-2\frac{x}{L})).$
(22)

Furthermore, we obtain

$\dot{x}=\frac{-v_{0}}{2\sqrt{\frac{x}{L}(1-\frac{x}{L})}}.$ (23)

It follows that $\dot{x}$ diverges as $x\rightarrow 0$, so do $v_{p}$ and
$v_{x}$. A plot of $x$ and $v_{p}$ as function of time are shown in Fig. 5.

(a) (b)

Figure 5: (a) The length of the top portion of the moving chain $x$ vs. $t$.
(b) velocity of the virtual pulley vs. $t$. Figure 6: A schematic showing the
periodic dynamics of a moving chain dynamics. The center of mass is
translating in the direction of arrow, and the system has a nonvanishing
constant angular momentum.

Since the initial momentum is zero, when the chain is straight, chain comes to
rest, except for the last piece; all the kinetic energy of the entire chain
ends up in the tip! The tip will rotate and move rapidly towards the left as
shown in Fig. 6. The part of the chain now at the bottom will move towards the
left and slow down, eventually reaching the conformation identical to the
initial one. In this simple model, the chain will continue this periodic
motion.

In this prelude, we have shown that a folded chain with moving segments can
concentrate all of its kinetic energy at one point - at one of the tips - with
or without an external force. The folded chain also appears to guide tensile
stress around the bend, and give rise to motion there opposite to the
direction of the force causing the motion. If allowed to continue, the
concentrated energy and momentum densities will spread out, and the process
will continue, without dissipation, indefinitely.

We believe these aspects constitute the dynamics normally associated with
whips. Essentially the same phenomena are exhibited by other flexible strips
of material; however, energy losses associated with dissipative processes
which are often present can mask the behavior.

##  3 Freely falling chain

###  3.1 Chain on a fixed pulley

We next consider the behavior of a chain of length $L$ and mass density $\rho$
on a fixed pulley in the presence of gravity. This is illustrated in Fig. 7.
The length of chain on the left is $x$, and at time $t=0$, it is moving up
with velocity $v_{0}$, while the length of the chain on the right is $y=L-x$,
and it is moving down. We assume the radius of the pulley is arbitrarily
small. The chain is pulled down by gravity, and the pulley is supported by an
upward force $F(t)$. The tension in the chain at the pulley is $T(t)$. The
chain is moving with speed $v(t)=-\dot{x}$. The net force on the pulley is

$F-2T=-2\rho\dot{x}^{2}=-2\rho\dot{y}^{2},$ (24)

where the right hand side is the rate of change of momentum of the chain on
the pulley.

Figure 7: Chain on a fixed pulley

Consider the chain on the left. The equation of motion is

$T-\rho xg=-\rho x\ddot{x}.$ (25)

Combining Eqs. (24) and (25) to get

$\rho x\ddot{x}+\rho\dot{x}^{2}-\rho xg=-\frac{1}{2}F.$ (26)

On the right side of the chain, we have

$\rho yg-T=\rho y\ddot{y},$ (27)

and proceeding similarly, we obtain

$-\rho(L-x)\ddot{x}+\rho\dot{x}^{2}-\rho(L-x)g=-\frac{1}{2}F.$ (28)

It is convenient, for simplicity, to use dimensionless variables. We measure
all lengths in units of $L$, time in units of $\sqrt{L/2g}$ and force in units
of $4\rho Lg$. We then have

$x\ddot{x}+\dot{x}^{2}-\frac{1}{2}x=-f,$ (29)

and

$-(1-x)\ddot{x}+\dot{x}^{2}-\frac{1}{2}(1-x)=-f.$ (30)

Subtracting these, we get

$\ddot{x}-x+\frac{1}{2}=0.$ (31)

Eq. (31) is the equation of motion when the chain is on the pulley. This can
be solved at once; the solution is

$x=\frac{1}{2}+Ae^{t}+Be^{-t},$ (32)

where $A$ and $B$ depend on initial conditions $x_{0}$ amd $v_{0}$ as

$A=\frac{1}{2}((x_{0}-\frac{1}{2})-v_{0}),$ (33)

and

$B=\frac{1}{2}((x_{0}-\frac{1}{2})+v_{0}).$ (34)

As expected, if $A<0$, $\dot{x}<0$, and the speed of the chain speed increases
nearly exponentially.

From Eq. (29) the force $f$ on the pulley can be readily calculated,

$f=\frac{1}{4}\left(1-8A^{2}e^{2t}-8B^{2}e^{-2t}\right).$ (35)

As the chain speeds up, the force $f$ exerted on the pulley by the chain
decreases with time and reaches zero at time $t^{\ast}$ given by

$t^{\ast}=\frac{1}{2}\ln\left(\frac{1}{16A^{2}}\left(\sqrt{1-256A^{2}B^{2}}+1\right)\right).$
(36)

At that time, the chain rises and lifts off the pulley! We again regard the
bend in the chain again as a ‘VP’.

A relevant question is the length of the upward moving part of the chain at
liftoff. A straightforward calculation gives

$x=\frac{1}{2}-\frac{1}{2}\sqrt{\frac{1}{2}+2(x_{0}-\frac{1}{2})^{2}-2v_{0}^{2}}.$
(37)

We note that the quantity under the square root needs to be nonnegative, and
that we must have $x<x_{0}$. These lead to

$(\frac{1}{2}-x_{0})^{2}+v_{0}^{2}<\frac{1}{4}.$ (38)

This must be satisfied for liftoff. In particular, starting from rest, with
$x_{0}<\frac{1}{2}$ should always work (with a frictionless pulley!); starting
at $x_{0}=1/2$ should always work if $v_{0}<1/2$.

###  3.2 After liftoff: the freely falling chain

After liftoff, the force from pulley vanishes, that is, $f=0$. Proceeding as
before, we now arrive at, for the chain on the left, in dimensionless units,

$x\ddot{x}+\dot{x}^{2}-x\dot{v}_{p}=0,$ (39)

and for the chain on the right,

$-(1-x)\ddot{x}+\dot{x}^{2}-(1-x)\dot{v}_{p}=0.$ (40)

We can eliminate $\dot{v}_{p}$, and in the presence of gravity, we have the
equations of motion

$2\ddot{x}x(1-x)+\dot{x}^{2}(1-2x)=0,$ (41)
$\dot{v}_{p}=\dot{x}^{2}\frac{1}{2x(1-x)}-\frac{1}{2}.$ (42)

Recall that without gravity, we have

$2\ddot{x}x(1-x)+\dot{x}^{2}(1-2x)=0,$ (43)
$\dot{v}_{p}=\dot{x}^{2}\frac{1}{2x(1-x)}.$ (44)

The only difference is the term $-1/2$ in the equation for $\dot{v}_{p}$, this
is just the dimensionless acceleration of gravity. So we should see
essentially the same phenomenon without gravity (say for the chain on a
frictionless table) except the falling of the center of mass. Without gravity,
we don’t have the initial conditions required for takeoff.

We also have the result that the height $h$ of the center of mass and the
height $h_{p}$ of the center of pulley are related;

$h=-x^{2}+x+h_{p}-\frac{1}{2}.$ (45)

Setting $h=0$, we have $h_{p}(t)$ once we have $x(t)$

$h_{p}=x^{2}-x+\frac{1}{2}.$ (46)

We also have

$v_{p}=(2x-1)\dot{x}.$ (47)

Therefore if $x<0.5$ and $\dot{x}<0$, the ‘VP’ rises.

###  3.3 Experimental Results

We carried out simple experiments to observe the dynamics of chains, initially
on a pulley. Our goal was to compare model predictions from Sects. 3.1 and 3.2
above with experimental observations. Due to the challenging requirements of
making reliable measurements on the length and timescales involved, we only
carried out simple observations at this preliminary stage of work.

Since the timescale of our experiments is $\sqrt{L/g}$, a length of $10m$
corresponds to $1s$. To observe the development of the dynamics and the
transition between the different regimes, we chose to work in an abandoned
stairwell of Kent State University Library. The building has 12 floors, as
well as a basement. The distance from the 10th floor, which was our station,
to the basement floor was $42.7m$. For our experiments, we used $3/16"$
nickel-plated steel ball-chains from McMaster-Carr, with mass density
$\rho=30g/m$ and maximum curvature $1.7cm^{-1}$. Instead of a pulley, we used
a polished steel pipe with $2cm$ o.d.

In the experiment described here, we used a $L=61m$ chain. The chain was
draped over the supporting steel pipe, with equal lengths on each side of the
bar. With the chain at rest, we pulled the chain on one side of the bar
downward, by hand, to make one side longer. At one point in this process,
chain started slip on the supporting pipe, with the longer side moving
downward, and the shorter side moving upward. We took a video of the chain
over the bar; we started the video once the chain was moving with appreciable
speed. The speed of the chain increased more and more rapidly in time, making
considerable noise as it passed over the pipe. As the chain speeded up, but
still in contact with the pipe, the curvature of the ‘VP’, that is, of the top
of the chain, decreased. About $2s$ after the video was started, the chain
lifted off the pipe, and the VP rose rapidly and accelerated upward, while its
curvature decreased significantly. We estimate that the in the first second
after liftoff, VP rose $1m$, in the second second it rose $2m$, and continued
to rise (off camera) for another $2$ seconds. At that time, we heard a sound,
which was consistent with the chain hitting the ceiling at the top of the
stairwell. We found that the chain had indeed hit the ceiling, leaving marks
and embedded pieces of chain in the fiber ceiling tiles. Subsequent
experiments produced similar results. We estimate that the VP in the video
rose in excess of $7m$ to hit the ceiling. Other experiments using 4” PVC pipe
for support, and $1/4"$ ball-chains we observed even greater heights. We
attach stills from our video in Fig. 8, which shows the chain supported by the
pipe, lifting off the pipe, the VP accelerating as it rises to the ceiling of
the stairwell. The last frame shows the chain falling back.

Figure 8: Sequence of stills from the video showing the chain dynamics. The
times are: a=0s,b=2s,c=3s,d=4s,e=5s,f=6s,g=7s and h=8s from start of the
video.

##  4 Summary

We have considered essential aspects of a lossless whip; these are the
concentration of the energy and momentum densities into a single point at the
tip, and the subsequent spreading out of these again throughout the chain. In
the absence of losses, this process repeats indefinitely. We have also argued
that if a straight chain at rest is folded in half, and if then one end is
pulled away from the fold, the VP, other end will move towards it.

We have considered a chain with free ends hanging from a pulley which is
supported from below. If the system is unbalanced, the chain speed grows
nearly exponentially, and for a given range of parameters, at a given time,
the chain lifts off the pulley. After this time, the dynamics of the chain is
the same as that of a whip; although the center of mass of the chain is in
free fall. Due to the whip dynamics, one end of the chain snaps upward, and if
the chain is lossless, this snapping of alternating ends is expected to
continue indefinitely.

We have carried out experiments where long slender ball-chain is supported by
a horizontal metal rod, so that one half of the chain hangs on one side, and
the other half hangs on the other side. If the chain is then positioned so
that there is more chain on one side than the other and the chain is brought
into motion, the chain will begin to slip with increasing speed, with the
longer part of the chain moving downward and the shorter part moving upward.
At a critical speed, the chain lifts off the supporting bar, with the curved
region where the chain changes direction, that is, with the ‘VP’, accelerating
rapidly upward. In our experiments, the top part of the chain rose 20’ above
the bar and came into collision with the ceiling tiles. Although the center of
mass of the chain is freely falling after lift-off, the chain exhibits whip
dynamics with the upward moving part snapping; in our case, the tip of the
chain disintegrated into numerous small fragments. The initiation of lift-off
of a chain from a pulley has been observed and examined by [8]; however, they
imposed a constant acceleration on the chain, and only a very small fraction
of the chain (less than the pulley circumference), at the end of the chain,
lifted off the pulley. Here, we see fully developed liftoff, far from the end
of the chain, with the upward and downward moving parts of the chain nearly
parallel. In addition to exhibiting fascinating dynamics, this mechanism may
be useful for imparting large upward velocities to objects at the end of the
chain, and thus for launching projectiles.

In this paper, using 1D dynamics, we have assumed that the curvature is a
constant at the fold, which is not the case, as experiments show. This will be
a topic for a future study.

A video showing the dynamics of the falling chain, together with additional
information, is available at http://mpalffy.kent.edu/chain.

## Acknowledgments

We are grateful for illuminating discussions with M. Warner (University of
Cambridge) and H. Yokoyama (Kent State University), and to R. Bammerlin and
Kent State Libraries for access to the building and stairwell. This work was
supported by the US Office of Naval Research [ONR 00014-18-1-2624]. Results
from this paper has been presented at the Minisymposium MS 125 “Liquid
Crystals, Elastomers and Beyond: In Memory of Mark Warner” at 2022 Annual
Meeting of Society for Industrial and Applied Mathematics.

## References

  * [1] Mould S., Investigating the “Mould effect”, TEDxNewcastle, (2014). https://www.youtube.com/watch?v=wmFi1xhz9OQ
  * [2] Biggins J. S., Warner M., Understanding the chain fountain, (2014). https://www.youtube.com/watch?v=-eEi7fO0_O0&t=472s
  * [3] Blundell J., Machado H., Fink T., Rope siphon, (2013). http://www.youtube.com/watch?v=X7CXzjFVUHQ
  * [4] Biggins, J.S. and M. Warner, Understanding the chain fountain. Proceedings of the Royal Society a-Mathematical Physical and Engineering Sciences, 2014. 470(2163).
  * [5] Virga, E.G., Dissipative shocks in a chain fountain. Physical Review E, 2014. 89(5).
  * [6] Virga, E.G., 2015. Chain paradoxes. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2173), p.20140657.
  * [7] Zheng, X. and Palffy-Muhoray, P., 2018. Acquiring momentum: simple strategies by simple objects. New Journal of Physics, 20(8), p.081001.
  * [8] Brun, P.T., Audoly, B., Goriely, A. and Vella, D., 2016. The surprising dynamics of a chain on a pulley: lift off and snapping. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472(2190), p.20160187.
  * [9] Yokoyama, H., 2018. Reexamining the chain fountain. arXiv preprint arXiv:1810.13008.
  * [10] Anghel, D.V., 2021, November. The theory of the chain fountain revisited. In Journal of Physics: Conference Series (Vol. 2090, No. 1, p. 012166). IOP Publishing.
  * [11] Corbin, N.A., Hanna, J.A., Royston, W.R., Singh, H. and Warner, R.B., 2018. Impact-induced acceleration by obstacles. New Journal of Physics, 20(5), p.053031.


backbone: An R package to extract network backbones

Zachary P. Neal1*

1 Psychology Department, Michigan State University, East Lansing, MI, USA

*<EMAIL_ADDRESS>
## Abstract

Networks are useful for representing phenomena in a broad range of domains.
Although their ability to represent complexity can be a virtue, it is
sometimes useful to focus on a simplified network that contains only the most
important edges: the backbone. This paper introduces and demonstrates a
substantially expanded version of the backbone package for R, which now
provides methods for extracting backbones from weighted networks, weighted
bipartite projections, and unweighted networks. For each type of network,
fully replicable code is presented first for small toy examples, then for
complete empirical examples using transportation, political, and social
networks. The paper also demonstrates the implications of several issues of
statistical inference that arise in backbone extraction. It concludes by
briefly reviewing existing applications of backbone extraction using the
backbone package, and future directions for research on network backbone
extraction.

## Introduction

Networks are useful for representing phenomena in a broad range of domains [1,
2]. Although their ability to represent complexity can be a virtue, in some
cases (e.g., computationally intensive analysis, presence of noise,
visualization) it is useful to simplify a network and focus instead on its
backbone. Given a complex network $\mathbf{N}$, which may be weighted or
unweighted, its backbone $\mathbf{N^{\prime}}$ is a sparse and unweighted
subgraph that aims to preserve or reveal important structural features. Many
backbone methods have been proposed for extracting a backbone
$\mathbf{N^{\prime}}$ from a network $\mathbf{N}$, with different methods
designed for different types of networks, or to preserve different types of
structural features. However, applying these methods in practice has been
challenging because they were either not implemented or implemented in
different software languages.

The backbone package for R aims to overcome this practical limitation of
network backbone extraction by providing an integrated implementation of
existing methods. The formal mathematical details of these methods, and
evidence of their performance as structure-preserving or structure-revealing
backbone models, is extensively documented elsewhere and referenced below.
Instead, the purpose of this paper is to provide a practical guide to using
the backbone package for backbone extraction. It is organized in six sections.
The first section provides an overview of the backbone package. The second,
third, and fourth sections illustrate how the package can be used to extract
the backbone from a weighted network, from a bipartite projection, and from an
unweighted network, respectively. Each of these sections begin with an
overview of the relevant backbone models, then provide a small toy example to
illustrate the model, followed by an empirical example to demonstrate its use
in practice. The fifth section reviews issues of statistical inference that
arise in backbone models. Finally, the sixth section concludes with a
discussion of past applications of the backbone package, its limitations, and
future directions for the implementation of backbone models.

## The backbone package

The most recent release of the backbone package can be installed in R [3] from
the Comprehensive R Archive Network (CRAN) using:

    > install.packages("backbone")
    trying URL ’https://cran.rstudio.com...’
    Content type ’application/x-gzip’ length 1758582 bytes (1.7 MB)
    ==================================================
    downloaded 1.7 MB

Once installed, the backbone package can be loaded using:

    > library(backbone)
     ____   backbone v2.1.0
    |  _ \  Cite: Neal, Z. P., (2022). backbone: An R package to extract network
    |#|_) |       backbones. PLOS ONE.
    |# _ <
    |#|_) | Help: type vignette("backbone"); email<EMAIL_ADDRESS>github zpneal/backbone
    |____/  Beta: type devtools::install_github("zpneal/backbone", ref = "devel")

The startup message displays the version of the backbone package that is
installed and has been loaded for use. It also displays the recommended
citation for the package, sources for help using the package, and the command
to install the beta release of the backbone package. Additional information
about the CRAN distribution is available at
https://CRAN.R-project.org/package=backbone, while additional materials
relating to backbone, including papers, presentations, workshop materials, and
data sets are available at http://www.zacharyneal.com/backbone. The code and
data necessary to replicate the examples shown in this paper are available at
https://osf.io/8tuc7/.

Fig 1 illustrates the typical workflow of the backbone package. A user begins
with source data, which can take the form of an R matrix object, sparse Matrix
object, an edgelist stored as a dataframe object, or an igraph object. The
source data may represent an unweighted bipartite network, a weighted
unipartite network, or an unweighted unipartite network. The relevant backbone
models depend on the type of network. For example, when starting with an
unweighted bipartite network, a user may extract the backbone from its
weighted projection using any of the following models: sdsm(), fdsm(),
fixedrow(), fixedcol(), fixedfill(), disparity(), or global(). Extracting the
backbone using one of these functions yields an unweighted unipartite network,
thus as the figure illustrates, the user may subsequently extract a second-
order backbone using the sparsify() function.

Fig 1: Workflow of the backbone package.

The backbone.suggest() function can be used to examine the source data,
identify the relevant backbone models, and suggest the most appropriate model.
For example:

    > N <- matrix(runif(100), 10, 10)
    > backbone.suggest(N)
    The disparity filter is suggested. Type "?disparity" for more information.

Here, the source data is a randomly generated $10\times 10$ matrix of values
between 0 and 1. Based on these characteristics, it is recognized as a
weighted and directed adjacency matrix representing a weighted unipartite
network. Therefore, the disparity filter model is recommended, for which
additional information is available by typing ?disparity in the R console.

In the sections below, we use randomly generated data like this for toy
examples, but also use real data for empirical examples. The empirical
examples use data on air transportation in the US [4], bill sponsorship in the
US Senate [5], and friendships among faculty at a UK university [6]. A random
number seed can be set to ensure that the randomly generated data are
replicable, and the empirical data can be loaded, using:

    > set.seed(1)
    > load("backbone2_tutorial.Rdata")

## Backbones of weighted networks

### Background

In a weighted network, each edge has a weight that captures the strength of
the relationship between the two nodes it connects. These weights can
represent capacities when larger values capture stronger relations, or can
represent costs when larger values capture weaker relations. The precise
meaning of these weights depends on the context of the network, but the models
implemented in the backbone package assume they represent capacities (i.e.,
stronger relations have larger weights). In a social network, the weights
might represent the strength of a friendship, or the frequency of
communication between two people. Alternatively, in a transportation network
the weights might represent the number of people flying from one city to
another, or the capacity of a road to carry vehicles. These weights provide
some information about the importance of edges in the network, and backbone
models designed for weighted networks use this information to determine which
edges should be preserved in the backbone.

The simplest backbone model applicable to weighted networks is the global
threshold. The global threshold model preserves all edges whose weight exceeds
a specified threshold. It can be classified as a ‘global’ model because it
applies the same criteria to all edges in the network, and can be classified
as a ‘structural’ model because it uses only structural information (here, the
edges’ weights). While it can be appealing for its simplicity, the global
threshold model can be problematic when applied to networks in which the edge
weights follow a long-tailed distribution (e.g., mostly weak edges, a few
strong edges) and networks in which different parts of the network have
different characteristic edge weights (e.g., multi-scale networks).

Many alternative backbone models have been developed for weighted networks
exhibiting these characteristics [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19, 20, 21], however among the most widely used is the disparity filter [4].
The disparity filter compares an edge’s observed weight to its expected weight
in a null model where a node’s total weight is uniformly distributed across
its edges. For example, it might compare the observed number of people
traveling between New York and Los Angeles (a large number) to the number that
would be expected if each of New York’s passengers randomly selected a
destination (likely a much smaller number). An edge is preserved in the
backbone if its observed weight is statistically significantly (at a specified
$\alpha$ significance level) stronger than expected under such a null model.
This model can be classified as a ‘local’ model because an edge’s importance
is evaluated using information in its neighborhood, and can be classified as a
‘statistical’ model because it makes inferences about edge importance by
reference to a statistical null model.

### Toy Example

To illustrate extracting the global threshold and disparity filter backbones
from weighted networks using the backbone package, we begin with a simple toy
example:

    > mat <- matrix(c(0,10,10,10,10,75,0,0,0,0,
                      10,0,1,1,1,0,0,0,0,0,
                      10,1,0,1,1,0,0,0,0,0,
                      10,1,1,0,1,0,0,0,0,0,
                      10,1,1,1,0,0,0,0,0,0,
                      75,0,0,0,0,0,100,100,100,100,
                      0,0,0,0,0,100,0,10,10,10,
                      0,0,0,0,0,100,10,0,10,10,
                      0,0,0,0,0,100,10,10,0,10,
                      0,0,0,0,0,100,10,10,10,0),10)

This network is shown in the left panel of Fig 2, where stronger edges are
drawn using thicker lines. The edge weights in this network clearly follow a
long-tailed distribution: a few edges are very strong (weights $\geq 75$),
while the majority of edges are fairly weak (weights $\leq 10$). Additionally,
the network has a clear multi-scale structure: the cluster in the lower left
is characterized by relatively strong edges, while the cluster in the upper
right is characterized by relatively weak edges.

Fig 2: Extracting the backbone of a weighted network (toy example).

We can extract a global threshold backbone from this network using the
global() function:

    > unweighted <- global(mat, upper = function(x)mean(x), class = "igraph")

In this example, the global() function takes three arguments: the source
network, the upper threshold for edges to discard, and the desired class of
the result. Here, the source network takes the form of an adjacency matrix
mat. The upper threshold for edges to discard is defined as the mean of all
edge weights, thereby yielding a backbone of stronger-than-average edges.
Finally, the result is returned as an igraph object, to facilitate subsequent
visualization and analysis. The global threshold backbone is shown in the
middle panel of Fig 2, which highlights its shortcomings when applied to a
multiscale network. By applying the same threshold to all edges, it preserves
edges in the high-weight cluster, but ignores the edges in the low-weight
cluster.

We can extract a disparity filter backbone from this network using the
disparity() function:

    > backbone <- disparity(mat, alpha = 0.05, class = "igraph")

The disparity() function also takes the adjacency matrix mat as input and
returns the backbone as an igraph object, however it does not use a global
threshold value. Instead, the user specifies alpha, the statistical
significance level. In this example, an edge is preserved if its weight is
larger than the weight of the corresponding edge it at least 95% of null model
networks. The disparity backbone is shown in the right panel of Fig 2, which
illustrates its ability to preserve the original network’s hub-and-spoke
structure despite its multi-scale edge weights.

### Empirical Example

The toy example illustrates the basic operation of the global() and
disparity() backbone extraction functions. To illustrate how these backbone
extraction functions might be used in practice, we apply them to extract the
backbone of the US air transportation network, which is known to have a hub-
and-spoke structure containing many low-degree nodes (e.g., most airports
serve few passengers) and a few very high-degree nodes (e.g., hub airports
that serve many passengers). These data closely resemble the airline data used
in the initial demonstration of the disparity filter [4], and were obtained
from the US Bureau of Transportation Statistics’ Airline Origin and
Destination Survey, which contains a 10% random sample of all domestic airline
tickets in 2019. The weighted, symmetric network records the number of
passengers traveling between 382 airports in the continental US. We can
inspect the entries in this network:

    > airport["JFK","LAX"]  #Passengers between New York and Los Angeles
    [1] 326756
    > airport["LAN","GRR"]  #Passengers between Lansing and Grand Rapids
    [1] 94
    > airport["LAN","LAX"]  #Passengers between Lansing and Los Angeles
    [1] 0

For example, we observe that many passengers fly between JFK (New York) and
LAX (Los Angeles), while very few fly between LAN (Lansing, Michigan) and GRR
(Grand Rapids, Michigan) and none fly between LAN and LAX because there is no
scheduled service.

The weighted original network is shown in the left panel of Fig 3, which
highlights that the network’s high density obscures any particular structure.
The network has a very large mean weighted degree $\langle k\rangle=374215$
because some airports such as Atlanta’s Hartsfield-Jackson serve very large
volumes of passengers. When a network has a hub-and-spoke structure, the
probability of observing a node with degree $k$ is approximately distributed
as $k^{\gamma}$, where typically $2<\gamma<3$ [22]. However, in this network
$\gamma=1.15$ [23], which falls outside this range and suggests the network
lacks the expected hub-and-spoke structure.

Fig 3: Extracting the backbone of a weighted network of continental US airline
traffic in 2019.

We can extract a global threshold backbone from this network using:

    > unweighted <- global(airport, upper = function(x)mean(x), class = "igraph")

The resulting backbone, which preserves only stronger-than-average edges, is
shown in the middle panel of Fig 3. It is clearly much sparser, with a smaller
mean degree $\langle k\rangle=14.81$. However, it remains fairly dense, and is
dominated by the high-volume routes on the East coast. Its degree scaling
exponent $\gamma=1.83$ remains inconsistent with a hub-and-spoke
transportation network [22, 23].

We can extract a disparity filter backbone from this network using:

    > backbone <- disparity(airport, alpha = 0.001, class = "igraph", narrative = TRUE)

    === Suggested manuscript text and citations ===
    We used the backbone package for R (v2.1.0; Neal, 2022) to extract the unweighted
    backbone of a weighted and undirected unipartite network containing 382 nodes.
    An edge was retained in the backbone if its weight was statistically significant
    (alpha = 0.001) using the disparity filter (Serrano et al., 2009). This reduced
    the number of edges by 91.4%, and reduced the number of connected nodes by 14.7%.

    Neal, Z. P. (2022). backbone: An R Package to Extract Network Backbones.
    arXiv:2203.11055 [cs.SI]. https://doi.org/10.48550/arXiv.2203.11055

    Serrano, M. A., Boguna, M., & Vespignani, A. (2009). Extracting the multiscale
    backbone of complex weighted networks. Proceedings of the National Academy of
    Sciences, 106(16), 6483-6488. https://doi.org/10.1073/pnas.0808904106

In this example, following [4] we use a conservative statistical significance
level of $\alpha=0.001$. We also include the narrative = TRUE argument, which
automatically generates text and citations describing what the disparity()
function has done, and that can be used or adapted for a manuscript. The
resulting backbone, which preserves edges whose weights are unusually large
compared to the null model, is shown in the right panel of Fig 3. It is quite
sparse: the narrative text indicates that the number of edges was reduced by
91.4%, and we observe it has a small mean degree $\langle k\rangle=4.34$. Its
degree scaling exponent $\gamma=2.34$ is now consistent with a transportation
infrastructure known to have a hub-and-spoke organization [22, 23]. Indeed,
the visualization clearly illustrates this organization, which is anchored by
the hub airports of ORD (Chicago), DEN (Denver), DFW (Dallas-Fort Worth), and
ATL (Atlanta).

## Backbones of bipartite projections

### Background

A bipartite projection is a special type of weighted network in which the edge
weights capture the number of artifacts shared in common by two agents, and
thus capture co-occurrence relations [24]. They arise in a wide range of
contexts, including co-author networks where authors (the agents) share papers
(the artifacts), co-sponsorship networks where legislators share bills, co-
attendance networks where people share events, and co-expression networks
where genes share expressed proteins. As a variety of weighted network, the
backbone of a bipartite projection can be extracted using global() or
disparity(). However, it is usually more appropriate to extract the backbone
of a bipartite projection using a backbone model designed specifically for
this purpose.

To understand why a specialized backbone model is needed for a bipartite
projection, it is helpful to consider how the edge weights in a bipartite
projection are defined. A bipartite network in which authors (the agents) are
connected to the papers (the artifacts) they have written can be represented
by an incidence matrix B, where where $B_{ik}=1$ if author $i$ wrote paper
$k$. We obtain the adjacency matrix of the bipartite projection $\mathbf{P}$
from the incidence matrix of the bipartite network $\mathbf{B}$ as
$\mathbf{P}=\mathbf{BB}^{\prime}$, where each entry $P_{ij}$ indicates the
number of papers co-authored by authors $i$ and $j$ (i.e., the edge weight).
Applying a conventional backbone model such as the disparity filter to this
bipartite projection would decide which edges to preserve based only on these
edge weights, but would fail to consider two important pieces of information
that are contained in the original bipartite network. First, the row sums of
$\mathbf{B}$ record how many papers the $i^{th}$ author wrote (i.e., the agent
degrees). This is important to consider because if two authors each wrote many
papers, it is likely they would have co-authored a few just by chance. Second,
the column sums of $\mathbf{B}$ record how many authors the $k^{th}$ paper has
(i.e., the artifact degrees). This is important to consider because if a paper
has many authors, then observing that two people are both on the author list
is not particularly noteworthy.

Backbone models designed specifically for bipartite projections are unique
because they are applied directly to the bipartite network, not to its
projection, so that they can incorporate this information. While there are
multiple such models [25, 26, 27, 24, 28], they all use a common approach: An
edge’s observed weight in the projection is compared to its expected weight in
a projection obtained from a random bipartite network generated by a null
model. An edge is preserved in the backbone if its observed weight is
statistically significantly (at a specified $\alpha$ significance level)
stronger than expected under the null model. These models differ in the
constraints they impose on the random bipartite networks. In this section, we
illustrate the Stochastic Degree Sequence Model (SDSM), which constrains the
random bipartite networks to have row sums and column sums that match those of
the observed bipartite network on average. The SDSM is faster than the Fixed
Degree Sequence Model, and more accurately recovers known structural
characteristics than the Fixed Row, Fixed Column, or Fixed Fill models [28].

### Toy Example

To illustrate extracting the SDSM backbone from a bipartite projection using
the backbone package, we begin with a simple toy example:

    > B <- rbind(cbind(matrix(rbinom(250,1,.8),10),
                       matrix(rbinom(250,1,.2),10),
                       matrix(rbinom(250,1,.2),10)),
                 cbind(matrix(rbinom(250,1,.2),10),
                       matrix(rbinom(250,1,.8),10),
                       matrix(rbinom(250,1,.2),10)),
                 cbind(matrix(rbinom(250,1,.2),10),
                       matrix(rbinom(250,1,.2),10),
                       matrix(rbinom(250,1,.8),10)))
    > sum(B[1,])  #Agent 1’s degree
    [1] 28
    > sum(B[,1])  #Artifact 1’s degree
    [1] 11

This generates a random bipartite network, represented as an incidence matrix
in B. This network consists of 30 agents and 75 artifacts, and is embedded
with a community structure. The agents are split into three groups of 10, and
the artifacts are split into three groups of 25. There is an 80% chance of an
edge between an agent and artifact belonging to the same group (i.e., people
are likely to attend their own group’s events), and only a 20% chance of an
edge between an agent and artifact belonging to different groups (i.e., people
are unlikely to attend another group’s events). Examining this bipartite
network, we observe that Agent 1 is associated with 28 artifacts (i.e., agent
degree), and that Artifact 1 is associated with 11 agents (i.e., artifact
degree).

Given this bipartite network, we can construct and examine its bipartite
projection:

    > P <- B %*% t(B)  #Construct bipartite projection
    > P[1,2]  #Edge weight for two agents in the same group
    [1] 17
    > P[1,20]  #Edge weight for two agents in different groups
    [1] 7
    > min(P)  #Smallest edge weight
    [1] 3

The projection is constructed by multiplying the bipartite network B by its
transpose t(B). In the projection, we observe that agents 1 and 2 shared 17
artifacts in common, a large number because they are members of the same
group. In contrast, agents 1 and 20 share only 7 artifacts in common, a small
number because they are members of different groups. Notably, however, the
smallest number of artifacts shared by any pair of agents is still 3. This
means that all agents share at least some artifacts in common, which makes the
bipartite projection not particularly useful as a network, and occurs because
the construction of a bipartite projection “induces an inflation in the number
of links” [29]. Indeed, as illustrated in the left panel of Fig 4, the
weighted bipartite projection is so dense that no particular structure is
visible, including the community structure that is known to exist.

Fig 4: Extracting the backbone of a bipartite projection (toy example).

Because the bipartite projection is a weighted network, we could extract its
backbone using the disparity filter:

    > disparity <- disparity(P, alpha = 0.05, class = "igraph")
    This object looks like it could be a bipartite projection. If so, consider extracting
    the backbone using a model designed for bipartite projections: sdsm, fdsm, fixedfill,
    fixedrow, or fixedcol.

The middle panel of Fig 4 shows the disparity backbone, which is empty. That
is, the disparity filter is unable to identify any significant edges to
preserve in the backbone. This occurs because the disparity filter is relying
only on information contained in the bipartite projection, but does not use
any of the information contained in the underlying bipartite network. Here,
the disparity() function detects that the source network appears to be a
bipartite projection, and recommends instead using a backbone model
specifically designed for such a network.

Following this advice, we can instead extract a backbone using the stochastic
degree sequence model and sdsm() function:

    > backbone <- sdsm(B, alpha = 0.05, class = "igraph")

Like the disparity() function, the sdsm() function also takes alpha and class
arguments. However, the source network is not the bipartite projection P, but
the original bipartite network B from which it was constructed. The right
panel of Fig 4 shows the SDSM backbone, which is sparse and clearly displays
the known three-community structure.

### Empirical Example

To illustrate the extraction of a bipartite projection backbone in practice,
we use data on bill sponsorship patterns in the US Senate’s $115^{th}$ session
(2017-2018) [30, 31]. In the US Senate, legislators can express support for a
bill by ‘sponsoring’ it. Therefore, bipartite projections of bill sponsorship
data yield co-sponsorship networks that are often used to study political
alliances among legislators. The bipartite network’s incidence matrix records
the bill sponsorships of 105 Senators on 3665 bills. We can inspect the
details of this bipartite network:

    > senate[1:2,1:2]  #First two rows and columns
                                 S.1 S.100
    Sen. Alexander, Lamar [R-TN]   0     1
    Sen. Baldwin, Tammy [D-WI]     0     0
    > sum(senate["Sen. Stabenow, Debbie [D-MI]",])  #Sponsorships by Sen. Stabenow
    [1] 316
    > sum(senate[,"S.1006"])  #Sponsors of Equality Act
    [1] 48

Looking at the first four entries of the incidence matrix, we see that neither
Senators Alexander nor Baldwin sponsored Senate Bill 1, and that Senator
Alexander sponsored Senate Bill 100 while Senator Baldwin did not. Row sums
indicate each Senator’s number of sponsorships, for example, Senator Stabenow
sponsored 316 bills in this session. Similarly, column sums indicate each
bill’s number of sponsors, for example, Senate Bill 1006 had 48 sponsors. This
bill, known as the Equality Act, would have amended the Civil Rights Act of
1964 to prohibit discrimination by sex, sexual orientation, and gender
identity, and was sponsored only by Democrats and Independents.

Given this bipartite network, we can construct and examine its bipartite
projection:

    > P <- senate %*% t(senate) #Construct bipartite projection
    > P["Sen. Stabenow, Debbie [D-MI]", "Sen. Peters, Gary C. [D-MI]"]
    [1] 151
    > P["Sen. Stabenow, Debbie [D-MI]", "Sen. Cruz, Ted [R-TX]"]
    [1] 14

We see that the Senators Stabenow and Peters co-sponsored 151 bills together,
reflecting their likely alliance as members of the same party and
representatives of the same state. In contrast, we see that Senators Stabenow
and Cruz sponsored only 14 bills together, reflecting their lack of
partnership as members of opposing parties and representatives of different
states. Notably, however, despite the sharp ideological differences between
Stabenow and Cruz, they still co-sponsored some bills together and are still
connected in the bipartite projection.

The weighted bipartite projection is shown in the left panel of Fig 5;
Democrats and Independents are shown in blue, while Republicans are shown in
red. It illustrates several problems with relying on a raw projection. First,
the network has a high density and large weighted mean degree $\langle
k\rangle=4688.48$ that obscures any structure. Second, the visualization fails
to capture the known partisan polarization of the US Senate [32, 30, 5], which
is confirmed by the network’s small modularity ($Q=0.18$) with respect to
political party affiliation. Finally, the network includes as connected two
nodes representing Senators Jeff Sessions (R-AL) and Jon Kyl (R-AZ). Although
these senators did co-sponsor some bills with their Republican colleagues,
they co-sponsored very few because their terms were unusually short: Sessions
served only 36 days before being appointed US Attorney General, while Kyl
served only 117 days to fill a vacancy left by the death of Senator John
McCain.

Fig 5: Extracting the backbone of a bipartite projection of bill co-
sponsorship in the 115th US Senate.

Because this bipartite projection is a weighted network, we could extract its
backbone using the disparity filter:

    > disparity <- disparity(P, alpha = 0.1, class = "igraph")
    This object looks like it could be a bipartite projection. If so, consider extracting
    the backbone using a model designed for bipartite projections: sdsm, fdsm, fixedfill,
    fixedrow, or fixedcol.

The disparity backbone is shown in the middle panel of Fig 5. Even using a
relatively liberal significance level ($\alpha=0.1$), this backbone is so
sparse ($\langle k\rangle=3.64$) that the majority of Senators are
disconnected. Additionally, although there are small clusters of Democrats and
Republicans, its modest modularity ($Q=0.24$) highlights that it fails to
capture the known partisan polarization. These problems arise because the
disparity filter is applied directly to the bipartite projection, and is
unable to incorporate the more detailed information contained in the original
bipartite network.

Finally, we can extract a backbone using the stochastic degree sequence model:

    > backbone <- sdsm(senate, alpha = 0.05, class = "igraph", narrative = TRUE)

    === Suggested manuscript text and citations ===
    We used the backbone package for R (v2.1.0; Neal, 2022) to extract the unweighted backbone
    of the weighted projection of an unweighted bipartite network containing 105 agents and
    3665 artifacts. An edge was retained in the backbone if its weight was statistically
    significant (alpha = 0.05) using the stochastic degree sequence model (SDSM; Neal,
    2014). This reduced the number of edges by 69.7%, and reduced the number of connected
    nodes by 1.9%.

    Neal, Z. P. (2022). backbone: An R Package to Extract Network Backbones. arXiv:
    2203.11055 [cs.SI]. https://doi.org/10.48550/arXiv.2203.11055

    Neal, Z. P. (2014). The backbone of bipartite projections: Inferring relationships
    from co-authorship, co-sponsorship, co-attendance and other co-behaviors. Social
    Networks, 39, 84-97. https://doi.org/10.1016/j.socnet.2014.06.001

By including the narrative = TRUE argument, the function generates text and
citations describing the backbone extraction. The SDSM backbone is shown in
the right panel of Fig 5. It is sparse ($\langle k\rangle=30.7$), but not as
sparse as the disparity backbone. Accordingly, with the desirable exceptions
of Sessions and Kyl, all the Senators are connected in a single component.
Moreover, this connected component clearly illustrates the known partisan
polarization of the US Senate, which is confirmed by the network’s high
modularity ($Q=0.47$).

## Backbones of unweighted networks

### Background

Backbone extraction from weighted networks and weighted bipartite projections
is facilitated by the availability of information about edge weights. A
different approach is required for extracting the backbone from unweighted
networks because there are no edge weights. Many backbone models for
unweighted networks have been proposed [33, 34, 35, 36, 37, 38], but they all
follow a common process: score, normalize, filter, connect:

  1. 1.

Because there are no weights associated with the edges, each edge is assigned
a score. These scores can be random [38], but more often are topologically-
derived. For example, an edge may be assigned a score that counts the number
of triangles [36] or quadrangles [33] it completes, or that captures the
amount of overlap in the neighborhoods of the two nodes it connects [37, 35].

  2. 2.

These edge scores can then be normalized. The normalization step is omitted by
some backbone models [38, 35], while other backbone models normalize through
rankings [37, 34] or more sophisticated transformations [36, 33].

  3. 3.

The (optionally normalized) edge scored are filtered according to a user-
specified sparsification parameter. Different backbone models interpret this
sparsification parameter in different ways. For example, this parameter may
specify the fraction of strongest scoring edges to retain [38], the threshold
edge score above which edges are retained [35, 36, 33], or a non-linear
parameter controlling the stringency of edge retention [37, 34].

  4. 4.

Because unweighted backbone models can frequently yield disconnected networks,
some models ensure the connectedness of the backbone by also including edges
that appear in the union of minimum spanning trees [33].

Each of these steps in the process represent a decision point in the
specification of a backbone model. The sparsify() function in the backbone
package allows users to explicitly specify how (or whether) each step is
performed, and therefore to customize a unique backbone model. However, some
combinations represent named backbone models that have been investigated in
the literature. In this paper, we focus on two such models that are fast and
often yield informative backbones: Local Sparsification (L-Spar) [37] and
Local Degree (LD) [34].

The L-Spar model [37] scores edges using the Jaccard coefficient, which
measures the amount of overlap in the neighbors of the two nodes it connects,
ranging from 0 (no overlap) to 1 (complete overlap). Here, the intuition is
that relationships are stronger between individuals who share many of the same
contacts, which makes this model most appropriate for extracting backbones
that preserve clustering structures. Next, the edge scores are normalized by
ranking them from the perspective of each node, such that a node’s top scoring
edge is ranked highest. Finally, for a node with degree $d$, the $d^{s}$
highest-ranked edges are preserved in the backbone, where $s$ is the
sparsification parameter ranging from 0 to 1. When $s=0$, which is the
smallest value yielding the sparsest backbone, every node’s one strongest edge
is retained because $d^{0}=1$ for all $d$. Alternatively, when $s=1$, which is
the largest value yielding the densest backbone, all of a node’s edges are
retained because $d^{1}=d$ for all $d$. Intermediate values of $s$ retain more
or fewer edges, but non-linearly with the goal of removing edges of “nodes
with higher degree more aggressively than nodes with lower degree” [37].

The LD model [34] scores edges, from the perspective of each node, using the
degree of the node at the other end. Here, the intuition is that the most
important edges are those that lead to hubs, which makes this model most
appropriate for extracting backbones that preserve branching, hub-and-spoke,
or hierarchical structures. After this scoring step, the LD model proceeds
identically to the L-Spar model: edges are normalized by rank, then the
$d^{s}$ highest-ranking edges for each node are preserved in the backbone.

### Toy Examples

To illustrate extracting backbones from unweighted networks using the backbone
package, we begin with two separate toy examples:

L-Spar and communities. The L-Spar backbone model is designed to extract
backbones that preserve hidden community structures. To illustrate, we begin
by generating a random unweighted network embedded with a hidden community
structure:

    > pref.matrix <- matrix(c(.75,.25,.25,.25,.75,.25,.25,.25,.75),3,3)
    > unweighted <- sbm.game(60, pref.matrix, c(20,20,20))

The sbm.game() function is from the igraph [39] package, and is used to
generate stochastic block models. In this case, it yields a 60-node
unweighted, undirected network composed of three 20-node communities, such
that there is a 75% chance of a within-community edge and a 25% change of a
between-community edge. The left panel of Fig 6 shows the resulting network,
which is so dense that the community structure is obscured.

Fig 6: Extracting the backbone of an unweighted network with embedded
communities (toy example).

We can extract the L-Spar backbone using either the flexible sparsify()
function or the simpler sparsify.with.lspar() function:

    > backbone <- sparsify(unweighted, escore = "jaccard", normalize = "rank",
                           filter = "degree", umst = FALSE, s = 0.5)
    > backbone <- sparsify.with.lspar(unweighted, s = 0.5)

The sparsify() function is highly customizable by allowing the user to specify
how to perform each of the four steps in an unweighted backbone model, while
the sparsify.with.lspar() function is a simplified wrapper that performs these
steps according to the L-Spar model. In either case, the sparsification
parameter s specifies the stringency of the backbone model in preserving
edges. The right panel of Fig 6 shows the L-Spar backbone, which clearly
reveals the known three-community structure of this network.

LD and hubs. The LD backbone model is designed to extract backbones that
preserve hidden hub-and-spoke or hierarchical structures. To illustrate, we
begin by generating a random unweighted netwokr embedded with hidden hubs:

    > unweighted <- as.undirected(sample_pa(60, m = 3), mode = "collapse")

The sample_pa() function is from the igraph [39] package, and is used to
generate random networks formed via preferential attachment. In this case, it
yields a 60-node unweighted, undirected network in which a few nodes occupy
high-degree hub positions. The left panel of Fig 7 shows the resulting
network, which is so dense that the hubs are obscured.

Fig 7: Extracting the backbone of an unweighted network with embedded
hub/spoke structure (toy example).

We can extract the L-Spar backbone using either the flexible sparsify()
function or the simpler sparsify.with.localdegree() function:

    > backbone <- sparsify(unweighted, escore = "degree", normalize = "rank",
                           filter = "degree", umst = FALSE, s = 0.1)
    > backbone <- sparsify.with.localdegree(unweighted, s = 0.1)

Here, the sparsify.with.localdegree() function is a wrapper for the more
flexible sparsify() function that performs each of the steps according to the
local degree model. The right panel of Fig 7 shows the LD backbone, which
clearly reveals the network’s hierarchical structure around a few high-degree
hubs.

### Empirical Examples

To illustrate the extraction of backbones from unweighted networks in
practice, we use two separate empirical examples:

L-Spar and friendship. To illustrate the use the L-Spar model to reveal a
hidden community structure, we use data on the undirected and unweighted
friendships among 79 faculty at a UK university [6]. Each individual is a
member of a single school, which provides an expectation for a community
structure. The left panel of Fig 8 shows the original network, with nodes
colored according to their school membership. For a small network, it is
relatively dense, with a mean degree $\langle k\rangle=14.25$, which can
hamper the identification of communities. In this case, the communities are
relatively obvious only because the nodes are already colored according to
their known community memberships. However, communities detected by applying a
fast-greedy algorithm to this network yields only a moderate recovery of the
true memberships as measured by the Adjusted Rand Index ($ARI_{fg}=0.76$)
[40].

Fig 8: Extracting the backbone of an unweighted network of faculty email
exchanges.

We can extract the network’s L-Spar backbone using:

    > backbone <- sparsify.with.lspar(faculty, s = 0.5, class = "igraph", narrative = TRUE)

    === Suggested manuscript text and citations ===
    We used the backbone package for R (v2.1.0; Neal, 2022) to extract the unweighted backbone
    of an unweighted and undirected unipartite network containing 79 nodes. Specifically,
    we used Satuluri et al.’s (2011) L-Spar model with a sparsification threshold of 0.5.
    This reduced the number of edges by 62%, and reduced the number of connected nodes by 0%.

    Neal, Z. P. (2022). backbone: An R Package to Extract Network Backbones. arXiv:
    2203.11055 [cs.SI]. https://doi.org/10.48550/arXiv.2203.11055

    Satuluri, V., Parthasarathy, S., & Ruan, Y. (2011, June). Local graph sparsification
    for scalable clustering. In Proceedings of the 2011 ACM SIGMOD International Conference
    on Management of data (pp. 721-732). https://doi.org/10.1145/1989323.1989399

By including the narrative = TRUE argument, the function generates text and
citations describing the backbone extraction. The resulting backbone is shown
in the right panel of Fig 8. It is substantially less dense, with a much
smaller mean degree $\langle k\rangle=5.32$, which more clearly reveals the
three-community structure. Indeed, detecting communities by applying a fast-
greedy algorithm to the backbone yields an excellent recovery of the true
memberships ($ARI_{fg}=0.92$)[40].

LD and air traffic. To illustrate the use the LD model to reveal hidden hubs,
we return to the US air traffic data discussed earlier. From the original
weighted network, we focus on the unweighted version:

    > unweighted <- simplify(graph_from_adjacency_matrix(airport, mode = "undirected",
                                                         diag = FALSE))
    > unweighted["JFK","LAX"]  #Are there flights between New York and Los Angeles?
    [1] 1
    > unweighted["LAN","LAX"]  #Are there flights between Lansing and Los Angeles
    [1] 0
    > sum(unweighted["LAN",])  #How many destinations are reachable from Lansing?
    [1] 40

The simplify() function from the igraph[39] package transforms the originally
weighted airline network into an unweighted network in which two airports are
connected if any number of passagers flew between them. Examining this
unweighted network, we observe that there are flights between JFK (New York)
and LAX (Los Angeles), but not between LAN (Lansing) and LAX (Los Angeles),
and indeed that it is possible to reach only 40 airports from LAN. The left
panel of Fig 9 shows this unweighted network, which is so dense and has such a
high mean degree ($\langle k\rangle=50.67$) that the underlying structure of
the US airline infrastructure is obscured. Additionally, its degree scaling
exponent $\gamma=6.84$ is inconsistent with the transportation network’s
expected in a hub-and-spoke transportation network [22, 23].

Fig 9: Extracting the backbone of an unweighted network of continental US
airline traffic in 2019.

We can extract the network’s LD backbone using:

    > backbone <- sparsify.with.localdegree(unweighted, s = 0.3, narrative = TRUE)

    === Suggested manuscript text and citations ===
    We used the backbone package for R (v2.1.0; Neal, 2022) to extract the unweighted backbone
    of an unweighted and undirected unipartite network containing 382 nodes. Specifically,
    we used Hamann et al.’s (2016) local degree model with a sparsification threshold of
    0.3. This reduced the number of edges by 90.7%, and reduced the number of connected
    nodes by 0%.

    Neal, Z. P. (2022). backbone: An R Package to Extract Network Backbones. arXiv:
    2203.11055 [cs.SI]. https://doi.org/10.48550/arXiv.2203.11055

    Hamann, M., Lindner, G., Meyerhenke, H., Staudt, C. L., & Wagner, D. (2016). Structure-
    preserving sparsification methods for social networks. Social Network Analysis and
    Mining, 6(1), 22. https:://10.1007/s13278-016-0332-2

By including the narrative = TRUE argument, the function generates text and
citations describing the backbone extraction. The resulting backbone is shown
in the right panel of Fig 9. It is much sparser, with a smaller mean degree
$\langle k\rangle=4.71$, which makes it easier to see the hub-and-spoke
structure of the network anchored by the hub airports of ORD (Chicago), DEN
(Denver), DFW (Dallas-Fort Worth), and ATL (Atlanta). Additionally, its degree
scaling exponent $\gamma=1.94$ is consistent with a transportation
infrastructure known to have a hub-and-spoke organization [22, 23].

## Issues in statistical inference

Some of the backbone models implemented in the backbone package are
‘structural’ models that choose edges to retain based on their topological
properties. For example, the global threshold model makes such decisions on
the basis of edge weights, while the L-Spar model makes such decisions on the
basis of an edge’s Jaccard coefficient, which functions as an imputed edge
weight. However, other models are ‘statistical’ in the sense that they choose
edges to retain by considering their probability relative to a statistical
null model. These statistical models compute a $p$-value for each edge, where
the $p$-value captures the probability of observing a larger weight associated
with the same edge in a random network. The statistical models are
identifiable from the presence of an alpha parameter; an edge is retained if
its $p$-value is smaller than the specified $\alpha$-value. This section
briefly reviews several issues related to statistical backbone models’
$p$-values.

### Exact versus approximate $p$-values

Most statistical models can compute edges’ $p$-values exactly using analytic
methods [28]. However, the fixed degree sequence model used for extracting the
backbone from bipartite projections and implemented in the fdsm() requires
estimating these $p$-values using Monte Carlo simulation. This function uses
efficient methods for simulating random networks under the FDSM null model
[41], however it is still necessary to perform a large number of simulations
to estimate $p$-values with sufficient confidence that they can be compared to
the $\alpha$-value. The trials parameter of the fdsm() function allows the
user to explicitly specify how many simulations are used to estimate the
$p$-value. When the trials parameter is not specified, the function
automatically determines the required number of simulations to estimate
$p$-value with sufficient confidence. Using the US Senate data, for example:

    > backbone <- fdsm(senate, alpha = 0.05, trials = 1000)
    Constructing empirical edgewise p-values using 1000 trials -
    |============================================================| 100%
    > backbone <- fdsm(senate, alpha = 0.05)
    Constructing empirical edgewise p-values using 80297 trials -
    |============================================================| 100%

When trials = 1000 is specified, the $p$-values are estimated from 1000 random
bipartite networks, which takes only a few seconds but yields potentially
unstable results. In contrast, when the trials parameter is omitted, fdsm()
automatically determines that testing the $p$-values against a significance
level of $\alpha=0.05$ will require 80,297 random bipartite networks, which
takes several minutes. This highlights why the FDSM backbone model is often
impractical for extracting the backbone from large bipartite projections, or
at conservative significance levels (e.g., $\alpha<0.05$).

### Obtaining $p$-values

Typically statistical backbone models evaluate edges’ p-values and return an
unweighted backbone that only contains the edges deemed statistically
significant. However, specifying alpha = NULL instead returns the $p$-values
themselves, rather than the backbone they imply. For example, in the case of
the US airport network:

    > bb.object <- disparity(airport, alpha = NULL)
    This matrix object is being treated as a weighted undirected network containing 382 nodes.

    > bb.object$G[1:3,1:3]
          ALB   ATL  AVP
    ALB     0 24349   13
    ATL 24349     0 5761
    AVP    13  5761    0

    > bb.object$Pupper[1:3,1:3]
                  ALB           ATL         AVP
    ALB 1.00000000000 0.00006619769 0.986723244
    ATL 0.00006619769 1.00000000000 0.001924381
    AVP 0.98672324370 0.00192438055 1.000000000

    > bb.object$Plower[1:3,1:3]
               ALB       ATL        AVP
    ALB 1.00000000 0.9999338 0.01327676
    ATL 0.99993380 1.0000000 0.99807562
    AVP 0.01327676 0.9980756 1.00000000

    > bb.object$model
    [1] "disparity"

By specifying alpha = NULL, the disparity() function returns a backbone object
that contains four items. First, it contains the original weighted network;
here we see that 24,349 passengers flew between Atlanta (ATL) and Albuquerque
(ALB). Second, it contains the upper-tail $p$-values, which capture the
probability of observing an edge weight as large or larger in a random
network. For example, we see that observing so many passengers flying between
ATL and ALB is very unlikely in a random network, which offers evidence that
the ATL-ALB edge is significant. Third, it contains the lower-tail $p$-values,
which capture the probability of observing an edge weight as small or smaller
in a random network. For example, we see that observing so few passengers
flying between ALB and Scranton (AVP) is very unlikely in a random network.
Finally, it contains a string indicating the backbone model that generated
these $p$-values.

Returning a backbone object containing $p$-values can be useful in cases where
computing the $p$-values is slow (e.g., using fdsm()), and backbones extracted
using multiple $\alpha$ significance levels are desired. In such cases,
backbone extraction is a two-step process: first generate a backbone object
containing the $p$-values, then use backbone.extract to extract a backbone at
the desired $\alpha$ significance level. For example,

    > bb.object <- disparity(airport, alpha = NULL) #Compute p-values
    > bb1 <- backbone.extract(bb.object, alpha = 0.001) #Backbone at alpha = 0.001
    > bb2 <- backbone.extract(bb.object, alpha = 0.005) #Backbone at alpha = 0.005

The first line computes the $p$-values for each edge in the airport network
using the disparity filter model. Then, the second line extracts the backbone
using a significance level of $\alpha=0.001$. The third line extracts a
different backbone using the more liberal significance level of
$\alpha=0.005$, which can be performed without needing to re-compute the
$p$-values.

### Correcting $p$-values

The chosen $\alpha$ significance level defines the Type-I error rate, which in
the backbone extraction context is the probability of deciding an edge has a
larger weight than would be expected in a random network, when in fact it does
not (i.e., a false positive). In practice, backbone extraction requires
evaluating every edge with a non-zero weight, which inflates the Type-I error
rate. For example, if $\alpha=0.05$, then the probability of incorrectly
deciding a single edge is significant is 0.05. However, the probability of
making at least one such incorrect decision across $m$ edges is
$1-(1-0.05)^{m}$, which can dramatically inflate the risk of false positives
even when extracting relatively small backbones. Performing a Multiple Test
Correction (MTC), which involves adjusting the computed $p$-values, is
necessary to control for this error inflation. The statistical backbone models
implemented in the backbone package can perform any of the MTC corrections
that are provided by R’s p.adjust() function via the mtc parameter.

We can use the empirical airport data to examine the effect of correcting for
multiple tests:

    > backbone_nomtc <- disparity(airport, alpha = 0.001, class = "igraph", mtc = "none")
    > backbone_mtc <- disparity(airport, alpha = 0.001, class = "igraph", mtc = "holm")

The left panel of Fig 10 shows the disparity filter backbone extracted from
the weighted network using an edgewise error rate of $\alpha=0.001$. That is,
each edge’s $p$-value is compared to 0.001 and retained if it is smaller. This
yields a backbone in which there is only a 0.001% chance that any given
retained edge is, in fact, not significant. However, because 9678 edges’
significance is tested, there is a $1-(1-0.001)^{9678}\approx 99.99\%$ chance
that at least one retained edge is not significant. The right panel of Fig 10
shows the disparity filter backbone extracted using a familywise error rate of
$\alpha=0.001$, achieved by applying a Holm-Bonferroni correction to the
$p$-values [42]. This yields a much sparser backbone in which there is only a
0.001% chance that one or more retained edge is, in fact, not significant.
Although the MTC-corrected backbone is much sparser, it still displays the
expected hub-and-spoke structure.

Fig 10: Disparity filter backbone of a weighted network of continental US
airline traffic in 2019 at a edgewise error rate (left) or familywise error
rate (right) of $\alpha=0.001$.

### Signed backbones

In most cases, backbone extraction involves identifying and retaining the
statistically significantly strong edges in the backbone. However, statistical
backbone models are also capable of identifying and retaining the
statistically significantly weak edges. Specifying signed = TRUE yields a
signed backbone in which significantly strong edges are retained in the
backbone as positive edges, while significantly weak edges are retained in the
backbone as negative edges. For example, returning to the US Senate data:

    > signed <- sdsm(senate, alpha = 0.05, signed = TRUE)

Fig 11 shows the signed backbone, with positive edges drawn in green and
negative edges drawn in red. This example illustrates that, as expected,
positive relations of cooperation exist primarily between legislators from the
same party, while negative relations of opposition exist primarily between
legislators from different parties [5].

Fig 11: Extracting the signed backbone of a bipartite projection of bill co-
sponsorship in the 115th US Senate.

Extracting an unweighted (i.e., not signed) backbone involves a one-tailed
significance test, that is, testing whether an edge’s observed weight is
larger than the weight of the corresponding edge in random networks generated
under the null model. In a one-tailed significance test, and edge is deemed
significant if $p<\alpha$. However, extracting a signed backbone involves a
two-tailed significance test because it requires evaluating whether an edge’s
observe weight it larger or smaller than expected in a random network. In such
two-tailed significance tests, an edge is deemed significant and negative if
its lower-tail $p<\frac{\alpha}{2}$, and is deemed significant and positive if
its upper-tail $p<\frac{\alpha}{2}$.

## Conclusion

### Applications

The models implemented in the backbone package are generic and can be used to
extract the backbone of networks observed in many different domains. To date,
the majority of applications have been in political, developmental, and
biological sciences. In political science, the backbone package has been used
to infer networks of political alliances among legislators from a bipartite
projection of their bill sponsorships [32, 30, 31, 43, 44, 45]. While two
legislators may be viewed as having an alliance when they are observed to
sponsor many of the same bills (i.e., the edge weight in a co-sponsorship
network), backbone models offer a way to identify pairs of legislators that
have sponsored significantly more bills together than expected. In
developmental science, the backbone package has been used to measure
childrens’ social networks from their peer-reported or observed group
memberships [46, 47]. This is useful because data on childrens’ social
relations are difficult to collect directly, but they can often be inferred
from bipartite projections. In the biological sciences, the backbone package
has been used to examine organism co-occurrences in both plants [48] and
animals [49], and is used by the genetonic package [50] to extract the
backbone of gene-geneset graphs. Beyond these scientific applications, the
backbone package has also been used to examine the collaboration relations
among the writers, pencilers, and editors responsible for the X-men comic book
series [51], and to recommend holiday movies [52]. These applications
highlight the flexibility of the backbone extraction functions implemented in
the backbone package and offer additional illustrations of its functionality.

### Limitations and future work

The backbone package remains under development, and the models it implements
are relatively new, which means that both the package and models are subject
to some limitations that highlight opportunities for future work. First, many
backbone models exist for both weighted and unweighted [7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 53, 54, 55] networks that are not yet implemented in the
backbone package. However, the backbone package provides a generic framework
for the implementation of additional models in future releases. Second,
relatively little research has focused on validating backbone models or
establishing their scope conditions, which leaves researchers with limited
guidance on model selection. The backbone.suggest() function already provides
rudimentary model selection assistance. However, by implementing multiple
models in a single framework, the backbone package can facilitate future work
comparing backbone models on data with a ground truth (for validation) or with
different characteristics (to establish scope conditions). Finally, because
the goal of backbone extraction is simplification, backbone models are most
useful in very large networks, however certain backbone models computationally
intensive. Thus, future work on backbone models and their implementations
should focus on scalability and efficiency.

## References

  * 1. Newman M.  Networks.  Oxford university press; 2018.
  * 2. Barabási AL.  Network Science.  Cambridge university press; 2016.
  * 3. R Core Team. R: A Language and Environment for Statistical Computing; 2018.  Available from: https://www.R-project.org.
  * 4. Serrano MÁ, Boguná M, Vespignani A.  Extracting the multiscale backbone of complex weighted networks.  Proceedings of the National Academy of Sciences. 2009;106(16):6483–6488.  doi:10.1073/pnas.0808904106.
  * 5. Neal ZP.  A sign of the times? Weak and strong polarization in the US Congress, 1973–2016.  Social Networks. 2020;60:103–112.  doi:10.1016/j.socnet.2018.07.007.
  * 6. Nepusz T, Petróczi A, Négyessy L, Bazsó F.  Fuzzy communities and the concept of bridgeness in complex networks.  Physical Review E. 2008;77(1):016107.  doi:10.1103/PhysRevE.77.016107.
  * 7. Tumminello M, Aste T, Di Matteo T, Mantegna RN.  A tool for filtering information in complex systems.  Proceedings of the National Academy of Sciences. 2005;102(30):10421–10426.  doi:10.1073/pnas.0500298102.
  * 8. Glattfelder JB, Battiston S.  Backbone of complex networks of corporations: The flow of control.  Physical Review E. 2009;80(3):036104.  doi:10.1103/PhysRevE.80.036104.
  * 9. Foti NJ, Hughes JM, Rockmore DN.  Nonparametric sparsification of complex multiscale networks.  PloS one. 2011;6(2):e16431.  doi:10.1371/journal.pone.0016431.
  * 10. Zhou F, Malher S, Toivonen H.  Network simplification with minimal loss of connectivity.  In: 2010 IEEE international conference on data mining. IEEE; 2010. p. 659–668.
  * 11. Dianati N.  Unwinding the hairball graph: Pruning algorithms for weighted complex networks.  Physical Review E. 2016;93(1):012304.  doi:10.1103/PhysRevE.93.012304.
  * 12. Zhang X, Zhang Z, Zhao H, Wang Q, Zhu J.  Extracting the Globally and Locally Adaptive Backbone of Complex Networks.  PLOS ONE. 2014;9(6):1–8.  doi:10.1371/journal.pone.0100428.
  * 13. Zhang RJ, Stanley HE, Ye FY.  Extracting h-Backbone as a Core Structure in Weighted Networks.  Scientific Reports. 2018;8:14356.  doi:10.1038/s41598-018-32430-1.
  * 14. Gursoy F, Badur B.  Extracting the signed backbone of intrinsically dense weighted networks.  Journal of Complex Networks. 2021;9(5):cnab019.  doi:10.1093/comnet/cnab019.
  * 15. Radicchi F, Ramasco JJ, Fortunato S.  Information filtering in complex weighted networks.  Phys Rev E. 2011;83:046101.  doi:10.1103/PhysRevE.83.046101.
  * 16. Rajeh S, Savonnet M, Leclercq E, Cherifi H.  Modularity-Based Backbone Extraction in Weighted Complex Networks.  In: Ribeiro P, Silva F, Mendes JF, Laureano R, editors. Network Science. Cham: Springer International Publishing; 2022. p. 67–79.
  * 17. Marcaccioli R, Livan G.  A Polya urn approach to information filtering in complex networks.  Nature Communications. 2019;10:745.  doi:10.1038/s41467-019-08667-3.
  * 18. Slater PB.  A two-stage algorithm for extracting the multiscale backbone of complex weighted networks.  Proceedings of the National Academy of Sciences. 2009;106(26):E66–E66.  doi:10.1073/pnas.0904725106.
  * 19. Gemmetto V, Cardillo A, Garlaschelli D.  Irreducible network backbones: unbiased graph filtering via maximum entropy.  arXiv preprint arXiv:170600230. 2017;doi:10.48550/arXiv.1706.00230.
  * 20. Coscia M, Neffke FM.  Network backboning with noisy data.  In: 2017 IEEE 33rd International Conference on Data Engineering (ICDE). IEEE; 2017. p. 425–436.
  * 21. Šubelj L.  Convex skeletons of complex networks.  Journal of The Royal Society Interface. 2018;15(145):20180422.  doi:10.1098/rsif.2018.0422.
  * 22. Barabási AL, Albert R.  Emergence of scaling in random networks.  science. 1999;286(5439):509–512.  doi:10.1126/science.286.5439.509.
  * 23. Clauset A, Shalizi CR, Newman ME.  Power-law distributions in empirical data.  SIAM review. 2009;51(4):661–703.  doi:10.1137/070710111.
  * 24. Neal ZP.  The backbone of bipartite projections: Inferring relationships from co-authorship, co-sponsorship, co-attendance and other co-behaviors.  Social Networks. 2014;39:84–97.  doi:10.1016/j.socnet.2014.06.001.
  * 25. Zweig KA, Kaufmann M.  A systematic approach to the one-mode projection of bipartite graphs.  Social Network Analysis and Mining. 2011;1(3):187–218.  doi:10.1007/s13278-011-0021-0.
  * 26. Saracco F, Di Clemente R, Gabrielli A, Squartini T.  Randomizing bipartite networks: the case of the World Trade Web.  Scientific reports. 2015;5(1):1–18.  doi:10.1038/srep10595.
  * 27. Tumminello M, Micciche S, Lillo F, Piilo J, Mantegna RN.  Statistically validated networks in bipartite complex systems.  PloS one. 2011;6(3):e17994.  doi:10.1371/journal.pone.0017994.
  * 28. Neal ZP, Domagalski R, Sagan B.  Comparing alternatives to the fixed degree sequence model for extracting the backbone of bipartite projections.  Scientific reports. 2021;11(1):1–13.  doi:10.1038/s41598-021-03238-3.
  * 29. Latapy M, Magnien C, Del Vecchio N.  Basic notions for the analysis of large two-mode networks.  Social networks. 2008;30(1):31–48.
  * 30. Aref S, Neal ZP.  Identifying hidden coalitions in the US House of Representatives by optimally partitioning signed networks based on generalized balance.  Scientific reports. 2021;11(1):1–9.  doi:10.1038/s41598-021-98139-w.
  * 31. Neal ZP, Domagalski R, Yan X.  Homophily in collaborations among US House Representatives, 1981–2018.  Social Networks. 2022;68:97–106.  doi:10.1016/j.socnet.2021.04.007.
  * 32. Aref S, Neal ZP.  Detecting coalitions by optimally partitioning signed networks of political collaboration.  Scientific reports. 2020;10(1):1–10.  doi:10.1038/s41598-020-58471-z.
  * 33. Nocaj A, Ortmann M, Brandes U.  Untangling the hairballs of multi-centered, small-world online social media networks.  Journal of Graph Algorithms and Applications: JGAA. 2015;19(2):595–618.  doi:10.7155/jgaa.00370.
  * 34. Hamann M, Lindner G, Meyerhenke H, Staudt CL, Wagner D.  Structure-preserving sparsification methods for social networks.  Social Network Analysis and Mining. 2016;6(1):1–22.  doi:10.1007/s13278-016-0332-2.
  * 35. Goldberg DS, Roth FP.  Assessing experimentally derived interactions in a small world.  Proceedings of the National Academy of Sciences. 2003;100(8):4372–4376.  doi:10.1073/pnas.0735871100.
  * 36. Nick B, Lee C, Cunningham P, Brandes U.  Simmelian backbones: Amplifying hidden homophily in facebook networks.  In: Proceedings of the 2013 IEEE/ACM international conference on advances in social networks analysis and mining; 2013. p. 525–532.
  * 37. Satuluri V, Parthasarathy S, Ruan Y.  Local graph sparsification for scalable clustering.  In: Proceedings of the 2011 ACM SIGMOD International Conference on Management of data; 2011. p. 721–732.
  * 38. Karger DR.  Random sampling in cut, flow, and network design problems.  Mathematics of Operations Research. 1999;24(2):383–413.  doi:10.1287/moor.24.2.383.
  * 39. Csardi G, Nepusz T. The igraph software package for complex network research; 2006.  Available from: http://igraph.org.
  * 40. Hubert L, Arabie P.  Comparing partitions.  Journal of classification. 1985;2(1):193–218.  doi:10.1007/BF01908075.
  * 41. Godard K, Neal ZP.  fastball: A fast algorithm to sample binary matrices with fixed marginals.  arXiv preprint arXiv:211204017. 2021;.
  * 42. Holm S.  A simple sequentially rejective multiple test procedure.  Scandinavian Journal of Statistics. 1979; p. 65–70.
  * 43. Schoch D, Brandes U.  Legislators’ roll-call voting behavior increasingly corresponds to intervals in the political spectrum.  Scientific reports. 2020;10(1):17369.  doi:10.1038/s41598-020-74175-w.
  * 44. Chen THY.  Informal Diplomacy Reinforces Formal International Cooperation: Evidence from Track Two Diplomacy.  SocArXiv. 2021;.
  * 45. de la Bruere J.  Uncovering Legislative Networks: A Socio-Economic Perspective.  Western Illinois University [PhD Thesis]. 2021;.
  * 46. Neal JW, Neal ZP, Durbin CE.  Inferring Signed Networks from Preschoolers’ Observed Parallel and Social Play.  PsyArXiv. 2022;.
  * 47. Neal ZP, Neal JW, Domagalski R.  False positives using social cognitive mapping to identify children’s peer groups.  Collabra: Psychology. 2021;7(1):17969.  doi:10.1525/collabra.17969.
  * 48. Custer GF, van Diepen LT, Stump W.  An Examination of Fungal and Bacterial Assemblages in Bulk and Rhizosphere Soils under Solanum tuberosum in Southeastern Wyoming, USA.  Applied Microbiology. 2021;1(2):162–176.  doi:10.3390/applmicrobiol1020013.
  * 49. Buerger AN, Dillon DT, Schmidt J, Yang T, Zubcevic J, Martyniuk CJ, et al.  Gastrointestinal dysbiosis following diethylhexyl phthalate exposure in zebrafish (Danio rerio): Altered microbial diversity, functionality, and network connectivity.  Environmental Pollution. 2020;265:114496.  doi:10.1016/j.envpol.2020.114496.
  * 50. Marini F, Ludt A, Linke J, Strauch K.  GeneTonic: an R/Bioconductor package for streamlining the interpretation of RNA-seq data.  BMC bioinformatics. 2021;22(1):610.  doi:10.1186/s12859-021-04461-5.
  * 51. Benton R. XMEN; 2022.  Available from: http://www.richardabenton.com/xmen.html.
  * 52. Neal ZP, Neal JW. Christmas Movie Recommendation Network; 2022.  Available from: https://twitter.com/zpneal/status/1474038138347610118.
  * 53. John E, Safro I.  Single-and multi-level network sparsification by algebraic distance.  Journal of Complex Networks. 2017;5(3):352–388.  doi:10.1093/comnet/cnw025.
  * 54. Auber D, Chiricota Y, Jourdan F, Melançon G.  Multiscale visualization of small world networks.  In: IEEE Symposium on Information Visualization 2003 (IEEE Cat. No. 03TH8714). IEEE; 2003. p. 75–81.
  * 55. Ghalmane Z, Cherifi C, Cherifi H, El Hassouni M.  Extracting modular-based backbones in weighted networks.  Information Sciences. 2021;576:454–474.  doi:10.1016/j.ins.2021.06.087.