We use topological data analysis to study “functional networks” that we construct from timeseries data from both experimental and synthetic sources. We use persistent homology with a weight rank clique filtration to gain insights into these functional networks, and we use persistence landscapes to interpret our results. Our first example uses timeseries output from networks of coupled Kuramoto oscillators. Our second example consists of biological data in the form of functional magnetic resonance imaging data that were acquired from human subjects during a simple motorlearning task in which subjects were monitored for three days during a fiveday period. With these examples, we demonstrate that (1) using persistent homology to study functional networks provides fascinating insights into their properties and (2) the position of the features in a filtration can sometimes play a more vital role than persistence in the interpretation of topological features, even though conventionally the latter is used to distinguish between signal and noise. We find that persistent homology can detect differences in synchronization patterns in our data sets over time, giving insight both on changes in community structure in the networks and on increased synchronization between brain regions that form loops in a functional network during motor learning. For the motorlearning data, persistence landscapes also reveal that on average the majority of changes in the network loops take place on the second of the three days of the learning process.
Computational topology is a family of methods that are based on topological ideas (e.g., they often arise from algebraic topology) and give insights into topological invariants, such as connectedness or holes in highdimensional data sets.^{1–3} Such efforts have come to be called topological data analysis, and a method known as persistent homology (PH) has been particularly helpful for understanding shapes and their persistence over multiple scales.^{4,5} Traditionally, PH has been applied to pointcloud data, though it has also been applied to networks in many applications, ranging from granular materials (see, e.g., Ref. 6) to functional brain networks.^{7,8} We employ these topological tools, which are designed to yield global, “higherorder” insights that go beyond revelations from pairwise connections (which are the norm in network science), in a study of functional networks constructed from both empirical and synthetic timeseries data. We use persistence landscapes to show that the topological tools can (1) capture dynamics of networks constructed from the data and (2) identify mesoscale features that we relate to community structure^{9,10} in the associated functional networks. To help readers optimally understand these insights, we also present an intuitive introduction to PH and how to apply it to networks.
I. INTRODUCTION
The human brain consists of many billions of neurons, whose major task is to receive, conduct, and transmit signals. Analysis of neuronal networks is crucial for understanding the human brain.^{11–16} Every neuron consists of a cell body and one long axon, which is responsible for propagating signals to other cells.^{17} Neurons or (on a larger scale) different brain regions can be construed as nodes of a network, whose edges represent either structural or functional connections between those nodes. Examining neuronal data using a networkbased approach allows one to use mathematical tools from subjects such as graph theory to better understand structural and functional aspects of neuronal interactions, identify key regions in the brain that are involved in physiological and pathological processes, and compare the structure of neuronal interactions with those of other complex systems. For example, data analysis using network theory has led to the insight that the brain has underlying modular structures, with small subunits that are able to carry out specific functions while minimally influencing other parts of the brain.^{11,14,18}
The standard methods from network theory are based on pairwise connections, which one can use to study the microscale, mesoscale, and macroscale structures.^{19} An alternative approach for studying networks^{20} is to use methods from computational topology, which explicitly incorporates “higherorder” structures beyond pairwise connections and includes algorithmic methods for understanding topological invariants, such as connectedness, loops, or holes in highdimensional data structures^{1–3}(see Section II B). Although one can also represent higherorder structures using formalisms such as hypergraphs^{21} (see, e.g., a recent paper^{22} by Bassett et al.), those other approaches may not be the most convenient means for optimally conveying information about the shape or scale of mesoscale structures in a network. Other recent work concerns clustering in networks using higherorder structures.^{23}
Methods from computational topology enable one to understand global lowdimensional structures in networks, and they have led to insights in an increasingly large number of applications^{5} in diverse topics, ranging from granular materials^{6} and contagions on networks^{24} to path planning^{25} and collective behavior in animals.^{26} In particular, persistent homology (PH), a mathematical formalism to explore the persistence of topological structures in data sets, has become increasingly prominent in neuroscience in the last few years.^{7,8} Among other applications, it has been used to determine differences in brain networks of children with hyperactivity disorders and autism spectrum compared to normal situations,^{27} study the effect of the psychoactive component of “magic mushrooms” (psilocybin mushrooms) on functional brain networks of humans,^{28} analyze covariates that influence neural spiketrain data,^{29} and study structural and functional organization of neural microcircuits.^{30} Other neuronal applications have included consideration of place cells in the hippocampus of rats during spatial navigation,^{31–33} analysis of mathematical models of transient hippocampal networks,^{34} and a demonstration that topological features of networks of brain arteries in humans are correlated with their age.^{35} We also note that PH is not the only topological method that has been used to study the human brain or time series. More than fifty years ago, for example, Zeeman^{36} used tolerance spaces and Vietoris homology theory to study aspects of visual perception. In the 1990s, Muldoon et al.^{37} developed a method to study the topology of manifolds that underlie timeseries data.
In the present investigation, we use timeseries data to construct socalled functional networks^{11,12,38,39} (but note that one can study coupled time series using a variety of different approaches^{40–42}). Functional brain networks consist of a set of nodes (e.g., brain regions) and a set of weighted edges between nodes, where the edge weights quantify the similarity of the associated time series according to a chosen measure. A functional network contrasts with a “structural network,” which refers to the underlying physical connections (e.g., anatomical connections) between nodes. For example, neurons are connected to each other in structural networks, but one can analyze the similarity in their firing patterns through functional networks. We use the term “functional network” in a more general way: by constructing a matrix of similarities between coupled time series using some measure (and enforcing the diagonal entries to be 0), one obtains a functional network whose weighted adjacency matrix (sometimes also called an “association matrix”) $ A \u0303 = ( a \u0303 i j ) i , j = 1 N $ has elements that indicate the similarity between the time series of entities i and j. Studying functional networks is common in neuroscience, and they are also used in a wealth of other applications (e.g., finance,^{43} voting among legislators,^{44} and climate^{45}). Importantly, the time series can come either from empirical data or from the output of a dynamical system (or stochastic process), and the latter can be helpful for validating methods for network analysis.^{46} In our paper, we will consider time series either from coupled oscillators (i.e., as the output of a dynamical system) or from a set of spatially distinct brain regions defined by a fixed anatomical atlas. In the context of functional brain networks, the adjacencymatrix element a_{ij} arises as a measure of “functional connectivity” (i.e., behavioral similarity) between the time series for nodes (i.e., brain regions) i and j. There are many different ways to measure similarity of time series,^{12,47,48} and that can be a major issue when it comes to interpreting results. Comparing the networks that arise from different similarity measures is beyond the scope of our work, so we will simply use two common measures (pairwise synchrony and wavelet coherence) of timeseries similarity. However, the methods that we employ can be applied to functional networks that are constructed using any measure of similarity between time series.
In many studies based on experimental data, functional networks are used to construct binary graphs (i.e., unweighted graphs).^{12} To do this, one typically applies a global threshold $ \xi \u2208 \mathbb{R} + $ to a weighted adjacency matrix to obtain a binary adjacency matrix $ A = ( a i j ) i , j = 1 N $ associated with an unweighted graph. The adjacencymatrix elements are then
The choice of threshold has a strong influence on the resulting matrix, and it thereby exerts a major influence on the structure of the associated graph.^{12} Some approaches to address this issue include determining a single “optimal” threshold, thresholding the weighted adjacency matrix at different values,^{49,50} examining the network properties as a function of threshold, or not thresholding at all and considering the weighted adjacency matrix itself.^{12,14} (One can also threshold a weighted adjacency matrix by setting sufficiently small entries to 0 but keeping the values of the other entries.) If one is thresholding and binarizing data, there is no guarantee that there exists an interval of thresholds that yield networks with qualitatively similar properties, and arbitrarily throwing away data can be problematic even when such intervals do exist. For example, parameters such as graph size (i.e., number of nodes) need to be taken into account when interpreting results on thresholded, binarized networks.^{51} An advantage of using persistent homology is that one can examine a graph “filtration” (see Section II C) generated by multiple—ideally all—possible global thresholds and systematically analyze the persistence of topological features across these thresholds. Such a filtration can also be created using decreasing local thresholds.
In our topological analysis, we focus on “loops” in a network. For our purposes, a loop in a graph is a set of at least four edges that are connected in a way that forms a topological circle.^{52} Loops are thus 1dimensional topological features. We choose to focus on loops rather than features with dimension 0, which correspond to connected components of a graph, are topologically simpler, and can be studied using many other approaches (e.g., through the number of 0 elements in the spectrum of the combinatorial graph Laplacian^{21}). It has been demonstrated in other applications (e.g., contagions on networks^{24}) that loops are important topological features of graphs, and a recent study^{53} demonstrated the importance of loops (and related higherdimensional objects) in structural neuronal networks. Structural and functional neuronal networks are related and share some network features,^{11} so we expect loops to provide interesting insights.
The remainder of our paper is organized as follows. In Section II, we give a brief and intuitive introduction to persistent homology, the weight rank clique filtration (WRCF), and persistence landscapes. In Section III, we introduce our first example, the Kuramoto model of nonlinearly coupled oscillators; and we present results from our application of persistent homology to timeseries data produced by coupled Kuramoto oscillators. In Section IV, we introduce and analyze our second example, which consists of timeseries functional magnetic resonance imaging (fMRI) data from experiments of humans performing a simple motorlearning task. We present our conclusions in Section V, and we provide a mathematical introduction to persistent homology in the Supplementary Material.
II. PERSISTENT HOMOLOGY
Persistent homology (PH)^{1–3} is a method from computational topology that quantifies global topological structures (e.g., connectedness and holes) in highdimensional data. One can think of PH as looking at the “shape” of data in a given dimension using a set of different lenses. Each lens conveys topological features inside data at a different resolution. One then construes structures that persist over a range of different lenses to represent one or more significant features of the data. Structures that are observed only through a small number of lenses are commonly construed as noise,^{54,55} especially in settings in which the data are sampled from a manifold. For empirical data, the relationship between low persistence of a feature and whether it constitutes noise in a data set rather than signal has not yet been verified statistically, but we will illustrate a situation in which some shortlived structures represent important features and possibly genuine geometrical (not just topological) features of data in Sections III and IV.
In this section, we provide an intuitive introduction to the mathematical concepts behind PH. In the Supplementary Material, we give a mathematically rigorous introduction (including precise definitions).
A. Simplicial complexes
One can study the properties of a topological space^{56,57} by partitioning it into smaller and topologically simpler pieces, which when reassembled include the same aggregate topological information as the original space. The most trivial topological space $ X = { \u2205 , x} $ consists of the empty set $\u2205$ and a single point x. If we want to simplify the description of the topological properties of X, we would simply choose a single node to represent it. However, a node or even a collection of nodes does not allow one to capture the topological properties of more complicated spaces, such as a 2sphere or the surface of the earth. In such cases, one needs a simple object that carries the information that the space is connected but also encloses a hole. For example, one could use a tetrahedron, which is an example of a mathematical object called a simplex.
The building blocks that one uses to approximate topological spaces are called ksimplices, where the parameter k indicates the dimension of the simplex. Every ksimplex includes k + 1 independent nodes: a point is a 0simplex, an edge is a 1simplex, a triangle is a 2simplex, and a tetrahedron is a 3simplex (see Fig. 1). Observe that the lowerdimensional simplices are contained in the higherdimensional simplices. This allows one to build higherdimensional simplices using lowerdimensional ones. The lowerdimensional simplices form the socalled faces of the associated higherdimensional objects.
One combines different simplices into a simplicial complex to capture different aspects of a topological space. For every simplex that is part of a simplicial complex, we demand that all of its faces are also contained in the simplicial complex. Additionally, two simplices that are part of a simplicial complex are allowed to intersect only in common faces. In Fig. 2, we show several examples of simplicial complexes and one example that is not a simplicial complex.
We take the dimension of a simplicial complex to be the dimension of its highestdimensional simplex. One can use simplicial complexes to represent topological spaces if and only if there exists a continuous deformation that can stretch and bend the simplicial complex into the topological space, and only then are topological properties of the topological space preserved by the simplicial complex.
B. Homology and Betti numbers
If one is interested in the nature of a simplicial complex of dimension k, one can either consider the full complex, which can be very large, or one can examine different subsets of simplices that are contained in the complex. For example, the set of all 1simplices consists of a collection of edges, some of which may be connected or even form a loop. However, one can consider a range of different topological features—e.g., in some cases, a collection of edges surrounding a hole or void could be more interesting than individual edges—and one typically seeks features that are invariant if one squeezes or stretches the edges. Homology is a formal way to quantitatively detect topological invariants in a given dimension to give insight into the nature of a topological space. By using homology, one can, for example, distinguish a 2sphere from a torus. For a simplicial complex of dimension k, one can define a vector space known as the pth homology group for every dimension $ p \u2208 { 0 , \u2026 , k} $. In dimension 1, for example, the elements of the homology group are called “loops.” The elements of the homology group can be divided into different homology classes, which each represent a hole in the topological space. For instance, in dimension 1, loops in the same homology class all surround the same 1dimensional hole. We give an example of two loops that surround the same hole in Fig. 3. The homology classes yield a family of vector spaces, whose dimensions are called Betti numbers, associated to a simplicial complex. One can interpret the first three Betti numbers, β_{0}, β_{1}, and β_{2}, to represent, respectively, the number of connected components, the number of 1dimensional holes, and the number of 2dimensional holes in a simplicial complex. As we pointed out in Section I, we focus on the number of loops (i.e., on β_{1}) in our network analysis rather than on β_{0}, which corresponds to the number of connected components in a graph. One can study connected components in graphs using many other approaches, such as by calculating the number of 0 eigenvalues in the spectrum of the combinatorial graph Laplacian.^{21}
C. Filtrations
Although homology gives information about a single simplicial complex, it typically is more relevant to study the topological features across sequences (called filtrations) of simplicial complexes. A filtration^{2,54,55} of a simplicial complex Σ is a sequence of embedded simplicial complexes,
starting with the empty complex and ending with the entire simplicial complex. One can use homology to study topological features (e.g., loops) in every step of the filtration and determine how persistent they are with respect to a given filtration. A topological feature h is born at Σ_{m} if the homology group of Σ_{m} is the first homology group to include the feature. Similarly, a topological feature dies in Σ_{n} if it is present in the homology group of $ \Sigma n \u2212 1 $ but not in that of Σ_{n}. One then defines the persistence p of the topological feature as
Persistence was first used as a measure to rank topological features by their lifetime^{1} in a filtration in $ \mathbb{R} 3 $.
There are many ways to define simplicial complexes and filtrations on weighted graphs, and the choice of filtration tends to be motivated either by the type of questions to be answered or by consideration of computation time.
1. Weight rank clique filtration
Although we focus on network data, we note that PH has been applied much more often to data in the form of point clouds.^{54,55} The simplest way to create a sequence of embedded graphs (e.g., a filtration) from a weighted network is to filter by weights.^{58} To do this, one creates a sequence of embedded (binary) graphs by ranking all edge weights ν_{t} in descending order. In filtration step t, one retains an edge if and only if its weight is at least ν_{t}. To construct the filtration, one repeats this procedure until the graph is complete in the last step. Using this method, one is able to study 0simplices (i.e., nodes) and 1simplices (i.e., edges). The weight rank clique filtration (WRCF),^{59} which we will use in our analysis and which has been applied previously for examining weighted neuronal networks,^{28,33,59} extends this definition to include higherdimensional simplices. One constructs a WRCF as follows:

