# A Geometric Chung Lu model and the Drosophila Medulla connectome

Susama Agarwala111Johns Hopkins University, Applied Physics Lab; Laurel, MD
20723<EMAIL_ADDRESS>
Franklin Kenter222Mathematics Department, United States Naval Academy;
Annapolis, MD 21402<EMAIL_ADDRESS>

###### Abstract

Many real world graphs have edges correlated to the distance between them,
but, in an inhomogeneous manner. While the Chung-Lu model and the geometric
random graph models both are elegant in their simplicity, they are
insufficient to capture the complexity of these networks. In this paper, we
develop a generalized geometric random graph model that preserves many graph
theoretic aspects of these real world networks. We test the validity of this
model on a graphical representation of the Drosophila Medulla connectome.

##  1 Introduction

In many networks, the presence of an edge between two nodes strongly
correlates with the spatial distance between them; this includes
infrastructure networks, transportation networks [2] and connectomes [17, 20].
Traditionally, scale free networks are used to model these systems, but there
is increasing evidence that this is too simplistic a paradigm [5]. For
instance, in large scale connectome networks of regions of the human brain, it
has been shown that the empirical probability that an edge exists between two
nodes is a function of distance which decays faster than any inverse power
decay[20].

To address is shortcoming, we build a composite model for that is a
combination of two prototypical random graph models: the random geometric
graph model and the Chung-Lu model. We validate this composite model on the
connectome formed by tracing the neurons in a small portion of the Drosophila
Medulla[26]. In doing so, we find that we are able to generate a family of
graphs that preserve many of the network properties observed in the original
connectome in a way that could not have been possible just using one or the
other of the two component models.

There is a lot of evidence in favor of spatial and geometric models for the
brain in meso- and macro-scale connectome work [6, 15, 22] that transcends
species. However, continuing this work on a micro, i.e. a cellular level has
been stymied by the difficulty and expense of obtaining the data. To address
this issue of sparse examples of mapped connectomes, we propose a means of
constructing a composite model that generates synthetic graphs which share
many of the graph characteristics of the graph defined by the Drosophilla
medulla [23].

In the random geometric graph model, vertices are randomly embedded into some
space with a metric (usually Euclidean); after which, two vertices $i,j$ are
adjacent with a probability connection function based upon the spatial
distance between them $d_{ij}$. The simplest case is when an edge is present
whenever $d_{ij}<R$ for some threshold value $R$ [18]. Though more complicated
connection functions can be used depending on the application [19, 10]. One of
the main shortcomings of random geometric graph models is the homogeneity of
the vertices which stands in contrast to most real-world networks [11]. The
Chung-Lu model [8] is a non-geometric heterogeneous model where each vertex is
randomly assigned a weight $w_{i}$ from a predetermined distribution. Then,
the probability that an edge exists between two vertices is a function of the
weights $w_{i}$ and $w_{j}$. Under mild conditions, the expected degrees of
the vertices are precisely $w_{i}$, allowing for heterogeneity [8]. This model
allows straightforward analysis of a general degree distribution, especially
those corresponding to “power-law” or “scale-free” networks [8]. However, it
has been observed that the Chung-Lu model tends to underestimate triangle
counts in real-world graphs, and miss estimates degrees for real world
networks that do not satisfy the Chung-Lu constraint [27]. The connectome
studied in this paper (and connectomes in general) has heterogeneous vertices,
as well as a geometrically defined connection function, making it an ideal
candidate for a composite of these two models.

There are several previously-used methods for combining geometric graphs with
the Chung-Lu model. One method is to use a hyperbolic embedding where the
connection function $f(d_{ij},w_{i},w_{j})$ is a function of the hyperbolic
distance as well as the two weights. The advantage of using hyperbolic space
is that the resulting clustering is similar to that found in social networks
[25]. In this model, as the number of vertices $\to\infty$, the limit of the
hyperbolic model is the classical Chung-Lu model [12]. The hyperbolic model
can be generalized using a toroidal embedding in different dimensions [4].
However, in all cases, the connection function is proportional to
$d_{i,j}^{-\beta}$ for some $\beta$ (i.e. an inverse power law) [4]. This
stands in contrast to the rate of decay for real-world graphs mentioned
previously.

To address these concerns, in Section 3.1, we introduce a simple yet robust
heterogeneous geometric random graph model that is able to consider both the
spatial position of the vertices as well as a randomly assigned weight. In the
hyperbolic model above, the probability that two nodes are connected, given
their distance is prescribed by a simple power law (namely
$d_{i,j}^{-\beta}$). In our model, we generalize this to allow for any
function that limits to the edge density of the graph as the distance goes to
$\infty$ and to the self loop number as the distance goes to $0$.

Then, in Section 3.2, we adapt this model to consider the boundary effects
when the spatial coordinates are restricted to a (usually non-convex) region.
Together, this allows us to build synthetic connectomes resulting in similar
local and global topological properties as the original as well as maintaining
the original spatial characteristics and spectral characteristics. In Section
4, we apply our model to part of the connectome of the Drosophila medulla
given in [23] wherein synaptic connections are detailed in their exact spatial
coordinates in three dimensions. In Section 5, we use a variety of network
measures, we show that synthetic connectomes generated under our model are
more similar to those created under a pure Chung-Lu model. Finally, Section 6,
we compare the accuracy of our model to other existing models.

##  2 Chung-Lu model and existing generalizations

For our purposes, the graphs we consider, $G=(V,E)$, are finite and
undirected; we do not allow for multiedges but we allow for graphs to have
self-loops. Throughout, we fix the number of vertices in the graph, $n=|V|$.
We write $i\sim j$ to denote that $i$ is adjacent to $j$. A self-loop occurs
whenever $i\sim i$. The degree of a vertex is the number edges incident to it,
denoted $\deg_{i}$; for our purposes, self-loops only contribute 1 toward the
degree of a vertex. The edge-density of a graph is
$\varepsilon=\frac{2|E|}{(n+1)n}$.

We work with graphs that are embedded in some geometric space, therefore each
vertex has an associated position. We denote the Euclidean distance between
two vertices, $i$ and $j$, as $d_{i,j}$. Note that this is the distance
between two vertices relative to the geometry the graph is embedded in, not
the graph-theoretic distance.

###  2.1 The Original Chung-Lu Model

The Chung-Lu model [8] is a heterogeneous random graph model whose parameter
is a single vector of expected degrees: $\ \vec{w}$ of length $n$. Given this
vector, a graph can be generated where the probability that two vertices (not
necessarily distinct) are adjacent is

$P(i\sim j|w_{i},w_{j})=\frac{w_{i}w_{j}}{\sum_{k}w_{k}}\;,$

and the event $i\sim j$ is independent from any other pair. This probability
is valid for all pairs $i,j$ if

$\displaystyle(\max_{\ell}w_{\ell})^{2}<\sum_{k}w_{k}$ (1)

We refer to equation 1 as the “Chung-Lu condition”.

###  2.2 Previous Geometric Chung-Lu Models

Bringmann, Keusch and Lengler introduced a Geometric Chung-Lu Model, called
Geometric Inhomogeneous Random Graphs (GIRG), [4] which generalizes a previous
geometrized Chung-Lu model, the hyperbolic model [12]. In this model, $n$
vertices are each assigned a random point chosen uniformly in the
$d$-dimensional unit torus, and each vertex is assigned a random weight from
$\ \vec{w}$ of length $n$. The connection function (i.e., the probability of
an edge $i\sim j$) is given by

$P(i\sim j|d_{ij},w_{i},w_{j})=\min\left\\{\frac{1}{d_{i,j}^{\alpha
d}}\cdot\left(\frac{w_{i}w_{j}}{W}\right)^{d},1\right\\}$

where $d_{i,j}$ is taken to be the $\infty$-norm on the torus and
$W=\sum_{k}w_{k}$. As before, this probability is independent for each pair.
Unlike the traditional Chung-Lu model, the vector $\vec{w}$ is not necessarily
the expected degree of each vertex. Indeed, under mild conditions for $d$,
$\alpha$ and $\vec{w}$, the probability that two randomly chosen vertices are
adjacent given their distance (but without knowing the $w_{i}$ or $w_{j}$) is

$P(i\sim j|d_{ij})=\frac{1}{d_{ij}^{\alpha d}}\cdot\frac{1}{n(n-1)}\cdot
W^{d}$

In particular, regardless of $\vec{w}$, $d$ and $\alpha$, the edge probability
decreases according to an inverse power law, $Cd_{ij}^{-\beta}$. However, this
is not necessarily the case for real-world networks, for instance in meso-
scale connectome networks [20, 15, 17]. In Section 6, we show that this is
also not true for the micro scale connectome of the Drosophilla Medulla
studied here.

In contrast, in the model presented in the paper, we relax this restriction by
allowing for a near arbitrary function with the potential of exponential
decay.

##  3 A Generic Geometric Chung-Lu Model

In this section, develop an generic geometric Chung-Lu Model. Unlike the
geometric inhomogeneous random graphs above, our connection function will
decay at an exponential rate. We call this a generic model because we relax
two fundamental assumptions of previous models to accommodate real world
structures. Namely, we do not require the nodes to either be distributed
uniformly over a space or have any periodic properties.

We do this in two stages. In section 3.1, we define a mathematically simple
model where vertices are embedded into a $d$-dimensional torus under some
distribution $\mu$. Most real-world graphs are not periodic in nature and have
boundary considerations. Therefore, in section 3.2, we restrict our graph to
one that is defined on a subset of the torus
$\mathcal{B}\subset\mathbb{R}^{d}/\mathbb{Z}^{d}$, with all vertices (and
edged between vertices) outside $\mathcal{B}$ removed, and all edges between
vertices inside and outside of $\mathcal{B}$ removed. Unlike the traditional
Chung-Lu model, or even the geometric inhomogeneous random graphs model, the
boundary condition implies that the vector of weights $\vec{\rho}$ does not
correspond directly with the expected degree of each node. Hence, the boundary
conditions requires a careful re-weighting of the vertices to find the
appropriate weights.

