Classifying Time-Varying Complex Networks on the Tensor Manifold
Romeil Sandhu
romeil.sandhu@stonybrook.edu

arXiv:1708.07904v1 [] 25 Aug 2017

Departments of Biomedical Informatics, Computer Science and Applied Mathematics/Statistics
Stony Brook University, Stony Brook, NY 11794

Abstract
At the core of understanding dynamical systems is the
ability to maintain and control the systems behavior that includes notions of robustness, heterogeneity, and/or regimeshift detection. Recently, to explore such functional properties, a convenient representation has been to model such
dynamical systems as a weighted graph consisting of a finite, but very large number of interacting agents. This said,
there exists very limited relevant statistical theory that is
able cope with real-life data, i.e., how does perform simple
analysis and/or statistics over a “family” of networks as
opposed to a specific network or network-to-network variation. Here, we are interested in the analysis of network
families whereby each network represents a “point” on an
underlying statistical manifold. From this, we explore the
Riemannian structure of the statistical (tensor) manifold in
order to define notions of “geodesics” or shortest distance
amongst such points as well as a statistical framework for
time-varying complex networks for which we can utilize in
higher order classification tasks.

1. Introduction
Notions of robustness, heterogeneity, and phase changes
are ubiquitous concepts employed in understanding complex dynamical systems with a variety of applications [1, 2,
3, 4]. This is seen in Figure 1. While recent advancements
in network analysis has arisen through the usage of spectral techniques [5], expander graphs [6], percolation theory
[7], information theoretic approaches [8, 9], scale-free networks [4], and a myriad of graph measures reliant on the
underlying discrete graph space (e.g., degree distribution
[10], shortest path [8], centrality [11, 12, 13]), such measures rarely incorporate the underlying dynamics in its construction. Limitations of such static measures can be more
aptly seen in multilayer networks whereby standard network measures are immediately ambiguous; shortest paths

Figure 1. This work focuses on the exploiting the underlying geometry in a statistical framework to elucidate regime-shifts and
bifurcation detection in dynamical systems. We apply our analysis
towards broader cell differentiation and signaling promiscuity in
biological systems with cancer applications (e.g., inducing terminal homogeneous state from a pluripotent stem-cell like state)

(geodesics) are non-trivial [14] while centrality measures
[15], clustering methods [16], and diffusion models [17]
require simplifications to make sense of real-world configurations. Moreover, static properties have been unable to
model time-varying nonlinear higher order effects (e.g., targeted therapy [18], cellular reprogramming [19]) which can
be seen pictorially in Figure 1 and Figure 2. This is akin to
waterbed theory which refers to the observation of when a
complex behavior is “pushed” down, the invariable effect of
complexity causes for that behavior to “pop-up” elsewhere
[20].
To this end, we have previously developed fundamental
relationships between network functionality [8, 21] and certain topological and geometric properties of the corresponding graph [22, 23]. In particular, recent work of ours has
shown the geometric notion of curvature (a measure of “flatness”) is positively correlated with a system’s robustness or
its ability to adapt to dynamic changes [24, 25]. This can
be seen in Figure 3A. In this regard, network curvature may
relate to anomaly detection [26], congestion in communication [27, 28], drug resistance [29], supply chain risk [30], to

Figure 2. Current modeling of static network properties leads to ill-posed “whack-a-mole”-based solutions whereby to alter functionality,
one “whacks” a vulnerabilities of a system at a given time t. However, due to heterogeneity, this leads to new vulnerabilities and further
(perhaps accelerated) breakdown via feedback.

even systemic risk in banking systems [31, 32]. Moreover,
this is particularly compelling since curvature can be readily
computed for large-scale networks using theory from optimal mass transport [33], and more importantly, highlights
the significance of employing geometry in network analysis. This said, there exists very limited relevant statistical
theory that is able cope with a family of networks and one
often focuses on a given network or is limited to measures
that exploit network-to-network variations as it pertains to
functionality. Here, we are interested in expounding upon
previous work [24] in order to exploit statistics over a family of networks such that higher order classification tasks
can be performed to better understand regime shifts.
To do so, we are interested in studying geometry of
networks as a dynamical system that evolves over time in
which each given network or observation may represent
a “point” (encoded by some matrix-based positive definite
model) on an underlying Riemannian manifold. From this,
we are interested in defining geodesics such that distances
between “points” (networks) in a given a family can be measured and proper statistical analysis can be achieved. For
a broad range of problems, such information may include,
but not limited to, vector-valued time series, power spectral
density [34], to density matrices from quantum mechanics
that aptly describes the mixed state of the system [35]. In
general, the set of Hermitian positive definite (PD) matrices
plays a fundamental role in a variety of engineering applications with the most recent attention renewed on developing
lower order statistics [36, 37, 38, 39]. In relationship to this
mathematical structure, significant related work on employing statistical analysis for Riemannian manifolds has focused on directional statistics, statistics on shape spaces, as
well as specific manifolds - we refer the reader [40, 41, 42]
and references therein. This being said, while these methods are mostly extrinsic whereby one embeds the manifold
in an ambient Euclidean space, we are interested in studying intrinsic Riemannian metrics that can be naturally applied to the network setting without restrictions on the network topology. For example, from a network geometric per-

