Can brain networks be computed with complex dynamics?

##  1 Introduction

Many natural systems can be described as collective behavior of an ensemble of
units part of a larger network (REF). To gain insight into the behavior of
these units, dynamical systems modeling has been widely used to understanding
nonlinear and chaotic phenomena in natural systems (REFs). Moreover, the
dynamics of each coupled unit aims to realistically model the evolution of the
actual natural system. However, it is challenging to address in a realistic
model of the network and are deemed analytically and computationally
intractable (REF). Rephrasing such questions in terms of a simplified model
using complex dynamics presents a mathematically approachable and canonical
way to address the limitation associated with current dynamical system
modeling (REF). In this paper, we describe the model, we summarize some of our
existing theoretical results associated with it, and we develop the first use
of this model to illustrate and classify ensemble dynamics in a brain
tractography data set.

###  1.1 Graph theory in networks

One recent direction of research has been focused on the graph theoretical
aspect of networks, particularly on investigating the relationship between the
network’s hardwired architecture and its function. By applying graph
theoretical measures of segregation (clustering coefficient, motifs,
modularity, rich clubs), integration (path length, efficiency) and influence
(node degree, centrality), various studies have been investigating the
sensitivity of a system’s temporal behavior to removing/adding nodes or edges
at different places in the network structure.

Graph theoretical approaches have been used in particular to investigate brain
networks, not only in the context of learning, memory formation and cognitive
performance, but also to understand more general organizational and functional
neural principles [7, sporns2011human, 9]. The brain can be viewed as a
“dynamic network,” self-interacting in a time-dependent fashion at multiple
spacial and temporal scales, to deliver an optimal range for biological
functioning.

Recent studies support certain generic topological properties of brain
architecture, such as modularity, small-worldness, the existence of hubs and
“rich clubs.” Such measures have been used on empirical data, to study whether
behavioral impairments are associated with specific abnormalities in
connectivity. However, no matter how well-designed or statistically powerful,
purely empirically-based analyses cannot explain in and of themselves the
mechanisms by which connectivity patterns actually act to change the system’s
dynamics, and thus the observed behavior. An underlying network model which
incorporates both graph structure and node-wise dynamics in one unified
framework, tractable theoretically or numerically, could be used _in
conjunction_ with data towards interpreting empirical results, and for making
predictions.

Along these lines, a lot of effort has been invested towards formal modeling
approaches that would explain how brain connectivity patterns may affect its
functional dynamics. Using combinations of analytical and numerical methods
from nonlinear dynamics and graph theory, such studies range from biophysical
models [gray2009stability] to simplified systems [siri2007effects]. These
analyses revealed a rich range of potential dynamic regimes and transitions
[brunel2000dynamics], shown to depend as much on the coupling parameters of
the network as on the arrangement of the excitatory and inhibitory connections
[gray2009stability]. The construction of a realistic, data-compatible
computational model has been subsequently found to present many difficulties
related to dimensionality, computational cost, addressing multiple scales,
etc. Such problems seem to transcend the existing methods from nonlinear
dynamics, and may in fact require: (1) new analysis and book-keeping methods
and (2) a new framework that would naturally encompass the rich phenomena
intrinsic to these systems.

Iterations in the complex plane are good candidates as building blocks for XXX
because dynamics on the complex plane are well understood. Further, they also
can be illustrated geometrically by sets with intricate, but topologically
quantifiable structure. This structure can be used to understand, describe,
and classify the long-term dynamics of the system.

Dynamic behavior of iterated complex quadratic maps refers to the behavior of
sequences in the complex plane obtained by picking an initial point z0 ∈ C,
and repeatedly applying a function in the family fc(z) = z2 + c (for a fixed
value of the parameter c ∈ C). This way, each complex number z0 has an
associated “orbit” under the function fc:

z →f (z )→f (f (z ))→…→f◦n(z )→… 0c0cc0 c0

For each value of the parameter c, the corresponding function fc presents with
different orbits, with potentially very different behaviors (which may vary
from trivial to very complicated, depending on both the parameter c and the
initial condition z0 itself). The study of these systems emerged historically
at the start of the 20th century, with the work of Fatou [1] and Julia [2].
These simple iterative rules describe rich, complex systems, and resulting
objects with intricate fractal geometry. Despite the seemingly elementary and
natural essence of the questions on orbit behavior for quadratic maps,
investigating them required the development of a new field of mathematics,
encompassing real and complex analysis, topology and computational methods.
Results span over a century of work, and some directions of inquiry remain
open.

Indeed, as already mentioned, complex quadratic iterations of individual
functions in the family fc = z2+c deliver incredibly rich dynamics, and
provide textbook recipes for creating fractal asymptotic sets, with history
going back more than a century [1, 2]. For a parameter value c ∈ C, the Julia
set of the corresponding map fc is defined as the boundary between initial
points z0 whose orbits remain bounded under iterations of fc (prisoners), and
those which escape to infinity under repeated iterations of the map
(escapees). The topological and fractal properties of Julia sets have been
well studied, with major results relating the geometry of the Julia set with
properties of the critical orbit of z0 = 0 [4, 5, 6]. For iterations of single
quadratic maps, it is known that the Julia set is either connected, if the
orbit of the critical point z = 0 is bounded, or

totally disconnected, if the orbit of the critical point is unbounded. The
postcritically bounded parameter locus is therefore the same as the Julia set
connectedness locus in the parameter complex plane, and is known as the
Mandelbrot set [3], the topology and properties of which have been extensively
studied since the 1960s.

Establishing relationships between quadratic dynamics and the topology of the
Mandelbrot set has been highly nontrivial; while many interesting properties
have been established (such as connectedness and full Hausdorff dimension),
there are open problems that persist to date (such as local connectedness and
density of hyperbolicity). The Mandelbrot set remains both a staple in
discrete dynamical systems, and a canonical representation for more general
phenomena.

However, application of complex dynamics has been limited in modeling natural
systems despite their potential to resolve outstanding issues. There are many
reasons for this lack of broad application but in part it pertains to complex
dynamics being considered a theoretical field of study and difficult to tie to
any natural system (REF). Specifically, one major impediment being view that
behavior of natural processes cannot truly be captured by repeated iterations
of one single quadratic function over time.