###  3.1 A Geometric Chung-Lu Model without uniformity

The following are the inputs for the generic geometric Chung-Lu model on a $d$
dimensional torus.

######  Definition 3.1 (Parameters of the Model).

  1. 1.

The dimension, $d$.

  2. 2.

The number of vertices, $n$.

  3. 3.

A probability measure. A measure on $\mathbb{R}^{d}/\mathbb{Z}^{d}$, $\mu$,
describing the distribution of nodes.

  4. 4.

One parameter family of connection probabilities. We define a function that
gives a joint probability that two vertices are adjacent and within a certain
distance from each other. We call this function

$\displaystyle F_{1}(x)=P(i\sim j,d_{i,j}\leq x)$ (2)

. We will place nontrivial constraints on $F_{1}$ that we describe later in
Definition 3.2.

  5. 5.

Intensity vector. A vector of intensities, $\vec{\rho}$, whose length is the
number of nodes $n$.

As with the Chung-Lu model, the vector $\vec{\rho}$ will, under mild
conditions, be the expected degree of the nodes. However, when we introduce
boundary considerations, this will no longer be the case. To emphasize this
distinction, we refer to the vector as a vector of “intensities” as opposed to
“weights.”

To generate a graph under this model, we first place $n$ vertices
independently on a $d$-dimensional torus, $\mathbb{R}^{d}$, according to the
distribution $\mu$.

Once the $n$ vertices have been assigned a position on the torus, one randomly
associates to each vertex an element of the intensity vector without
replacement. (Throughout, we will mildly abuse notation and represent this by
writing $P(\rho_{i})$ which is $\frac{1}{n}$).

Given these inputs, we next construct the connection function for the model.
In order to do so, we need to impose some additional assumptions and
constraints. In particular, we need to ensure that the expected edge densities
from the Chung-Lu and geometric aspects of the model agree. This, and other
necessary assumptions and constraints are given below.

######  Definition 3.2 (Additional Constraints and Assumptions).

  1. 1.

Chung-Lu generalization We impose that the probability that two nodes are
connected given their intensity similar to the Chung-Lu model

$\displaystyle P(i\sim
j|\rho_{i},\rho_{j})=\min\left\\{\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}},1\right\\}\;.$
(3)

  2. 2.

Compatible expected edge densities The expected edge density of the graph from
the geometric model is given by
$\varepsilon_{geom}=\lim_{x\to\infty}F_{1}(x)$. To get the expected edge
density from the Chung-Lu generalization note that

$\displaystyle\mathbb{E}(\deg_{i})=\sum_{j}P(i\sim
j|\rho_{i},\rho_{k})=\sum_{j}\min\\{\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}},1\\}\;.$

Therefore, the edge density of the graph is

$\displaystyle\varepsilon_{CL}=\frac{\sum_{k}\mathbb{E}(\deg_{k})}{n^{2}}\;.$

The expected edge densities computed by the geometric and by the Chung-Lu
models must coincide:

$\displaystyle\varepsilon_{geom}=\varepsilon_{CL}=\varepsilon\;.$

We denote by $\varepsilon$ edge density of the model.

  3. 3.

Independence The distances between two vertices are independent of their
intensities:

$\displaystyle
P(d_{i,j}<x,\rho_{i},\rho_{j})=P(d_{i,j}<x)P(\rho_{i},\rho_{j}).$ (4)

Furthermore, we require that the distances between two vertices and the
expected degrees of the graph remain independent, given that they are
adjacent:

$\displaystyle P(d_{i,j}<x,\rho_{i},\rho_{j}|i\sim j)=P(d_{i,j}<x|i\sim
j)P(\rho_{i},\rho_{j}|i\sim j).$ (5)

Note that the Chung-Lu generalization in Definition 3.2 is a generalization of
the Chung-Lu condition in Section (1). This is the same generalization of this
condition adopted by [4]. Namely, we allow the intensities to be such that the
value of $\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}}$ may be greater than 1.
When this occurs, we set the probability that the two nodes are connected to
$1$.

Under the conditions laid out in Definition 3.2 we may calculate the
probability $P(d_{i,j}<x)$ for any two randomly chosen points $i,j$ given
their positions in $\mathbb{R}^{d}/\mathbb{Z}^{d}$. Call this function
$F_{2}(x)$ given by

$\displaystyle
F_{2}(x)=P(d_{ij}<x)=\idotsint_{\|\mathbf{y}-\mathbf{z}\|_{2}<x}(\mu\times\mu)(\mathbf{y},\mathbf{z})~{}dy_{1}\ldots
dy_{d}~{}dz_{1}\ldots dz_{d}$ (6)

where $(\mu\times\mu)$ is the product measure. Note that the derivative of
$F_{2}(x)$ can be interpreted as a probability density function, while the
derivative of $F_{1}(x)$ cannot. Namely, the derivative $F_{1}^{\prime}(x)$
gives the probability that two nodes are connected and distanced exactly $x$
apart: $F_{1}^{\prime}(x)=P(i\sim j,d_{i,j}=x)$.

We are now ready to build the connection function for this model.

######  Lemma 3.3.

We may find the probability that two nodes are connected given that they are a
specific distance apart as

$\displaystyle P(i\sim
j|d_{i,j}=x)=\frac{F^{\prime}_{1}(x)}{F^{\prime}_{2}(x)}\;.$

###### Proof.

First note that by Bayes rule, (2) and (6), we may write the probability that
two nodes are connected given that they are within a given distance as

$\displaystyle P(i\sim j|d_{i,j}<x)=\frac{F_{1}(x)}{F_{2}(x)}\;.$ (7)

One may calculate the probability that two nodes are connected given that they
are exactly a fixed difference apart by taking the derivatives of the
numerator and denominator of the above display:

$\displaystyle\frac{F_{1}^{\prime}(x)}{F_{2}^{\prime}(x)}$
$\displaystyle=\lim_{h\to 0}\frac{P(i\sim j,x<d_{ij}<x+h)}{P(x<d_{ij}<x+h)}$
(8) $\displaystyle=\lim_{h\to 0}P(i\sim j|x<d_{i,j}<x+h)=P(i\sim j|d_{i,j}=x)$
(9)

as desired. ∎

######  Remark 3.4.

Note that there is some delicacy in how we define the function $F_{1}(x)$ and
its derivatives. Namely, we set it uniformly to $0$ for $x<0$. However, if
there are self loops, this function is no longer continuous. We have

$\displaystyle F_{1}(x)=\begin{cases}0&x<0\\\
\frac{\mathbb{E}(\\#\;\textrm{self loops})}{n^{2}}&x=0\\\ P(i\sim
j,d_{i,j}<x)&x>0\end{cases}$

When $x=0$, we have

$\displaystyle F_{1}(0)=\frac{P(i\sim
i)}{n}=\frac{\mathbb{E}(\\#\;\textrm{self loops})}{n^{2}}\;.$

As a result, under mild conditions on the model, the function $F_{1}(x)$ is
continuous on and differentiable on $(0,\infty)$ with
$\lim_{x\to\infty}F_{1}(x)=\varepsilon.$ We define

$\displaystyle F^{\prime}_{1}(0)=\frac{\mathbb{E}(\\#\textrm{self
loops})}{n^{2}}\delta(0)$

where $\delta(\cdot)$ is the Dirac-delta function. Thus,
$\int_{0}^{\infty}F^{\prime}_{1}(x)~{}{dx}=\varepsilon$.

We may now use Lemma 3.3 and the fact that the distance between the nodes is
not only independent of their intensity, but that both are conditionally
independent on whether the nodes are adjacent to define the connection
function.

######  Theorem 3.5.

The connection function for this model is

$\displaystyle P(i\sim
j|d_{i,j}=x,\rho_{i},\rho_{j})=\min\\{\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}},1\\}\frac{F_{1}^{\prime}(x)}{F^{\prime}_{2}(x)}\frac{1}{\varepsilon}\;,$

###### Proof.

By Bayes rule, and assumption (5), we may write

$\displaystyle P(i\sim
j|d_{i,j}=x,\rho_{i},\rho_{j})P(d_{i,j}=x,\rho_{i},\rho_{j})$
$\displaystyle=P(d_{i,j}=x,\rho_{i},\rho_{j}|i\sim j)P(i\sim j)$
$\displaystyle=P(d_{i,j}=x|i\sim j),P(\rho_{i},\rho_{j}|i\sim j)P(i\sim j)$
$\displaystyle=\frac{P(i\sim j|d_{i,j}=x)P(i\sim
j|\rho_{i},\rho_{j})P(d_{i,j}=x)P(\rho_{i},\rho_{j})}{P(i\sim j)}\;,$

or by (4)

$\displaystyle P(i\sim j|d_{i,j}=x,\rho_{i},\rho_{j})=\frac{P(i\sim
j|d_{i,j}=x)P(i\sim j|\rho_{i},\rho_{j})}{P(i\sim j)}\;.$

Equations (3) and (9) give expressions for $P(i\sim j|d_{i,j}=x)$ and $P(i\sim
j|\rho_{i},\rho_{j})$ respectively. By construction $P(i\sim j)=\varepsilon$.
Combining these, we have

$\displaystyle P(i\sim
j|d_{i,j}=x,\rho_{i},\rho_{j})=\min\\{\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}},1\\}\frac{F_{1}^{\prime}(x)}{F^{\prime}_{2}(x)}\frac{1}{\varepsilon},$

