Persistent homology of time-dependent functional networks constructed from coupled time series

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.¹⁷ 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 network-based 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.¹⁹ An alternative approach for studying networks²⁰ is to use methods from computational topology, which explicitly incorporates “higher-order” structures beyond pairwise connections and includes algorithmic methods for understanding topological invariants, such as connectedness, loops, or holes in high-dimensional data structures^1–3(see Section II B). Although one can also represent higher-order structures using formalisms such as hypergraphs²¹ (see, e.g., a recent paper²² 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 higher-order structures.²³ 

Methods from computational topology enable one to understand global low-dimensional structures in networks, and they have led to insights in an increasingly large number of applications⁵ in diverse topics, ranging from granular materials⁶ and contagions on networks²⁴ to path planning²⁵ and collective behavior in animals.²⁶ 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,²⁷ study the effect of the psychoactive component of “magic mushrooms” (psilocybin mushrooms) on functional brain networks of humans,²⁸ analyze covariates that influence neural spike-train data,²⁹ and study structural and functional organization of neural microcircuits.³⁰ 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,³⁴ and a demonstration that topological features of networks of brain arteries in humans are correlated with their age.³⁵ 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³⁶ used tolerance spaces and Vietoris homology theory to study aspects of visual perception. In the 1990s, Muldoon et al.³⁷ developed a method to study the topology of manifolds that underlie time-series data.

In the present investigation, we use time-series data to construct so-called 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 ̃ = ( a ̃ 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,⁴³ voting among legislators,⁴⁴ and climate⁴⁵). 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.⁴⁶ 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 adjacency-matrix 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 time-series 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).¹² To do this, one typically applies a global threshold $ξ ∈ ℝ +$ 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 adjacency-matrix 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.¹² 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.⁵¹ 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.⁵² Loops are thus 1-dimensional 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²¹). It has been demonstrated in other applications (e.g., contagions on networks²⁴) that loops are important topological features of graphs, and a recent study⁵³ demonstrated the importance of loops (and related higher-dimensional objects) in structural neuronal networks. Structural and functional neuronal networks are related and share some network features,¹¹ 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 time-series data produced by coupled Kuramoto oscillators. In Section IV, we introduce and analyze our second example, which consists of time-series functional magnetic resonance imaging (fMRI) data from experiments of humans performing a simple motor-learning task. We present our conclusions in Section V, and we provide a mathematical introduction to persistent homology in the Supplementary Material.