spective, there has been very interesting work modeling the
space of tree-like structures in understanding phylogenetic
trees for cancer and evolutionary biology [43, 44]. Recently,
phylogenetic trees were shown to live in a negatively curved
space [43] resulting in a new push to develop an appropriate statistical framework (i.e., if a manifold is negatively
curved, then one can define a unique geodesic). This has
led to methods in hierarchical learning and clustering [45].
Nevertheless, as seen in Figure 4, non tree-like networks
still “live” upon some Riemannian (statistical) manifold for
which we need to construct a framework capable of resolving unique geodesics (shortest path) amongst such points.
In turn, we can then perform statistics on this space allowing
us to classify certain behaviors and operating regimes which
will be key to maintaining control and avoiding shotgunbased solutions during “black swan” events [46] whereby
the continuous failing of interacting agents may result in
catastrophic failure.
The remainder of the present paper is outlined as follows: In the next section, we first provide important preliminaries in motivating the theoretical need of understanding
geometry by revisiting curvature and entropy as it pertains
to functionality. Here, we see that changes in the geometry provides intrinsic information about the dynamic system and thus, its exploration fits towards a larger program
developed in this paper and current related work [24, 25].
From this, Section 3 lays the foundation of a Riemannian
framework previously employed for Diffusion Tensor Imaging whereby we now adapt such notions towards complex
networks. Then, Section 4 presents preliminary results on
synthetic and biological data. We conclude with a summary
and future work in Section 5.

2. Geometry and System Behavior
To illustrate how geometry elucidates the functional behavior of a dynamical system and how such ingredients are
ever-important the analysis and control of time-varying networks, let us revisit optimal mass transport (OMT) [33].

Figure 3. Visual illustration of curvature. (A) The thematic and motivating vision of this work is based on the above-illustrated result
[24, 25]. (B) Other forms of curvature that play important roles in (projective) geometry and machine learning.

The first notion of OMT was proposed by Gaspar Monge
in 1781 with the concern of finding the minimal transportation cost for moving a pile of soil from one site to another.
The modern formulation, given by Kantorovich, has been
ubiquitously used in fields of econometrics, fluid dynamics, to shape analysis [33, 47, 48] and recently, has received
a renewed mathematical interest [49, 50]. More formally,
let (X, µ0 ) and (Y, µ1 ) be two probability spaces and let
π(µ0 , µ1 ) denote the set of all couplings on X × Y whose
marginals are µ0 and
R µ1 . As such, the Kantorovich costs
seeks to minimize c(x, y)dπ(x, y)∀π ∈ π(µ0 , µ1 ) where
c(x, y) is the cost for transporting one unit of mass from x to
y. The cost originally defined in a distance form on a metric
space leads to the Lp Wasserstein distance as follows:

Wp (µ1 , µ2 ) :=

Z Z
inf

 p1
d(x, y) dµ(x, y))
p

µ∈π(µ0 ,µ1 )

(1)
From this, let us begin considering M to be a Riemannian
manifold for which “points” are now defined
Z

P := µ ≥ 0 : µ dvol(M ) = 1
Z
(2)

Tµ P := η : η dvol(M ) = 0
as the space of probability densities and the tangent space at
a given point µ, respectively. Due to the work of Benamou
and Brenier [48], one can naturally compute the geodesic
(in the Wasserstein sense) between two densities µ0 , µ1 ∈
P as the below optimal control problem:
Z Z 1
inf
µ(t, x)k∇g(t, x)kdtdvol(M )
µ,g

0

subject to

∂u
+ div(µ∇g) = 0
∂t


µ(0, .) = µ0 ,

(3)

µ(1, .) = µ1

which leads us to give P a Riemannian structure due to
some very nice work of Jordan et. al [51]. Armed with
this, we are now able to see how changes to the underlying geometry (namely curvature) may serve as a proxy for

system behavior (namely functional robustness). As such,
let
Z
H(µt ) :=
log µt dvol(M )
(4)
M

represent Boltzmann entropy where the dependency on x
has been dropped for convenience and we consider a family
of densities evolving over time. Taking the second variation
with respect to time t in the Wasserstein sense (i.e., our gradient steps are taken to minimize or maximize the Wasserstein distance rather than the classical Euclidean norm) and
noting that, by construction, η := ∂µ
∂t |t = 0, we have that
d2
H(µt )|t=0 = hHess(H)(η), ηiW
dt2 Z