as desired.

∎

Note that if we remove the Chung-Lu generalization, i.e. demand that
$\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}}\leq 1$, we see the the intensities
are exactly the expected degree of a vertex, given the associated intensity.

######  Lemma 3.6.

Restricting to the case where the Chung-Lu condition holds, i.e. to when
$\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}}\leq 1$ for all $i,j$, we have that
for any vertex $i$, the expected degree of the vertex is equal to its given
intensity:

$\mathbb{E}(deg_{i}|\rho_{i})=\rho_{i}$

###### Proof.

Note that by definition

$\displaystyle\mathbb{E}(deg_{i}\;|\;\rho_{i})=\sum_{j}{P}(i\sim j|\rho_{i})$

which, by the law of total probability can be rewritten

$\displaystyle\mathbb{E}(deg_{i}\;|\;\rho_{i})=\sum_{j}\sum_{\ell}{P}(\rho_{j}=\rho_{\ell}){P}(i\sim
j\;|\;\rho_{i},\rho_{\ell})\;.$

Similarly, we may write ${P}(\rho_{j}=\rho_{\ell}){P}(i\sim
j\;|\;\rho_{i},\rho_{\ell})=\int_{0}^{\infty}F_{2}^{\prime}(x){P}(\rho_{j}=\rho_{\ell}\;|\;d_{ij}=x){P}(i\sim
j\;|\;\rho_{i},\rho_{\ell},d_{ij}=x)~{}{dx}$. Therefore, we may write

$\displaystyle\mathbb{E}(deg_{i}\;|\;\rho_{i})=\sum_{j}\sum_{\ell}\int_{0}^{\infty}F_{2}^{\prime}(x){P}(\rho_{j}=\rho_{\ell}\;|\;d_{ij}=x){P}(i\sim
j\;|\;\rho_{i},\rho_{\ell},d_{ij}=x)~{}{dx}\;.$

Note that since the $\rho_{j}$ are distributed randomly amongst the vertices,
${P}(\rho_{j}=\rho_{\ell}\;|\;d_{ij}=x)={P}(\rho_{j}=\rho_{\ell})=\frac{1}{n}$.
Combined with Theorem 3.5 we rewrite this

$\displaystyle\mathbb{E}(deg_{i}\;|\;\rho_{i})=\sum_{j}\sum_{\ell}\int_{0}^{\infty}F_{2}^{\prime}(x)\frac{1}{n}\min\\{\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}},1\\}\frac{F_{1}^{\prime}(x)}{F^{\prime}_{2}(x)}\frac{1}{\varepsilon}\;.$

Since there is not $\ell$ in this statement anymore, by the assumption that
$\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}}\leq 1$, this simplifies to

$\displaystyle\mathbb{E}(deg_{i}\;|\;\rho_{i})=\frac{1}{\varepsilon}\sum_{j}\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}}\int_{0}^{\infty}F_{1}^{\prime}(x)~{}{dx}\;.$

By Remark 3.4, we see that
$\int_{0}^{\infty}F^{\prime}_{1}(x)~{}{dx}=\varepsilon$. In other words,

$\displaystyle\mathbb{E}(deg_{i}\;|\;\rho_{i})=\sum_{j}\frac{\rho_{i}\rho_{j}}{\sum_{k}\rho_{k}}=\rho_{i}$

as desired.

∎

###  3.2 The generic geometric Chung-Lu model: without uniformity and with
boundaries

While the generic model in Section 3.1 is mathematically convenient, the
toroidal nature is not present in real-world networks. To address this, we
modify the original model by additionally prescribing an addition parameters.

######  Definition 3.7.

In addition to the parameters given in Definition 3.1, we add an additional
input, a region within the torus:

  1. 6.

A region $\mathcal{B}\subset\mathbb{R}^{d}/\mathbb{Z}^{d}$.

This restriction to a region induces several changes to the parameters given
in Definition 3.1.

First, we need to change the graph itself, as well as the measure on this
space. To generate a random graph, we select $n$ points within $\mathcal{B}$
using the measure $\frac{1}{\mu(\mathcal{B})}\mu$ and generate edges according
to the probabilities Theorem 3.5. We then, consider the resulting induced
subgraph of the vertices in $\mathcal{B}$. We call the resulting network
$G_{\mathcal{B}}$.

Second, imposing a region $\mathcal{B}$ changes the previously defined
probabilities $F_{1}(x)$ and $F_{2}(x)$.

######  Definition 3.8.

Let $F_{\mathcal{B},1}(x)=P(i\sim j,d_{i,j}\leq x|i,j\in\mathcal{B})$ and
$F_{\mathcal{B},2}(x)d_{i,j}\leq x|i,j\in\mathcal{B})$.

As in Lemma 3.3, we write

$\displaystyle\frac{F_{\mathcal{B},1}(x)}{F_{\mathcal{B},2}(x)}=P(i\sim
j|d_{ij}<x\;;i,j\in\mathcal{B})\;.$

Then we may relate this conditional probability to the conditional probability
in the generic model, $P(i\sim j|d_{ij}<x)$.

To understand the bias between $F_{1}(x)$ and $F_{\mathcal{B},1}(x)$ and
$F_{2}(x)$ and $F_{\mathcal{B},2}(x)$, we define two functions,
$G_{1}(x)=P(i,j\in\mathcal{B}\;|\;i\sim j,d_{ij}<x)$, and
$G_{2}(x)=P(i,j\in\mathcal{B}\;|\;d_{ij}<x)$

######  Theorem 3.9.

$\frac{F_{\mathcal{B},1}^{\prime}(x)}{F_{\mathcal{B},1}^{\prime}(x)}=\frac{G_{1}^{\prime}(x)}{G_{2}^{\prime}(x)}\frac{F_{1}^{\prime}(x)}{F_{2}^{\prime}(x)}.$

###### Proof.

This proof follows from the arguments presented in Lemma 3.3. Namely

$\displaystyle\frac{F_{\mathcal{B},1}^{\prime}(x)}{F_{\mathcal{B},2}^{\prime}(x)}$
$\displaystyle=P(i\sim j|d_{ij}=x\;;i,j\in\mathcal{B})$ (10)

while by a similar calculation,

$\displaystyle\frac{G_{1}^{\prime}(x)}{G_{2}^{\prime}(x)}$
$\displaystyle=\frac{P(i,j\in\mathcal{B}\;|\;i\sim
j,d_{ij}=x)}{P(i,j\in\mathcal{B}\;|d_{ij}=x)}\;.$

The result follows from Bayes rule.

∎

The ratio $G_{1}^{\prime}(x)/G_{2}^{\prime}(x)$ measures bias introduced to
the conditional probability
$\frac{F^{\prime}_{1}(x)}{F^{\prime}_{2}(x)}=P(i\sim j|d_{i,j}=x)$ caused by
the shape of $\mathcal{B}$. If the ratio
$\frac{G_{1}^{\prime}(x)}{G_{2}^{\prime}(x)}$ is approximately 1 for all $x$,
that is, if the probability that two nodes are connected and the probability
that they are both in $\mathcal{B}$ are conditionally independent on their
distance, then

$\frac{F_{\mathcal{B},1}^{\prime}(x)}{F_{\mathcal{B},2}^{\prime}(x)}\approx\frac{F_{1}^{\prime}(x)}{F_{2}^{\prime}(x)}.$

Finally, by restricting the region, one must also change the intensities of
the vertices accordingly.

######  Definition 3.10.

Let $\vec{\rho}_{\mathcal{B}}$ be the vector of intensities on the graph
$G_{\mathcal{B}}$. Namely, as a generalization of Lemma 3.6,
$\vec{\rho}_{\mathcal{B}}=\mathbb{E}(deg_{i}^{\mathcal{B}}|\rho_{i},d_{i,j})$.

We may relate this quantity to the original vector of intensities. However, in
order to do so, we first need to define an expected valency of each vertex
given its position in the graph, using Lemma 3.3.

######  Definition 3.11.

As above, write $d_{i,k}$ to indicate the distance between the vertex $i$ and
$k$. Let $\omega_{i}=\sum_{j}P(i\sim j|d_{i,j},i,j\in\mathcal{B})$ be the
geometric weight of vertex $i$, or the expected degree of the vertex, given
its position relative to all other vertices.

Note that $\omega_{i}$ is a real number associated to each vertex, and not a
function of $x$.

######  Lemma 3.12.

Let $d_{i,j}$ be the distance between vertex $i$ and vertex $j$, while
$\deg_{j}$ is the degree of the vertex $j$. Then the geometric weight of a
vertex given its position can be expressed

$\displaystyle\omega_{i}=\sum_{j}\frac{F_{\mathcal{B},1}^{\prime}(d_{i,j})}{F_{\mathcal{B},2}^{\prime}(d_{i,j})}\;.$

###### Proof.

Summing both sides of equation (10), we see that

$\displaystyle\sum_{i}P(i\sim
j|d_{i,j}=x;\;i,j\in\mathcal{B})=\sum_{i}\frac{F_{\mathcal{B},1}^{\prime}(x)}{F_{\mathcal{B},2}^{\prime}(x)}\;.$

Evaluating this at the distances between the nodes as they have been
distributed in $\mathbb{R}^{d}/\mathbb{Z}^{d}$, we may rewrite this as

$\displaystyle\omega_{i}=\sum_{j}P(i\sim
j|d_{i,j};\;i,j\in\mathcal{B})=\sum_{j}\frac{F_{\mathcal{B},1}^{\prime}(d_{i,j})}{F_{\mathcal{B},2}^{\prime}(d_{i,j})}$