Define filtration step 0 as the set of all nodes.

Rank all edge weights $ { \nu 1 , \u2026 , \nu end } $, with $ \nu 1 = \nu max $ and $ \nu end = \nu min $. (We will use τ to denote the number of distinct weights in a graph.)

In filtration step t, threshold the graph at weight ν_{t} to create a binary graph.

Find all maximal ccliques for $ c \u2208 \mathbb{N} $, and define them to be csimplices.
This is a valid simplicial complex: every $ ( c + 1 ) $clique in the graph guarantees the existence of a cface on that clique, because cliques are closed under both intersection and taking subsets. Consequently, they satisfy the requirements for a simplicial complex. This type of simplicial complex on a graph is called a clique complex.
One can visualize the persistence of homology classes of a filtration of a simplicial complex using barcodes.^{54} A barcode for a given dimension is a collection $ { b l , d l } l = 1 m $ of interval endpoints, where every interval (b_{l}, d_{l}) represents a topological feature l of the given dimension (examples of such features include connected components and loops), b_{l} denotes the birth time of feature l with respect to the filtration step, and d_{l} denotes its death time. The length $ d l \u2212 b l $ of the bar measures the persistence of the feature. In Fig. 4, we show an example of a WRCF and its corresponding barcode.
D. Persistence landscapes
As an alternative topological summary to barcodes, one can use persistence landscapes,^{60,61} which consist of piecewiselinear functions in a separable Banach space. For a given barcode interval (b, d), one defines the function
For a barcode $ { b l , d l } l = 1 m $ and $ q \u2265 0 $, the qth persistence landscape is given by the set of functions
If the qthlargest value does not exist, then $ \lambda q ( x ) = 0 $. One can think of the 0th persistence landscape as being the outline of the collection of peaks created by the images of the collection of functions f associated to a barcode. To obtain the 1st persistence landscape, one peels away this topmost “layer” of peaks and then considers the outline of the remaining collection of peaks. This gives the 1st persistence landscape, and one continues in this manner to obtain subsequent persistence landscapes. The persistence landscape λ of the barcode $ { b l , d l } l = 1 m $ is then defined as the sequence $ { \lambda q } $ of functions λ_{q} (where λ_{q} is the qth layer of the landscape).
Even though persistence landscapes visualize the same information as barcodes for an individual filtration and one can construct a bijective correspondence between the two objects, the former have distinct advantages over the latter. For example, one can calculate a unique “average landscape” for a set of persistence landscapes from several filtrations by taking the mean over the function values for every landscape layer (see Fig. 5). This is not possible for barcodes, as they are not elements of a Banach space. For an average landscape, it is thus not possible to find a corresponding average barcode. We show a schematic illustration of how to obtain an average persistence landscape in Fig. 5.
One can also define an L^{p} distance between two landscapes (either individual ones or average ones) and thereby use a variety of statistical tools.^{60} This allows one to compare multiple groups of landscapes (and thereby also compare groups of barcodes) by calculating a measure of pairwise similarity between landscapes. Persistence landscapes have been used to study conformational changes in protein binding sites,^{62} phase separation in binary metal alloys,^{63} and music audio signals.^{64}
E. Computational tools
For our PH calculations, we use Matlab code that we construct using Javaplex,^{65,66} a software package for persistent homology. For the WRCFs, we also use a maximal cliquefinding algorithm from the Mathworks library^{67} that is based on the Bron–Kerbosch algorithm, which is the most efficient algorithm known for this problem. For statistical analysis and interpretation of our barcodes, we apply the Persistence Landscapes Toolbox.^{61}
III. EXAMPLE I: COUPLED KURAMOTO OSCILLATORS
A. The Kuramoto model
The Kuramoto model^{68–72} is a wellstudied model for a set of coupled phase oscillators whose natural frequencies are drawn from a prescribed distribution. The model was developed in the 1970s to understand synchronization in a large system of oscillators. It has subsequently been used as a toy model by many neuroscientists (as well as scholars in many other areas), as some of the characteristics of its synchronization patterns resemble some of the ones in neuronal communities.^{73–76} The Kuramoto model and its generalizations have also been applied to numerous other applications in chemistry, biology, and other disciplines.^{70,71,77}
When all oscillators are coupled to each other, the Kuramoto model is most commonly written as^{69,71}
where θ_{i} denotes the phase of oscillator i, the parameter ω_{i} is its natural frequency, $ K \u2265 0 $ parametrizes the coupling strength between different oscillators, and N is the number of oscillators in the model. The normalization factor $ 1 N $ ensures that the righthand side of Eq. (5) is bounded as $ N \u2192 \u221e $. The distribution from which the frequencies ω_{i} are drawn is usually assumed to be unimodal and symmetric about its mean frequency, which can be set to 0 due to the rotational symmetry of the model (because Eq. (5) is invariant under translation of θ_{i}). The parameter ω_{i} then denotes the deviation from the mean frequency.
We also adapt Eq. (5) to examine a network of N oscillators with uniform coupling between the oscillators.^{40,46,70,71,78} We consider the following generalized version of Eq. (5):
where $ \kappa \u2265 0 $ denotes the normalized coupling strength and the entries of the coupling matrix $ A = ( A i j ) i , j = 1 N $ indicate whether oscillators i and j are coupled. That is, A is an unweighted adjacency matrix, and A_{ij} = 1 for coupled oscillators and A_{ij} = 0 for uncoupled oscillators. The coupling matrix A thereby imposes a “structural network” between the oscillators. One can further generalize Eq. (6) by using heterogeneous coupling strengths κ_{ij} or by considering functions other than sine on the righthand side.
We divide the oscillators into 8 separate communities^{79} of 16 distinct oscillators each, and we suppose that every oscillator has exactly 14 connections, 13 of which are with oscillators in the same community and 1 of which is to an oscillator outside the community. As in Bassett et al.,^{46} we choose a coupling strength of $ \kappa = 0.2 $, consider a network with N = 128 oscillators, and suppose that the ith natural frequency is $ \omega i \u223c N ( 0 , 1 ) $. (That is, we draw natural frequencies from a Gaussian distribution with mean 0 and standard deviation 1.) However, our network architecture differs somewhat from that in Bassett et al.,^{46} where every oscillator had at least 13 connections inside its community and at least 1 connection outside its community.
We simulate the Kuramoto model [Eq. (6)] using the Runge–Kutta Matlab solver ODE45 (with an integration time interval of $ [ 0 , T max ] $, where $ T max = 10 $).^{80} We observe the system for M = 500 time steps in total (including the initial time step) and obtain time series $ T ( i ) = ( \theta i ( t 0 ) , \u2026 , \theta i ( t 499 ) ) $ as the output of the model for each oscillator θ_{i}. Kuramoto oscillators with a similar imposed community structure were demonstrated previously to initially synchronize rapidly within their communities, followed by a phase of global synchronization in an entire network.^{46} (There have also been other studies of community structure via synchronization of Kuramoto oscillators.^{78,81}) To study the dynamics of the coupled Kuramoto oscillators, we follow the approach of Bassett et al.^{46} and partition the time series into two time regimes, which we denote by $ k \u0302 = 1 $ and $ k \u0302 = 2 $. In our example, these time regimes each consist of 250 time steps. (The authors of Ref. 46 also split their time series into two equal parts, but their time series consist of 100 time steps in total rather than 500.)
To quantify the pairwise synchrony of two oscillators i and j, we use the local measure^{46,78}
where the angular brackets indicate that we take a mean over 20 simulations. We use the absolute value both to facilitate comparison with Arenas et al.^{78} and Bassett et al.^{46} (by making the same choice that they made) and to avoid negative values, which can complicate interpretation and pose other difficulties in network analysis.^{43,47,82}
In each simulation, we choose the initial values for the phases θ_{i} from a uniform distribution on $ ( 0 , 2 \pi ) $ and draw the natural frequencies ω_{i} from $ N ( 0 , 1 ) $. We apply the same underlying coupling matrix $ A = ( A i j ) i , j = 1 N $ for all 20 simulations and then use the values $ \varphi i j $ to define the edge weights in the fully connected, weighted network of Kuramoto oscillators for each time regime. We also study a network based on one full time regime that consists of 500 time steps. In analogy to neuronal networks, we call these networks “functional networks.” In Fig. 6, we illustrate our pipeline for creating a functional network from the output of a simulation of the Kuramoto model.
B. Null models for the Kuramoto data
To assess whether our observations illustrate meaningful dynamics of the Kuramoto model or whether they can be explained by a random process, we consider two different null models based on the timeseries output. In the first null model, which we call the “simple null model,” we independently reassign the order of the time series of each oscillator according to a uniform distribution before computing the similarity measure with Eq. (7). The second null model, which we call the “Fourier null model,” is based on creating surrogate data using a discrete Fourier transformation. This approach^{83} has the advantage of preserving not only the mean and the variance of the original time series but also the linear autocorrelations and cross correlations between the different time series.
To construct the Fourier null model, we start by taking the discrete Fourier transform
of a timeseries vector $ T $ (with components $ T m $) of length μ. In our case, μ = 250 or μ = 500, depending on whether we are examining two different time regimes or just one. We then construct surrogate data by multiplying the Fourier transform $ T \u0302 n $ by phases a_{n} chosen uniformly at random from the interval $ ( 0 , 2 \pi ) $, aside from the constraint that they must satisfy the following symmetry property: for every $ n \u2264 \mu $, there exists $ n \u0303 $ such that $ a n = \u2212 a n \u0303 $. This symmetry ensures that the inverse Fourier transform yields real values. The surrogate data $ \sigma = ( \sigma 1 , \u2026 , \sigma \mu ) $ are thus given by
Both the simple null model and the Fourier null model were used previously on timeseries output of coupled Kuramoto oscillators, and they exhibit different dynamics from those of the oscillators.^{22,46}
C. Persistent homology applied to the Kuramoto model and null models
We apply the WRCF to functional networks created from the output of two time regimes of the Kuramoto model, one time regime (i.e., the full timeseries output) for the Kuramoto model, the simple null model, and the Fourier null model. We run the filtrations up to filtration step 1800 for the first time regime and up to 2000 for the second; we go up to filtration step 1100 for cases in which we only consider one time regime. The total number of edges in the network, and thus the total number of possible filtration steps, is 8128. The number of filtration steps thereby corresponds to respective edge densities of 0.22, 0.25, and 0.14 for the three examples above; in each case, this amounts to a threshold that is approximately in the middle of the range of the edgeweight values. The Masters thesis of Stolz,^{84} which is a precursor to the present paper, also applied PH to networks created from the Kuramoto model, and such an example was subsequently also studied using Betti curves by other authors.^{85}
As we described in Section I, we focus our analysis on topological features in dimension 1, so we examine loops in the network. In the first row of Fig. 7, we show the 1dimensional barcodes for the networks constructed from time regime 1 (i.e., the first 250 time steps of the dynamics) and time regime 2 (i.e., time steps 251–500 of the dynamics) for the WRCF of the Kuramoto model. The barcode for each time regime includes several very shortlived bars between filtration steps 50 and 300. For the second time regime, we find more short bars for a longer filtration range at the beginning of the barcode. The loops that correspond to these short bars are all formed within the strongly synchronized communities. In fact, in time regime 1, the first 44 bars in the barcodes represent intracommunity loops; in time regime 2, only 2 of the first 28 bars represent intracommunity loops. As strong intracommunity edges are added to the simplicial complexes, they start to cover the loops with triangles (i.e., 2simplices), and the loops disappear from the filtration.
In the second row of Fig. 7, we show the persistence landscapes that we construct from the 1dimensional barcodes. We ignore infinitelypersisting bars in the barcode. (We also studied persistence landscapes that include the infinite bars as features with a death time that corresponds to the maximum filtration value but did not obtain any additional insights that way.)
As expected, the landscapes have a group of small peaks early in the filtration for both time regimes. This feature occurs in a longer filtration range in the second time regime before morepersistent loops appear. In the second time regime, some of the peaks that occur early in the filtration appear to almost double their heights to values of about 100. In contrast, in the first time regime, peaks at a similar location are about half as high (i.e., they are less persistent).
The persistence landscapes reveal more persistent loops in the second time regime (i.e., between time steps 251 and 500) than in the first (i.e., between time steps 1 and 250), and the second time regime also appears to reveal a clearer separation between a group of very early short peaks and a group of mediumsized peaks towards the end of the filtration. For this second group of mediumsized peaks, we observe a larger absolute increase in persistence in the second time regime than for the shorter peaks early in the filtration. These observations reflect the dynamics of the two time regimes in the Kuramoto model.^{46} In time regime 1, there is strong synchronization within the communities, and such dynamics are reflected by the appearance of shortlived intracommunity loops (corresponding to the short peaks in the persistence landscapes) early in the filtration. In the second time regime, the amount of global synchronization is more prominent than in the first time regime. Moreover, in addition to intracommunity loops, some of the peaks early in the filtration now represent intercommunity loops, which are more persistent than the loops within communities. Additionally, as some of the peaks that correspond to intercommunity loops have shifted to earlier parts of the filtration, there is an increase in the gap between the initial group of peaks and the group of mediumsized peaks at the end of the filtration. In general, we observe increased persistence of the peaks in the landscapes due to the stronger synchronization between the communities. These observations are much easier to visualize using persistence landscapes than using barcodes.
We calculate pairwise L^{2} distances between all dimension1 persistence landscapes, and we note that L^{2} distance has been used previously to compare persistence landscapes in an application to protein binding.^{62} The L^{2} distance between the two time regimes is 27078. (Here, and in subsequent discussions, we round the L^{2} distances to the nearest integer.) Given the length of the support of the landscapes and the function values that the landscapes attain, this is a large distance, which captures the aforementioned visible differences between the landscapes. The L^{2} distance is unable to capture the fact that the peaks that appear early in the filtration in the first time regime correspond to loops between nodes within one community, whereas they correspond to loops that form between nodes of different communities in the second time regime. Consequently, this feature does not contribute to the value of the distance.
We also compare the Kuramoto model with the two null models that we discussed in Section III B. To do this, we construct a functional network by considering a single time regime that consists of 500 time steps. In Fig. 8, we show the weighted adjacency matrices of the three functional networks, and we also show their corresponding persistence landscapes based on WRCFs of the functional networks. One can observe clearly that there is stronger intracommunity synchronization for the Kuramoto time series than for the null models, as there is a very distinct group of short peaks early in the filtration (which, as we discussed above, is also the case for the Kuramoto model when performing separate calculations in the two time regimes).
Again, the corresponding loops occur within communities. The peaks in the Kuramoto landscape appear to be separated from a second group of short peaks further along in the filtration. Between the two groups of peaks, there are two strikingly higher peaks that correspond to persistent loops that appear to be formed by connections between different communities. For both null models, we also observe groups of short peaks early in the filtration, but these are less persistent and separated less clearly from other peaks than for the Kuramoto model. Indeed, we do not see any separation at all for the Fourier null model, which exhibits a much weaker intracommunity synchronization than the simple null model. Moreover, the persistence landscape for the Fourier null model appears to be “noisier,” as the majority of the peaks in the landscape have similar persistences and appear in similar regions of the filtration.
The peaks in the landscapes of the null models appear to have a very different distribution along the filtration than is the case for the Kuramoto model. They also possess more mediumsized and long persisting features than we observe in the Kuramoto data. These features occur in parts of the filtration in which the Kuramoto data have a smaller number of peaks. They consist of intercommunity loops and are a symptom of the weaker intracommunity and stronger intercommunity synchronization. The null models thus appear to have more topological features in the form of loops than is the case for the Kuramoto data. This is consistent with previous examinations of null models in other studies.^{33,53,59} The fact that there are fewer persistent loops in the Kuramoto model than in the null models implies that there are more highdimensional simplices (e.g., triangles and tetrahedra) in the corresponding network than in the networks constructed from the null models.
To distinguish between the three landscapes, we calculate the L^{2} distances between them. The L^{2} distance between the Kuramoto landscape and the Fourier nullmodel landscape is 13540; the L^{2} distance between the two nullmodel landscapes is 13263; and the L^{2} distance between the Kuramoto landscape and the simple nullmodel landscape is 11703. Again considering the support of the landscapes and the magnitudes of the function values, we see that three distances can be construed as large.
We find that PH can detect Kuramotomodel dynamics and that the persistence landscapes are rather different for the Kuramoto model and the null models. The L^{2} distances between landscapes underscore these differences. We are also able to distinguish between the two null models using persistence landscapes. In contrast to conventional wisdom,^{54,55} we do not find for our examples that only the persistence of topological features distinguishes signal from noise. In fact, the short bars early in the filtration of the Kuramoto model carry important information about the dynamics, and the mediumsized persistent peaks in the Fourier null model are a symptom of the weaker intracommunity and stronger intercommunity synchronization in that model. We therefore assert that the position of features in the barcode is as important as persistence length for their interpretation in our examples; this provides an important point to consider for future studies. Note that persistence landscapes alone do not provide enough information to assess a system's dynamics. It is only by combining them with information about nodes that form loops (which are represented by certain groups of peaks) that we are able to obtain conclusions about intracommunity and intercommunity synchronization.
IV. EXAMPLE II: TASKBASED fMRI DATA
A. Human brain networks during learning of a simple motor task
We use a data set of functional brain networks from experiments that were first analyzed by Bassett et al.^{86} The data set was collected to study human subjects during learning of a simple motor task, and a full description of the experiments conducted is available in Ref. 86. We apply a WRCF to the functional networks, and we compare our findings with previous studies on these and similar networks.^{86–88} The functional networks are based on functional magnetic resonance imaging (fMRI) time series^{89,90} from 20 healthy subjects who undertook a motorlearning task on three days (during a fiveday period). During the imaging of the subjects, an “atlas” of 112 brain areas was monitored while they were performing a simple motorlearning task (similar to a musical sequence), which they executed using four fingers of their nondominant hand. For each subject and for each day of the study, the fMRI images are interpreted as 2000 time points for each monitored brain region. The brain regions and their time series were used subsequently to construct functional networks based on a functional connectivity measure known as the coherence of the wavelet scale2 coefficients. This measure was applied to the time series to determine edge weights between each pair of brain regions in the network. The weighted adjacency matrices for the functional networks were then corrected for a falsediscovery rate, as matrix elements under a certain threshold (which represents a coherence amount that one expects to occur at random) were set to 0. The other matrix elements were retained.
The functional networks that we just described were studied previously using community detection by Bassett et al.,^{86} whose results suggest that there is a significant segregation of the nodes in the functional networks into a small number of different communities with stronglyweighted connections inside the communities and weaklyweighted connections to nodes in other communities. Within these communities, certain nodes appeared to remain in the same community, whereas others (the “flexible” ones) often switched between different communities.
There have also been studies of networks from a similar experiment that examined mediumterm learning (instead of shortterm learning) and included training sessions.^{87,88} These networks have a noticeable core–periphery organization, with the sensimotor and visual regions of the brain grouped into a temporally “stiff” core of nodes, whose community memberships (in contrast to flexible, peripheral nodes) do not change much over the course of the learning task.^{87} It was also shown subsequently that the interaction of the primary and secondary sensorimotor regions with the primary visual cortex decreases as the regions (presumably) become more autonomous with task practice.^{88}
Because we observed shortlived loops early in the filtrations for the Kuramoto model in a simulated setting with community structure in oscillator connectivity, we will explore whether the fMRI data exhibit similar features during the three observation days.
B. Persistent homology applied to the taskbased fMRI data
We run the WRCF until filtration step 2600, which is when $ 42 % $ of the edges are present in the network. (Note that using more filtration steps leads to very long computational times.) We again focus our analysis on topological features in dimension 1. We construct persistence landscapes for dimension 1 (omitting infinitelypersisting loops). In Fig. 9, we summarize our results for one particular subject and for the whole data set. We use this subject to illustrate a representative example of the particular landscape features that we observe in the data.
Similar to the Kuramoto oscillators in Section III, we find a group of small peaks early in the filtration (between filtration steps 1 and 200). We can see this group very clearly by magnifying both the landscape of individual subjects and the average landscape, whose peak heights are slightly smaller than the heights of the peaks in the individual landscape that we show. This feature of the heights indicates that a group of short peaks arises early in the filtration in the majority of the barcodes. We also consider the standard deviation from the average landscapes in the first 200 filtration steps. For all three days, it is very small: it is 127 for the first day, 167 for the second day, and 126 for the third day.
We expect the observed short peaks early in the filtration to be associated with network communities, which have been observed previously using other methods.^{86} We observe, in particular, that these short peaks undergo changes on day 2: during filtration steps 20 to 60, some of the peaks that are present in the landscapes for days 1 and 3 vanish, and a larger number of persistent peaks occur for day 3 than for the other two days between filtration steps 80 and 200. This appears to suggest that there is a change in community structure that takes place on day 2, with either (1) very strong synchronization in some of the communities, leading to very shortlived loops; or (2) very strong individual differences between the subjects, leading to the vanishing of peaks in the average landscapes for the first 50 filtration steps. The particularly persistent peaks on day 2 could represent either persistent intercommunity loops or loops that occur due to sparse intracommunity connections.
We calculate pairwise L^{2} distances between each pair of dimension1 persistence landscapes. We create distance vectors, which we use as an input for kmeans clustering (with k = 3) and average linkage clustering (until we obtain 3 groups), and we obtain the same qualitative result for both methods. We find that 9 of the 20 distance vectors that correspond to persistence landscapes from day 1 are assigned to a common group (together with a small number of landscapes from days 2 and 3), whereas 11 and 10 landscapes from days 2 and 3, respectively, are assigned together to a separate group. We summarize our results in Table I.
.  Cluster 1 .  Cluster 2 .  Cluster 3 . 