1
=−
h∇gη , ∇∆gη i + ∆ k∇gη k2 µ0 dvol(M ) (5)
2
M
where µ0 and gη satisfy equation (3). Using the Bochner
formula [54], which relates harmonic functions on a Riemannian manifold to Ricci curvature (a geometric measure
of “flatness” and herein denoted as “Ric”), we can further
assume Ric ≥ kI as quadratic forms where k is a constant
and I is the identity matrix. Then, due to a beautiful result of Sturm [53] as well as Lott and Villani [33], one
can show that the Hess(H) is k-convex with respect to the
Wasserstein norm:
H(µt )≤ tH(µ0 )+(1 − t)H(µ1 ) − φ(k, t, µ0 , µ1 )∀t ∈ [0, 1]
(6)
where the right hand portion of the above equation denoted as φ(.) can be shown to be φ(k, t, µ0 , µ1 ) = k2 t(1 −
t)W2 (µ0 , µ1 )2 allowing for k-convexity. That is, changes in
entropy and curvature are positively correlated, i.e., ∆H ×
∆Ric ≥ 0. Now, let p (t) denote that probability that the
mean of a given observable deviates by more than  from
the original (unperturbed) value at time t. Then from large
deviations theory, one can define a rate function as:


1
R := lim
− log p (t)
(7)
t→∞,→0
t
Therefore, a large R means a fast return to the original state
(robustness), and a small R corresponds to a slow return

(fragility). In thermodynamics, it is well-known that entropy and rate functions from large deviations are closely
related [21]; the Fluctuation Theorem [8] is an expression
of this fact for complex networks and may be expressed as
∆H × ∆R ≥ 0. With the result of equation (6), we have
∆Ric × ∆R ≥ 0

(8)

which states that changes in Ricci curvature are positively
correlated to changes in robustness. Perhaps more importantly, the manner in which the underlying geometry
changes elucidates certain behaviors of a complex dynamical system and this stylized fact, sets the motivation
for this note in which we attempt to generalize statistics to
a family of networks that live on a certain manifold.

2.1. Scalar Model: Network-to-Network Analysis
As highlighted in the previous section, geometry of networks plays an integral role in understanding functionality. In particular, in the scalar case with the structure given
by equation (2), one can discretely defined Ricci Curvature
(due to Ollivier [22]) between any two “points” x and y as:
κ(x, y) := 1 −

W1 (µx , µy )
d(x, y)

(9)

This definition, motivated by coarse geometry, is applicable
to the graph setting whereby the geodesic distance d(x, y)
is given by the hop metric and W1 can be computed simply
via linear programming [55]. However, this measure is limited to local network-to-network variation as opposed to a
family of networks and as such, we require a generalization
to the matrix valued setting; a non-trivial problem that has
been recently studied [50]. Therefore, this note and latter
(future) work will seek to construct matrix-valued curvature
based measures. To do so, we first need notion of matrixvalued distance dM (x, y).

3. Manifold of Time-Varying Networks
This section will introduce a framework capable of defining statistics on a family of networks as seen in Figure 4.

3.1. Geodesics on the Network Manifold
Consider a dynamical system that evolves over time and
whose information is encapsulated by a positive definite
matrix (tensor). In this regard, it is well-known that the
space of tensors is not a vector space, but instead forms a
Riemannian manifold M [56]. In particular, we will leverage the theory of symmetric spaces which has been extensively studied since the seminal work of Nomizu [57]; a
comprehensive work on tensor manifolds can be found in
[58, 59]. Riemannian manifolds are endowed with a metric
that smoothly assign to each point ρ ∈ M an inner product

on Tρ , the tangent space to M at ρ. More formally, let us
denote Υ and Υ+ and as the set of all Hermitian matrices
and the cone of positive-definite matrices, respectively. We
may then define an analogous structure to equation (2) as:

Λ+ := ρ ∈ Υ+ | tr(ρ) = 1

Tρ := χ ∈ Υ | tr(χ) = 0

(10)

where Λ+ represents our space of networks and Tρ is
the corresponding tangent space. From this, we let
the observations of our dynamic system be recorded as
{ρ0 , ρ1 , ρ2 , ..., ρt } and where each “point” (network) ρi ∈
Υ+ is given by a tensor. As seen in Figure 4, we need a
sensible notion of distance that is dependent on the underlying geometry to define statistics for classification. To do,
one defines the exponential map as the function that maps
→
+
to each vector −
ρ−
0 ρ1 ∈ Tρ0 , the point ρ1 ∈ Λ on the manifold M that is reached after unit time by the geodesic starting at ρ0 with this tangent vector. The exponential map,
herein denoted as expρ : Tρ 7→ M at point ρ is defined on
the whole tangent space and for the particular space of tensors, this map is also one-to-one. Moreover, one can also
define a unique inverse map denoted as the logarithm map
Logρ0 : M 7→ Tρ0 that maps a point ρ1 ∈ M to the unique
tangent vector χ ∈ Tρ0 at ρ0 whose initial velocity is that of
the unique geodesic γ with γ(0) = ρ0 and γ(1) = ρ1 .
Thus, the problem of classifying time-varying networks
amounts to being able to properly computing geodesics on
this space and for which, we must now consider a family of
curves γ(t) on the manifold and its speed vector γ̇(t). Then,
in the general setting, our geodesic will be a curve that realizes the minimum distance (length) between any two points,
e.g., ρ0 and ρ1 . That is, we want to compute:
Lρρ10 = min
γ