∎

We can use this to calculate $\vec{\rho}_{\mathcal{B}}$, the intensity vector
in the bounded region.

######  Theorem 3.13.

The expected degree of a vertex $i$ on the induced subgraph in the region
$\mathcal{B}$, given $\rho_{i}$ and its distance from all points, $d_{i,j}$,
is

$\rho_{\mathcal{B},i}=\rho_{i}\frac{\omega_{i}}{n\varepsilon}$

###### Proof.

If its intensity is known, as well as the position of the other nodes in the
graph. Therefore, we may write

$\displaystyle\mathbb{E}(deg_{i}^{\mathcal{B}}|\rho_{i},d_{i,j})$
$\displaystyle=$ $\displaystyle\sum_{j\in\mathcal{B}}P(i\sim
j|\rho_{i},d_{i,j})$ (11) $\displaystyle=$
$\displaystyle\sum_{j\in\mathcal{B}}\sum_{\ell}P(\rho_{\ell}|\rho_{i},d_{i,j})P(i\sim
j|\rho_{i},\rho_{j}=\rho_{\ell},d_{i,j})$ (12)

where the second equality comes from the marginal law of probabilities.

By construction, $P(\rho_{\ell}|\rho_{i},d_{i,j})=\frac{1}{n}$. Therefore,

$\displaystyle\mathbb{E}(deg_{i}^{\mathcal{B}}|\rho_{i},d_{i,j})=\frac{1}{n}\sum_{j\in\mathcal{B}}\sum_{\ell}P(i\sim
j|\rho_{i},\rho_{j}=\rho_{\ell},d_{i,j})$

By substituting the restricted functions $F_{\mathcal{B},1}(x)$ and
$F_{\mathcal{B},1}(x)$ into Theorem 3.5 and evaluating at the distances
between nodes, the right hand side becomes

$\displaystyle\mathbb{E}(deg_{i}^{\mathcal{B}}|\rho_{i},d_{i,j})=\frac{1}{n}\sum_{j\in\mathcal{B}}\frac{F_{\mathcal{B},1}^{\prime}(d_{ij})}{\hat{F}_{2}^{\prime}(d_{ij})}\frac{1}{\varepsilon}\cdot\min\\{\sum_{\ell}\frac{\rho_{i}\rho_{\ell}}{\sum_{k}\rho_{k}},1\\}$

Note that $\sum_{\ell}\frac{\rho_{i}\rho_{\ell}}{\sum_{k}\rho_{k}}=\rho_{i}$.
Therefore, by Lemma 3.12, we get

$\displaystyle\mathbb{E}(deg_{i}^{\mathcal{B}}|\rho_{i},d_{i,j})=\omega_{i}\frac{1}{n\varepsilon}\cdot\rho_{i}\;.$

Definition 3.10 gives the desired statement. ∎

###  3.3 Estimating the parameters given a network

In many cases of real world data, one does not actually have the parameters
needed to fit the model described above. Rather, they must be estimated from
the data observed in a few given networks.

In this section, given a fixed network $G$, we estimate the parameters given
in Definitions 3.1 and 3.7. In Section 4 we used these estimated values to
synthesize networks similar to $G$. We call $G$ the reference graph.

######  Definition 3.14 (Inputs from the data).

  1. 1.

The reference graph $G$ with a list of vertices and edges $G=(V,E)$ where $E$
may include self-loops. (Note that this input includes the degree of each
vertex, ${deg}_{i}^{G}$.)

  2. 2.

A set of coordinates $\mathcal{P}$ for the vertices in the reference graph
(from which all $d_{i,j}$ can be calculated).

######  Definition 3.15 (Estimated quantities).

  1. 1.

A function $\hat{F}_{1}(x)$ fit from the location and adjacency data in the
reference graph: $F_{\mathcal{B},1}(x)=P(i\sim
j\;,d_{i,j}<x|i,j\in\mathcal{B})$.

  2. 2.

A function $\hat{F}_{2}(x)$ fit from the empirical cumulative probability
function of the geometric distribution of the vertices of the reference graph:
$F_{2}(x)=P(d_{i,j}<x|i,j\in\mathcal{B})$.

  3. 3.

A vector of intensities for each vertex, $\hat{\rho}_{i}$ to be calculated
according to Theorem 3.13.

To calculate $\hat{F}_{1}(x)$, $\hat{F}_{2}(x)$, we simply fit the
corresponding equally-spaced quantiles based upon the inputs of $G$ and
$\mathcal{P}$. To calculate $\hat{\rho}$, we may write $\deg_{i}^{G}$ as a
point estimate for $\mathbb{E}(deg_{i}^{\mathcal{B}}|\rho_{i},d_{i,j})$. By
rearranging Theorem 3.13 may take the estimate the vector of intensities,
$\rho$, as

$\displaystyle\hat{\rho}_{i}=\frac{deg_{i}^{G}\cdot n\varepsilon}{w_{i}}.$
(13)

In other words, the estimator for $\rho_{i}$ is proportional to the ratio
between the degree of the references graph and the degree one would expect
based on the geometries of the nodes.

There are several remarks to note about these estimates.

First, $\rho_{i}$ (and $\hat{\rho}_{i}$) represent the expected degree in the
full unbounded model; hence, one should expect that $\hat{\rho}_{i}$ is larger
than the corresponding degree in the reference graph, $\deg^{G}_{i}$.

Second, note that the numerator of $\hat{F}_{1}$ is $\varepsilon$. From
Definition 3.1, one expects the limit

$\displaystyle\lim_{x\to\infty}F_{1}(x)=\lim_{x\to\infty}P(i\sim
j,d_{i,j}<x|i,j\in\mathcal{B})=\varepsilon$

to be the edge density of the graph. Similarly, the numerator of $\hat{F}_{2}$
is $1$. From Definition 3.1, one expects the limit

$\displaystyle\lim_{x\to\infty}F_{1}(x)=\lim_{x\to\infty}P(d_{i,j}<x|i,j\in\mathcal{B})=1$

to be 1, as the volume of $\mathcal{B}$ is finite.

Further, one may notice that we omit several of the more trivial parameters of
the model as $n$ and $d$ are subsumed within the data given by $\mathcal{P}$.
On the other hand, $\mu$, is not taken as an input. Instead, given a reference
graph $G$, we take $\mu$ to be the point-masses induced by the locations of
the vertices given by $\mathcal{P}$. That is, we fix the locations of the
vertices of the input data.

##  4 Simulated Drosophila Medulla connectome

In this section, we consider two graphs arising from the connectome of the
Drosophila Medulla and fit the parameters of the model described in section 3.

###  4.1 Drosophila Medulla Connectome

As an example graph to validate this model, we consider a connectome of the
visual system of the Drosophila medulla consisting of 1781 different neurons
with 33,641 synapses [26]. Note that this is not the full Drosophila medulla,
which has approximate 40,000 neurons [14]. Instead, this is a mapping of a
physically extracted sample of full medulla. Briefly, the authors of [23]
consider seven neuron columns important for motion detection. The neurons in
the central column are traced, that is, all of their connections to neurons in
neighboring six columns are recorded. Each synapse is given a pair of
locations, a pre-synaptic and and post-synaptic position, corresponding to the
two halves of the connection corresponding to the two neurons involved.

Biologically speaking, synapses take up small volumes, while the neurons
themselves can be large with many long, spindly dendrites that branch across a
large volume and intertwine with dendrites of other neurons. Each synapse
represents a directed edge between the neurons (from the pre-synaptic neuron
to the post synaptic neuron), where a pair of neurons may share multiple
directed connections between each other, in either direction. Occasionally
neurons have synapses with themselves.

In order to validate this model, we represent the neurons as the nodes of the
graph while the synapses are represented by edges. We further simplify the
connectome by ignoring the directions and multiplicities of the edges between
the nodes. That is, we do not keep track of how many synapses connect one
neuron to another, or in which direction the synapse allows for current flow.
We only track whether or not two neurons have a connection between them. We
do, however, keep track of when a neuron is connected to itself. That is, we
allow for self loops in the graph.

Finally, as we are studying a geometric model, we assign positions to the
neurons. Each synapse in the connectome is given a pre-synaptic and a post-
synaptic coordinate, as well as a pre-synaptic neuron and a post-synaptic
neuron. On a directed multi-graph, one may think of this as each directed edge
being decorated by a pair of coordinates (pre and post coordinates). To assign
a single coordinate to a neuron, we take the geometric centroid of the
synaptic coordinates associated to it. I.e. if we were to model this as a
directed multi-graph, we take the centroid of the pre-synaptic coordinates of
all the edges from a node, and all the post-synaptic coordinates of all the
edges to the node. In this manner, we define a point position to a
fundamentally three-dimensional, non-convex volume filling object. Note that
we have no reason to believe that the coordinate thus assigned to the neuron
lies in the physical neuron in the medulla, let alone reason to believe that
it corresponds to any meaningful part of the neuron, such as the location of
the soma, or the center of mass of the neuron. However, given the data at
hand, and the intrinsic difficulty of mapping the geometry of all the neurons
in a connectome, we use this as an approximation of the “location” of each
neuron, see Figure 1. Figure 2 is a graphic of the connectome described here.

Figure 1: Three views of the graph of the connectome of the Drosophila Medula,
with the vertex locations computed as the centroid of the synapse locations.
Note that the dimensions of the bounding box of the connectome is $7290\;\mu
m\times 8685\;\mu m\times 1212\;\mu m$. Note that $98\%$ of the data lies in a
bounding box of dimensions $4300\;\mu m\times 4005\;\mu m\times 1130\;\mu m$.
Figure 2: The connectome of the of the Drosophila Medula connectome