Day 1  9  6  5 
Day 2  5  4  11 
Day 3  5  5  10 
.  Cluster 1 .  Cluster 2 .  Cluster 3 . 

Day 1  9  6  5 
Day 2  5  4  11 
Day 3  5  5  10 
We also consider the average dimension1 landscapes for WRCF steps 1–2600 and calculate the L^{2} distances between them. We show the results of these calculations in Fig. 10.
The distances between the average landscape for day 1 and the subsequent days of the experiment indicate that the WRCFs on average are able to detect changes in the functional networks across the filtration range. Based on the distances, we observe that most of these changes occur between the first and second days. However, the standard deviations from the average landscapes are a factor of about 4 larger than the distances between the landscapes, and one therefore needs to be cautious about interpreting the results of these calculations. In a permutation test with 10000 regroupings of the landscapes, we do not find the distances to be statistically significant. We obtain pvalues of about 0.4 for the distance between the average landscapes of day 1 and day 2, about 0.85 for the distance between the average landscapes of day 2 and day 3, and about 0.6 for the distance between the average landscapes of day 1 and day 3.
We also find that the primary peak in the average landscapes in Fig. 10 shifts to the left over the course of the three days. This implies that the edge weights (between the brain regions) that give rise to persistent loops increase on average over the three days (presumably due to stronger synchronization). This can imply either that loops present on the first day consist of edges with higher edge weights on days 2 and 3 (i.e., the same brain regions synchronize more on the latter two days) or that new loops that appear on days 2 and 3 consist of edges with higher weights than those of day 1 but involve different brain regions than on day 1 (i.e., different brain regions exhibit stronger synchronization than what occurs on day 1). Brain regions that synchronize in a loop in a network may be an indication of an interesting neurobiological communication pattern that in this case also becomes stronger over the course of the learning process. To analyze the most frequently occurring edges involved in these loops, we extract “representatives” for all loops in dimension 1 across all subjects and days. (See Fig. 3 for an illustration of two different representatives of a loop in a network.) For each day, we construct a network, which we call an “occurrence network,” using the same nodes (i.e., brain regions) that we used before and assign every edge an edge weight that is equal to the number of occurrences of that edge in loops in the subjects on the given day. We then perform a WRCF on the three occurrence networks and study the representative loops that we obtain. In Table II in the Appendix, we list the brain regions that we find in loops consisting of edges that occur at least 50 times in functional networks in the subjects. We now examine loops in the occurrence networks. These particular loops may not correspond exactly to loops in associated functional networks. For example, individual edges with high edge weights that are part of a loop in the occurrence network may be part of a variety of different loops in associated functional networks, rather than part of one specific loop that occurs in many of the functional networks. Nevertheless, it is very likely that such loops are also loops in a functional network. One also needs to consider that the representative loops given by the software JavaPlex are not necessarily chosen optimally or are “geometrically nice”^{91} representatives of the loop.^{66} We address the issue of JavaPlex's choice of representatives to some extent by using PH on occurrence networks, but even then we cannot rule out possible artifacts. There exist loops in the occurrence networks that remain stable across the three days, although other loops occur on only one or two days. There also seem to be more loops that occur at least 50 times in the functional networks on days 2 and 3 than on day 1. It would be useful to study the brain regions (see Table II in the Appendix) involved in these loops to investigate their biological role(s) in motorlearning tasks.
Loop .  Day 1 .  Day 2 .  Day 3 . 