Here we develop a mathematical framework using complex quadratic dynamics to
address outstanding questions in applied dynamical systems where the network
aspect is crucial to the function of the system. Specifically, we utilize
modeled and structural brain networks derived from diffusion Magnetic
Resonance Imaging (dMRI). To preempt our results, we show that complex
quadratic dynamics reveal XXX, YYY, and ZZZ

Figure 1: _Example asymptotic sets for single map iterations. The left panels
illustrate the Julia set corresponding to the following parameter values: A.
$c=0$ (unit circle); B. $c=-0.65$; C. $c=-0.8+0.156i$; D. $c=-0.13+0.77i$; E.
$c=-0.117+0.856i$; F. $c=i$ (dendritic set). The prisoner set is shown in
black; the colors represent how fast (in how many iterations) escapees fall
outside the disc of radius 2. G. Traditional Mandelbrot set, illustrating (in
black) the values of $c$ for which the Julia set is connected, and for which
critical point $z_{0}=0$ is a prisoner. The colors outside of the set
represent how fast $z_{0}=0$ iterates out of the disc of radius 2._

###  1.2 Network iterations, and applications to dynomics

While studying iterations of one single map at the time delivers tractability
to the study of rich dynamics, realistic systems (as previously mentioned) are
generally formed of many coupled units, and cannot be encoded by the dynamics
of a single iterated map. Most natural systems, such as the brain, are
organized as networks of interacting nodes. The study of networks has
particularly been the subject of great interest. In this study, we will focus
our attention on a working example from network neuroscience. One particular
direction of investigation (currently identified as _dynomics_) addresses the
question of how the hardwired _structure_ of a brain circuit affects its
_function_ [7]. A recent trend in computational neuroscience research is
focused on studying rigorously (through a combined graph theoretical and
dynamic approach) how configuration patterns affect the efficiency of network
function.

Graph theoretical approaches have been used in particular to investigate brain
networks, and better understand organizational and functional neural
principles [7, 9, 8]. The brain is viewed as a “dynamic network,” self-
interacting in a time-dependent fashion at multiple spatial and temporal
scales, to deliver an optimal range for biological functioning. By applying
graph theoretical measures, various studies have been investigating how
different properties in the network architecture impact the ensemble dynamics,
and have been studying the sensitivity of a system’s temporal behavior to
removing/adding nodes or edges at different places in the network structure.
Understanding the way in which various parts of the brain (from the micro-
scale of neurons to the macro-scale of functional regions) are wired together
and control complex cognition and behavior is one of the great scientific
challenges of the 21st century, currently being addressed by large-scale
research collaborations, such as the Human Connectome Project [10]. The
construction of a realistic, data-compatible computational model has been
subsequently found to present many difficulties related to dimensionality,
computational cost, addressing multiple scales, etc. Such problems seem to
transcend the existing methods from nonlinear dynamics, and may in fact
require: (1) new analysis and book-keeping methods and (2) a new framework
that would naturally encompass the rich phenomena intrinsic to these systems –
both of which aspects are central to our work.

In our work, we use quadratic iterations to develop a new framework for
understanding how dynamic behavior emerges in large networks of nodes, and how
it depends on the network structure. Each node is viewed symbolically as an
integrator of internal and external inputs. We use the asymptotic behavior of
multi-dimensional orbits (via the topological and fractal structure of Julia
and Mandelbrot sets) to quantify dynamic behavior under perturbations of the
network architecture. Topological landmarks of Mandelbrot sets provide
valuable means of classification and comparison between different systems’
behavior, as will be further summarized in the following section.

In our CQN network model, the node-wise dynamics is considered to be complex
quadratic dynamics, within the family $f_{c}\colon\mathbb{C}\to\mathbb{C}$,
$f_{c}(z)=z^{2}+c$. More specifically, in this framework, each network node
adds all weighted inputs from adjacent nodes, and integrates the sum of
inputs, in discrete time, as a complex quadratic map. Then the system takes
the form of an iteration in $\mathbb{C}^{n}$:

$\displaystyle z_{j}(t)\longrightarrow z_{j}(t+1)$ $\displaystyle=$
$\displaystyle f_{j}\left(\sum_{k=1}^{n}{g_{jk}A_{jk}z_{k}}\right)$

where $n$ is the size of the network, $\displaystyle A=(A_{jk})_{j,k=1}^{n}$
is the binary adjacency matrix of the oriented underlying graph, that is
$A_{jk}=1$ if there is an edge from the node $k$ to the node $j$, and
$A_{jk}=0$ otherwise. The coefficients $g_{jk}$ are the signed weights along
the adjacency edges (in particular, $g_{jk}=0$ if there is no edge connecting
$k$ to $j$, that is if $A_{jk}=0$). In isolation, each node $z_{j}(t)\to
z_{j}(t+1)$, $1\leq j\leq n$, iterates as a quadratic function
$f_{j}(z)=z^{2}+c_{j}$. When coupled as a network with adjacency $A$, each
node will act as a quadratic modulation on the sum of the inputs received
along the incoming edges (as specified by the values of $A_{jk}$, for $1\leq
k\leq n$).

This is a simplified framework, which turned out to preserve, often in a
weaker form, some of the properties and results determined in the traditional
case of a single iterated quadratic map (e.g., the existence of an escape
radius is guaranteed in some types of networks, but may fail in others). We
used this to our advantage, and proceeded to study theoretically the effects
of network architecture on its long-term dynamics. We have been generally
investigating whether one can use properties of multi-dimensional orbits in
$\mathbb{C}^{n}$, in particular their asymptotic behavior (via the topological
and fractal structure of Julia and Mandelbrot multi-sets), to classify dynamic
behavior for different network architectures. To this aim, we defined network
asymptotic sets as follows:

######  Definition 1.1.

We define the multi-Mandelbrot set (or the multi-M set) of the network the
parameter locus of ${\bf c}=(c_{1},...,c_{n})\in\mathbb{C}^{n}$ for which the
multi-orbit of the critical point $(0,...,0)$ is bounded in $\mathbb{C}^{n}$.
For fixed parameter ${\bf c}=(c_{1},...,c_{n})\in\mathbb{C}^{n}$, the multi-
prisoner set of the network is the locus of initial points
$(z_{1},...,z_{n})\in\mathbb{C}^{n}$ which have a bounded multi-orbit in
$\mathbb{C}^{n}$. The muti-Julia set of the network is the boundary (in
$\mathbb{C}^{n}$) of the network’s prisoner set.