In summary, we model the connectome found in [26] as a simple graph with self
loops with $1781$ and $9016$ edges, $105$ of which are self loops. However,
there is a further wrinkle that must be explored. As stated in [23], the
authors trace the $379$ neurons in the central neural column extracted from
the Drosophila medulla. That is, all of neighborhoods of these central $379$
vertices are traced, but connections between other neurons are not necessarily
traced. In other words, this connectome may be better modeled if seen as a
union of stars of the fully traced neurons.

In the dataset present in [26], there are $358$ “named” neurons (i.e. these
are the neurons with alpha-numeric labels) that correspond roughly to the
central $379$ nodes.333This correspondence has been verified by the authors of
[26].. In aggregate, the location (as determined above) of the named neurons
are slightly closer to the centroid than all vertices. Furthermore, the mean
degree of the named neurons is 30.05 versus 10.01 for all neurons.

Therefore, in this paper, we consider two graphs: the _full connectome_ and
the _named subgraph_. The full connectome is the graphical representation of
the entire data set (as described above). The named subgraph is the induced
subgraph of the full connectome defined by the $358$ named vertices. The named
subgraph has 3254 edges, 72 of which are self loops. For our purposes, each of
these networks has a select few nodes with an overwhelmingly high degree,
violating the Chung-Lu condition. For simplicity, we study the networks with
these outlier nodes removed which we denote $G_{\textrm{named}}^{(-2)}$ and
$G_{\textrm{full}}^{(-1)}$.

###  4.2 Justification for a geometric model

The generic geometric Chung-Lu model presented in this paper is particularly
well suited for creating synthetic models for a connectome. In particular,
neurons are more likely to connect to those closer than farther from it,
indicating that there is a geometric component to the formation of this
network. On a mesoscale (studying the connections between brain regions
roughly the size of the entire graph studied in this paper) there is
successful work modeling the strength of connections between brain regions
using a continuous spatial propagation model [21]. Furthermore, different
types of neurons have different numbers of synapses, but this number is by no
means fixed. In fact, there is evidence that the these connections change over
an organisms lifespan [3]. Furthermore, different neurons of the same type
having different valencies, for reasons that are not fully biologically
understood. This can be most easily seen in the small Caenorhabditis elegans
connectome [24]. Therefore, it is natural to work with a Chung-Lu type model
where one associates an expected degree or intensity to each vertex.

In our specific context, we make a few more observations. First, for the full
connectome minus one vertex, $G_{\textrm{full}}^{(-1)}$, the mean distance
between any two nodes is 2000.40 $\mu$m; in contrast, the mean distance for
any edge 963.25$\mu$m. This substantial discrepancy indicates a strong
relationship between the physical distance between two neurons and the
likelihood of a connection between them. Also, as illustrated in Figure 4, the
distributions of edge lengths appear to be similar but shifted based upon the
valency of the node. This suggests that both the distance between two nodes as
well as their degrees (i.e., intensities) should play a role in the connection
function. Figure 3 shows the neighborhood of the first, hundredth and two
hundredth most valent vertex in the connectome of the Drosophila medulla
graph.

(a) Neighorbohood of the vertex with the highest valence.

(b) Neighorbohood of the vertex with the $100^{th}$ highest valence.

(c) Neighorbohood of the vertex with the $200^{th}$ highest valence.

Figure 3: Note that while the diameter of the largest connected component in
Figure 2 has diameter 6, the neighborhoods of the vertices with the highest,
$100^{th}$ highest and $200^{th}$ highest valencies are 927, 37 and 23
respectively, suggesting a correlation between the valency of a vertex and is
average distance from neighbors.

(a) Empirical cumulative distribution of edge distances for the highest
valency node in the full connectome minus a vertex.

(b) Empirical cumulative distribution of edge distances for the 50-$th$
highest valency node in the full connectome minus a vertex.

(c) Empirical cumulative distribution of edge distances for the 100-$th$
highest valency in the full connectome minus a vertex.

Figure 4: The empirical cumulative distribution of edge distances have similar
‘S’ shaped in spite of its valence.

Finally, it is worth noting that the last assumption of our model, that the
distances and intensities of our vertices are independent is not a biological
observation. It is an assumption that is needed for mathematical simplicity.
In fact, there is biological evidence that this assumption may not hold, For
the Wistar rat, the number of synaptic connections vary greatly between
different layers of the brain [1]. It is also well known that the neurons in a
brain have non-uniform physical structure. For example, this is explored in
detail for the Drosophila Medulla in [10].

###  4.3 Estimated Parameters for the Drosophila Medulla connectome

In this section, we fit the parameters presented in section 3.3 for the two
graphs $G_{\textrm{named}}^{(-2)}$ and $G_{\textrm{full}}^{(-1)}$ defined in
section 4.1.

First, we note that after fitting $\hat{F}_{1}$ and $\hat{F}_{2}$ on the named
subgraph $G_{\textrm{named}}$ and the full connectome, $G_{\textrm{full}}$ the
Chung-Lu condition $\hat{\rho}_{i}^{2}<\sum_{k}\hat{\rho}_{k}$ does not hold
for the two highest values of $\hat{\rho}_{i}$ for the named subgraph, and the
highest value for the full graph. One potential solution to this is to follow
[4], and set any probability greater than 1 generated in this manner to 1.
Another solution (the one adopted in this paper) is to omit the vertex with
highest degree from our calculations. For the named subgraph, we disregard the
2 most valent vertices and call the graph $G_{\textrm{named}}^{(-2)}$. For
full graph, we omit the single most valent vertex and call the graph
$G_{\textrm{full}}^{(-1)}$. After these omissions, both
$G_{\textrm{named}}^{(-2)}$ and $G_{\textrm{full}}^{(-1)}$ satisfy the Chung-
Lu condition.

Since our goal is to generate a new simplified model that empirically
reproduces many features of the underlying data based on a single instance, it
is not appropriate to formally estimate parameters with confidence intervals
nor test a null hypothesis. Instead, the parameters learned from the data are
a necessary step towards creating the connection function in the model
presented in section 3.3.

Since there is no posited true simple model, there is no null hypothesis to
reject and no substantive interpretation of any “standard errors” that could
be calculated, and thus they are not reported. Instead, to validate our
parameters, we present the mean squared error of the observed data from the
fit functions.

In particular, we wish to estimate the probability $P(i\sim
j,d_{i,j}|\mathcal{B})$ logistically with the function

$\displaystyle\hat{F}_{1}=\frac{\varepsilon}{1+\exp(\alpha_{1}+\beta_{1}x)}$

and the probability $P(d_{i,j}|\mathcal{B})$ with the logistic function

$\displaystyle\hat{F}_{2}=\frac{1}{1+\exp(\alpha_{2}+\beta_{2}x)}\;.$

Note that we write the logistic functions above as
$\frac{L}{1+\exp(\alpha+\beta x)}$, and not $\frac{L}{1+\exp(-k(x-x_{0}))}$
which is not the standard notation. In the notation of this paper, $-\beta$
give the growth rate of the function. In particular, $\beta$ needs to be a
negative number, while the logistic function is centered at
$x=-\frac{\alpha}{\beta}$.

In Table 1, we report the estimated values for the constants $\alpha_{1}$ and
$\beta_{1}$ that define $\hat{F}_{1}(x)$ and the constants $\alpha_{2}$ and
$\beta_{2}$ that define $\hat{F}_{2}(x)$.

| $G_{\textrm{named}}^{(-2)}$ | $G_{\textrm{full}}^{(-1)}$
---|---|---
$\alpha_{1}$ | 2.94589 | 2.76129
$\beta_{1}$ | -0.00293112 | -0.00300996
$\alpha_{2}$ | 3.16284 | 3.31531
$\beta_{2}$ | -0.0018741 | -0.00170436

Table 1: Estimated parameters defining $\hat{F}_{1}$ and $\hat{F}_{2}$ using
the reference graphs $G_{\textrm{named}}^{(-2)}$ and
$G_{\textrm{full}}^{(-1)}$.

Next we consider the estimated and observed values at $x=0$ in Table 2. Recall
from Remark 3.4 that that $\lim_{x\to 0^{+}}F_{1}(0)=\frac{|\textrm{self
loops}|}{n^{2}}$. This leads to the observed values of 0.00056022 for
$G_{\textrm{named}}^{(-2)}$ and 0.0000331398 for $G_{\textrm{full}}^{(-1)}$.
However, the estimated values from the model are an order of magnitude higher:
0.0020432 for $G_{\textrm{named}}^{(-2)}$ and 0.000269171 for
$G_{\textrm{full}}^{(-1)}$. Furthermore, the observed value for
$F_{2}(0)=1/n$, or the probability that the node $i$ is the same as the node
$j$. This leads to the observed values of 0.00280899 for
$G_{\textrm{named}}^{(-2)}$ and 0.000561798 for $G_{\textrm{full}}^{(-1)}$.
Again, the estimated values are orders of magnitude larger than the observed
values.

| $G_{\textrm{named}}^{(-2)}$ | $G_{\textrm{full}}^{(-1)}$
---|---|---
| estimated | observed | estimated | observed
$\lim_{x\rightarrow 0}\hat{F}_{1}(x)$ | 0.00224016 | 0.00056022 | 0.000303395 | 0.0000331398
$\lim_{x\rightarrow 0}\hat{F}_{2}(x)$ | 0.0405882 | 0.00280899 | 0.0350498 | 0.000561798
Table 2: Model-estimated and observed limiting values of the functions
$\hat{F}_{1}$ and $\hat{F}_{2}$ for loops and edges of small length. The
impact of these discrepancies are shown in Figure 6.