–lSuppMA–rSuppMA–rPreGy–lPreGy–  x  x  x 
–lOFG–lOP–rOP–rOFG–  x  x  x 
–lSupTempGy ant–lPP–lHG–lPT–lSupTemGy post–  x  x  x 
–lIC–rIC–rPP–lPP–  x  Variant: –rIC–rPP–lPP–lHG–lCOC–lIC–  Variant: –rIC–rPP–rHG–lPP–lIC– 
–rIC–lIC–lPut–rPut–  x  x  x 
–lIntCalC–lLinGy–lOFG–rOFG–rLinGy–rIntCalC–rSupCalC–lSupCalC–  x  Variant: –lIntCalC–lLinGy–rLinGy–rIntCalC–  
–lFP–lPaCinGy–rPaCinGy–rFP–  x  x  
–rPP–lPP–lHG–lCOC–lIC–lFOC–lInfFGyPO–lInfFGyPT–lFP–LSuppMA–lPreGy–RPreGy–rPostGy–rSupMargGy ant–rParOpC–rPT–rSupTempGy post–rSupTempGy ant–  x  Variant: –rPP–rHG–rPT–rSupTempGy post–rSupTempGy ant–  
–rOFG–lOFG–lLinGy–rLinGy–  x  
–lPaCinGy–rPaCinGy–rCinGy ant–lCinGy ant–  x  
–lFP–lFMedC–rFMedC–rFP–  x  
–lIC–lCOC–lHG–lPP–  x  
–lPHGy ant–lPHGy–rPHGy–rPHGy ant–  x  
–lInfFGyPT–lInfFGyPO–lFOC–lIC–lPP–lSupTempGy ant–lSupTemGy post–lMTGy post–lMTGy ant–lFP–lOFC–  x  
–rSupPL–rSupMargGy post–rSupMargGy ant–rPostGy–RPreGy–  x  
–rPostGy–rSupPL–lSupPL–lSupMargGy ant–rSupMargGy ant–  x  
–lIntCalC–lLinGy–rLinGy–rIntCalC–  x  
–lSupMargGy post–lAnGy–rAnGy–rSupMargGy post–rSupMargGy ant–lSupMargGy ant–  x  
–lPT–lHG–lPP–lSupTempGy ant–lSupTemGy post–  x  
Total number of loops  7  12  13 
Loop .  Day 1 .  Day 2 .  Day 3 . 