Z

ρ1

ρ0
Z ρ1

= min
γ

kγ̇(t)kγ(t) dt


hγ̇(t), γ̇(t)iγ(t)

(11)

ρ0

where γ(0) = ρ0 and γ(1) = ρ1 . Given this, we seek to
exploit the Riemannian structure of the tensor space proposed by Pennec [59] towards the network setting; most recently used in Diffusion Tensor Imaging [60, 61]. Here,
authors utilize a result from differential geometry regarding
geodesics for invariant metrics on affine symmetric spaces
[62, 63] together with the Sylvester equation from control,
to propose an affine invariant metric whereby the unique
geodesic curve from point ρ (and at the origin) on our network manifold with a tangent vector χ can be shown to be
1

1

1

1

γ(t)(ρ,χ) = ρ 2 exp(tρ− 2 χρ− 2 )ρ 2
γ(t)(I,χ) = exp(tχ)

(12)

Figure 4. (A) To understand more global functionality, we treat time-varying networks as “points” on a manifold. (B) The study geodesics
is needed in order to compute a statistical framework for a family of networks. (C) Manifold of Positive Definite Matrices

P+∞
k
where exp(χ) = k=0 = χk! is the usual matrix exponential. From this, the speed vector γ̇(t) can be computed:
dγ(t)
d
=
exp(tχ)
dt
dt
 tχ
e 1,1
d

0
Q
=
dt
0

0

etχ2,2
0

= QDIAG χi exp(tχi ) QT
1


!
...
0
T

...
0
Q
... etχn,n

1

Σ=

= exp(tχ) 2 χ exp(tχ) 2 )
1
2
= γ(t)(I,χ)
?χ

(13)
T

where the A ? B = ABA . Noting the geodesic γ(t)(ρ,χ)
between any two points ρ0 , ρ1 ∈ M where ρ1 = expρ0 (χ)
at t = 1, the logarithm map can be seen as
1

−1

−1

1

γ(1)(ρ0 ,χ) = ρ1 = ρ02 exp(ρ0 2 χρ0 2 )ρ02
−1

−1

−1

It has been noted that although no closed form expression
for the mean exists, the above converges towards a steadystate solution in a few iterations. As in the traditional Euclidean distance and similar to the Karcher or Frechet mean,
the above mean is constructed on the premise that we seek
a “point”
the sum of squared distance.e.,
PN that minimizes
2
ρ =
dist
(ρ,
ρ
).
In
a similar fashion, the variance
i
i=1
can be computed as follows:

− 12

log(ρ0 2 ρ1 ρ0 2 ) = ρ0 2 χρ0

1
1
−1
−1
→
ρ02 log(ρ0 2 ρ1 ρ0 2 )ρ02 = χ = −
ρ−
0 ρ1 = Logρ0 (ρ1 )

(14)

From this, the distance between networks can finally be
computed as
dist2 (ρ0 , ρ1 ) = kLogρ0 (ρ1 )k2ρ0
−1

−1

= klog(ρ0 2 ρ1 ρ0 2 )k22 .

(15)

We now have a natural distance that measures the length of
the geodesic curve that connects any two “points” on our
network manifold. This is particularly important as we can
now begin to perform statistics over a family of networks.

3.2. Mean, Variance, Mahalanobis Distance
Given the above geodesic distance, one can formulate necessary statistics on our network manifold. Using
the framework established by Pennec [59] with the above
affine invariant metric, the unique mean ρ̄ of N “points”
{ρ0 , ρ1 , ρ2 , ..., ρN } on the tensor manifold can be computed
via gradient descent:
 X
N
 1
 1
1
1
−
−1
ρ̄t+1 = ρ̄t2 exp
log ρ̄t 2 ρi ρ̄t 2
ρ̄t2
(16)
N i=1

N

T
1 X
Vecρ̄ Logρ̄ (ρi ) Vecρ̄ Logρ̄ (ρi )
(17)
N − 1 i=1


where Vecρ̄ (ρi ) = VecI log(ρ̄ ? ρi ) and where the operation VecI (A) is a vectorized projection of the independent coefficients of our
√ tensor in√the above√formulation, i.e.,
VecI (A) = (a1,1 , 2a1,2 , ... 2a1,n , ... 2an−1,n , an,n ).
From this, one obtains a generalized form of the Gaussian
distribution on the tensor manifold as:


1
N(ρ̄,Σ) (ρi ) = k exp
ξ(ρ̄,Σ) (ρi )
(18)
2
where ξ(ρ̄,Σ) represents the classic Mahalanobis distance
now defined in the Riemannian setting as seen below
T