In spite of these inconsistencies, we show in Section 5.1 that the synthetic
graphs generated by this model match the graphical characteristics of the
reference graphs $G_{\textrm{named}}^{(-2)}$ and $G_{\textrm{full}}^{(-1)}$.
Furthermore, we show that the estimated values of $\hat{F}_{1}(x)$ and
$\hat{F}_{2}(x)$ overall do not differ greatly from the observed value.
Namely, in Table 3 give the Mean Squared Error (MSE) of the model as well as
the mean percent error for $\hat{F}_{1}(x)$ and $\hat{F}_{2}(x)$ over the
entire range of distances, and away from $x=0$ where the model fit is the
worst.

For $\hat{F}_{1}(x)$ the mean percent error is less than $5\%$ overall, and
less than $1\%$ for values of $x>500$. The mean percent error for $F_{2}(x)$
is greater, as indicated by large difference in the estimated vs observed
values near $x=0$ (Table 2). However, we see that away from this region,
namely for $x>500$, the mean percent error falls to less than $5\%$. There is
a similar pattern with the MSE: the overall MSE is greater for F2 than it is
for F1, but for both fits, this is dominated by the behavior for low $x$. For
$x>0$ the MSE drops by at least an order of magnitude.

|  | $G_{\textrm{named}}^{(-2)}$ | $G_{\textrm{full}}^{(-1)}$
---|---|---|---
|  | all $x$ | $x>500$ | all $x$ | $x>500$
$\hat{F}_{1}(x)$ | MSE | 0.0138996 | 0.000175405 | 0.0201616 | 0.000250721
Mean Percent error | 4.24187 % | 0.822747 % | 4.68995 % | 0.901061 %
$\hat{F}_{2}(x)$ | MSE | 0.0538964 | 0.00856598 | 0.0694437 | 0.0127325
Mean Percent error | 12.215 % | 5.32331 % | 12.655 % | 4.9259 %
Table 3: Mean square errors and Mean percentage errors of $F_{1}$ and $F_{2}$
for the graphs $G_{\textrm{named}}^{(-2)}$ and $G_{\textrm{full}}^{(-1)}$.
These fits are good for larger edges.

This behavior is clear from Figure 5, which shows the logistic fits of the
cumulative functions $\hat{F}_{1}$ and $\hat{F}_{2}$ for the graph
$G_{\textrm{named}}^{(-2)}$ and $G_{\textrm{full}}^{(-1)}$.

(a) Fits for $F_{1}(x)$ for $G_{\textrm{named}}^{(-2)}$ (left) and
$G_{\textrm{full}}^{(-1)}$ (right).

(b) Fits for $F_{2}(x)$ for $G_{\textrm{named}}^{(-2)}$ (left) and
$G_{\textrm{full}}^{(-1)}$ (right).

Figure 5: The empirical probabilities (blue) for $F_{1}$ (top) and $F_{2}$
(bottom) against their modeled probabilities (orange) ($\hat{F}_{1}$ and
$\hat{F}_{2}$ respectively.). The fits are good, particularly for longer
edges.

Figure 6 shows the empirical and estimated ratio
$\hat{F}_{1}^{\prime}(x)/\hat{F}^{\prime}_{2}(x)$, which is central to the
connection function defined in Theorem 3.5. As with the estimated versus
observed values for $\hat{F}_{1}(x)$ and $\hat{F}_{2}(x)$, we see that the
estimated value is a much better fit for large values of $x$ for both
$G_{\textrm{full}}^{(-1)}$ and $G_{\textrm{named}}^{(-2)}$. In other words,
the synthetic models generated by these graphs will have fewer edges between
nodes that are physically close together than observed in the reference
graphs.

Figure 6: Plot of $\hat{F}_{1}^{\prime}(x)/\hat{F}^{\prime}_{2}(x)$ (blue)
against the empirical probability of an edge at a fixed distance ($P(i\sim
j~{}|d_{i,j}=x)$) (orange). The estimated value is a much better fit for large
values of $x$ for both $G_{\textrm{full}}^{(-1)}$ and
$G_{\textrm{named}}^{(-2)}$

Finally, Figure 7 addresses the parameters $\hat{\rho_{i}}$ as defined in
display (13). The figure shows that there is a linear relationship between the
observed degree and the estimated intensity near the centroid of the neurons,
i.e. where the boundary effects laid out in Section 3.2 have the smallest
influence on the parameters. Furthermore, this proportionality constant near
the centroid is less than one. As one moves further from the centroid, i.e.
when the boundary of the graph plays a greater role in determining the
observed degree of the graph, the estimated intensity increases for any given
degree. This is exactly as is expected: the closer a given vertex is to the
boundary (i.e. the further it is from the centroid) the more likely it is to
have connections to other vertices that are not within the boundary, and thus
are cut off.

Figure 7: Scatter plots of the degree of each node in the reference graph
against its associated intensity, $\rho_{i}$. $G_{\textrm{named}}^{(-2)}$ is
on the left and $G_{\textrm{full}}^{(-1)}$ is on the right. Colors indicate
the distance of that vertex to the centroid of all node. There is a linear
relationship between the observed degree and the estimated intensity near the
centroid of the neurons.

##  5 Simulation results

In this section, we provide validation for the model that we have developed in
Section 3.3 and fit in Section 4.3. We do this by generating 200 simulated
graphs from the connection function (Theorem 3.5) and plot the histograms of
several key network properties on this sample of 200 graphs. We then show that
the reference graphs $G_{\textrm{named}}^{(-2)}$ and
$G_{\textrm{full}}^{(-1)}$ fall well within the range of values from these
simulated graphs.

For our simulations we use the estimated parameters in Section 4.3 and
construct 200 graphs for $G_{\textrm{named}}^{(-2)}$ and 100 graphs for
$G_{\textrm{full}}^{(-1)}$. We emphasize that, the positions of the nodes are
fixed; however, the assignment of the $\rho_{i}$ are randomly permuted among
them.

###  5.1 Edges, vertices and connected components

If the connection function produces a good model for the graph, then the
simulated graphs would have similar basic graph characteristics, such as the
number of edges and number of connected components, etc. In Table 4, we report
the mean and standard deviation of various graph characteristics for models
built from $G_{\textrm{named}}^{(-2)}$ and $G_{\textrm{full}}^{(-1)}$ with the
corresponding values from the original graph. For many of the characteristics
including the maximum degree, the number of triangles and the number of closed
4-walks, the values for the reference graph are within two standard deviations
of the mean of the simulated graphs. Figure 8 provides a comparison of the
simulated distribution of these characteristics. In all of these cases, the
reference graph has characteristic counts that fall within the range of the
simulations. We emphasize that the random assignment of the $\rho_{i}$ among
the vertices appears induce a large variance for each of these
characteristics; for instance, if nodes closer to the boundary are assigned
larger $\rho_{i}$, then the resulting network will have fewer edges (and hence
triangles and 4-walks, etc.). This variance is realized in the simulations for
both $G_{\textrm{named}}^{(-2)}$ and $G_{\textrm{full}}^{(-1)}$, as many of
the parameters for the reference graph fall in the lower of range of the
simulations. None the less, the reference graphs are is consistent with the
positive correlation among these parameters as illustrated in Figure 8(g) and
8(h).

While the simulations produce appropriate numbers of triangles and closed
4-walks as well as a comparable maximum degree, there are significant
discrepancies between the model and the reference graph are in the number of
self-loops and the number of connected components. Both of these short-comings
can be attributed to the effect of the model on node(s) that are a short
distance apart (i.e. the mismatch of the expected and observed values for
$F_{1}(x)$ and $F_{2}(x)$ near $x=0$ discussed in Section 4.3.)

In particular, from Table 2 and Figure 5, at $x=0$, $F_{2}$ is significantly
more overestimated than $F_{1}$. As a result, the ratio
$\hat{F}_{1}/\hat{F}_{2}$ underestimates the conditional probability $P(i\sim
j|d_{i,j})$. However, the connection function given in Theorem 3.5 depends on
the ratio $\hat{F}_{1}/\hat{F}_{2}$. From Figure 13, we see that the estimated
ratios in that define the connection function are also lower than the
empirically observed ratios. This leads to an underestimation of the number of
self loops and connections at short length in the model than is observed. In
particular, if one only considered Remark 3.4, were $F_{1}(0)$ gives the
number of self loops, one would expect the simulated graphs to have, on
average an larger number of self-loops than the reference graphs, which is in
contradiction with the data presented in Table 4. However, the observed
difference in the estimated ratio $\hat{F}_{1}/\hat{F}_{2}$ versus the
observed value of $P(i\sim j|d_{i,j}=0)$ counteracts this overestimation.

In addition, the biophyisical nature of the reference graph makes it nearly a
priori connected with few isolated neurons. For instance, neurons that fail to
make appropriate connection often die (see for example [7]); such neurons
would not be realized in our reference graph. As a result, the discrepancy
concerning the number of self-loops and the number connected components
appears can be attributed to the nuances of the specific application of the
model (i.e., connectome networks) as opposed to the overall model itself.

| $G_{\textrm{named}}^{(-2)}$ | $G_{\textrm{full}}^{(-1)}$
---|---|---
| Ref. Graph | mean | St. Dev. | Ref. Graph | Mean | St.Dev
Self-loops | 71 | 47.265 | 4.60552 | 105 | 69.76 | 7.16913
Connected Components | 8 | 20.415 | 3.06098 | 26 | 234.81 | 10.4134
Non-Iso. Conn. Comp. | 1 | 1.17 | 0.40238 | 7 | 3.46 | 1.77195
Max Valency | 71 | 100.725 | 13.0237 | 150 | 263.72 | 53.0845
Number of Edges | 2851 | 2979.1 | 94.9285 | 8089 | 8194.41 | 221.677
Number of Triangles | 4400 | 5546.18 | 795.223 | 10224 | 9621.73 | 1402.54
Number of Closed 4-Walks | 981091 | 1.37876 $\times 10^{6}$ | 236890 | 2.95063 $\times 10^{6}$ | 3.00237 $\times 10^{6}$ | 549348.
Table 4: Graph statistics for 200 graph generated by the model built from the
Named graph, as well the values for the Named graph. Many of values for the
reference graph are within two standard deviations of the means for the
simulated graphs.