Broadly speaking, one can interpret iterated orbits as describing temporal
trajectories of an evolving system. Along these lines, we view a parameter
${\bf c}$ in the Mandelbrot set (for which the critical orbit that escapes to
$\infty$) as representing a system with unsustainable long term dynamics when
initiated from rest, while the ${\bf c}$ range for which the critical orbit
remains bounded can be viewed as the “sustainable dynamics locus.”

Since topological properties of $n$-dimensional complex sets are both
difficult to visualize and hard to investigate, we started by considering
Mandelbrot sets for networks with identical nodes ${\bf c}=(c,...c)$. These
are object in$\mathbb{C}$, which we called equi-M sets.

######  Definition 1.2.

We call the equi-Mandelbrot set (or the equi-M set) of the network, the locus
of $c\in\mathbb{C}$ for which the critical multi-orbit is bounded for equi-
parameter ${\bf c}=(c,c,...c)\in\mathbb{C}^{n}$.

In order to theoretically test such hypotheses, and to tie connectivity-based
theories to the dynamic patterns we had observed in empirical brain signals,
we first centered our work around simple, low-dimensional networks, which are
both analytically tractable and allow easier visualization and interpretation
of the results. We documented the effects that excitation versus inhibition in
the system have upon the size, shape and fractality of the uni-Julia sets. We
showed that weak inhibitory coupling is more efficient in breaking down the
uni-Julia set into connected components than introducing even excessive levels
of excitation. A complex network (e.g., the brain) may have to pre-compute
which hardwired structure is most effective to use in order to obtain the
desired effect, and how to plastically modify this structure, constantly
adapting to new behavioral requirements.

Figure 2: _Examples of uni-Julia sets for three different small networks._
Figure 3: _Examples of disconnected equi-M sets for three networks with N = 3
nodes, interacting according to $z_{1}(n+1)=z_{1}^{2}(n)+c$,
$z_{2}(n+1)=(a_{21}z_{1}(n)+z_{2}(n))^{2}+c$,
$z_{3}(n+1)=(z_{1}(n)+z_{2}(n)+a_{33}z_{3}(n))^{2}+c$, where the connectivity
weight parameters are as follows: A. $a_{21}=−1$, $b_{33}=−1$; B.
$a_{21}=−1/3$, $b_{33}=−1/3$; C. $a_{21}=−2/3$, $b_{33}=−1/3$._ The equi-M set
is shows in black. The colors outside the set illustrate how fast the multi-
critical orbit iterates outside of the disc of radius $R=20$.

In our previous work, we showed that the topology of equi-Mandelbrot sets
captures crucial information on both network architecture and connection
strengths. For example, inserting or deleting just one connection may
significantly change the geometry of the equi-M set of the network; this
effect depends tightly, however, on the position of the inserted/deleted edge
in the network (better off, on the graph-theoretical contribution of the edge
to the network as a whole), as well as on the weight of the edge.

In general, we showed that the changes to the architecture and the weights of
the connectome can induce the type of changes in network equi-M sets that are
not only visually striking, but also topologically quantifiable. In [11,
radulescu2019asymptotic] we illustrated how variations in the connectome even
in small networks of 3 nodes can produce equi-M sets with a wide range of
topological properties. While equi-M sets still have a main cardioid, the cusp
of which is robust to network variations, they don’t necessarily preserve
connectedness, or the hyperbolic bulb structure in the traditional sense, or
even the position of the cusp. One important theoretical task is to tease
apart which type of connectome changes trigger which type of topological
perturbations (e.g., what graph-theoretical conditions lead to 3-node
structure with connected versus disconnected equi-M sets). One equally
important (and not unrelated) practical question to ask is whether natural
network are more prone to produce some equi-M topologies versus others. While
the first question is the subject of a different study, the second issue will
be central to this paper.

Our previous work hints to the potential of equi-M sets as statistical
descriptors and classifiers. This paper aims to relate for the first time our
technique with its potential applications to assessing brain dynamics and
behavior in human subjects. The application stems from the potential to use
our modeling framework to produce a fractal, topological signature of the
equi-M set for any given subject’s brain network. This signature (1) is based
on the connectivity architecture of the brain, since the node-wise quadratic
dynamics are used in conjunction with the specific connectome for each
individual; (2) is designed to encode the long-term, global dynamic outcome
for the network, without having to biophysically model each node. In this
paper, we investigate to what extent this signature encapsulates, on one hand,
basic patterns of the connectome, and, on the other hand, basic physiological
information (gender, age) of the corresponding subjects. This represents a
preliminary exploration of the potential for using this signature as a
classification and prediction tool. In prior work, we also developed a
stochastic version of the equi-M set, that can be computed to describe cross-
sectional properties over a whole collection of connectomes (which we will use
in this paper).

##  2 Methods

###  2.1 Subjects and DTI data

For a number of $N=197$ subjects, DTI scans provide connectivity information
within a network consisting of $n=116$ brain regions, delimited by masks. Info
on the source of the data, masks etc needs to be added here.