Persistent homology (PH)^1–3 is a method from computational topology that quantifies global topological structures (e.g., connectedness and holes) in high-dimensional 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 short-lived 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).

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 = { ∅ , x }$ consists of the empty set $∅$ 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 2-sphere 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 k-simplices, where the parameter k indicates the dimension of the simplex. Every k-simplex includes k + 1 independent nodes: a point is a 0-simplex, an edge is a 1-simplex, a triangle is a 2-simplex, and a tetrahedron is a 3-simplex (see Fig. 1). Observe that the lower-dimensional simplices are contained in the higher-dimensional simplices. This allows one to build higher-dimensional simplices using lower-dimensional ones. The lower-dimensional simplices form the so-called faces of the associated higher-dimensional 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 highest-dimensional 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.

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 1-simplices 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 2-sphere 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 ∈ { 0 , … , 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 1-dimensional 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, β₀, β₁, and β₂, to represent, respectively, the number of connected components, the number of 1-dimensional holes, and the number of 2-dimensional holes in a simplicial complex. As we pointed out in Section I, we focus on the number of loops (i.e., on β₁) in our network analysis rather than on β₀, 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.²¹ 

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 $Σ n − 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¹ in a filtration in $ℝ 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.

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.⁵⁸ 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 0-simplices (i.e., nodes) and 1-simplices (i.e., edges). The weight rank clique filtration (WRCF),⁵⁹ 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 higher-dimensional simplices. One constructs a WRCF as follows:

1. Define filtration step 0 as the set of all nodes.

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

3. In filtration step t, threshold the graph at weight ν_t to create a binary graph.

4. Find all maximal c-cliques for $c ∈ ℕ$⁠, and define them to be c-simplices.

This is a valid simplicial complex: every $( c + 1 )$-clique in the graph guarantees the existence of a c-face 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.⁵⁴ 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 − 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.

As an alternative topological summary to barcodes, one can use persistence landscapes,^60,61 which consist of piecewise-linear 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 ≥ 0$⁠, the q-th persistence landscape is given by the set of functions

If the qth-largest value does not exist, then $λ 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 ${ λ 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.⁶⁰ 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,⁶² phase separation in binary metal alloys,⁶³ and music audio signals.⁶⁴ 

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 clique-finding algorithm from the Mathworks library⁶⁷ 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.⁶¹ 

The Kuramoto model^68–72 is a well-studied 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 ≥ 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 right-hand side of Eq. (5) is bounded as $N → ∞$⁠. 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 $κ ≥ 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 right-hand side.

We divide the oscillators into 8 separate communities⁷⁹ 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.,⁴⁶ we choose a coupling strength of $κ = 0.2$⁠, consider a network with N = 128 oscillators, and suppose that the ith natural frequency is $ω i ∼ 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.,⁴⁶ 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$⁠).⁸⁰ We observe the system for M = 500 time steps in total (including the initial time step) and obtain time series $T ( i ) = ( θ i ( t 0 ) , … , θ 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.⁴⁶ (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.⁴⁶ and partition the time series into two time regimes, which we denote by $k ̂ = 1$ and $k ̂ = 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.⁷⁸ and Bassett et al.⁴⁶ (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 π )$ 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 $ϕ 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.

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 time-series 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⁸³ 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 time-series 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 ̂ n$ by phases a_n chosen uniformly at random from the interval $( 0 , 2 π )$⁠, aside from the constraint that they must satisfy the following symmetry property: for every $n ≤ μ$⁠, there exists $n ̃$ such that $a n = − a n ̃$⁠. This symmetry ensures that the inverse Fourier transform yields real values. The surrogate data $σ = ( σ 1 , … , σ μ )$ are thus given by

Both the simple null model and the Fourier null model were used previously on time-series output of coupled Kuramoto oscillators, and they exhibit different dynamics from those of the oscillators.^22,46

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 time-series 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 edge-weight values. The Masters thesis of Stolz,⁸⁴ 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.⁸⁵ 

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 1-dimensional 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 short-lived 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 intra-community loops; in time regime 2, only 2 of the first 28 bars represent intra-community loops. As strong intra-community edges are added to the simplicial complexes, they start to cover the loops with triangles (i.e., 2-simplices), and the loops disappear from the filtration.

In the second row of Fig. 7, we show the persistence landscapes that we construct from the 1-dimensional barcodes. We ignore infinitely-persisting 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 more-persistent 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 medium-sized peaks towards the end of the filtration. For this second group of medium-sized 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.⁴⁶ In time regime 1, there is strong synchronization within the communities, and such dynamics are reflected by the appearance of short-lived intra-community 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 intra-community loops, some of the peaks early in the filtration now represent inter-community loops, which are more persistent than the loops within communities. Additionally, as some of the peaks that correspond to inter-community 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 medium-sized 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² distances between all dimension-1 persistence landscapes, and we note that L² distance has been used previously to compare persistence landscapes in an application to protein binding.⁶² The L² distance between the two time regimes is 27078. (Here, and in subsequent discussions, we round the L² 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² 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 intra-community 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 intra-community 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 medium-sized 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 inter-community loops and are a symptom of the weaker intra-community and stronger inter-community 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 high-dimensional 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² distances between them. The L² distance between the Kuramoto landscape and the Fourier null-model landscape is 13540; the L² distance between the two null-model landscapes is 13263; and the L² distance between the Kuramoto landscape and the simple null-model 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 Kuramoto-model dynamics and that the persistence landscapes are rather different for the Kuramoto model and the null models. The L² 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 medium-sized persistent peaks in the Fourier null model are a symptom of the weaker intra-community and stronger inter-community 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 intra-community and inter-community synchronization.

We use a data set of functional brain networks from experiments that were first analyzed by Bassett et al.⁸⁶ 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 motor-learning task on three days (during a five-day period). During the imaging of the subjects, an “atlas” of 112 brain areas was monitored while they were performing a simple motor-learning task (similar to a musical sequence), which they executed using four fingers of their non-dominant 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 scale-2 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 false-discovery 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.,⁸⁶ whose results suggest that there is a significant segregation of the nodes in the functional networks into a small number of different communities with strongly-weighted connections inside the communities and weakly-weighted 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 medium-term learning (instead of short-term 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.⁸⁷ 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.⁸⁸ 

Because we observed short-lived 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.

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 infinitely-persisting 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.⁸⁶ 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 short-lived 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 inter-community loops or loops that occur due to sparse intra-community connections.

We calculate pairwise L² distances between each pair of dimension-1 persistence landscapes. We create distance vectors, which we use as an input for k-means 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.

We also consider the average dimension-1 landscapes for WRCF steps 1–2600 and calculate the L² 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 p-values 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”⁹¹ representatives of the loop.⁶⁶ 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 motor-learning tasks.

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 edge-weight 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 more-persistent 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 day-2 landscape is strikingly different visually from the other two landscapes. When calculating L² 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 day-1 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.

We have illustrated applications of persistent homology to functional networks constructed from time-series output of the Kuramoto model, null models constructed from the Kuramoto time series, and task-based fMRI data from human subjects. In all cases, we observed that non-persistent loops occur early in the filtrations. Although such non-persistent 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 densely-connected 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 intra-community 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 non-persistent 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.⁹² 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.⁹³ 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⁹² of data.⁹⁴ 

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 infinitely-persisting features and only included filtration steps 1–2600 when creating the landscapes, our result also suggests that the medium-lived (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.,³⁵ who observed in their study that medium-scale 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 medium-lived 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 task-based 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,⁹⁵ 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 medium-sized persistent features in the persistence landscape for the Fourier null model is a symptom of the weak synchronization in the communities, whereas the medium-sized persistent bars capture increasing synchronization in loops for the task-based 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,⁹⁶ 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 time-series 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.

See Supplementary Material for additional mathematical background, definitions, examples, and theorems.

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.

In Table II, we indicate the brain regions that occur often in loops.