(a) Histogram of maximum degree over 200 model-generated graphs for
$G_{\textrm{named}}^{(-2)}$.

(b) Histogram of maximum degree over 100 model generated
$G_{\textrm{full}}^{(-1)}$.

(c) Histogram of the number of closed 4-walks over 200 model-generated graphs
for $G_{\textrm{named}}^{(-2)}$.

(d) Histogram of the number of closed 4-walks over 100 model-generated graphs
for $G_{\textrm{full}}^{(-1)}$.

(e) Histogram of three cycles over 200 model-generated graphs for
$G_{\textrm{named}}^{(-2)}$.

(f) Histogram of three cycles over 100 model-generated graphs for
$G_{\textrm{full}}^{(-1)}$.

(g) Point plot of the maximum degree versus the number of triangles over 200
model-generated graphs for $G_{\textrm{named}}^{(-2)}$.

(h) Point plot of the maximum degree versus the number of triangles over 100
model-generated graphs for $G_{\textrm{full}}^{(-1)}$.

.

Figure 8: Histograms for various graph parameters of randomly-generated graphs
from the model built from based on $G_{\textrm{named}}^{(-2)}$ (left) and
$G_{\textrm{full}}^{(-1)}$ (right) with the value for the reference graph
indicated in red. In all these cases, the reference graph observation lies
within the support of the simulation observations.

###  5.2 Eigenvalues

One method for comparing networks is to compare the eigenvalues of various
graph-theoretic matrices of the two networks, such as the adjacency matrix or
Laplacian matrix. The eigenvalues of these matrices strongly correspond to the
local and global properties of the network such as degree, clustering
coefficient, motif counts and so on [9]. It possible for two non-isomorphic
networks to have the exact same spectrum (i.e., set of eigenvalues) for these
graph-theoretic matrices [13]; however, such cospectral networks have similar
properties. Hence, a network model used to generate networks similar to a
desired target network, such as in our case, should produce networks with a
similar spectrum.

Figure 9 shows the histograms for the eigenvalues of the adjacency matrices
for both $G_{\textrm{named}}^{(-2)}$ and $G_{\textrm{full}}^{(-1)}$, overlaid
with eigenvalues for a corresponding simulated graph. In general, the randomly
generated network has a similar distribution of eigenvalues to the appropriate
reference graph. One noticeable discrepancy is that the randomly generated
network has distinctly more eigenvalues near 0. In fact, the multiplicity of
the exact eigenvalue of 0 for the randomly generated network is significantly
higher that the fixed network. We attribute this phenomena to the fact that
the randomly generated network is likely to create isolated nodes with no
incident edges whereas the original network is a priori taken to be connected.
Each isolated node corresponds exactly to an additional instance of 0 as an
eigenvalue of the adjacency matrix. Hence, since our model produces more
isolated nodes than the reference graph, one should expect substantially
higher multiplicity of 0 as an eigenvalue in the simulations.

(a) A histogram of the eigenvalue densities for the adjacency matrix for the
reference graph $G_{\textrm{named}}^{(-2)}$ (yellow) compared that of a
randomly-generated model graph (blue).

(b) A histogram of the eigenvalue densities for the adjacency matrix for the
reference graph $G_{\textrm{full}}^{(-1)}$ (yellow) compared that of a
randomly-generated model graph (blue).

Figure 9: The distribution of eigenvalues for the simulated graphs is close to
the distribution of eigenvalues for the reference graphs.

###  5.3 Centrality Measures

Next, we compare the centrality measure of the graphs generated by the two
models and to those of the original graphs. Centrality measures are various
heuristic measures to determine how important a vertex is based upon the
structure of the network. There are many different centrality measures; here,
we focus on three traditional measures: betweenness, closeness, and
eigenvector centrality. The betweenness of a vertex $i$ is
$\sum_{u,v}f_{u,v}(i)$ where $f_{u,v}(i)$ is the proportion of shortest paths
between $u$ and $v$ that contain $i$. The closeness of a vertex $i$ is
$N/\sum_{j\neq i}d_{i,j}$ where $1/\infty$ is taken to be 0. And the
eigenvector centrality of vertex $i$ is given by the principle eigenvector of
the adjacency matrix of the network corresponding to the maximum eigenvalue;
Perron-Frobenius theory guarantees such an eigenvalue exists and its
eigenvector is nonnegative.

We focus on these centrality measures, as they generally capture distinct
aspects of the graph including local motifs (eigenvector centrality), distance
(closeness), and flow (betweenness). As such, these measures are typically
hard to replicate jointly. Figure 10 shows the plots for eigenvector
centrality, betweenness, and closeness. For $G_{\textrm{named}}^{(-2)}$, the
simulated network emulate the centralities of the reference graph very well.
On the other hand, for $G_{\textrm{full}}^{(-1)}$, the centralities of the
simulated network do not match as well. We attribute this to the fact that the
network data originates from focusing on select nodes (neurons), and that
$G_{\textrm{full}}^{(-1)}$ is mostly-star like. On the other hand, the model
we present here does not indicate any focused nodes. This discrepancy effects
the local structure of the simulated graph which in turn, affects the
centralities measures of the vertices. In contrast, the named graph, is
effectively an induced subgraph on the selected nodes.

(a) Eigenvector centrality of the Named graph compared to a randomly generated
graph generated by the corresponding model.

(b) Eigenvector centrality of the Full connectome compared to a randomly
generated graph generated by the corresponding model.

(c) Betweenness measure of vertices in the Named graph compared to a randomly
generated graph generated by the corresponding model.

(d) Betweenness measure of vertices in the Full connectome compared to a
randomly generated graph generated by the corresponding model.

(e) Closeness of vertices in the Named graph compared to a randomly generated
graph generated by the corresponding model.

(f) Closeness of vertices in the Full connectome compared to a randomly
generated graph generated by the corresponding model.

Figure 10: Centrality comparisons of a representative sample from the model
built from the Named graph (left) or the Full graph (right) and the reference
graph. In all these cases, the centrality measures, sorted by rank, mimic that
of the reference graph.

##  6 Comparison to other models

The model built in this paper is built to accommodate the biological structure
of neurons, which are more likely to form connections with near by neurons
than ones further away [15, 22, 6]. While much of this work has been conducted
on a meso scale, in this paper we have shown a correlation between the
distance between two nodes and the presence of an edge between them (Figures 6
and 13). Even so, it is plausible that a simpler network model with fewer
parameters generates the network which is later embedded spatially in a more
optimized way (e.g., like a spring network).

The logistic fit for $\hat{F}_{1}$ and $\hat{F}_{2}$ in this generic geometric
Chung-Lu model is much more complicated than existing models such as the
classical Chung-Lu model and the polynomial decay models such as GIRG model
(as mentioned in Section 2). However when we compare our logistically fit
generic geometric model to both the Chung-Lu model and the polynomial decay
models, we see that the our model does indeed have a better fit.

The classical Chung-Lu model has no consideration for the physical geometry of
the vertices. Hence, the nature of clustering should be distinct in the
classical model over geometric ones. The Chung-Lu model generates a network
whose nodes are substantially closer together than the logistic-geometric
model, resulting in a substantially higher closeness centrality. Figure 11
illustrates that the closeness centrality is nearly-uniformly higher in the
classical Chung-Lu model. This stands in direct contrast to the closeness of
the logistic-geometric model (see Figure 10) which far better mimics the
reference graphs.

Figure 12 compares the number of triangles each vertex is involved in for a
single sample graph from each model compared to the reference graph, and to a
random graph generated by a simple Chung-Lu model.

(a) The closeness centrality of a random graph based on the classical Chung-Lu
network using $G_{\textrm{named}}^{(-2)}$ as reference.

(b) The closeness centrality of a random graph based on the classical Chung-Lu
network using $G_{\textrm{full}}^{(-1)}$ as reference.

Figure 11: Closeness centrality of the classical Chung-Lu model for the
reference graphs $G_{\textrm{named}}^{(-2)}$ and $G_{\textrm{full}}^{(-1)}$.
The closeness of the Chung-Lu model is significantly higher than the reference
graph compared to the generalized geometric model as seen in 10.

(a) Number of triangles each vertex is involved in for the model generated by
the Named graph, a Chung-Lu model of the Named graph, and the Named graph.

(b) Number of triangles each vertex is involved in for the model generated by
the Full connectome, a Chung-Lu model of the Full connectome, and the Full
connectome.

Figure 12: The clustering coefficients (or “triangle degree”), sorted by rank,
of the original connectome, minus two vertices with high degree (blue), a
simulated graph based on the generalized geometric (orange) and a simulated
one based on the original Chung-Lu Model (green). This is done for the Named
graph (a) and the Full connectome (b). Observe the triangle degrees for the
generalized geometric model are closer to that for the reference graph than
the Chung-Lu model.

Previous geometric generalizations of the Chung-Lu model (i.e., GIRG. etc.)
have utilized connection functions with polynomial decay. In Figure 13
illustrates that the best fits for a connection function with polynomial decay
$(e.g.,kx^{-\beta})$ are unable to capture the appropriate decay rate for
longer distance functions.