The data for each subject included a large collection of additional measures
related to physiological, behavioral and cognitive data. For the purposes of
this study, we used only a few of these measures. In our comparisons, we use
the gender and age of the subjects ($N=90$ female, mean age $\mu\sim 28$
years, and $N=107$ male, mean age $\mu\sim 25.8$ years (where the subjects’
age was recorded in terms of four age compartments: 22-25, 26-30, 31-35, $>36$
years). Out of the cognitive tasks, we focused on the ones that had complete
data for all subjects, as listed in Table LABEL:cog_measures.

###  2.2 Data processing

The network regions/nodes were labeled as a vector $(x_{1},...x_{n})$ (the
listing order was established arbitrarily, but was preserved consistently
between subjects). With this notation, the DTI connectivity information for
each subject was stored into the subject’s $n\times n$ connectivity matrix
$C$, with entries $C_{ij}\in[0,1]$ representing the normalized thickness of
white matter projections between nodes $x_{i}$ and $x_{j}$. This is always a
symmetric matrix, with positive entries. Many entries are zero (signifying the
virtual lack of connection between the respective nodes), hence the matrix is
sparse. The strictly positive entries (corresponding to the presence of
projections between the node pair) are generally very small, which is a
characteristic of brain connectivity.

For each data, we computed two different types of measures. First, we
estimated basic topological measures of each network M set: position of the
cusp and tail, area, horizontal and vertical diameter. Second, we computed
graph theoretical measures for each connectome $C_{ij}$: average degree,
average connection strength, eigenspectrum, clustering coefficient,
centrality, betweenness and local efficiency, as well as number of 3-node and
4-node motifs. In the next sections, we provide additional information on the
definition of these measures, their computation and significance.

Figure 4: _Examples of M sets for three connectomes in our data set. Subtle
differences can be distinguished in the equi-M set geometries, as well as
between the corresponding connectivity matrices._

###  2.3 Topological measures of asymptotic dynamics

The network equi-M set was computed for each individual connectome, using 100
iteration steps and 400 by 400 spacial resolution. These parameters were
chosen to allow a first round of estimates, explored in this paper as a proof
of principle; these can be easily improved in future iterations of this
analysis, at the cost of higher computational power. Figure 4 illustrates
three examples of equi-M sets obtained this way, for three distinct
connectomes in our data set. Notice how differences in position, size, shape
of the tail and boundary between the equi-M sets of the three networks (top
panels) correspond to subtle differences between the brain connectomes (bottom
panels). It is interesting to note that, while significant differences exist
between equi-M sets, these are not as qualitatively distinct as the general
theory allows them to be (e.g., all are connected). This was more generally
true for all the 197 equi-M sets computed for the connectomes in our data set.
An interesting question to pursue is to understand not just the source of the
differences between sets, but also the source of the common properties. These
may be the specific of brain connectomes on equi-M sets, and may have
potential correlations with actual brain function.

To this end, we estimated, from the equi-M sets obtained under in our
computational setup, a few descriptors of their size and position, which are
indirectly related to their common topological landmarks (e.g., the existence
of a cusp and a tail, connectivity, symmetry of the set in the complex plane).
Additional and more computationally advanced measures, likely better
candidates for capturing more subtle topological details, are proposed for
future analyses in Section 4.2.

_Position $\gamma$ of the cusp and $\tau$ of the tail._ Our prior work showed
that the presence of the cusp is a robust, network-independent feature of
equi-M sets – unlike the tail structure, which is very fragile to network
modifications. However, the position $\gamma$ of the cusp along the real axis
varies with the network properties, as illustrated in Figure 5. For a fixed
equi-M set, $\gamma$ is easy to compute (as the largest coordinate reached by
the boundary of the equi-M set along the real axis), and relatively robust to
choosing the spacial resolution (in fact, it preserves quite faithfully the
elephant valley formation from the traditional Mandelbrot set). The accuracy
of the tail position $\tau$, computed as the leftmost point of the equi-M set
along the real axis, is more resolution-dependent, since the tail is connected
by thin filaments.

_Area ${\cal A}$ of the set in the complex plane_. Our computation algorithm
provides an overestimate of the exact value of ${\cal A}$, due to the finite-
iteration based approximation of the asymptotic equi-M set (some of the equi-M
points retained in our set representation may in reality escape upon further
iteration). Note that spacial resolution of the computation may limit
detection of thin filaments; these, however, have negligible area, hence they
are not expected to have a major contribution to the value of ${\cal A}$.

_Diameter $d_{h}$ along the real axis_, that is the distance between the cusp
and the end of the tail.

_Vertical diameter $d_{v}$_. We obtain rough estimates of this measure. This
may be in general problematic (due to the presence of thin filaments around
the boundary).

Figure 5: _Landmarks for Equi-M sets. Each panel shows one DTI-generated
equi-M set, with labels that facilitate computation of the measures proposed
in the text: A. $\gamma=1.12$, $\tau=-7.13$,$d_{h}=8.25$; B. $\gamma=0.82$,
$\tau=-6.2$, $d_{h}=7.02$ ; C. $\gamma=0.94$, $\tau=-6.2$, $d_{h}=7.14$._

###  2.4 Graph theoretical measures of the network architecture

For each subject’s network, we computed basic graph theoretical measures, in
particular on measures which can capture global features of the network
architecture into one or a few values. For this study, we focus in particular
on the following measures (applied to the adjacency matrix, and to the
weighted connectivity matrix).

_Degree centrality_ is the number of edges incident on each node. We will use
the mean degree centrality $\delta$ as a global network measure, describing
the information “flow,” where the mean is considered over all nodes in the
network. While in this project our networks are not oriented (hence we do not
differentiate between in and out degree centrality), this measure is otherwise
equivalent to the edge density considered in our previous studies on oriented
networks.

_Eigenspectrum_

_Clustering Coefficient_ is a measure of the degree to which nodes in a graph
tend to cluster together. For each node, its (local) clustering coefficient is
given by a proportion of the number of links between the vertices within its
neighborhood divided by the theoretical maximum number of links that can exist
between them. We used the mean of the node-wise clustering coefficients, as a
measure of overall clustering in the network.

_Betweenness centrality_ quantifies for each node the number of times a node
acts as a bridge along the shortest path between two other nodes
[freeman1977set]. For each connectome, we computed the mean betweenness
centrality $\beta$ over all nodes in the network.

_Eigenvector Centrality_

_Local Efficiency_

_Number of 3-node motifs_

_Number of 4-node motifs_

###  2.5 Correlation analysis

First, we aim to determine the strength of correlations within our set of
topological measures, and independently within our set of graph theoretical
measures. In either set, while it seems likely that not all the measures are
independent from one another, we are interested to see in what way and to what
extent these measures are inter-related.

Second, we investigate inter-correlations between the toplogical and the graph
theoretical measures. Since the topology of the M set can be viewed as a
symbolic representation of the network’s asymptotic dynamics, and since the
network architecture (underlying graph) is a factor in determining these
dynamics, correlations between properties of one and the other may help
understand which graph theoretical properties lead to which types of long-term
dynamics. Understanding this relationship at the level of simple quadratic
dynamics may be of crucial importance when extrapolating to the more complex,
neural dynamics that occurs in reality in these brain networks.

###  2.6 Gender based statistics

We investigated whether the hardwiring of the network (as captured by the set
of graph theoretical measures), and the quadratic network dynamics (as
captured by the collection of topological measures of the Mandelbrot set)
reflect gender differences in our subject set.

To study the relationship with gender (which was assessed for each subject on
a binary, M/F, basis), we use a two-way Kolmogorov-Smirnoff nonparametric
test. This allows us to assess whether each specific measure differs
significantly between the male and female subjects (the null hypothesis being
that they were extracted from the same overall distribution).

###  2.7 Generative modeling

When visually inspecting the equi-M sets for tractography-based graphs, one
immediate point of interest was the elaborate intertwining of unifying and
differentiating topological properties. On one hand, the sets all had common
basic features which are not generically assumed by all equi-M sets. For
example, although equi-M sets don’t have to be connected in general (in our
prior theoretical work we showed that examples of disconnected sets can occur
even in small, 3-node networks [mccrae2004contemplated]), our tractography-
generated equi-M sets turned out to be connected for all the 197 subjects in
our data-set. On the other hand, there are many subtle differences between
equi-M sets for different brain connectomes, as described in Section 2.3. We
aim to further investigate which of the graph theoretical properties of the
connectome are those that create the unifying trends, and which type of
variability between connectomes create the subtle differences in asymptotic
dynamics. To do this, we create surrogate networks, replicating specific graph
theoretical properties of the actual connectomes, and we investigate what
properties may be required to replicate the shape of tractography-based M sets
in a theoretical setup.

##  3 Results

###  3.1 Topological measures

As mentioned in the Methods section, we computed for each subject a set of
topological measures consisting of $(\gamma,\tau,{\cal A},d_{h},d_{v})$, that
use known topological landmarks of the M set to estimate its position and
size. Let us note first that these measures are not independent in the context
of our data. Figure 6a suggests that they are strongly correlated with each
other (e.g., a position of the cusp more to the right also corresponds to a
tail shifted to the left, larger diameters and a larger area).

The statistics (mean and standard deviation) of the topological measures
calculated for the M sets are described in Table 1, not only for the overall
data set, but also separately by the subjects’ gender (male versus female).
One can easily notice the significant gender differences (the M sets
corresponding to the male subgroup showing statistically larger diameters and
area than those corresponding to the female subgroup, with the cusp situated
more to the right, and the tail more to the left). Figure 6, panels b and c,
present a visual illustration of the differences between the mean M sets
computed over the male and respectively female subgroups.

| Cusp pos. $\gamma$ | Tail pos. $\tau$ | Area ${\cal A}$ | Real diam $d_{h}$ | Vert. diam $d_{v}$
---|---|---|---|---|---
| Mean | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std
Overall stats | 0.77 | 0.36 | -4.88 | 1.87 | 18.52 | 15.29 | 5.65 | 2.22 | 5.10 | 2.45
Male stats | 0.88 | 0.36 | -5.52 | 1.76 | 23.08 | 16.19 | 6.40 | 2.10 | 5.85 | 2.42
Female stats | 0.64 | 0.33 | -4.12 | 1.72 | 13.11 | 12.19 | 4.76 | 2.03 | 4.21 | 2.19
KS p-value | $<10^{-7}$ | $<10^{-8}$ | $<10^{-7}$ | $<10^{-7}$ | $<10^{-6}$
Table 1: Topological measures statistics and comparison between the Male and
Female subjects. For each measure (from left to right: the cusp position
$\gamma$, the tail position $\tau$, the total area, the diameter along the
real axis $d_{h}$, the maximum vertical diameter $d_{v}$) the table includes
the mean and standard deviation over the whole subject set (first row), as
well as mean and standard deviation over the Male group (second row) and over
the Female group, (third row, respectively). The last row represents the
$p$-value for the Kolmogorov-Smirnoff two-way statistics between the Male and
Female groups.

###  3.2 Graph theoretical measures

For each connectome, we computed the set of graph theoretical measures
described in the Methods section. As one would expect, not all these measures
represent independent assessments of the network structure, and some showed
significant correlations with others (as shown in Figure 6a). For example, it
is not surprising that the mean Degree Centrality is strongly correlated with
the mean eigenspectrum ($R=0.95$, $p<10^{-10}$), and the Clustering
Coefficient is strongly correlated with the Local Efficiency ($R=0.99$,
$p<10^{-19}$). Interestingly, there is also a consistent positive correlation
for mean number of motifs, between different types of 3 and 4-node motifs (the
red upper right 8 by 8 corner of the correlation square).

| Degree centr. | Eigensp. | Clustering Coeff. | Between Centr. | Eigen Centr. | Local Eff.
---|---|---|---|---|---|---
| Mean | Std | Mean | Std | Man | Std | Man | Std | Mean | Std | Mean | Std
Overall | 26.81 | 3.36 | 9.03 | 1.15 | 0.35 | 0.06 | 3.68 | 0.38 | 4.28 | 1.14 | 0.47 | 0.07
Male | 24.84 | 2.20 | 8.40 | 0.84 | 0.32 | 0.04 | 3.66 | 0.39 | 4.28 | 1.08 | 0.44 | 0.05
Female | 29.14 | 2.98 | 9.79 | 1.01 | 0.38 | 0.05 | 3.70 | 0.36 | 4.28 | 1.21 | 0.52 | 0.07
p-value | $<10^{-18}$ | $<10^{-16}$ | $<10^{-10}$ | – | – | $<10^{-11}$

| Motifs 3-9 | Motifs 3-13 | Motifs 4-113 | Motifs 4-126
---|---|---|---|---
| Mean | Std | Mean | Std | Man | Std | Man | Std
Overall stats | 2.8947 | 0.3860 | 1.4552 | 0.2584 | 37.9429 | 5.1461 | 11.8501 | 2.2029
Male stats | 2.6826 | 0.2719 | 1.3323 | 0.1964 | 35.3522 | 4.0317 | 10.7585 | 1.5965
Female stats | 3.1469 | 0.3481 | 1.6013 | 0.2469 | 41.0229 | 4.6018 | 13.1480 | 2.1191
KS p-value | $<10^{-12}$ | $<10^{-14}$ | $<10^{-10}$ | $<10^{-9}$
| Motifs 4-128 | Motifs 4-187 | Motifs 4-193 | Motifs 4-199
| Mean | Std | Mean | Std | Mean | Std | Mean | Std
Overall stats | 39.6685 | 7.3643 | 1.6916 | 0.3480 | 13.4973 | 3.1462 | 5.8630 | 1.8147
Male stats | 36.3657 | 5.8735 | 1.5757 | 0.2989 | 12.3780 | 2.7253 | 5.2784 | 1.5164
Female stats | 43.5951 | 7.0221 | 1.8294 | 0.3519 | 14.8281 | 3.0945 | 6.5582 | 1.8934
KS p-value | $<10^{-9}$ | $<10^{-4}$ | $<10^{-7}$ | $<10^{-8}$

Table 2: Graph theoretical measures statistics and comparison between the Male
and Female subjects. For each measure, the table includes the mean and
standard deviation over the whole subject set (first row), as well as mean and
standard deviation over the Male group (second row) and over the Female group
(third row, respectively). The last row represents the $p$-value of the
Kolmogorov-Smirnoff two-way statistics between the Male and Female groups.
Figure 6: _Correlations within and between graph theoretical and topological
measures. A. The panel represents in color the correlation values for all
pairs of measures for which the correlation value was significant (with a
p-value threshold $\alpha=0.05$). The color coding for sign and strength of
correlations are specified by the color bar. The pairs for which correlations
were nor significant are shown in white. On both axes, topological measures
are first labeled 1-5 (representing the following order: 1 = cusp position
$\gamma$, 2 = absolute value of the tail position $\tau$, 3 = area ${\cal A}$,
4 = diameter along real axis $d_{h}$, 5 = vertical diameter); then graph
theoretical measures are labeled 1-13 (representing the following sequence: 1
= degree centrality; 2 = eigenspectrum; 3 = clustering coefficient; 4 =
betweenness centrality; 5 = eigenvalue centrality; 6 = local efficiency; 7-8 =
numbers of 3-motifs; 9-13 = numbers of 4-motifs._

###  3.3 Inter-correlations between graph and topological measures

When looking to study the relationship between the topological and the graph
theoretical measures we also found strong correlations for most mutual pairs.
Higher Degree Centrality, Eigenspectrum and Local Efficiency correlated
significantly with smaller equi-M sets (landmarks closer to the origin,
smaller diameters and area). In contrast, higher Eigencentrality correlated
positively with the size of the equi-M set. Betweenness Centrality did not
show significant correlations. Interestingly, the number of motifs showed
consistent significant negative correlations with the size of the equi-M set
(the more motifs, the smaller the topological measures).

One of the limiting factors of this latter analysis is that all the
topological measures considered relate to the size of the equi-M set. More
subtle measures (representing properties such as fractal dimension of the
boundary, number of pinch points and pseudo-bulbs) may better and more
thoroughly underline the distinctions between the effects of different graph
structures on the dynamics. These concepts will be further defined and
explored in the Discussion section.

###  3.4 Gender based statistics

Table 2 shows that the values of some of the graph theoretical measures differ
significantly between males and females. This reflects into gender differences
in the geometry of the equi M set between the male and female groups (as
described in Table 1). Figure 7 additionally illustrates gender differences,
in the case of the area and horizontal diameter of the M set, by means of
histograms comparing these values between the male and female groups. Notice
the significantly lower mean area and diameter in the case of females (red)
versus their male counterparts (blue), and the relatively small overlap
between the two distributions (reflecting the small $p$ value associated with
the two-way Kolmogorov-Sinai test).

Figure 7: _Histograms of male versus female values, for two example
topological measures: A. The area ${\cal A}$ of the equi M set; B. the
horizontal diameter $d_{h}$ of the equi M set._

Below are two additional approaches to visualizing the statistical gender
differences between equi M sets, each highlighting different aspects at play.
First, we computed two group “mean” (or “prototypical”) connectomes, by
averaging the individual connectomes over each subject group (Males and
Females, respectively). The Male and Female prototype connectomes are shown in
Figure 8 (top row). We then computed the prototypical Male and respectively
Female equi-M sets from these two connectomes, as illusrtated in Figure 8
(bottom row). Notice that the two prototypical connectome are hardly
distinguishable with the naked eye; without the statistical instument at play,
some of the individual detail in connectivity patterns were averaged out,
resulting in two almost identical connectomes. However, these remaining,
seemingly unnoticeable differences were enough to produce visibly different
equi M prototypical sets for Males and Females (with a deeper cusp, and detail
around the cusp wider and more pronounced in the set corresponding to the
Female group than in the set corresponding to Males).

Figure 8: _Prototypical connectomes and equi M sets for the Male and Female
groups. Top. Prototypical connectomes for Male subjects (left) and Female
subjects (right), obtained by averaging the individual connectomes for the
respective subject groups. Bottom. Prototypical equi-M sets for Male (left)
and Female subjects (right), obtained from the prototypical connectomes for
the respective groups._

As a second approach, we took a more inclusive, statistical approach, by
creating stochastic equi-M sets for each subject group (Figure LABEL:). To
each point in the complex plane, we associated the number of subjects in the
group for which that point is in the equi M set.

Figure 9: _Correlations within and between graph theoretical and topological
measures. A. “Mean” equi-M set for male subjects. B. “Mean” equi-M set for
female subjects._

###  3.5 Generative modeling

We want to further investigate the effect of different patterns in the
connectome onto the resulting ensemble dynamics of the network. We want to
understand how much of the variability in dynamics can be easily tracked down
to measurable differences in the connectome, and to what extent the diversity
in the shapes of the M sets is the result of small, seemingly undetectable
variations in the connectome (as suggested by the results in the previous
section).

To illustrate the importance of each of these two aspects, we used a
generative model to create connectomes with graph theoretical properties
similar to those observed in the human tractography-derived connectomes in our
data-set.

Add results here.

Altogether, this section reinforces an important point: the geometry of the
equi M set is to a large extent encoded in measurable graph theoretical
properties of the connectome. Finer information on not just the size, but also
the distribution of the weights throughout the connectome is crucial to
producing the topology of the connectomes observed empirically. However, some
of the unifying as well as some of the distinguishing features between the
tractography-generated equi-M sets within out data seem to be driven by detail
in the connectome that is not traditionally included in a graph theoretical
analysis.

To obtain a better understanding of how these two aspects are tied, one may
have to (1) refine the connectome measures (to include specific information on
the distribution of weights) and (2) refine the topological assessment of
equi-M set, to include more specific measures of its topology. Both of these
directions are further discussed in our last section on Future owrk.

##  4 Discussion

###  4.1 Specific comments on the model

This first analysis of the topological trends in assymptotic sets of
structural brain networks is very promising. The equi-Mandelbot sets
corresponding to brain networks (1) are significantly different than the ones
obtained from surrogate data with similar basic connectivity properties such
as number of nodes, number of edges and average connection strengths; (2)
exhibit a set of common topological properties, some of them not generally
found in randomly generated equi-M sets (size, connectivity, robustness of the
cusp and tail structure); (3) exhibit vast detail variability, stemming in
graph configurational and connectivity differences (see Figure). A LOT MORE TO
ADD HERE.

###  4.2 Limitations

One intrinsic limitation comes from the fact that the methods are still
developing (hence their applicability and our understanding of their
significance will evolve with the theoretical and methodological progress).
For example, it is not yet don’t know whether an escape radius exists for all
networks. Thus far, our existing proofs and computations of escape radius only
apply to a few restricted network classes, which may not be comprise all our
data-generated networks. What this implies is that our numerically-generated
plots of equi-M sets may not be a reliable representation of the actual sets.
With knowledge of escape radius, a traditional computation of the Mandelbrot
set operates with positive errors (i.e., shows a larger set of points,
narrowing down to the true Mandelbrot set with increased spatial resolution
and with a higher number of iterates used). Without knowledge of escape radius
appears the potential for negative errors (i.e., missing a point which escaped
the test radius, but which may return close to the origin with later
iterates). As our proofs develop, together with optimization of the cost of
our numerical algorithms, equi-M sets will be computed more reliably and
efficiently, in higher resolution, permitting better estimation of the
topological invariants.

Finally, we briefly discuss an important question related to the inherent
limitation of tractography data, which only produces non-oriented, unsigned
connectomes. In reality, the direction of information propagation through the
network, and the signature of connections (excitatory versus inhibitory
projections) are very important to the emerging dynamics, and down to line to
the brain’s function and observed behavior. While sign information is missing
from our data, we focus on extending our analysis with a theoretical
exploration of potential effects of introducing different degrees and
distributions of signed projections. In our future work, we will also include
generalizations related to edge orientation. The study also had to abide
within data-related limitations. While convenient as a first step, using DTI
data will always present the issues related to the fact that the brain
networks are not oriented (making impossible to infer projection
directionality and strength), and connections are not signed (the thickness of
the track does not differentiate between excitatory and inhibitory
connections). Figure 10 illustrates the importance of this aspect, which will
be further addressed in future work, in conjunction with different imaging
modalities.

Figure 10: _Exploration of one random brain-derived network, artificially
enhanced to include inhibitory connections (negative signs were assigned
randomly to a fixed fraction of the nonzero edges). The M-sets exhibit
consistent topological differences related to the fraction of positive versus
negative signs, but also subtler distinctions in response to changing their
positions in the network architecture. A. All positive or all negative
weights. B. Roughly 1/4 positive weights. C. Roughly half positive weights. D.
Roughly 3/4 positive weights. _

###  4.3 Future directions

A primary direction in our future work will be to expand and refine this
analysis by including assitional assessment measures of the equi-M set
topology. Based on our first observations, we propose below two representative
candidates, which are expected to be able to efficiently differentiate between
equi-M shapes. However, substantial work needs to be invested in the direction
of better understanding and computing these measures, as well as the original
ones, with the accuracy needed for making valuable assessments,
classifications and predictions.

_Vertical narrowing and pinching._ The traditional Mandelbrot set has a
hyperbolic bulb structure. In particular, the bulbs that are centered along
the real axis are tangent to each other at the Feigenbaum points along the
axis, which we also call “pinch points” (because their removal would break the
set into connected components, changing the homeomorphism class). While for
equi-M sets the bulb structure is not typically preserved, some of the pinch
points persist – either as actual pinch points or just as a narrowing of the
set (which we call “narrow bridge”), delimiting “pseudo-bulbs.” We will use
the _number of pinch points $p$_ as a measure of bulb structure preservation.
We expect these to be easier to compute than the maximal vertical diameter in
general, since the narrowing points of the sets are not typically associated
with filaments (similarly with the seahorse valley in the traditional
Mandelbrot set).

Figure 11: _Counting pinch points for a DTI-generated equi-M set. The
algorithm detects the pinch points (purple) points where the M set narrows out
vertically to a point on the x-axis (virtually separating the connected
component obtained for the equi-M set after removing the x-axis)._

_Approximate entropy along the boundary._ A one-value measure which is related
to the presence of narrow bridges is the approximate entropy of the top
boundary, as described in Figure 11 (since it captures the regularity of the
vertical variations between high and low points along the top boundary of the
set). We compute the top of the equi-M set as follows:, for each real
coordinate along the x-axis reached by the equi-M set, we consider the highest
point such that its vertical projection to the axis is contained in the M set.
We then compute the approximate entropy of this function, using a simple
algorithm by Pincus , effective for short time series.

A desirable and feasible future direction is to unify the results into a more
statistically comprehensive approach. We will aim to select a set of
topological measures that reflect the differences in the brain-generated
M-sets, as well as their relationship with the graph architecture. We will
investigate statistically which combination of graph theoretical measures and
topological measures delivers the most consistent clusters across the two. The
optimal combination of measures will be henceforth used as our _topological
profiling toolbox_ (TPT) for the data-set at hand. We will then further
validate our TPT by creating additional surrogate networks within the same
graph-theoretical properties as in the data, and verifying whether they fall
into the TPT-predicted clusters. If they consistently fail to do so, the
previous step may need to be revisited, to make modifications to the TPT. We
will then further validate the procedure of selecting a TPT and of using it
for classification purposes with a separate human DTI data set, with larger
networks (hence finer brain region masks).

To address the DTI related limitation in our future work implementing our TPT
framework, we propose expanding it to imaging data that provides information
on the connection orientation and sign. For example, Granger Causality, or
Dynamic Causal Modeling can be used on functional MRI data to inform a model.
Once validated at this level, our framework can be used to produce a more
accurate “brain profile,” encoding the efficiency in each individual’s brain
regulation.This is the long-term aim of or proposed research. The steps,
including troubleshooting, can be accomplished within the five years of the
award, addressing the PI’s original question, which lead her down this
research path in the first place. Our overarching goal is working on
developing TPT-like lookup chart towards a neurobiologically-based diagnostic
instrument that can be used by clinicians to optimize classification,
prediction and intervention techniques, as stated in Broader Impacts.

## References

  * [1] Fatou P, 2019.  Sur les équations fonctionnelles.  Bulletin de la Société mathématique de France. 47: 161-271.
  * [2] Julia G, 1918.  Mémoire sur l’itération des fonctions rationnelles.  J. Math. Pures. Appl. 47–246.
  * [3] Branner B, 1989.  The mandelbrot set.  In Proc. Symp. Appl. Math. 39: 75-105.
  * [4] Branner B, Hubbard JH, 1992.  The iteration of cubic polynomials Part II: Patterns and parapatterns.  Acta mathematica. 169(1): 229-325.
  * [5] Qiu WY, Tin YC, 2009.  Proof of the Branner-Hubbard conjecture on Cantor Julia sets.  Science in China Series A: Mathematics. 52(1): 45.
  * [6] Devaney RL, Look DM, 2006.  A criterion for Sierpinski curve Julia sets.  Topology Proceedings. 30: 163-179.
  * [7] Bullmore E, Sporns O, 2009.  Complex brain networks: graph theoretical analysis of structural and functional systems.  Nature Reviews Neuroscience. 10(3): 186–198.
  * [8] Sporns O, 2002.  Graph theory methods for the analysis of neural connectivity patterns.  Neuroscience databases: A practical guide. 171–186.
  * [9] Sporns O, 2011.  The non-random brain: efficiency, economy, and complex dynamics.  Frontiers in Computational Neuroscience. 5: 5.
  * [10] Van Essen DC, Smith SM, Barch DM, Behrens TEJ, Yacoub E, Ugurbil K, 2013.  The WU-Minn Human Connectome Project: An overview.  NeuroImage 80; 62-79.
  * [11] Rǎdulescu A, Pignatelli A, 2016.  Real and complex behavior for networks of coupled logistic maps.  Nonlinear Dynamics. 87(2); 1295–1313.
  * [12] Rǎdulescu A, Evans S, 2018  Asymptotic sets in networks of coupled quadratic nodes.  Journal of Complex Networks. 7(3): 315–345.

## Appendix A

### Statistics for topological measures

### Statistics for graph theoretical measures

| Degree centr. | Eigensp. | Clustering Coeff. | Between Centr. | Eigen Centr. | Local Eff.
---|---|---|---|---|---|---
| Mean | Std | Mean | Std | Man | Std | Man | Std | Mean | Std | Mean | Std
Overall stats | 0.2681 | 0.0336 | 0.0903 | 0.0115 | 0.0035 | 0.0006 | 0.0368 | 0.0038 | 0.0428 | 0.0114 | 0.0047 | 0.0007
Male stats | 0.2484 | 0.0220 | 0.0840 | 0.0084 | 0.0032 | 0.0004 | 0.0366 | 0.0039 | 0.0428 | 0.0108 | 0.0044 | 0.0005
Female stats | 0.2914 | 0.0298 | 0.0979 | 0.0101 | 0.0038 | 0.0005 | 0.0370 | 0.0036 | 0.0428 | 0.0121 | 0.0052 | 0.0007
KS p-value | $<10^{-18}$ | $<10^{-16}$ | $<10^{-10}$ | – | – | $<10^{-11}$
Corr. with age | R-val | p-val | R-val | p-val | R-val | p-val | R-val | P-val | R-val | p-val | R-val | p-val
| 0.2198 | 0.0019 | 0.1561 | 0.0285 | – | – | -0.1492 | 0.0364 | 0.1935 | 0.0064 | – | –
Table 3: Title | Motifs 3-9 | Motifs 3-13 | Motifs 4-113 | Motifs 4-126 | Motifs 4-128 | Motifs 4-187 | Motifs 4-193 | Motifs 4-199
---|---|---|---|---|---|---|---|---
| Mean | Std | Mean | Std | Man | Std | Man | Std | Mean | Std | Mean | Std | Mean | Std | Mean | Std
Overall stats | 2.8947 | 0.3860 | 1.4552 | 0.2584 | 37.9429 | 5.1461 | 11.8501 | 2.2029 | 39.6685 | 7.3643 | 1.6916 | 0.3480 | 13.4973 | 3.1462 | 5.8630 | 1.8147
Male stats | 2.6826 | 0.2719 | 1.3323 | 0.1964 | 35.3522 | 4.0317 | 10.7585 | 1.5965 | 36.3657 | 5.8735 | 1.5757 | 0.2989 | 12.3780 | 2.7253 | 5.2784 | 1.5164
Female stats | 3.1469 | 0.3481 | 1.6013 | 0.2469 | 41.0229 | 4.6018 | 13.1480 | 2.1191 | 43.5951 | 7.0221 | 1.8294 | 0.3519 | 14.8281 | 3.0945 | 6.5582 | 1.8934
KS p-value | $<10^{-12}$ | $<10^{-14}$ | $<10^{-10}$ | $<10^{-9}$ | $<10^{-9}$ | $<10^{-4}$ | $<10^{-7}$ | $<10^{-8}$
Corr. with age | R-val | p-val | R-val | p-val | R-val | p-val | R-val | P-val | R-val | p-val | R-val | p-val | R-val | p-val | R-val | p-val
| 0.2783 | $<10^{-4}$ | 0.2626 | $<10^{-3}$ | 0.2914 | $<10^{-4}$ | 0.2400 | $<10^{-3}$ | 0.2801 | $<10^{-4}$ | 0.2029 | 0.0042 | 0.2367 | $<10^{-3}$ | 0.2325 | 0.0010
Table 4: Title