ξ(ρ̄,Σ) (ρi ) = Vecρ̄ Logρ̄ (ρi ) Σ Vecρ̄ Logρ̄ (ρi ) . (19)
While the above provides a few basic tools, we have omitted, for the sake of brevity varying other advantages that
can be employed towards the network setting through the
framework established above. This will be a subject of future work in which powerful tools, previously employed for
DTI imaging, will now be used to exploit time-varying networks. Of course, to do so, we need to define our matrixbased model that will allow for networks to be represented
on the tensor manifold. This is discussed next.

3.3. Matrix-Based Model: Graph Laplacian
Given that we are focused network control, we employ an approximated graph Laplacian as our matrix-based
model whose structure provides a global representation of a
given network with many useful properties that are intrinsically tied to system functionality. While a complete review

Figure 5. (A) Graphical view of the three synthetic “toy” graphs that can be representative of varying types of communication networks,
each with 200 nodes and 400edges (i.e., only topology changes). (B) Distances of network-to-network under varying weights or “noise.”
We see that the proposed framework is able to properly quantify difference as compared to Euclidean measure, i.e., the star-like network is
“closest” to the random network as supported in [8, 24].

is beyond the scope of this note (we refer the reader to [5]),
a few areas that have garnered attention have included construction of expander graphs [6] to Cheegers inequality with
increasing attention on connections between spectral graph
theory and its respective geometry [64]. Mathematically,
one can define the Laplacian
P operator for a specific network
as ∆f (x) = f (x) − y f (y)µx (y) with f being a realvalued function which coincides with the usual normalized
graph Laplacian given by:
1

1

L = I − D− 2 AD− 2

(20)

It is interesting to note in this connection that if k ≤ κ(x, y)
is a lower bound for Ricci curvature defined in equation (9),
then the eigenvalues of Laplacian is bounded k ≤ λ2 ≤
... ≤ λN ≤ 2 − k [64]. This relationship is important
since 2 − λN measures the deviation of the graph from being bipartite, i.e., a graph whose vertices can be divided into
two disjoints sets U and V such that every edge connects a
vertex in U to one in V . As such ideas appear in resource
allocation and control, the study of such structures in the
time-varying setting motivate our focus on the Laplacian
matrix. Lastly, while one can form a normalized Laplacian
with unitary trace, we note that the Laplacian matrix by
definition is a positive semi-definite matrix. Therefore, to
enforce positive definiteness, we utilized an approximated
Laplacian L̂ := L + L I where 0 < L << 1 such that
xT L̂x > 0 for all nonzero x ∈ Rn .

4. Results
We now present preliminary results. While these experiments are by no means exhaustive (and will be extended),
these results highlight potential of the above framework for
time-varying network analysis.

4.1. Toy Network Configurations
To provide initial intuition of the proposed framework,
we constructed three toy “communication” based networks.

As one can see in Figure 5, the networks are composed of
200 nodes and 400 edges with varying topology that represents a chain, a random, and a star-like network. Similar
to our earlier work [24], we are interested in using the proposed framework to measure distances amongst these structures under varying “noise”. As highlighted by network
curvature [24] and network entropy[8], one should expect
from a communication perspective the star-like network is
“closer” to a “random” network. That is, given that their
functional properties such as robustness (a measure of “flatness”) are similar compared towards the chain-like network,
we should expect that such networks on the manifold lie
closer to one another. This is exactly what is shown in Figure 5B. Moreover, as we move from an unweighted graph
to a weighted graph in which weights are chosen randomly
in a uniform manner, this relationship still holds. On the
other hand, if we utilize the classical Euclidean distance
(Frobenius norm) to measure network-to-network distance,
we find that the chain-like network is in fact “closer” to the
star-like network. From a communications perspective, this
does not make sense and highlights the proposed framework. Lastly, it is interestingly to see that as we increase
the support of weight selection, the distances between such
networks increases as to be expected.

4.2. Scale-Free and Random Networks
This section presents clustering results and the importance of exploiting the underlying geometry as it relates
scale-free networks. In particular, we generate 50 scale-free
networks (i.e., degree distribution follows power law) [1]
and 50 random (Erdos-Renyi) networks with each network
composed 200 nodes and 1164 edges (differences are only
in topology). Such scale-free networks stylistically could
represent “business as usual” behavior while such random
networks may represent a particular crisis in which structure (correlation of agents) breaks down. From this, one
can generate an approximated normalized graph Laplacian

Figure 6. Heat map of distances between scale-free and random networks measured by the proposed Riemannian metric to that of the
Frobenius norm. (B) Average distance of network-to-network under varying weights or “noise.” We see that the proposed framework is
able to properly quantify difference as compared to Euclidean measure. Note: Each random and scale-free network possess same number
of nodes and edges, only topology changed.