(a) Plot of the empirical probability $y=P(i\sim j|d_{i,j}=x)$ (blue), the
modelled function $F_{1}^{\prime}(x)/F_{2}^{\prime}(x)$ (orange) and the best
inverse power fit $y=kx^{-\beta}=18.0418x^{-0.8088}$ (green) for
$G_{\textrm{named}}^{(-2)}$.

(b) Plot of the empirical probability $y=P(i\sim j|d_{i,j}=x)$ (blue), the
modelled function $F_{1}^{\prime}(x)/F_{2}^{\prime}(x)$ (orange) and the best
inverse power fit $y=kx^{-\beta}=1.40837x^{-0.732098}$ (green) for
$G_{\textrm{full}}^{(-1)}$.

Figure 13: The empirical probabilities of edges given distance (blue) versus
modeled probabilities using both the logistic model (orange) and an inverse
power model (green). Observe that the inverse power model is unable to
effectively model the edge-distance probabilities.

Finally, we end with a note on the nonmonotonicity of the connection function
developed in this paper (see Figures 6 and Figure 13). A nonmonotonic
connection function may appear to be unintuitive. However, Figure 6 shows that
both graphs $G_{\textrm{full}}^{(-1)}$ and $G_{\textrm{named}}^{(-2)}$ have a
non-monotonic $P(i\sim|d_{i,j},\;\hat{\rho}_{i},\;\hat{\rho}_{j})$.
Specifically, Tables 5 and 6 shows that for both $G_{\textrm{named}}^{(-2)}$
and $G_{\textrm{full}}^{(-1)}$ there is an increase in the probabilities for
connections in the short range while this probability decreases as the
distances get large. For another example, probabilities in endemic networks
have been observed to be non-monotonic and logistic-like [16].

Mathematically, the ratio of the derivative of two logistic functions (e.g.,
$\frac{F_{1}^{\prime}(x)}{F_{2}^{\prime}(x)}$) has the form

$\left.K~{}\frac{exp(\alpha_{1}+\beta_{1}x)}{(1+exp(\alpha_{1}+\beta_{1}x))^{2}}\right/\frac{exp(\alpha_{2}+\beta_{2}x)}{(1+exp(\alpha_{2}+\beta_{2}x))^{2}}$

which will be nonmonotonic and have a local maximum whenever
$\beta_{1}<\beta_{2}<0$ as is the case here. Upon considering this, we should
expect the connection function to be nonmonotonic for small distances.

$P(i\;\text{self-loop})$ | .199
---|---
$P(i\sim j\;|\>0<d_{ij}<200)$ | .270
Table 5: Empirical probabilities for connections of short length in $G_{\textrm{named}}^{(-2)}$. Observe the non-monotonic nature of the empirical probabilities. $P(i\;\text{self-loop})$ | .0.0727
---|---
$P(i\sim j\;|\>0<d_{ij}<100)$ | 0.107
$P(i\sim j\;|\>100<d_{ij}<200)$ | 0.0435
Table 6: Empirical probabilities for connections of short length in
$G_{\textrm{full}}^{(-1)}$. Observe the non-monotonic nature of the empirical
probabilities.

##  7 Conclusion

In this paper, we have developed a hybrid model for spatial networks with
heterogeneous nodes, which we call the generic geometric Chung-Lu model. The
connection function for this model is given in Theorem 3.5. This model is
designed specifically to model real world situation, where key assumptions
from the graph theory literature, namely, uniform distribution of the nodes on
a torus, is relaxed.

After relaxing these constraints, we use this model to generate networks that
are similar to the connectome Drosophila Medulla. Connectomes are a natural
candidate for this type of geometric model with heterogeneous nodes because of
the biological evidence both that different types of neurons have different
number of connections to other neurons, and that, due to energy conservation,
neurons are more likely to connect to other neurons that are physically
closer. We find that the model developed in this paper does a particularly
good job of modeling the features of the graph that control the information
flow across the network (i.e. eigenvalues, and centrality measures).

The approach taken in this paper is particularly well suited for generating
graphs with similar information flow properties to a given handful of
reference graphs. By relaxing the distributional and periodic constraints on
the model, we allow for a paradigm where one can learn the distribution of
nodes and the connections function from the reference data. In this manner,
the work in this paper can be adapted to a large number of scenarios where
there is reason to believe that the network has heterogeneous nodes and has a
spatial component to its connection function. For instance, ground based and
air transport networks are two other classes of data sets that include
distances between the nodes. Both data sets have key features that
differentiated it from the assumptions in this paper. (For instance, the there
are non-planarity costs in all rail networks, and air networks are both more
significantly and much more explicitly non-uniform than the connectome studied
here.) However, we believe that, given the assumptions made in this paper, the
model developed here can synthesize these real world networks with greater
accuracy than currently existing models.

## References

  * [1] L. Anton-Sanchez, C. Bielza, Merchán-Pérez A., J.R. Rodríguez, J. DeFelipe, and P. Larrañaga.  Three-dimensional distribution of cortical synapses: a replicated point pattern-based analysis.  Frontiers in Neuroanatomy, 85, 2014.
  * [2] Marc Barthélemy.  Spatial networks.  Physics Reports, 499(1-3):1–101, 2011.
  * [3] Sophie H. Bennett, Alastair J. Kirby, and Gerald T. Finnerty.  Rewiring the connectome: Evidence and effects.  Neuroscience and Biobehavioral Reviews, 88:51–62, 2018.
  * [4] Karl Bringmann, Ralph Keusch, and Johannes Lengler.  Geometric inhomogeneous random graphs.  Theoretical Computer Science, 760:35–54, 2019.
  * [5] Anna D Broido and Aaron Clauset.  Scale-free networks are rare.  Nature communications, 10(1):1–10, 2019.
  * [6] Ed Bullmore and Olaf Sporns.  Complex brain networks: graph theoretical analysis of structural and functional systems.  Nature reviews neuroscience, 10(3):186–198, 2009.
  * [7] Bruce M Carlson.  Human Embryology and Developmental biology E-book.  Elsevier Health Sciences, 2018.
  * [8] Fan Chung and Linyuan Lu.  Complex graphs and networks (CMBS regional conference series in mathematics).  2006\.
  * [9] Fan RK Chung.  Spectral graph theory.  Number 92. American Mathematical Soc., 1997.
  * [10] Carl P Dettmann and Orestis Georgiou.  Random geometric graphs with general connection functions.  Physical Review E, 93(3):032313, 2016.
  * [11] Ernesto Estrada.  Quantifying network heterogeneity.  Physical Review E, 82(6):066102, 2010.
  * [12] Nikolaos Fountoulakis.  On a geometrization of the chung–lu model for complex networks.  Journal of Complex Networks, 3(3):361–387, 2015.
  * [13] Chris D Godsil and Brendan D McKay.  Constructing cospectral graphs.  Aequationes Mathematicae, 25(1):257–268, 1982.
  * [14] Eri Hasegawa, Yusuke Kitada, Masako Kaido, Rie Takayama, Takeshi Awasaki, Tetsuya Tabata, and Makoto Sato.  Concentric zones, cell migration and neuronal circuits in the drosophila visual center.  Development, 138:983–993, 2011.
  * [15] Marcus Kaiser.  A tutorial in connectome analysis: topological and spatial features of brain networks.  Neuroimage, 57(3):892–907, 2011.
  * [16] John C Lang, Hans De Sterck, Jamieson L Kaiser, and Joel C Miller.  Analytic models for sir disease spread on random spatial networks.  Journal of Complex Networks, 6(6):948–970, 2018.
  * [17] Zhen-Qi Liu, Ying-Qiu Zheng, and Bratislav Misic.  Network topology of the marmoset connectome.  Network Neuroscience, 4(4):1181–1196, 2020.
  * [18] Mathew Penrose et al.  Random geometric graphs, volume 5.  Oxford university press, 2003.
  * [19] Mathew D Penrose et al.  Connectivity of soft random geometric graphs.  The Annals of Applied Probability, 26(2):986–1028, 2016.
  * [20] James A Roberts, Alistair Perry, Anton R Lord, Gloria Roberts, Philip B Mitchell, Robert E Smith, Fernando Calamante, and Michael Breakspear.  The contribution of geometry to the human connectome.  Neuroimage, 124:379–393, 2016.
  * [21] P. A. Robimson.  Physical brain connectomics.  Physical Review E, 99, 2019.
  * [22] H Francis Song, Henry Kennedy, and Xiao-Jing Wang.  Spatial embedding of structural similarity in the cerebral cortex.  Proceedings of the National Academy of Sciences, 111(46):16580–16585, 2014.
  * [23] Shin-ya Takemura, Arjun Bharioke, Zhiyuan Lu, Aljoscha Nern, Shiv Vitaladevuni, Patricia K Rivlin, William T Katz, Donald J Olbris, Stephen M Plaza, Philip Winston, et al.  A visual motion detection circuit suggested by drosophila connectomics.  Nature, 500(7461):175–181, 2013.
  * [24] L.R. Varshney, B.L. Chen, E. Paniagua, D.H. Hall, and D.B. Chklovskii.  Structural properties of the Caenorhabditis elegans neuronal network.  PLoS Computational Biology, 7, 2011.
  * [25] Kevin Verbeek and Subhash Suri.  Metric embedding, hyperbolic space, and social networks.  In Proceedings of the thirtieth annual symposium on Computational geometry, pages 501–510, 2014.
  * [26] Joshua Vogelstein, Randal Burns, and Ben Falk.  Neurodata: Open connectome project.
  * [27] M Winlaw, H DeSterck, and G Sanders.  An in-depth analysis of the chung-lu model.  Technical report, Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States), 2015.