–lSuppMA–rSuppMA–rPreGy–lPreGy–  x  x  x 
–lOFG–lOP–rOP–rOFG–  x  x  x 
–lSupTempGy ant–lPP–lHG–lPT–lSupTemGy post–  x  x  x 
–lIC–rIC–rPP–lPP–  x  Variant: –rIC–rPP–lPP–lHG–lCOC–lIC–  Variant: –rIC–rPP–rHG–lPP–lIC– 
–rIC–lIC–lPut–rPut–  x  x  x 
–lIntCalC–lLinGy–lOFG–rOFG–rLinGy–rIntCalC–rSupCalC–lSupCalC–  x  Variant: –lIntCalC–lLinGy–rLinGy–rIntCalC–  
–lFP–lPaCinGy–rPaCinGy–rFP–  x  x  
–rPP–lPP–lHG–lCOC–lIC–lFOC–lInfFGyPO–lInfFGyPT–lFP–LSuppMA–lPreGy–RPreGy–rPostGy–rSupMargGy ant–rParOpC–rPT–rSupTempGy post–rSupTempGy ant–  x  Variant: –rPP–rHG–rPT–rSupTempGy post–rSupTempGy ant–  
–rOFG–lOFG–lLinGy–rLinGy–  x  
–lPaCinGy–rPaCinGy–rCinGy ant–lCinGy ant–  x  
–lFP–lFMedC–rFMedC–rFP–  x  
–lIC–lCOC–lHG–lPP–  x  
–lPHGy ant–lPHGy–rPHGy–rPHGy ant–  x  
–lInfFGyPT–lInfFGyPO–lFOC–lIC–lPP–lSupTempGy ant–lSupTemGy post–lMTGy post–lMTGy ant–lFP–lOFC–  x  
–rSupPL–rSupMargGy post–rSupMargGy ant–rPostGy–RPreGy–  x  
–rPostGy–rSupPL–lSupPL–lSupMargGy ant–rSupMargGy ant–  x  
–lIntCalC–lLinGy–rLinGy–rIntCalC–  x  
–lSupMargGy post–lAnGy–rAnGy–rSupMargGy post–rSupMargGy ant–lSupMargGy ant–  x  
–lPT–lHG–lPP–lSupTempGy ant–lSupTemGy post–  x  
Total number of loops  7  12  13 
Finally, we also apply WRCF to the average networks for each of the three days. To create an average network, we take the mean of the edgeweight values over all 20 subjects for each day separately and study the resulting network. We show the corresponding landscapes in Fig. 11.
As with the average landscapes, we find that the landscapes for the average networks have very short peaks early in the filtration. There are morepersistent features (e.g., larger peaks) on day 1 and day 3 than on day 2, and we even find (as in the average landscapes) that the larger peaks appear earlier (at about filtration step 400) in the filtration on day 3 than on day 1 (where they appear at about step 900). Additionally, on day 2, we observe many short peaks, especially in the later stages of the filtration. This is not the case for day 1 and day 3, so the day2 landscape is strikingly different visually from the other two landscapes. When calculating L^{2} distances, we again find that the landscape distance between days 1 and 2 and that between days 1 and 3 are larger than the landscape distance between days 2 and 3. From visual inspection, we see that this arises from the fact that the day1 landscape appears to have a clearer separation of short and high peaks than the landscapes for the later days. Taken together, the results for the landscapes of the average networks mirror our prior results for the average landscapes.
V. CONCLUSION AND DISCUSSION
We have illustrated applications of persistent homology to functional networks constructed from timeseries output of the Kuramoto model, null models constructed from the Kuramoto time series, and taskbased fMRI data from human subjects. In all cases, we observed that nonpersistent loops occur early in the filtrations. Although such nonpersistent features are commonly construed as noise in topological data analysis,^{54,55} we observed that these features appear to be consistent with prior segregations of the studied networks into communities of denselyconnected nodes. In one case (the Fourier null model), we even found that particularly persistent features appear to be linked to a network with a weak intracommunity synchronization. These very persistent features in the null model may thus represent noise. In other studies of PH using (different) null models,^{33,53,59} it was also observed that the null models often exhibit a richer topological structure than the empirical data. One could thus perhaps interpret the persistent features in the Fourier null model as features of the null model rather than as noise. Our results on the importance of nonpersistent features resemble previous observations for synthetic examples with barcodes that consist of short intervals (which are commonly construed as noise), but the differences between the corresponding persistence landscapes for the various spaces are nevertheless statistically significant.^{92} Our results are also consistent with the findings of a study on protein structure using PH that reported that bars of any length in the barcodes were equally important.^{93} For weighted networks, we suggest that when using a filtration based on edge weights, one needs to consider the actual birth and death times of filtration features (such as loops) in addition to their persistence to be able to determine whether they should be construed as noise or part of a signal. In particular, in the present paper, we observed that the early appearance of loops in a filtration is an important distinguishing feature of these data. They may also yield important insights on the geometry^{92} of data.^{94}
We also found—both by calculating average persistence landscapes and studying landscapes of average networks—that persistence landscapes for dimension 1 of the weight rank clique filtration (WRCF) are able to capture changes in the studied functional brain networks during the process of learning a simple motor task. Because we did not consider infinitelypersisting features and only included filtration steps 1–2600 when creating the landscapes, our result also suggests that the mediumlived (when compared with the full filtration length) persistent loops are able to capture changes in the networks, so it is not always necessary to consider a full WRCF to study the dynamics of a system. This observation is similar to a finding of Bendich et al.,^{35} who observed in their study that mediumscale barcode features were able to distinguish human brain artery networks from different age groups. This again suggests that persistence length should not be the only measure of signal versus noise when applying PH. We also found that the persistent features that dominate the middle part of the filtrations appear in earlier filtration steps on days 2 and 3 of the experiment than they do on day 1, which suggests that interesting dynamics in synchronization patterns are captured by mediumlived bars in the middle of a barcode.
As in other biological contexts, where PH has been applied successfully and has lead to new insights,^{28,30–33,35} we find that PH can lead to fascinating insights about the dynamics of a system. We were able not only to detect symptoms of previously observed community segregation, but we also found notable differences between a setup with strong community structure (in the coupled Kuramoto oscillators) and weakly synchronized communities (in the associated null models). For the taskbased fMRI data, we found that we can detect symptoms of community structure over the three days (in the short peaks early in the landscapes) of the data as well as changes in the loops that we observe in the average persistence landscapes of the functional networks. On average, most of these changes appear to take place on day 2 of the learning task. In particular, brain regions that yield loops in the functional networks on days 2 and 3 seem to exhibit stronger synchronization on average than those that yield loops on day 1. We obtained this observation both by calculating average persistence landscapes of the WRCF performed on individual functional networks and by calculating persistence landscapes based on the WRCF performed on average networks for each day. Although the landscape distances between the average landscapes are not statistically significant, our similar observations in both of our approaches suggest that our observations indeed reflect the average dynamics of the system. Our findings on loops thereby provide novel insights that complement previous studies of synchronization in functional brain networks. Of course, it is desirable to repeat our study using larger data sets.
There is a known relation between homology and graph Laplacians,^{95} and an interesting possible direction for future research would be to study the possible connections between graph Laplacians (and, more generally, spectral graph theory) and our results on barcodes and persistence landscapes.
Using methods from topological data analysis for studying networks has the important benefit of being both mathematically principled and generalizable. However, for biological interpretation, it is necessary to include information on the specific nodes that are part of topological features such as loops. Moreover, the interpretation of the results and importance of persistence versus position of a topological feature in a barcode can differ depending on which type of filtration is employed. Different topological features can also have different levels of relevance for different dynamical systems. For example, the occurrence of many mediumsized persistent features in the persistence landscape for the Fourier null model is a symptom of the weak synchronization in the communities, whereas the mediumsized persistent bars capture increasing synchronization in loops for the taskbased fMRI data. It would be interesting to apply WRCF (and other types of filtrations) to different synthetic networks with underlying communities (e.g., using stochastic block models) to investigate such ideas further. Importantly, one should include both the persistence and the position of topological features in analysis of PH. It would also be beneficial to combine topological tools with additional methods, such as persistence images,^{96} to determine the exact topological features that are responsible for the detected differences between the persistence landscapes of the different networks.
In conclusion, we have shown that persistent homology and persistence landscapes can be applied successfully to functional networks (from either experimental data or timeseries output of models), and that they can lead to fascinating insights, such as segregation of a network into communities and changes of network structure over time.
VI. SUPPLEMENTARY MATERIAL
See Supplementary Material for additional mathematical background, definitions, examples, and theorems.
ACKNOWLEDGMENTS
The experimental data were collected originally by Nicholas F. Wymbs and Scott T. Grafton through funding from the Public Health Service Grant No. NS44393, and we thank Nicholas and Scott for access to the data. We thank Danielle S. Bassett for help in providing the data, use of her Matlab code when we were debugging our code, and helpful discussions. We also thank Pawel Dłotko for useful discussions, his help with the Persistence Landscapes toolbox, and providing us with new versions of his code during our work. We also thank Alex Arenas and Nina Otter for helpful comments. B.J.S. thanks the Berrow foundation for funding during her M.Sc. degree, and she also gratefully acknowledges the EPSRC, MRC (Grant No. EP/G037280/1), and F. Hoffmann–La Roche AG for funding her doctoral studies. H.A.H. acknowledges funding from EPSRC Fellowship No. EP/K041096/1.
APPENDIX: TABLE OF OFTENOCCURRING BRAIN REGIONS IN LOOPS
In Table II, we indicate the brain regions that occur often in loops.
References
In the present paper, we use the terms “network” and “graphs” synonymously, although the former is often used in a way that includes structures that are more complicated than ordinary graphs.
Conventionally, loops (other than selfloops) in undirected graphs must have at least 3 edges, and loops in directed graphs must have at least 2 edges. In our paper, we adapt the terminology to represent the topological features that we detect in our simplicial complexes.
We use an input time step of $\Delta t=0.02$, but we note that ODE45 uses an adaptive step size.
See Ref. 90 for a recent discussion of fMRI inference and potential perils in the statistical methods in use in neuroimaging.
For example, a loop can be represented by a double loop.