as defined in Section 3.3 whereby each network now represents a point on our tensor manifold. In doing so, we compute pairwise distances amongst all potential network pairs
based on the Euclidean distance, defined by the Frobenius
norm, as well as the Riemannian distance, defined by equation (15). Figure 6A presents a “heat map” of the distances.
As one can see, Euclidean distance fails to cluster networks
as the distance between scale-free to random networks are
“closer” than random to random networks. The framework proposed here is able to properly cluster the aforementioned networks highlighting the thematic importance
of using Riemannian geometry when exploiting distances
and/or similarity measures. Figure 6B also presents average distances amongst scale-free-to-scale-free, randomto-scale-free, and random-to-random networks. Lastly, we
note that if the above experiments were repeated and/or the
number of nodes/edges were changed, the clustering result
would still hold.

4.3. Waddingtons Landscape: Cell Differentiation
In this section, we now turn our attention in utilizing the
proposed framework towards a regime shifts of cellular differentiation of Waddington’s Landscape as depicted in Figure 1. This epigenetic landscape provides a stylistic hierarchal view cellular differentiation [65], i.e., landscape elevation associates to varying degrees of pluripotency or the
ability for a cell to differentiate to multiple function specific tissues. That is, pluripotent and progenitor cells are
“higher” compared to a further differentiated state whether
that be another function specific stem-cell or tissue. Nevertheless, at a high level, one can consider the intricate and
complex processes of cellular differentiation as a interconnected dynamical system (network) whereby such bifurcations or regime-shift relate to differentiated states that lie on
Waddington’s landscape [66].
To illustrate this and from a classification perspective, we

accessed the Gene Expression Omnibus (GEO). In particular, we selected two public datasets, GSE11508 and GSE
30652, that possess 59 pluripotent / 160 non-pluripotent and
159 pluripotent / 32 non-pluripotent cell samples, respectively. Within each pluripotent samples of both datasets,
there contains varying types that include human embryonic
stem cells (hESCs) and induced pluripotent stem cells (iPSCs). This results in a potential mix of differentiation,
i.e., iPSCs pluripotency may be much lower than hESCs
yet are deemed pluripotent. Here, we are primarily interested to examine how “far” apart differentiated samples
are from their counterparts and would be key in devising
a necessary quantitive framework characterizing Waddingtons Landscape. As such, through the use of the KEGG
database downloaded from Pathway Commons, we construct a single pluripotent network whereby the weights are
the gene-to-gene correlation of twenty random pluripotent
samples selected from each of the respective data sets. This
was done similarly for a single non-pluripotent network.
Together, this was repeated to generate 50 pluripotent and
non-pluripotent networks from each GSE11508 and GSE
30652. Figure 7 presents heat map results as well as the average network distance between such networks. As one can
see, we see a marked increase in distance between pluro
and non-pluro samples. We note the Frobenius norm follows this same pattern and can primarily be attributed to the
fact that we fixed the topology between pluro/non-pluro networks, the small number of samples for correlation, as well
as the dependency of our underlying interactome KEGG
network. This said, we have presented such results here
to motivate the need for future work.

5. Summary & Conclusion
This note introduces a statistical geometric framework
[59] towards time-varying network analysis. While such

Figure 7. Heat map of Riemannian distances between pluripotent and non-pluripotent networks under two public datasets. (B) Average
numerical distance between the pluro and non-pluro networks. We see the proposed measure is able to classify such samples as well as the
Euclidean measure. We note the topology here is fixed across each networks unlike previous synthetic results.

mathematical concepts have been successfully employed
for DTI Imaging [61], the underlying structure is equally
relevant in the classification of regime-shifts in networks.
To illustrate the importance of geometry in this setting, we
presented preliminary results that show classical Euclidean
distances are unable to properly differentiate between not
only toy network configurations, but that of scale-free and
random networks. We concluded with an example of how
the proposed framework may be utilized on real-life biological stem-cell data. Putting this all together, these preliminary results and methodology set the foundation for a
more expansive study on global notions of robustness, heterogeneity, and/or regime-shift detection.

References
[1] A. Barabasi. “The Network Takeover,” Nature
Physics. 2012. 1, 6
[2] H. Kitano. “Cancer as a robust system: implications
for anticancer therapy.” Nature Reviews Cancer. 2004.
1
[3] M. Csete and J. Doyle. “Reverse engineering of biological complexity.” Science. 2002. 1
[4] A. Barbasi and R. Albert. “Emergence of Scaling Random Networks.” Science. 1999. 1
[5] R.K. Fan Chung. “Spectral Graph Theory.” American
Mathematical Society. 1997. 1, 6
[6] S. Hoory, N. Linial, and A. Widgerson. “Expander
Graphs and Their Applications.” American Mathematical Society. 2006. 1, 6

[9] J. West, G. Bianconi, S. Severini, and A. Teschendorff. “Differential Network Entropy Reveals Cancer
System Hallmarks.” Nature (Scientific reports) 2012.
1
[10] I. Taylor, R. Lindling, D. Warde-Farley, Y. Liu, C.
Pesquita, D. Faria, S. Bull, T. Pawson, Q. Morris,
and J. Wrana. “Dynamic Modularity in Protein Interactions Networks Predicts Breast Cancer Outcome.”
Nature Biotechnology. 2009 1
[11] S. Borgatti. “Centrality and Network Flow.” Social
Networks. (2005). 1
[12] S. Borgatti and M. Everett. “A Graph-Theoretic Perspective on Centrality.” Social Networks. (2006). 1
[13] G. Ghoshal and A. Barabasi. “Ranking Stability and
Super-Stable Nodes in Complex Networks.” Nature
Communications. 2011. 1
[14] P. Brodka, P. Stawiak, and P. Kazienko. “Shortest Path
Discovery in the Multilayered Social Network.? International Conference on Advances in Social Networks
Analysis. 2011. 1
[15] M. Domenico, A. Ribalta, E. Omodei, S. Gomez,
and A. Arenas. “Ranking in Interconnected Multilayer
Networks Reveals Versatile Nodes.” Nature Communications. 2015. 1
[16] P. Mucha, T. Richardson, K. Macon, M. Porter, and
J.P. Onnela. “Community Structure in Time- Dependent, Multiscale, and Multiplex Networks.” Science.
2010. 1

[7] D. Stauffer and A. Aharony. “Introduction to Percolation Theory.” Taylor and Francis. 1994 1

[17] M. Domenico, C. Granell, M. Porter, A. Arenas. “The
Physics of Spreading Processes in Multilayer Networks,” Nature Physics. 2016. 1

[8] L. Demetrius and T. Manke. “Robustness and Network
Evolution: Entropic principle.” Physica A: Statistical
Mechanics and its Applications. 2005. 1, 4, 6

[18] L. Kwong, T. Heffernan, L. Chin. “A systems biology
approach to personalizing therapeutic combinations.”
Cancer Discovery. 2013. 1

[19] J. Utikal, J. Polo, M. Stadtfeld, N. Maherali, W.
Kulalert, R. Walsh, A. Khalil, J. Rheinwald, and
K Hochedlinger. “Immortalization eliminates a roadblock during cellular reprogramming into iPS cells.”
Nature. 2009. 1

[33] C. Villani. Optimal Transport, Old and New. SpringerVerlag. 2008. 2, 3

[20] E. Nicola. “When Bode meets Shannon: Controloriented feedback communication schemes.” IEEE
transactions on Automatic Control. 2004. 1

[35] S. Gustafson and I. M. Sigal. Mathematical Concepts
of Quantum Mechanics. Springer. 2011. 2

[21] S. R. S. Varadhan. “Large Deviations and Applications.” SIAM. 1984. 1, 4

[36] Z. Chebb and M. Moahker. “Means of Hermitian positive-definite matricies based on the logdeterminant alpha-divergence function.” Linear Algebra and Applications. 2012. 2

[22] Y. Ollivier. “Ricci Curvature of Markov Chains on
Metric spaces.” Journal of Functional Analysis. 2009.
1, 4
[23] Y. Ollivier. “Ricci Curvature of Metric Spaces.” C. R.
Math. Acad. Sci. Paris. 2007. 1
[24] R. Sandhu, T. Georgiou, E. Reznik, L. Zhu, I. Kolesov,
Y. Senbabaoglu, and A. Tannenbaum. “Graph Curvature for Differentiating Cancer Networks.” Nature
(Scientific Reports). 2015. 1, 2, 3, 6
[25] R. Sandhu, T. Georgiou, and A. Tannenbaum.
“Ricci Curvature: An Economic Indicator for Market Fragility and Systemic Risk.” Science (Science Advances). 2016. 1, 2, 3
[26] B. Anderson, D. Quist, J. Neil, C. Storlie, and T.
Laine. “Graph-Based Malware Detection Using Dynamic Analysis.” Journal in Computer Virology. 2011.
1
[27] C. Wang, E. Jonckheere, and R. Banirazi. “Wireless
network capacity versus Ollivier-Ricci curvature under Heat Diffusion (HD) protocol.” IEEE American
Control Conference. 2014. 1
[28] P. Ranjan and E.H. Abed. “Nonlinear dynamics of
TCP with RED control.” Workshop on Large-Scale
Communication Networks: Topology, Routing, Traffic,
and Control. Institute for Pure and Applied Mathematics (IPAM). 2002. 1
[29] R. Sandhu, S. Lamhamedi-Cherradi, S. Tannenbaum,
J. Ludwig, and A. Tannenbaum. “An Analytical Approach for Insulin-like Growth Factor Receptor 1 and
Mammalian Target of Rapamycin Blockades in Ewing
Sarcoma,” arXiv. 2016. 1
[30] S. Wagner and N. Neshat. “Assessing the Vulnerability
of Supply Chains Using Graph Theory.” International
Journal of Production Economics. 2010. 1

[34] T. Georgiou. “An Instrinisc metric for Power Spectral
Density Functions.” IEEE Signal Processing Letters.
2007. 2

[37] F. Hai and D. Petz. “Riemannian geometry on positive
definite matricies related to means.” Linear Algebra
and Applications. 2009. 2
[38] S. Bonnabel and R. Sepulchre. “Geometric distance
and mean for positive semi-definite matrices of fixed
rank.” SIAM Journal of Matrix Analysis and Applications. 2009. 2
[39] H. Karcher. “Riemannian center of mass and mollifer smoothing.” Communications on Pure and Applied Mathematics.1977. 2
[40] P. Jupp and K. Mardia. “A unfied view of the theory of
directional statistics.” International Statistical Review.
1989. 2
[41] M. Moakher. “Means and averaging in the group of
rotations.” SIAM Journal of Matrix Analysis and Applications. 2002. 2
[42] A. Edelman, T. Arias, and S Smith. “The Geometry
of Algorithms with Orthonganlity Constraints.” SIAM
Journal of Matrix Analysis and Applications. 2002. 2
[43] L. Billera, S. Holmes, K. Vogtmann. “Geometry of
the space of phylogenetic trees.” Advances in Applied
Mathematics. 2001. 2
[44] C. Whidden and F. Matsen. “Ricci-Ollivier curvature of the rooted phylogenetic subtree-prune-regraft
graph.” Theoretical Computer Science. 2017. 2
[45] S. Zairis, H. Khiabanian, A. Blumberg, R. Rabadan.
“Moduli spaces of phylogenetic trees describing tumor evolutionary patterns.” Lectures Notes in Computer Science. 2014. 2
[46] N. Taleb. “The Black Swan: The Impact of the Highly
Improbable.” Random House. 2007. 2

[31] A. Haldane and R. May. “Systemic risk in banking
ecosystems.” Nature 2011. 1

[47] S. Haker, L. Zhu, and A. Tannenbaum. “Optimal
Mass Transport for Registration and Warping,” International Journal of Computer Vision. 2004. 3

[32] S. Battiston, G. Caldarelli, C. Georg, R. May, and J.
Stiglitz. “Complex derivatives.” Nature Physics. 2013.
1

[48] J. Benamou and Y. Brenier. “A Computational Fluid
Mechanics Solution to the Monge-Kantrovich Mass
Transfer Problem.” Numerische Mathematik. 2000. 3

[49] W. Li, S. Osher, and W. Gangbo. “A Fast Algorithm
for Earth Mover?s Distance based on Optimal Transport and L1 Type Regularization.” arXiv. 2016. 3
[50] Y. Chen, T. Georgiou, and A. Tannenbaum. “Matrix
Optimal Mass Transport: A Quantum Mechanical Approach.” arXiv. 2016. 3, 4
[51] R. Jordan, D. Kinderlehrer, and F. Otto. “The Variational Formulation of the Fokker-Planck Equation.”
SIAM Journal on Mathematical Analysis. 1998. 3
[52] J. Lott and C. Villani. “Ricci Curvature for MetricMeasure Spaces via Optimal Transport.” Annals of
Mathematics. 2009.
[53] K. Sturm. “On the geometry of metric measure
spaces.” Acta Math. 2006. 3
[54] B. Chow, P. Lu, L. Ni. “Hamilton?s Ricci Flow.”
American Mathematical Society. 2006. 3
[55] Y. Rubner, C. Tomasi, and L. Guibas. “The Earth
Mover’s Distance as a Metric for Image Retrieval.” International Journal on Computer Vision. 2000. 4
[56] M. DoCarmo. Riemannian Geometry. Birkhauser.
1992. 4
[57] K. Nomizu. “Invariant affine connections on homogeneous spaces.” American Journal of Math. 1954. 4
[58] P. Fletcher and S. Joshi. “Principal geodesic analysis
on symmetric spaces: Statistics of diffusion tensors.”
ECCV. 2004 4
[59] X. Pennec, P. Fillard, and N. Ayache. “A Riemannian
Framework for Tensor Computing.” IJCV. 2005. 4, 5,
7
[60] V. Arsingy, P. Fillard, X. Pennec, and N. Ayache.
“Log-Euclidean metrics for fast and simple calculus on diffusion tensors.” Magnetic Resonance in
Medicine. 2006. 4
[61] Y. Rathi, A. Tannenbaum, O. Michailovich. “Segmenting Images on the Tensor Manifold.” CVPR. 2007
4, 8
[62] S. Helgason. Differential Geometry, Lie Groups, and
Symmetric Spaces. Academic Press. 1978 4
[63] S. Kobayashi and K. Nomizu. Foundations of Differential Geometry. Interscience Tracts in Pure and Applied Mathematics. 1969 4
[64] F. Bauer, J. Jost, and S. Liu. “Ollivier-Ricci Curvature
and the Spectrum of the Normalized Graph Laplace
Operator.” arXiv. 2013. 6
[65] G Keller. “Embryonic stem cell differentiation: emergence of a new era in biology and medicine.” Genes
Dev 19, 1129?1155 2005. 7

[66] C. Banerji, D. Saavedra, S. Severini, M. Widshwendter, T. Enver, J. Zhou, and A. Teschendorff.
“Cellular network entropy as the energy potential in
Waddington differentiation landscape.” Nature (Scientific Reports). 2013. 7

