A novel statistical methodology for quantifying the spatial arrangements of axons in peripheral nerves

A thorough understanding of the neuroanatomy of peripheral nerves is required for a better insight into their function and the development of neuromodulation tools and strategies. In biophysical modeling, it is commonly assumed that the complex spatial arrangement of myelinated and unmyelinated axons in peripheral nerves is random, however, in reality the axonal organization is inhomogeneous and anisotropic. Present quantitative neuroanatomy methods analyze peripheral nerves in terms of the number of axons and the morphometric characteristics of the axons, such as area and diameter. In this study, we employed spatial statistics and point process models to describe the spatial arrangement of axons and Sinkhorn distances to compute the similarities between these arrangements (in terms of first- and second-order statistics) in various vagus and pelvic nerve cross-sections. We utilized high-resolution transmission electron microscopy (TEM) images that have been segmented using a custom-built high-throughput deep learning system based on a highly modified U-Net architecture. Our findings show a novel and innovative approach to quantifying similarities between spatial point patterns using metrics derived from the solution to the optimal transport problem. We also present a generalizable pipeline for quantitative analysis of peripheral nerve architecture. Our data demonstrate differences between male- and female-originating samples and similarities between the pelvic and abdominal vagus nerves.

1. Introduction

Understanding the functionalities of the peripheral nerves and developing neuromodulation tools require an in-depth quantitative characterization of the anatomy of the nerves. A large portion of the quantitative neuroanatomical studies focus on counting and comparing the number of myelinated and unmyelinated axons in the peripheral nerves in different animals (Hoffman and Schnitzlein, 1961; Krous et al., 1985; Asala and Bower, 1986; Prechtl and Powley, 1990; Pereyra et al., 1992; Soltanpour and Santer, 1996; Safi et al., 2016). There are studies on analyzing the changes in the number of myelinated and unmyelinated axons as the function of animals' age (Krous et al., 1985; Pereyra et al., 1992; Soltanpour and Santer, 1996). The morphometric characteristics of the axons, such as area of axon cross-section, diameter, myelin thickness, are also well-developed and helpful for estimating electrode distances for neuromodulation purposes (Asala and Bower, 1986; Prechtl and Powley, 1990; Walter and Tsiberidou, 2019; Pelot et al., 2020; Havton et al., 2021; Settell et al., 2021).

The vagus is a complex, multi-functional peripheral nerve of the autonomic nervous system, containing both sensory and motor axons that regulate a wide variety of functions (Câmara and Griessenauer, 2015; Breit et al., 2018). These include regulation of the heart, respiratory tract, and many areas of the gastrointestinal system, influencing motility, secretions, and communication with the immune system. This breadth of activity and its bidirectional connectivity with the central nervous system have led to the vagus becoming a promising target for bioelectric medicine, through development of specific protocols for vagal nerve stimulation (VNS) (Bonaz et al., 2017a,b; Horn et al., 2019).

In addition to providing an alternative therapeutic approach for drug-resistant clinical conditions within organs, the vagal afferent connections to the brain provide opportunities for novel therapies directed to various psychiatric disorders. To improve the efficacy and specificity of VNS for each type of clinical condition, a greater understanding of the intra-vagal neural elements relating to each organ system is required (Howland, 2014; Thompson et al., 2019, 2023), as demonstrated by a recent study showing that fascicle-selective stimulation can reduce off-target effects of VNS (Thompson et al., 2023). This includes understanding the spatial organization of different functional classes of axons within and between fascicles. This spatial organization has not been investigated in depth within the visceral nervous system.

We have begun to address this knowledge gap using an extensive dataset of transmission electron microscopy (TEM) images derived from multiple cross-sections of the rat vagus. We have included in our study additional TEM images from the rat pelvic nerve, another multi-functional major nerve of the autonomic nervous system that supplies sensory and motor axons to the urogenital organs and lower bowel. Both TEM data sets have been published through the SPARC Portal RRID:SCR_017041 under a CC-BY 4.0 license (Plebani et al., 2022).

Cross-sections of large peripheral nerves reveal a variety of components (myelinated and unmyelinated axons, Schwann cells), high-order structures (fascicles and Remak bundles), and raise questions regarding the spatial arrangement of these components, similarities between multiple arrangements, and their relationship to various biological factors including age, sex, and diseases. This study focuses on the unmyelinated axons segmented utilizing our high-throughput deep learning model (Plebani et al., 2022). We aim to define a notion of similarity (or dissimilarity) between the spatial arrangements of the unmyelinated axons and quantify the distances between them. We resort to spatial point patterns to represent the image data conveniently for analyzing the axons' spatial organization. We use the centroids of the segmented axons to construct spatial point patterns. We consider spatial inhomogeneity and anisotropy to be the spatial features to represent the spatial arrangement of the axons, and be used for quantification. A known way to get an intuitive sense of spatial inhomogeneity (inhibition and/or attraction between points) and anisotropy of a point pattern is to investigate its second-order statistics (Ripley, 1976, 1977; Sengupta et al., 2013; Dixon, 2014). Since Ripley's summary of spatial statistical methods in 1977 the techniques for spatial pattern analysis have been occasionally employed in neuroscience, often by statisticians who saw the extraordinary complexity of neuroanatomical patterns to be a perfect demonstration of the spatial statistics inference ability (Bjaalie et al., 1991; Diggle et al., 1991; Prodanov et al., 2007; Jafari-Mamaghani et al., 2010; Waller et al., 2011).

Although the standard spatial statistical measures can quantify overall global differences between point interactions, they are not well suited for calculating distances between complex non-random patterns with multiple distinct local interactions. Therefore, we propose a method that involves computing the local second-order spatial statistics for the nerve fascicles to capture their spatial arrangement and utilizing a revised optimal transport distance (Sinkhorn distance) to measure similarities between the second-order spatial statistics of every pair of nerve cross-sections. We visualize the resulting Sinkhorn distance matrix in a new metric space using multi-dimensional scaling that helps interpret the similarity (dissimilarity) of the spatial features in the nerve cross-sections. Our concept of Sinkorn distance embedding was influenced by related work on optimal transport-based morphometry applications in cell biology (Wang et al., 2013; Basu et al., 2014).

In addition to addressing a neuroscientific problem of quantifying the vagus nerve anatomy with computer science tools, we intend to bring together a variety of approaches from various computer science domains to extend the toolkit for point-pattern comparisons in biology. Utilizing the optimal transport framework to establish a quantitative measure for spatial point patterns and advance the workflow of quantitative analysis of peripheral nerve architecture, our method could be used beyond neuroscience. We believe that our findings contribute to the establishment of spatially selective stimulation of nerve axons to improve the efficacy of VNS.

Figure 1 depicts a high-level overview of the steps in the quantitative analysis of the spatial arrangement of axons in the peripheral nerve architecture presented in this paper. We explain each component of the pipeline in the following sections. The biological data and the data preprocessing steps are described in Section 2.1. The basics of spatial point pattern, spatial statistics, and optimal transport framework are introduced in Section 2.2 and Section 2.3, respectively. We provide details on the experimental setup and results in Section 3, which cover steps 2, 3, and 4 of the computational pipeline. We discuss the results of the empirical study in Section 4 before concluding.

FIGURE 1

Figure 1. The pipeline for the quantitative analysis of the spatial arrangement of axons in the peripheral nerve cross-sections.

2. Materials and methods

2.1. Biological data and automated segmentation

Although the vagus nerve anatomy was the primary motivation behind this study, we use the TEM images of the vagus and pelvic nerve cross-sections in rats for comparisons. A list of the TEM images used in this study is shown in Table 1. The protocols and techniques followed for nerve sample collection, processing, and imaging are documented in our previous work (Plebani et al., 2022). The data is publicly available via NIH-supported SPARC Pennsieve database (Havton et al., 2022). Briefly, the unmyelinated axons in some of these TEM images were manually annotated and used as labeled data to train, validate, test, and evaluate an automated segmentation model based on the U-Net architecture (Ronneberger et al., 2015; Plebani et al., 2022). The segmentation model is a U-Net with four stages: the convolutional layers have a batch normalization layer followed by a ReLU activation layer, and the bottleneck stage has extra dropout layers between convolutions (Plebani et al., 2022). The model classifies the TEM image pixels as one of the three following classes: (a) fiber if it is inside an unmyelinated axon, (b) border if it is in a boundary region between an axon and the rest of the image defined by the outer edge of each axon, and (c) background. An updated version of the model¹ was used here to segment the unmyelinated axons. The resulting axon counts are listed in Table 1. We used the open-source image processing package Fiji (Schindelin et al., 2012) to extract the centroid coordinates of the segmented unmyelinated axons and the functions in R packages to establish the outer boundaries and inner void spaces of the nerve cross-sections, to construct the spatial point patterns. Images 15 (vagus) and 29 (pelvic) listed in Table 1 are shown in Figures 2A, D, respectively, along with their corresponding automated segmentations and spatial point patterns.

TABLE 1

Table 1. The list of the TEM images of vagus and pelvic nerve cross-sections in rats used in this study.

FIGURE 2

Figure 2. (A, D) The transmission electron microscopy (TEM) images of the nerve cross-section of Image 15 (vagus) and Image 29 (pelvic) listed in Table 1, respectively. The visible void spaces in the nerve cross-sections are blood vessels. The tiny light gray regions without any border are the unmyelinated axons. The myelinated axons have slightly darker gray borders. (B, E) The automated segmentation of the unmyelinated axons (the white regions) in the nerve cross-sections. (C, F) The spatial point patterns constructed with the centroid locations (the black circles) of the segmented unmyelinated axons.

Before reporting the experimental details, we provide a brief overview of spatial point patterns, spatial statistics concepts, and optimal transport framework in the following two subsections.

2.2. Spatial point patterns and spatial statistics

A spatial point pattern (SPP) is a set of spatial locations associated with entities of interest in 2-D or 3-D space, encompassed by an observation window (Møller and Waagepetersen, 2003; Stoyan, 2006; Jafari-Mamaghani et al., 2010; Baddeley et al., 2015). The objective of SPP analysis is to examine the spatial arrangement of the points in an SPP and recognize trends that define the point pattern. Two fundamental descriptive characteristics of an SPP are intensity and interaction. The intensity ρ or ρ(u) of a point pattern, a first-order statistics, is the average number of points per unit area, and it can be uniform across the observation window (homogeneous), or it can vary according to an intensity function (inhomogeneous). A point pattern's intensity is usually denoted by λ in literature, but we use ρ to avoid confusion with another notation related to the optimal transport problem. Three SPPs of different average intensity, with 20 (sample 1), 100 (sample 2), and 200 (sample 3) points per unit area are illustrated in Figure 3.

FIGURE 3

Figure 3. An illustration of spatial point patterns with different spatial intensities. Samples 1, 2, and 3 have 20, 100, and 200 points per unit area, respectively.

The interaction is associated with a distance r and describes the influence the points have on their neighbors within r radius. The interaction is termed complete spatial randomness (CSR) if the points are independent. The points can exhibit positive interaction (spatial attraction), negative interaction (spatial inhibition), or a combination of both.

It is common practice to use second-order statistics such as Besag's centered L-function, which is a transformation of Ripley's K-function (Ripley, 1976, 1977; Besag, 1977), to investigate the interaction in point patterns. Let X be a point pattern and t(u, r, X) be the number of points in X which lie within distance r of the location u. Assuming X is a homogeneous point pattern with intensity ρ, the number of points within distance r of a specific point is represented by ρK(r) (Ripley, 1976).

$\begin{array}{l}\hfill \text{Ripley's }K\text{-function}:K\left(r\right)\text{=}\frac{𝔼\left[t\left(u,r,X\right)\text{|}u\in X\right]}{\rho },\hfill \\ \hfill \text{Besag's centered }L\text{-function}:L(r)\text{=}\sqrt{\frac{K\left(r\right)}{\pi }}-r\hfill & (1)\end{array}$

An estimator for the empirical K-function $\widehat{K}(r)$ is formulated in Equation (2), and represents the cumulative average number of neighbors within r radius of a typical point, standardized by the intensity and corrected for edge effects.

$\begin{array}{ll}\widehat{K}(r)=\frac{W}{n(n-1)}{\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{j=1,j\ne i}^n}\nu (d_{ij}⩽r)e_{ij}(r)& (2)\end{array}$

where ν(.) is an indicator function that equals 1 if the argument is true and otherwise is 0. Here n is the number of points; W is the area of the observation window; r is the interaction distance; d_ij is the Euclidean distance between x_i and x_j; and e_ij denotes weights for edge correction (Baddeley et al., 2015). The K- and L-functions can illustrate the non-random spatial arrangement of the points if compared with CSR. They are invariant to the intensity of a point pattern and to missing random points (Ripley, 1976; Baddeley et al., 2000), which allows these second-order statistics to be compared when the number of points and observation window vary in the point patterns under consideration. Positive values of the centered L-function depict spatial attraction, and negative values describe spatial inhibition in a point pattern.

When dealing with inhomogeneous point patterns, an inhomogeneous L-function, based on the inhomogeneous K-function K_inhom(r) can be evaluated. The estimator for the inhomogeneous K-function ${\widehat{K}}_{\mathrm{\text{inhom}}}(r)$ is formulated in Equation (3) below.

$\begin{array}{l}{\widehat{K}}_{\text{inhom}}\left(r\right)=\frac{1}{D^{p}W}{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1,j\ne i}^n\frac{\nu (d_{ij}⩽r)}{\widehat{\rho }\left(x_i\right)\widehat{\rho }\left(x_j\right)}}}{e}_{ij}(r),\\ D^p={\left(\frac{1}{W}{\displaystyle \sum _{i=1}^n\frac{1}{\widehat{\rho }\left(x_i\right)}}\right)}^p,p\in \left\{1,2\right\}& (3)\end{array}$

where $\widehat{\rho }(u)$ is an estimator of the intensity function ρ(u), obtained using a kernel-smoothed (described later in the paper) intensity estimator (Baddeley et al., 2015). Figure 4A shows three SPPs portraying complete spatial randomness (Sample 1), spatial inhibition (Sample 2), and spatial attraction (Sample 3). These point patterns were formed by the Poisson process, the hardcore process, and the Matern cluster process (Baddeley et al., 2015). The L-function illustrates the differences in spatial interaction of these SPPs, as shown in Figure 4B. The blue curve's fluctuation about zero within its significance band indicates that the point pattern in Sample 1 is most likely random, but the green and orange curves' prominent negative and positive peaks imply spatial inhibition and attraction in Samples 2 and 3, respectively.

FIGURE 4

Figure 4. An illustration of spatial point patterns with different spatial interactions. (A) SPPs portraying complete spatial randomness (sample 1), spatial inhibition (sample 2), and spatial attraction (sample 3). (B) Besag's centered L-function computed for the patterns shown in (A). The solid lines illustrate the spatial interaction of the patterns compared to CSR. The shaded area around the solid lines shows the boundaries of 95-percentile confidence interval.

Similar analysis of spatial interaction can be done for the biological point patterns as well. Figure 5 shows the Besag's centered inhomogeneous L-function computed for Image 15 (vagus) and Image 29 (pelvic), listed in Table 1 and displayed in Figure 2. For both samples, this demonstrates spatial inhibition over a narrow interaction distance range, followed by spatial attraction. The clustering tendency (spatial attraction) is more pronounced in the pelvic sample.

FIGURE 5

Figure 5. Besag's centered inhomogeneous L-function computed for Images 15 and 29. The shaded area around L(r) = 0 shows the significance bands of complete spatial randomness (CSR). The solid lines illustrate the non-random spatial arrangement of the point patterns compared to CSR. The shaded area around the solid lines shows the boundaries of 95-percentile confidence interval.

An SPP is anisotropic if any of its statistical characteristics change when the point pattern is rotated about any axis in 2-D or 3-D space. The K- and L-functions can be modified in various ways to estimate anisotropy (Ohser and Stoyan, 1981; Chiu et al., 2013). Computing the cumulative distribution of the neighbors within a section of the disc of r radius between two directional preferences θ₁ and θ₂, instead of the entire disc of r radius, gives the sector K- and L-functions (Baddeley et al., 2015). Figure 6A depicts three random point patterns with no preferential direction (Sample 1) as well as horizontal (Sample 2) and vertical anisotropies (Sample 3). We compute sector K-functions for sectors forming a 15° segment around the horizontal (0°) and vertical (90°) axes, and if these two functions are not approximately equal, we conclude that the point patterns are anisotropic (Baddeley et al., 2015). Figure 6B shows the differences between the horizontal and vertical K-functions for the three aforementioned point patterns. We observe that the difference between the sector K-functions fluctuates around zero for Sample 1, indicating the absence of anisotropy, holds positive values for Sample 2 and negative values for Sample 3, indicating horizontal and vertical preferences, respectively, up to a certain interaction distance.

FIGURE 6

Figure 6. (A) Inhomogeneous random spatial point patterns with no directional preference (Sample 1), horizontal (Sample 2), and vertical (Sample 3) directional preferences. (B) Differences computed between the horizontal and vertical K-functions for the samples in (A).

When an SPP exhibits different interactions in different places, it is beneficial to look into the spatial statistics locally by decomposing them into contributions from individual points (Baddeley et al., 2015). If the K-function estimator $\widehat{K}(r)$ shown in Equation (2) is decomposed, the contributions from individual points are referred to as local K-functions and can be formulated as follows:

$\begin{array}{ll}\widehat{K}(r,x_i)=\frac{W}{n-1}{\displaystyle \sum _{j=1,j\ne i}^n}\nu (d_{ij}⩽r)e_{ij}(r),\text{ for }i=1,\dots ,n& (4)\end{array}$

The estimator $\widehat{K}(r)$ is simply the average of all the $\widehat{K}(r,x_i)$s for i = 1, …, n. The centered local L-functions are formulated as $\widehat{L}(r,x_i)={(\widehat{K}(r,x_i)/\pi )}^{1/2}-r$. This notion of decomposition is applicable for the inhomogeneous and anisotropic K- and L-functions as well.

The points in an SPP may be of different types (multitype point pattern). The additional information attached to each point in point patterns is called a mark and can hold categorical or continuous-valued, physical or statistical characteristics. The points can carry additional attributes (forming marked point patterns) or be linked to the space of interest (covariates). It is often helpful to apply spatial smoothing to the marks of a point pattern for visualization and various post-processing purposes. The result of the kernel-smoothing (usually with Gaussian kernel) at a location u is a spatially weighted average of the marks attached to the points in the neighborhood of u (Baddeley et al., 2015). This is also known as the Nadaraya-Watson smoother (Nadaraya, 1964, 1989; Watson, 1964).

There are several partially similar, but differently interpreted spatial statistical functions employed to describe the dependence (K-function, L-function, pair correlation function) and the spacing (nearest-neighbor function G, empty-space function F, and their combination called the J-function) between points in an SPP (Baddeley et al., 2015). Although we used local inhomogeneous and anisotropic L-functions to describe the spatial arrangement of unmyelinated axon point patterns in nerve cross-sections, our method could be easily reimplemented utilizing pair correlation function (PCF) in place of the L-functions. PCF, which is related to K and L-functions, is of particular interest because it has previously been used in the context of biological microscopy (Sengupta et al., 2011, 2013; Veatch et al., 2012).

In addition intensity (or spatial density) and interaction, regionality is another aspect of the spatial organization that influences spatial descriptors and contributes to the concept of spatial similarity (or location). Regionality is the most intuitively understandable feature of spatial organization and denotes the absolute position of structures of interest within an object or region of interest (ROI). Regionality is not translation- and rotation-invariant, whereas interaction and intensity are (if not anisotropic). All of these characteristics are captured when spatial point patterns are processed using the tools described in this report.

2.3. Optimal transport framework and Sinkhorn distance

The transport problem distributes a certain amount of mass from a set of sources to a set of destinations at minimum cost. There are two major factors in a transport problem: the cost function and the transportation plan. The cost function defines a fixed, non-negative effort required to transport unit mass from a source to a destination. This cost may only depend on the distance between the source and the destination or on other additional factors; in the former case a Euclidean distance matrix between the sources and the destinations is a reasonable representation of effort.

Once the cost of transportation is represented, the remaining part of the problem involves transporting a non-negative amount of mass between sources and destinations, as described by a transportation plan. Various transportation plans result in different total costs, and the optimal transport problem (OT) aims to minimize this cost. The OT problem is balanced if the total mass at the sources equals the total mass at the destinations, and unbalanced otherwise (Peyré and Cuturi, 2019).

Let r and c be two d dimensional vectors representing the amount of mass at the d sources and the d destinations, respectively. The number of sources and destinations could differ, but they can be considered equal without loss of generality. Let U(r, c) be the set of all non-negative d × d matrices with row and column summing to r and c, respectively. Any matrix P ∈ U(r, c) describes a transportation plan that transports the mass in r to c. Given a d × d cost matrix M, the total cost of mapping r to c using the transportation plan P is $\sum _{i,j}{P}_{ij}{M}_{ij}$. Thus the OT problem between r and c given cost M can be formulated by Equation (5), where D_M(r, c) is the optimal transport distance:

$\begin{array}{l}\hfill D_M(r,c)=\underset{P\in U(r,c)}{\min }{\displaystyle \sum _{i,j}{P}_{ij}}{M}_{ij},\hfill \\ \hfill \text{subject to }{\displaystyle \sum _j{P}_{ij}}=r_i,{\displaystyle \sum _i{P}_{ij}}=c_j,P_{ij}\ge 0,\forall i,j\le d.\hfill & (5)\end{array}$

The masses in r and c could be normalized to sum to one, and then both r and c can be interpreted as probability distributions.

For D_M(r, c) to be a metric, the cost matrix M has to be a metric matrix (Avis, 1980; Brickell et al., 2008; Villani, 2009) satisfying the conditions shown in Equation (6).

$\begin{array}{l}\hfill \text{Non-negativity: }{M}_{ij}\ge 0,\hfill \\ \hfill \text{Identity: }{M}_{ii}=0,\hfill \\ \hfill \text{Symmetry: }{M}_{ij}=M_{ji},\hfill \\ \hfill \text{Triangle inequality: }{M}_{ij}\le M_{ik}+M_{kj},\forall i,j,k\le d.\hfill & (6)\end{array}$

The OT is a convex optimization problem that can be solved using various approaches (Ahuja et al., 1993; Orlin, 1993). For a general cost matrix the computational cost scales as ${O}\left[d^3\mathrm{\text{log}}(d)\right]$ (Pele and Werman, 2009), which prevents scaling the solution to large problem sizes. Earlier approximate solutions obtained by putting constraints on the cost matrix could result in a loss of applicability and performance (Grauman and Darrell, 2004). A later approximation to the original OT problem using an entropic regularization scheme was proposed by Cuturi (2013) to reduce the computational complexity. The scheme employs the Sinkhorn-Knopp matrix scaling algorithm (Sinkhorn and Knopp, 1967; Knight, 2008), and hence the name Sinkhorn distance for its objective function.

2.3.1. Sinkhorn distance

A straightforward way of thinking about a transportation plan is by noticing that if a source contains more mass, it should originate more, and if a destination requires more mass, it should receive proportionally more. Such a transportation plan is represented by rc^T, and the optimal plan P should be somewhere around the distribution rc^T. Simply speaking, the idea of the entropic regularization scheme by Cuturi (2013) is to choose P from a smaller set near rc^T, instead of the entire set U(r, c).

To capture these ideas, Cuturi (2013) imposes an additional constraint of Kullback-Leibler (KL) divergence on the OT formulation, as shown in Equation (7), and computes the Sinkhorn distance $D_{M,\alpha }^{*}(r,c)$. This constraint introduces a set U_α(r, c) ⊂ U(r, c) from which an optimal transportation plan P is selected. The KL divergence distance between P and rc^T is set to be smaller than a predefined parameter α. In other words, P should belong to a distribution near rc^T.

$\begin{array}{l}\hfill D_{M,\alpha }^{*}(r,c)=\underset{P\in U_{\alpha }(r,c)}{\min }{\displaystyle \sum _{i,j}{P}_{ij}}{M}_{ij},\hfill \\ \hfill \text{subject to }KL(P|r{c}^T)\le \alpha ,{\displaystyle \sum _j{P}_{ij}}=r_i,{\displaystyle \sum _i{P}_{ij}}=c_j,\forall i,j\le d.\hfill & (7)\end{array}$

The entropy (h) of the transportation plan (P) and the mass vectors (r and c) are given in Equation (8):

$\begin{array}{l}\hfill h(P)=-{\displaystyle \sum _{ij}{P}_{ij}}\log P_{ij},\hfill \\ \hfill h(r)=-{\displaystyle \sum _i{r}_i}\log r_i,h(c)=-{\displaystyle \sum _j{c}_j}\log c_j.\hfill & (8)\end{array}$

We proceed to express the KL divergence constraint in terms of the entropy:

$\begin{array}{l}KL(P|r{c}^T)={\displaystyle \sum _{ij}{P}_{ij}}\log \frac{P_{ij}}{r_i{c}_j}\\ ={\displaystyle \sum _{ij}{P}_{ij}}\log P_{ij}-{\displaystyle \sum _{ij}{P}_{ij}}\log r_i-{\displaystyle \sum _{ij}{P}_{ij}}\log c_j\\ ={\displaystyle \sum _{ij}{P}_{ij}}\log P_{ij}-{\displaystyle \sum _i{r}_i}\log r_i-{\displaystyle \sum _j{c}_j}\log c_j\\ [\because {\displaystyle \sum _j{P}_{ij}}=r_i,{\displaystyle \sum _i{P}_{ij}}=c_j]\\ =-h(P)+h(r)+h(c)\le \alpha .& (9)\end{array}$

Thus the new constraint states that the entropy of P should be large enough to satisfy

$h(P)\ge h(r)+h(c)-\alpha ,$

which constrains P to be chosen from the Kullback-Leibler ball of level α centered about rc^T (see Figure 1 in Cuturi, 2013).

This interpretation makes the OT problem non-convex, and an alternative formulation of Sinkhorn distance is required for ease of optimization. For every pair (r, c), each α corresponds to a Lagrange multiplier λ ∈ [0, ∞) such that $D_{M,\alpha }^{*}(r,c)=D_M^{\lambda }(r,c)$. The distance $D_M^{\lambda }$, shown in Equation (10), is called the dual-Sinkhorn divergence by Cuturi (2013).

$\begin{array}{l}\hfill D_M^{\lambda }(r,c)={\displaystyle \sum _{i,j}{P}_{ij}^{\lambda }}{M}_{ij},\text{ where }{P}^{\lambda }=\underset{P\in U(r,c)}{\mathrm{argmin}}{\displaystyle \sum _{i,j}{P}_{ij}}{M}_{ij}-\lambda h(P),\hfill \\ \hfill \text{subject to }{\displaystyle \sum _j{P}_{ij}}=r_i,{\displaystyle \sum _i{P}_{ij}}=c_j,\forall i,j\le d.\hfill & (10)\end{array}$

By introducing two dual variables ϕ and ψ for each of the two equality constraints of Equation (10), the Lagrangian of the objective function can be written as Equation (11).

$\begin{array}{l}{L}(P,\varphi ,\psi )={\displaystyle \sum _{i,j}{P}_{ij}}{M}_{ij}-\lambda h(P)+{\displaystyle \sum _i{\varphi }_i}({\displaystyle \sum _j{P}_{ij}}-r_i)\\ +{\displaystyle \sum _j{\psi }_j}({\displaystyle \sum _i{P}_{ij}}-c_j).& (11)\end{array}$

The derivative of the Lagrangian objective function with respect to P_ij, for any pair (i, j), can be set to zero to obtain an extremum; the second derivative of the Lagrangian, ($\frac{\lambda }{P_{ij}}$), is positive since both the numerator and the denominator are positive, and thus we have obtained a minimizer of the Lagrangian.

$\begin{array}{l}\frac{\partial {L}}{\partial P_{ij}}=M_{ij}+\lambda +\lambda \log P_{ij}+{\varphi }_i+{\psi }_j=0.\\ \Rightarrow P_{ij}=e^{-\frac{{\varphi }_i}{\lambda }-\frac{1}{2}}.e^{-\frac{M_{ij}}{\lambda }}.e^{-\frac{{\psi }_j}{\lambda }-\frac{1}{2}}\\ \equiv u_i{K}_{ij}{v}_j[u_i=e^{-\frac{{\varphi }_i}{\lambda }-\frac{1}{2}},v_j=e^{-\frac{{\psi }_j}{\lambda }-\frac{1}{2}},K=e^{-\frac{M}{\lambda }}]& (12)\end{array}$

Given K, r, and c, the Sinkhorn-Knopp matrix scaling algorithm converges to a solution P^λ of the following form:

$\begin{array}{ll}\exists u,v:P^{\lambda }=\text{diag}(u)K\text{diag}(v).& (13)\end{array}$

P^λ should have the correct row and column sums, as shown in Equation (10). We deduce the update rule for the Sinkhorn-Knopp algorithms from those constraints in the following manner:

$\begin{array}{l}{\displaystyle \sum _j{P}_{ij}^{\lambda }}=r_i,{\displaystyle \sum _i{P}_{ij}^{\lambda }}=c_j\\ \Rightarrow {\displaystyle \sum _j{u}_i}{K}_{ij}{v}_j=r_i[\text{Equation 12}]\Rightarrow {\displaystyle \sum _i{u}_i}{K}_{ij}{v}_j=c_j\\ \text{ [Equation 12}]\\ \Rightarrow u_i{\displaystyle \sum _j{K}_{ij}}{v}_j=r_i\Rightarrow v_j{\displaystyle \sum _i{u}_i}{K}_{ij}=c_j\\ \Rightarrow u_i=r_i/{\displaystyle \sum _j{K}_{ij}}{v}_j\Rightarrow v_j=c_j/{\displaystyle \sum _i{u}_i}{K}_{ij}\end{array}$

Thus the update rule for the Sinkhorn-Knopp algorithm can be written as Equation (14), where v can be initialized randomly.

$\begin{array}{l}u=r./(Kv),\\ v=c./(K^{T}u).& (14)\end{array}$

Cuturi (2013) observes that the number of iterations in the Sinkhorn-Knopp algorithm is bounded independent of d. Thus, the cost of computing $D_M^{\lambda }$ is ${O}(d^2)$, which is an improvement over ${O}\left[d^3\mathrm{\text{log}}(d)\right]$. Cuturi (2013) describes an approach to compute the Sinkhorn distance $D_{M,\alpha }^{*}(r,c)$ through the dual-Sinkhorn divergence $D_M^{\lambda }(r,c)$, and also reports that the dual-Sinkhorn divergence does not perform worse than the classic optimal transport distances. Therefore, we use the dual-Sinkhorn divergence to measure the distance between the spatial statistics of the point patterns in our experiments and refer to as the Sinkhorn distance. We utilize the R packages T4transport (You, 2022) and Barycenter (Klatt, 2018) for computing the dual-Sinkhorn divergences, and spatstat (Baddeley et al., 2015) for spatial point pattern analysis.

3. Experiments and results

We represent the unmyelinated axonal arrangements in the vagus and pelvic nerve cross-sections as spatial point patterns. We intend to quantify (using the Sinkhorn distance) similarities between the point patterns in terms of the following spatial features:

1. spatial intensity,

2. local inhomogeneous L-function,

3. local inhomogeneous anisotropic L-function with

(a) horizontal and

(b) vertical sectors.

For the horizontal and vertical cases, we choose sectors forming 15° segment around the horizontal (0°) and vertical (90°) axes, respectively. We attach the above-mentioned spatial features to the point patterns as marks (described in Section 2.2). We compute the Sinkhorn distance between every pair of point patterns, for the four spatial features, in two different manners:

1. using the spatial point patterns directly (in Section 3.4) and

2. using the map of the spatial features constructed by kernel-smoothing (in Section 3.5).

The Sinkhorn distance between every pair of nerve cross-sections is then used to construct a symmetric Sinkhorn distance matrix and visualized in an embedded space via multi-dimensional scaling. In the following three subsections, we describe a few preprocessing and parameter selection tasks required for configuring the spatial features, before going into the experimental details.

3.1. Interaction distance configuration

A critical issue regarding the computation of local inhomogeneous L-functions is determining the interaction distance (r), described in Section 2.2. The point patterns constructed from the nerve cross-sections differ in size, as do their interaction ranges. The preferred choice for the interaction distance is the one that can reasonably separate the spatial features of interest present in the point patterns. We compute a range of interactions that are common for all the point patterns and configure the interaction distance using two approaches: (a) based on the standard deviations of the inhomogeneous L-function of all the point patterns and (b) based on the F-ratio (analysis of variance) of the inhomogeneous L-function of the point patterns grouped as vagus vs. pelvic, within the expected range. Other strategies for choosing r or a linear combination of multiple r values are also possible (See Section 3.4.1).

3.2. Translation and rotation normalization

The optimal transport distance is not invariant under translations and rotations (Wang et al., 2013). This is critical in the case of analyzing the point patterns because spatial inhomogeneity and anisotropy depend largely on the placement and orientation of the point patterns. To provide translation invariance, we scale the point patterns maintaining proportionality, align the center of mass to the origin, and apply the necessary 0-padding around the biological structures.

Ensuring rotation invariance is non-trivial. In an ideal setting, the information regarding the orientation of the biological structures would be available directly to the analyst. Unfortunately, the experimental and instrumental setting may not always allow the orientation of the samples to be maintained during the specimen preparation and the imaging process. Therefore, we implemented a post-hoc minimization process as a workaround. While computing the Sinkhorn distance between the spatial intensities of a pair of point patterns, we keep the orientation of one of them unchanged and rotate the other one about the origin by multiple θ values (θ= 45° in our experiments). We compute the Sinkhorn distance for all possible values of θ and keep the orientation that provides the smallest Sinkhorn distance result. We use this identified orientation for the computation of other spatial features. This method provides a reproducible procedure in the absence of known anatomical orientation data.

3.3. The entropic regularization parameter

We refer to the coefficient of the entropy of the transportation plan h(P) in the dual-Sinkhorn divergence formulation shown in Equation (10), λ, as the entropic regularization parameter. As λ → 0, Sinkhorn distance approaches the optimal transport distance (Wasserstein distance, provided that the cost is Euclidean distance). As λ increases, the computation results in different approximations of the optimal transport distance, i.e., the Sinkhorn distances. Tuning the appropriate entropic regularization parameter is an important task. We can consider two scenarios: (a) selecting an entropic regularization parameter that provides better separation between analyzed instances and (b) selecting an entropic regularization parameter that makes the Sinkhorn distance a more accurate approximation of the exact optimal transport distance. Thus, parameter tuning is a trade-off between favoring the utility of the method (and lower computational cost) and the accuracy of the approximation.

Smaller values of the λ parameter produce Sinkhorn distances that more closely approximate Wasserstein distances (Cuturi, 2013). Therefore, if computational resources are not limited, smaller λ values are preferable to larger values. To estimate the computational cost, we tried λ = 0.01, 0.05, 0.1, 0.5, 1.0, 2.0, 5.0, but we settled on λ = 0.01.

3.4. Sinkhorn distance between spatial point patterns

In this section, we compute Sinkhorn distance between the spatial point patterns (directly) of the nerve cross-sections, for the four spatial features mentioned earlier. Let S₁ and S₂ be two spatial point patterns, with n₁ and n₂ number of points respectively. Considering an optimal transport problem between S₁ and S₂, we assume that each point in S₁ contains $\frac{1}{n_1}$ amount of mass (r), therefore the total mass $=\frac{1}{n_1}\times n_1=1$. Similarly, each point in S₂ requires $\frac{1}{n_2}$ amount of mass (c), so the total mass $=\frac{1}{n_2}\times n_2=1$. Thus the problem is to transport the mass from S₁ to S₂ (balanced). The inputs for the computations are the spatial locations of the points in S₁ and S₂. Therefore, we can compute the Euclidean distance matrix between them to be the cost matrix M, of dimension n₁ × n₂. Here, the cost matrix M captures the spatial intensity of the point patterns. Further, we compute the transportation plan P, of dimension n₁ × n₂, using the formulation described in Section 2.3.1 and obtain the Sinkhorn distance, which provides a measure of similarity of the spatial intensity between S₁ and S₂.

In the cases of the three other spatial features, the cost matrix (M) remains the same (the spatial location of the points are unchanged), but the amount of mass produced (r) or required (c) at each point changes. The local inhomogeneous L-function attaches a numerical value to each point in a point pattern that captures its local spatial interaction within a certain interaction distance. Instead of a uniform mass amount, we may assume that each point in S₁ and S₂ is assigned an amount of mass that equals its local inhomogeneous L-function value. We normalize the values assigned to each point pattern to sum to one. Then we can compute the transportation plan P in the same manner as described above. The resultant Sinkhorn distance gives us a measure of similarity of the local spatial interaction between S₁ and S₂. The Sinkhorn distances for the anisotropic spatial features, the local inhomogeneous anisotropic L-function with horizontal and vertical sectors, are also computed similarly.

We have 29 nerve cross-sections in the dataset and once we compute the Sinkhorn distance between every pair, we can construct Sinkhorn distance matrices of dimension 29 × 29, for each of the four spatial features. We use multi-dimensional scaling to embed the Sinkhorn distance matrices in 2-D to illustrate the computed Sinkhorn distance between the spatial features and interpret the notion of similarity (or dissimilarity) between the nerve cross-sections. We denote the new embedded 2-D space as the Sinkhorn space.

The application of the optimal transport (OT) solution to define the similarity of spatial point patterns is dependent on all of the spatial organization characteristics described in Section 2.2. Intensity differences are conveyed by the number of points at a given location relative to other sites where mass for transportation is present. As point weights in OT, the values of the local (inhomogeneous and anisotropic) L-functions represent information about point-point interactions. Lastly, the movement of mass from one area in ROI to another directly captures the sense of regionality. Although it is possible to simulate pattern arrangements that emphasize only one aspect of spatial organization while diminishing the influence of the others, this will not lead to the development of an intuitive sense of computed distance in a real-world setting. As with numerous other mathematical concepts (and distance in particular), the interplay between the various aspects of spatial organization manifested in biological samples defies simple models.

The architecture of the unmyelinated axons in the vagus and pelvic nerve cross-sections are examples of spatial complexity arising from the combination of multiple aspects of organization. Therefore, it might be difficult to develop an intuitive sense of the Sinkhorn distance between the axonal organization of these structures and their placement in the Sinkhorn space. Before discussing biological point patterns, we present a number of simulated point patterns to illustrate only a few particular scenarios emerging as realizations of pre-defined spatial point-pattern processes (they do not represent the entire landscape of possible phenotypic manifestations).

3.4.1. Sinkhorn distance for simulated point patterns

Figures 7A–C show 18 simulated spatial point patterns with different spatial interactions—inhibition, randomness, and clustering. These point patterns were formed by inhomogeneous point processes, namely the hardcore process, the Poisson process, and the Matern cluster process (Baddeley et al., 2015). For the simulation, we used a hardcore process with 400 points per unit area and a hardcore distance of 0.04 (the points are not allowed to be within 0.04 unit of distance from each other, ensuring inhibition). We used a Poisson process with inhomogeneous intensity function ρ_r(x, y), shown in Equation (15), and a Matern cluster process with inhomogeneous intensity function ρ_c(x, y) for the cluster centers, shown in Equation (16), with cluster radius 0.10 and 25 points per cluster. The intensity functions ρ_r(x, y) and ρ_c(x, y) introduce some anisotropy in the corresponding point patterns. The embedding of the local inhomogeneous L-function of the simulated point patterns in the Sinkhorn space is depicted in Figure 7D, where point patterns with similar spatial interaction are mapped close to each other and a clear separation can be seen between samples with different spatial interaction.

$\begin{array}{ll}{\rho }_r(x,y)=100\times exp(-5x)& (15)\end{array}$

$\begin{array}{ll}{\rho }_c(x,y)=10\times exp(2\left|x\right|-1)& (16)\end{array}$

FIGURE 7

Figure 7. An illustration of spatial point patterns with different inhomogeneous spatial interactions. (A) Spatial inhibition. (B) Spatial randomness. (C) Spatial clustering. The intensity functions of the random and clustered patterns introduce some anisotropy. (D) An embedding of the point patterns in the Sinkhorn space. The inhibited, random, and clustered point patterns are shown in green, blue, and orange, respectively.

Similar simulations can be used for any type of spatial point pattern, providing an explainable, semi-mechanistic rationale for the emergence of the patterns and an interpretable representation of their properties and reasons for separation.

Figure 8A depicts a set of spatial point patterns with different interactions and regionality. The points are concentrated in the upper right (examples 1, 2, 5, 6, 9, 10, 13, and 14) and the lower left (examples 3, 4, 7, 8, 11, 12, 15, and 16) corners to demonstrate difference in regionality. The points are organized randomly in the odd-numbered examples, and clustered in the even-numbered examples, within their corresponding regions. Figures 8B–D illustrate the notion of similarity capturing interaction and regionality as the two key aspects of spatial organization.

FIGURE 8

Figure 8. An illustration of spatial point patterns with different interactions and regionality. (A) Simulated examples of spatial point patterns: concentrated in the upper right (1, 2, 5, 6, 9, 10, 13, 14) and the lower left (3, 4, 7, 8, 11, 12, 15, 16) corners demonstrate regionality. The points are organized randomly in the odd-numbered examples and clustered in the even-numbered instances. (B) Euclidean distance between the inhomogeneous L-functions of the simulated models. (C) Sinkhorn distance between the local inhomogeneous L-functions of the simulated examples at a large interaction distance (r = 0.679). (D) Sinkhorn distance between the first principal component (PC) of local inhomogeneous L-functions of the simulated examples over a set of interaction distances. The point patterns of different interactions and regionality are shown in different colors.

When averaged for all the points, the second-order spatial statistics, such as the inhomogeneous L-functions, describe the overall spatial interaction of a point pattern (as shown in Figures 4B, 5) but lose the ability to capture the notion of regionality. Therefore, by computing pairwise Euclidean distances between the simulated examples' inhomogenous L-functions and embedding the results in the Euclidean space, we communicate only the distance between point-pattern interactions while excluding the regionality contribution entirely. See Figure 8B, where the random patterns are primarily present on the right side of the plot, whereas the clustered patterns lie on the left.

Figure 8C illustrates another valid concept of similarity but constructed with a different emphasis. In this example, the local inhomogeneous L-functions were computed with a very large interaction distance. Therefore, the effect of regionality dominates. The resulting Sinkhorn embedding shows patterns with points concentrated on one side of ROI separating from the patterns in which points were concentrated on the other side of ROIs. Changing the interaction distance r at which the local L-functions are computed provides flexibility on how much influence of regionality is incorporated into the statistics and, correspondingly, to which extent the notion of similarity is shaped by regionality vs. point-point interactions.

Alternatively, one can employ singular value decomposition to compress all the information regarding interactions at various r values. In this case, the Sikhorn distance would operate on the linear combination of r distances that convey the largest variance. This case is illustrated in Figure 8D, where both aspects of the spatial organization are employed to compute point-pattern similarity. Consequently, the visualized distances between point patterns are influenced by the interaction and the regional organization. The degree to which different aspects of spatial architecture and interaction distance should shape the metric is a choice that must be made based on domain knowledge regarding the anatomy and significance of varying levels of structural organization. Finally, it is worth mentioning that the Sinkhorn distance between two point patterns with similar interaction but different regionality might be approximately equal to the Sinkhorn distance between two point patterns with vastly different interaction characteristics but similar regionality.

3.4.2. Sinkhorn distance for biological point patterns

Now that we have an idea of how different aspects of the spatial organization can be captured with Sinkhorn distance and visualized in the Sinkhorn space, we describe the experimental results of the biological point patterns. Figure 9B shows an embedding of the spatial intensity of the spatial point patterns (used directly) of the nerve cross-sections in the Sinkhorn space for entropic regularization parameter λ = 0.01. The vagus and the pelvic samples are shown in cyan and orange, respectively, and labeled with the Image ID listed in Table 1. Figure 9H shows Image 13 (vagus), which is positioned far from the other samples in the Sinkhorn space. It is the largest sample in our dataset regarding image size and the number of segmented unmyelinated axons. It is also the only sample from a left cervical trunk. Figures 9A, C display Image 12 and Image 3 (both vagus) respectively. They are collected from the right cervical trunks. Figure 9I shows Image 7 (vagus) from abdominal vagus posterior trunk. Considering the spatial intensity, these four vagus samples are embedded far apart and are visually different. The rest of the samples are positioned in proximity, yet we can see a rightward tendency in the vagus samples than the pelvic ones. Figures 9D, E show Image 1 and Image 16 (both vagus), respectively. They look spatially different from the vagus samples discussed so far and are embedded at the leftmost part of the Sinkhorn space, far from those samples. However, they have a similar spatial organization as Image 28 and Image 22 (both pelvic) shown in Figures 9F, G and are embedded closer in the Sinkhorn space.

FIGURE 9

Figure 9. (A, C–I) A set of images of the segmented unmyelinated axons in the nerve cross-sections, labeled with the Image ID. (B) An embedding of the spatial intensity of the spatial point patterns in the Sinkhorn space for entropic regularization parameter λ = 0.01. The vagus and the pelvic samples are shown in cyan and orange [circles for female (F) and triangles for male (M)], respectively, and labeled with the Image ID listed in Table 1.

Although one can intuitively understand the global differences in spatial intensity of the point patterns by looking at the images of segmented unmyelinated axons in the nerve cross-sections, our approach can quantify and visualize the differences with the Sinkhorn distance between every pair of samples, resulting in a map of patterns. For instance, in Figure 9B, the Sinkhorn distances between the spatial intensity of Image 1, and Image 3 and Image 12 are 0.257 and 0.254, respectively, whereas the distance between Image 3 and Image 12 is 0.107. Again, Image 16 and Image 28 have respectively distances 0.054 and 0.048 from Image 1.

Interpreting spatial statistics, such as local inhomogeneous and anisotropic L-functions, can be more challenging than understanding raw spatial intensity. Figures 10C–E show three embeddings of the spatial features in the Sinkhorn space for λ = 0.01: the local inhomogeneous L-function, the local inhomogeneous L-function with horizontal sector and vertical sector, respectively. With a few exceptions, the overall landscape in the embeddings is similar to the one for the spatial intensity shown in Figure 9B. Figure 10A showing the segmented unmyelinated axons for Image 18 (vagus) contains several elongated axons. The elongated axons make the spatial arrangement of centroids in the corresponding point pattern (see Figure 10B) quite sparse and direction-oriented (anisotropic) in certain regions. The different positioning of Image 18 in the embeddings can reflect these characteristics.

FIGURE 10

Figure 10. (A, B) The segmented unmyelinated axons and the spatial point pattern of Image 18 (vagus) listed in Table 1. (C–E) The embeddings of the local inhomogeneous and anisotropic L-functions (no sector, horizontal sector, and vertical sector) of the spatial point patterns in the Sinkhorn space for entropic regularization parameter λ = 0.01. The vagus and the pelvic samples are shown in cyan and orange [circles for female (F) and triangles for male (M)], respectively, and labeled with the Image ID listed in Table 1.

3.5. Sinkhorn distance between maps of spatial features

Here, we compute the Sinkhorn distance between every pair of point patterns using the map of the spatial features constructed by kernel-smoothing. When we consider spatial point patterns directly (as in Section 3.4), the mass (corresponding to the spatial intensity or any other spatial feature attached as marks) to be transported is concentrated at the exact location of a point. As we apply kernel-smoothing to the point pattern, the concentrated mass at any point diffuses into its neighborhood. This step can help capture the notion of regionality in the kernel-smoothed maps while computing Sinkhorn distances. Notably, the kernel-smoothing reduces the influence of the differing number of points in the compared point patterns on the resulting Sinkhorn distances. The transportation-based metrics are well-suited for quantifying differences between bitmaps in which pixel values can be interpreted as transportable mass without strict geometric constraints (Rubner et al., 2000; Grauman and Darrell, 2004; Haker et al., 2004; Chefd'Hotel and Bousquet, 2007; Wang et al., 2011, 2013).

The kernel-smoothed spatial intensities, as well as marks attached to a point pattern, can be depicted as bitmaps, where the pixel values represent kernel-smoothed intensity values or other quantities derived from the marks (e.g., local inhomogeneous and anisotropic K- and L-functions). The kernel-smoothed spatial features of Images 3 (vagus) and 29 (pelvic) listed in Table 1 are shown in Figure 11. The bitmaps for the local inhomogeneous L-function, demonstrated in Figures 11B, F, have higher values of the spatial feature compared to their anisotropic counterparts shown in Figures 11C, D, G, H and slight shifts in values at certain locations are observed between the bitmaps of the horizontal and vertical sectors. Quantifying similarities between the kernel-smoothed bitmaps can be performed using Sinkhorn distance just like quantifying similarities between the spatial features of the original point patterns.

FIGURE 11

Figure 11. Visualizing the kernel-smoothed spatial features of Image 3 (vagus) and Image 29 (pelvic) listed in Table 1. (A, E) Spatial intensity. (B, F) Local inhomogeneous L-function. (C, G) Local inhomogeneous L-function with the horizontal sector. (D, H) Local inhomogeneous L-function with the vertical sector. The scale bars show the range of values for each spatial feature separately (column-wise). The kernel-smoothed bitmaps were downsampled for reasonable runtime and memory requirements.

Let I₁ and I₂ be the centered (0-padded as necessary) kernel-smoothed maps of the spatial intensity of the point patterns S₁ and S₂, respectively. The pixel values in I₁ and I₂ are normalized to sum to one, and the value at each pixel is considered the amount of mass contained (r) or required (c) at that pixel. The location of the pixels is not known beforehand, so we construct a unique grid [0, 1]² over which the pixel locations of I₁ and I₂ are defined. The cost matrix M is the Euclidean distance matrix computed from the [0, 1]² grid. The transportation plan P and the Sinkhorn distance between I₁ and I₂ are computed in the previously described manner. The Sinkhorn distances between the maps representing the other three spatial features are also calculated in the same fashion. Constructing the Sinkhorn distance matrix and visualizing embeddings in the Sinkhorn space are also done in the same way as described in Section 3.4.

Figure 12A shows an embedding of the kernel-smoothed maps of the spatial intensity of the point patterns in the Sinkhorn space of λ = 0.01. The vagus and the pelvic samples are shown in cyan and orange, respectively, and labeled with the Image ID listed in Table 1. The vagus samples 3, 7, 12, and 13 are embedded at a distance from the rest of the samples, and this trend was also observed in Figure 9B, when we processed the point patterns directly. However, vagus samples 6 and 15, and pelvic samples 19 and 26 (see Figures 12E–H), which were positioned close to the rest of the samples in Figure 9B, are located far apart in the right-most region of the embedding in Figure 12A. Therefore, some characteristics of the spatial intensities that were not captured during the processing of the raw point patterns became apparent when spatial feature maps were employed. The embeddings of kernel-smoothed local inhomogeneous and anisotropic L-function are illustrated in Figures 12B–D, portraying a similar trend overall, where the perceptually comparable vagus and pelvic samples are positioned in proximity, the rest are far apart, and the vagus samples are more spread out.

FIGURE 12

Figure 12. (A–D) The embeddings of the kernel-smoothed spatial intensity and the local inhomogeneous and anisotropic L-functions (no sector, horizontal sector, and vertical sector) of the spatial point patterns in the Sinkhorn space for entropic regularization parameter λ = 0.01. The vagus and the pelvic samples are shown in cyan and orange [circles for female (F) and triangles for male (M)], respectively, and labeled with the Image ID listed in Table 1. (E–H) The segmented unmyelinated axons of Image 15 and 6 (vagus) and Image 19 and 26 (pelvic) listed in Table 1.

3.6. Insights regarding the spatial architecture

We computed the Sinkhorn distances between every pair of nerve cross-sections for the four spatial features using spatial point patterns directly (data shown in Section 3.4) and kernel-smoothed bitmaps representation of the spatial features (data shown in Section 3.5). The resulting Sinkhorn embeddings are displayed in Figures 9, 10, 12. The created Sinkhorn space allows us to observe the similarities (or dissimilarities) of spatial intensities and second-order spatial properties.

The secondary statistical analysis performed on the embedded patterns generated by kernel-smoothing to mitigate the effects of the unequal number of axons revealed that the difference in spatial architecture between the vagus and pelvic nerves is relatively small (Mahalanobis distance Δ = 0.91). However, the sample size is insufficient to determine whether this observed difference reflects biological reality or results from random chance. With n_pelvic = 11 and n_vagus = 18, the achieved power (1-β) is only 0.6. In order to confirm the spatial architectural difference between vagus and pelvic nerve cross-sections, the required data set size should be at least n = 26 per class for 1-β = 0.8 and α = 0.05 in 2-D embedding, according to the collected preliminary results. In other words, any future research on the potential architectural difference between peripheral nerves' axonal organization (or modulation of a such organization due to pathology or treatment) must use these preliminary effect sizes as a reasonable basis for necessary power analysis needed for experimental design.

On the other hand, there is a substantial difference between the nerve cross-sections of males and females (Δ=1.246, Hotelling T² test p-value = 0.013). However, this effect must be confirmed and replicated with an unconfounded sample set in which the correlation between sex and cross-section origin (pelvic vs. vagus) is absent. The result reported here is based on the assumption that there is indeed no statistically significant difference between pelvic and vagus.

Regarding intraclass variability, vagus samples exhibit a significantly greater diversity of spatial architecture than pelvic samples when the raw spatial patterns are directly compared (Figure 9B). However, this difference disappears when the kernel-smoothed spatial patterns are compared, indicating that it was driven mainly by the difference in the number of axons rather than the spatial architecture (Figures 12, 13).

FIGURE 13

Figure 13. Boxplots displaying the ranges of Sinkhorn distance between the spatial intensity of the nerve cross-sections. Pelv: pelvic; Vag: vagus. (A) Analysis of the spatial point patterns directly. (B) Analysis of the kernel-smoothed maps of the spatial features. The points show the individual Sinkhorn distances, revealing the hidden distribution.

The vagus samples in our dataset are collected from the abdominal and cervical regions (see Table 1). Thus, looking into the degree of variability of the Sinkhorn distance within the sub-categories of the vagus samples and between the pelvic samples is helpful. Figure 14 illustrates the range of Sinkhorn distance between the spatial intensity of every pair of nerve cross-sections (for both types of analysis), categorized as the following:

1. within vagus samples

(a) intra-class measurements (i) (abdominal vagus vs. abdominal vagus)

(b) intra-class measurements (ii) (cervical vagus vs. cervical vagus)

(c) inter-class measurements (abdominal vagus vs. cervical vagus)

2. between pelvic and vagus samples

(a) inter-class measurements (i) (pelvic vs. abdominal vagus)

(b) inter-class measurements (iii) (pelvic vs. cervical vagus)

FIGURE 14

Figure 14. Boxplots displaying the ranges of Sinkhorn distance between the spatial intensity of the sub-categories of the nerve cross-sections. Abd: abdominal vagus; Cerv: cervical vagus, and pelvic. (A, B) Analysis of the spatial point patterns directly. (C, D) Analysis of the kernel-smoothed maps of the spatial features. The points show the individual Sinkhorn distances, revealing the hidden distribution.

In Figure 14, we see the ranges of the Sinkhorn distances for the abovementioned sub-categories. The degree of variability is more prominent in the analysis of the raw spatial point patterns (Figures 14A, B) compared to the analysis of the kernel-smoothed bitmaps representing spatial features (Figures 14C, D). This observation applies not just to spatial intensity but also statistical second-order spatial statistics. Figure 14 shows that variability within the abdominal vagus samples is substantially lower than the spatial variability within the cervical vagus cross-sections (p-value < 0.001). However, this notion has not been confirmed when looking at Figure 14C, which is based on pre-processing that eliminates the effect of axons' number. As before, this is most likely due to the low statistical power of the available sample size. The number of abdominal (55) and cervical (21) pairwise measurements is too small to confidently demonstrate the observed standardized effect size of Δ = 0.59. The required number of measurements for such effect size should be at least 45 per class to achieve 1 = β−0.8 with α = 0.05.

The much smaller standardized difference (Δ = 0.3) between two sites of vagus nerves sampling and pelvic nerves shown in Figure 14D can be confidently demonstrated due to the considerably larger number of available data points (121 pelvic vs. abd. vagus and 77 pelvic vs. cervical vagus pairwise measurements). Therefore, we can state that the dissimilarity between pelvic and abdominal vagus spatial architectures is much smaller than between pelvic and cervical vagus nerves (p-value = 0.0352). In other words, abdominal vagus samples resemble pelvic cross-sections to a higher degree than cervical vagus cross-sections.

4. Discussion

While many modern feature learning methods can directly classify biological images based on structural differences, the critical issue is the ability to quantify the specific architectural aspects of biological structures in order to relate them to a function or pathology. Neuroanatomy is one of the fields in which black-box image classifiers are undesirable, as the objective of the research is to link the image attributes to actual anatomical and physiological knowledge regarding cell and tissue organization, as opposed to simply sorting the images into predetermined categories. The approach presented here adds another module to our multi-step sample and data processing pipeline, which also includes the data acquisition and image segmentation modules described previously (Havton et al., 2021; Plebani et al., 2022).

We pursued the representation of the segmented unmyelinated axons in the TEM images of the peripheral nerve cross-sections as spatial point patterns not only to gain a better understanding of their neuroanatomy but also to express the observed differences in a quantitative manner, which would enable a variety of automated analysis tasks in the future, including automated image queries, image database retrieval, and biological image comparisons. While visual inspection of segmented images and their associated point patterns might provide some basic intuition regarding the spatial intensity, it is impossible to rely on the investigator's visual perception and judgments when examining more complex pattern characteristics such as local heterogeneous and anisotropic spatial features. Although global bulk measures of second-order spatial statistics, such as the K- and L-functions, help to represent and explore spatial interactions (randomness, inhibition, or clustering), they fail to capture local variations within biological structures. On the other hand, the local form of these spatial statistics functions generates yet another complex spatial pattern, leaving scientists with an equally tricky quantification problem. In this context, our analytical approach that captures differences between arrangements of any spatial distributions to form a visualizable embedding that enables straightforward comparison between complex structures provides a simple-to-use tool for neuroanatomists and computational neuroscientists.

There are at least three notable limitations associated with the demonstrated methods and their specific implementation. As previously stated, the claimed differences are relatively small and, despite being statistically significant, may be biologically unimportant. There is no reason to anticipate that the spatial arrangement would be dramatically altered in samples that do not represent a recognized disease. We realize that the value of the method would be more clearly demonstrated if the detected differences were associated with a specific biological mechanistic model, especially one associated with a disease or an abnormality. Although we lack such examples, we hope that researchers working on projects involving anatomical pathologies will be able to easily reproduce our methodology for quantifying observable differences in a biomedical context.

The second concern stems from the first: because there are no established alternative methods to quantify the spatial organization of axons, there is no way to validate the results by relating them to known physiology. This would indeed be a valuable exercise if the relevant nerves' physiology was already defined with sufficient precision. Unfortunately, this is not the case, so in the absence of sufficient data of this type, we have considered a number of well-established anatomical characteristics of these nerves that are consistent with our new study. For instance, our research aimed to distinguish between the cervical and abdominal vagus, and it revealed that the abdominal vagus and pelvic nerve share some similarities. This correlates with the higher prevalence of myelinated axons in the cervical vagus (Hulsebosch and Coggeshall, 1982; Prechtl and Powley, 1990), but it does not necessarily imply that the overall patterning of myelinated axons within each nerve type will be distinct. Physiological evidence of the type required to validate the current findings regarding sex differences in the vagus is also lacking. As many of the motor and sensory pathways supply sexually dimorphic targets, a sex difference may be anticipated for the pelvic nerve. However, in the present study, only samples from male rats were available.

The third concern relates to the orientation of biological structures. The anatomical rotational positions (orientations) of the analyzed fascicles were unknown (they were not recorded during sample processing), and we employed the post-hoc method described in Section 3.2 to identify the preferable orientation of the specimens. Although this method finds the rotations representing the smallest discrepancy between specimens, there is no guarantee that the identified orientations are biologically relevant. Sample alignment and orientation labeling is a broader problem in microscopy, not only affecting our analysis but also other techniques, such as multimodal imaging.

Despite these limitations, the presented method is an important contribution to the microscopy analytical toolbox. Notably, the availability of analytical tools is essential for the collection and evaluation of a large, comprehensive set of neuroanatomy and neuropathology data. As a result, the current scarcity of labeled and segmented images is attributable in part to the absence of an established analytic framework, casting doubt on the systematic value of acquiring comprehensive peripheral nerve data. We certainly hope that the conception and presentation of our method will inspire neuroanatomists to collect more data on peripheral nerves, resulting in broader quantitative anatomical studies. Importantly, our approach is simple to reproduce because it employs existing libraries for a popular statistical prototyping language.

Although we focused here on unmyelinated axons, the computational pipeline is applicable to multi-type point patterns and spatial research outside of neuroscience. This work demonstrates that the spatial architecture of unmyelinated axons in peripheral nerve cross-sections is neither uniform nor random but forms complex and rich arrangements. In order to simulate such a complicated spatial form, hybrid point processes are required. In the future, we plan to focus on spatial modeling and further classification of peripheral nerve cross-sections. The similarity (or dissimilarity) measure we established in this study will be the foundation for these modeling tasks.

5. Conclusions

In this report, we examined one of the key research problems in neuroanatomy, the fundamental description, measurement, and quantification of the spatial arrangement of axons in peripheral nerves such as the vagus and pelvic nerves (Hulsebosch and Coggeshall, 1982; Prechtl and Powley, 1990). This topic is significant not only from the basic neuroanatomical standpoint, but also due to the growing importance of peripheral nerve electrostimulation approaches, which rely on a precise understanding of the peripheral nerve architecture during the modeling and development phases (Pelot et al., 2020; Eiber et al., 2021). We believe that quantitative analysis, comparisons, and visualization of spatial arrangement can provide valuable insight to neuroanatomists, computational neuroscientists, and engineers working in the field of electrostimulation. We also believe that the presented method can be easily adapted to other biological fields, including spatial proteomics and genomics (Ji et al., 2020; Hickey et al., 2022).

Data availability statement

The microscopy data associated with this study were collected as part of the Stimulating Peripheral Activity to Relieve Conditions (SPARC) program and are available at the SPARC Portal (RRID: SCR_017041) under CC-BY 4.0 license (Havton et al., 2022).

Author contributions

BR conceived, planned, and supervised the spatial analysis study. AP contributed to the mathematical models. EP and MD preprocessed the microscopy data and designed the segmentation pipeline. TP and JK collected and processed the biological samples. TP provided neuroscience expertise and envisioned the quantitative peripheral nerve research. LH and NB collected and annotated the microscopy images. DJ curated the data and performed image preprocessing and mosaicing. AS and BR executed the study and co-wrote the manuscript with input from all the researchers. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Institutes of Health, Office of the Director, Stimulating Peripheral Activity to Relieve Conditions (SPARC) Program under Award Number OT2OD023847 (TP, BR, DJ, MD, EP, and AS), OT2OD023872 (JK), and OT2OD026585 (LH and NB); and by the U.S. Department of Energy's Advanced Scientific Computing Research program through grant DE-SC-0022260 (AP and AS).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Author disclaimer

The content of this work is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or the Department of Energy.

Footnotes

1. ^At https://github.com/Banus/umf_unet.

References

Ahuja, R. K., Magnanti, T. L., and Orlin, J. B. (1993). Network Flows - Theory, Algorithms and Applications. Prentice Hall.

Google Scholar

Asala, S. A., and Bower, A. J. (1986). An electron microscope study of vagus nerve composition in the ferret. Anat. Embryol. 175, 247–253. doi: 10.1007/BF00389602

PubMed Abstract | CrossRef Full Text | Google Scholar

Avis, D. (1980). On the extreme rays of the metric cone. Can. J. Math. 32, 126–144. doi: 10.4153/CJM-1980-010-0

CrossRef Full Text | Google Scholar

Baddeley, A., Rubak, E., and Turner, R. (2015). Spatial Point Patterns: Methodology and Applications With R. London: Chapman and Hall; CRC Press. doi: 10.1201/b19708

CrossRef Full Text | Google Scholar

Baddeley, A. J., Møller, J., and Waagepetersen, R. (2000). Non- and semi-parametric estimation of interaction in inhomogeneous point patterns. Stat. Neerland. 54, 329–350. doi: 10.1111/1467-9574.00144

CrossRef Full Text | Google Scholar

Basu, S., Kolouri, S., and Rohde, G. K. (2014). Detecting and visualizing cell phenotype differences from microscopy images using transport-based morphometry. Proc. Natl. Acad. Sci. U.S.A. 111, 3448–3453. doi: 10.1073/pnas.1319779111

PubMed Abstract | CrossRef Full Text | Google Scholar

Besag, J. (1977). Comment on “Modelling spatial patterns” by B.D. Ripley. J. R. Stat. Soc. Ser. B 39, 193–195. doi: 10.1111/j.2517-6161.1977.tb01615.x

CrossRef Full Text | Google Scholar

Bjaalie, J. G., Diggle, P. J., Nikundiwe, A., Karagulle, T., and Brodal, P. (1991). Spatial segregation between populations of ponto-cerebellar neurons: statistical analysis of multivariate spatial interactions. Anat. Rec. 231, 510–523. doi: 10.1002/ar.1092310413

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonaz, B., Sinniger, V., and Pellissier, S. (2017a). The Vagus nerve in the neuro-immune axis: implications in the pathology of the gastrointestinal tract. Front. Immunol. 8:1452. doi: 10.3389/fimmu.2017.01452

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonaz, B., Sinniger, V., and Pellissier, S. (2017b). Vagus nerve stimulation: a new promising therapeutic tool in inflammatory bowel disease. J. Intern. Med. 282, 46–63. doi: 10.1111/joim.12611

PubMed Abstract | CrossRef Full Text | Google Scholar

Breit, S., Kupferberg, A., Rogler, G., and Hasler, G. (2018). Vagus nerve as modulator of the brain-gut axis in psychiatric and inflammatory disorders. Front. Psychiatry 9:44. doi: 10.3389/fpsyt.2018.00044

PubMed Abstract | CrossRef Full Text | Google Scholar

Brickell, J., Dhillon, I. S., Sra, S., and Tropp, J. A. (2008). The metric nearness problem. SIAM J. Matrix Anal. Appl. 30, 375–396. doi: 10.1137/060653391

CrossRef Full Text | Google Scholar

Câmara, R., and Griessenauer, C. J. (2015). “Chapter 27: Anatomy of the vagus nerve,” in Nerves and Nerve Injuries, eds R. S. Tubbs, E. Rizk, M. M. Shoja, M. Loukas, N. Barbaro, and R. J. Spinner (San Diego, CA: Academic Press), 385–397.

Google Scholar

Chefd'hotel, C., and Bousquet, G. (2007). “Intensity-based image registration using Earth Mover's Distance,” in Medical Imaging 2007: Image Processing, Vol. 6512 (San Diego, CA: SPIE), 801–808. doi: 10.1117/12.709490

CrossRef Full Text | Google Scholar

Chiu, S. N., Stoyan, D., Kendall, W. S., and Mecke, J. (2013). “Point processes II – General theory,” in Stochastic Geometry and its Applications (Chichester; West Sussex: John Wiley & Sons, Ltd), 108–157. doi: 10.1002/9781118658222.ch04

CrossRef Full Text

Cuturi, M. (2013). “Sinkhorn distances: lightspeed computation of optimal transport,” in Advances in Neural Information Processing Systems, Vol. 26, eds C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Weinberger (Curran Associates, Inc.).

Google Scholar

Diggle, P. J., Lange, N., and Benes, F. M. (1991). Analysis of variance for replicated spatial point patterns in clinical neuroanatomy. J. Am. Stat. Assoc. 86, 618–625. doi: 10.1080/01621459.1991.10475087

CrossRef Full Text | Google Scholar

Dixon, P. M. (2014). “Ripley's K Function,” in Wiley Stats Ref: Statistics Reference Online (John Wiley & Sons, Ltd). doi: 10.1002/9781118445112.stat07751

CrossRef Full Text

Eiber, C. D., Payne, S. C., Biscola, N. P., Havton, L. A., Keast, J. R., Osborne, P. B., and Fallon, J. B. (2021). Computational modelling of nerve stimulation and recording with peripheral visceral neural interfaces. J. Neural Eng. 18:066020. doi: 10.1088/1741-2552/ac36e2

PubMed Abstract | CrossRef Full Text | Google Scholar

Grauman, K., and Darrell, T. (2004). “Fast contour matching using approximate Earth Mover's Distance,” in Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004, Vol. 1 (Washington, DC: IEEE), 1–8. doi: 10.1109/CVPR.2004.1315035

CrossRef Full Text | Google Scholar

Haker, S., Zhu, L., Tannenbaum, A., and Angenent, S. (2004). Optimal mass transport for registration and warping. Int. J. Comput. Vis. 60, 225–240. doi: 10.1023/B:VISI.0000036836.66311.97

CrossRef Full Text | Google Scholar

Havton, L. A., Biscola, N. P., Plebani, E., Rajwa, B., Shemonti, A., Jaffey, D., et al. (2022). High-Throughput Segmentation of Rat Unmyelinated Axons by Deep Learning (Version 1). SPARC Consortium. doi: 10.26275/K0MX-JCTH

CrossRef Full Text | Google Scholar

Havton, L. A., Biscola, N. P., Stern, E., Mihaylov, P. V., Kubal, C. A., Wo, J. M., et al. (2021). Human organ donor-derived Vagus nerve biopsies allow for well-preserved ultrastructure and high-resolution mapping of myelinated and unmyelinated fibers. Sci. Rep. 11:23831. doi: 10.1038/s41598-021-03248-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Hickey, J. W., Neumann, E. K., Radtke, A. J., Camarillo, J. M., Beuschel, R. T., Albanese, A., et al. (2022). Spatial mapping of protein composition and tissue organization: a primer for multiplexed antibody-based imaging. Nat. Methods 19, 284–295. doi: 10.1038/s41592-021-01316-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoffman, H. H. and Schnitzlein, H. N. (1961). The numbers of nerve fibers in the Vagus nerve of man. Anat. Rec. 139, 429–435. doi: 10.1002/ar.1091390312

PubMed Abstract | CrossRef Full Text | Google Scholar

Horn, C. C., Ardell, J. L., and Fisher, L. E. (2019). Electroceutical targeting of the autonomic nervous system. Physiology 34, 150–162. doi: 10.1152/physiol.00030.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Howland, R. H. (2014). Vagus nerve stimulation. Curr. Behav. Neurosci. Rep. 1, 64–73. doi: 10.1007/s40473-014-0010-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hulsebosch, C. E., and Coggeshall, R. E. (1982). An analysis of the axon populations in the nerves to the pelvic viscera in the rat. J. Comp. Neurol. 211, 1–10. doi: 10.1002/cne.902110102

PubMed Abstract | CrossRef Full Text | Google Scholar

Jafari-Mamaghani, M., Andersson, M., and Krieger, P. (2010). Spatial point pattern analysis of neurons using Ripley's K-function in 3D. Front. Neuroinformatics 4:9. doi: 10.3389/fninf.2010.00009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, A. L., Rubin, A. J., Thrane, K., Jiang, S., Reynolds, D. L., Meyers, R. M., et al. (2020). Multimodal analysis of composition and spatial architecture in human squamous cell carcinoma. Cell 182, 497–514.e22. doi: 10.1016/j.cell.2020.05.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Klatt, M. (2018). Barycenter: Regularized Wasserstein Distances and Barycenters. R Package Version 1.3.1. Available online at;: https://CRAN.R-project.org/package=Barycenter

Knight, P. A. (2008). The Sinkhorn–Knopp algorithm: convergence and applications. SIAM J. Matrix Anal. Appl. 30, 261–275. doi: 10.1137/060659624

CrossRef Full Text | Google Scholar

Krous, H. F., Jordan, J., Wen, J., and Farber, J. P. (1985). Developmental morphometry of the Vagus nerve in the opossum. Dev. Brain Res. 20, 155–159. doi: 10.1016/0165-3806(85)90100-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Møller, J., and Waagepetersen, R. P. (2003). Statistical Inference and Simulation for Spatial Point Processes. Boca Raton, FL: CRC Press. doi: 10.1201/9780203496930

CrossRef Full Text

Nadaraya, E. A. (1964). On estimating regression. Theory Probabil. Appl. 9, 141–142. doi: 10.1137/1109020

CrossRef Full Text | Google Scholar

Nadaraya, E. A. (1989). Nonparametric Estimation of Probability Densities and Regression Curves. Dordrecht: Springer Netherlands. doi: 10.1007/978-94-009-2583-0

CrossRef Full Text | Google Scholar

Ohser, J., and Stoyan, D. (1981). On the second-order and orientation analysis of planar stationary point processes. Biometr. J. 23, 523–533. doi: 10.1002/bimj.4710230602

CrossRef Full Text | Google Scholar

Orlin, J. B. (1993). A faster strongly polynomial minimum cost flow algorithm. Oper. Res. 41, 338–350. doi: 10.1287/opre.41.2.338

CrossRef Full Text | Google Scholar

Pele, O., and Werman, M. (2009). “Fast and robust Earth Mover's Distances,” in 2009 IEEE 12th International Conference on Computer Vision (Kyoto: IEEE), 460–467. doi: 10.1109/ICCV.2009.5459199

CrossRef Full Text | Google Scholar

Pelot, N. A., Goldhagen, G. B., Cariello, J. E., Musselman, E. D., Clissold, K. A., Ezzell, J. A., et al. (2020). Quantified morphology of the cervical and subdiaphragmatic vagus nerves of human, pig, and rat. Front. Neurosci. 14:601479. doi: 10.3389/fnins.2020.601479

PubMed Abstract | CrossRef Full Text | Google Scholar

Pereyra, P. M., Zhang, W., Schmidt, M., and Becker, L. E. (1992). Development of myelinated and unmyelinated fibers of human Vagus nerve during the first year of life. J. Neurol. Sci. 110, 107–113. doi: 10.1016/0022-510X(92)90016-E

PubMed Abstract | CrossRef Full Text | Google Scholar

Peyré, G., and Cuturi, M. (2019). Computational optimal transport: with applications to data science. Found. Trends Mach. Learn. 11, 355–607. doi: 10.1561/9781680835519

CrossRef Full Text | Google Scholar

Plebani, E., Biscola, N. P., Havton, L. A., Rajwa, B., Shemonti, A. S., Jaffey, D., et al. (2022). High-throughput segmentation of unmyelinated axons by deep learning. Sci. Rep. 12:1198. doi: 10.1038/s41598-022-04854-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Prechtl, J. C., and Powley, T. L. (1990). The fiber composition of the abdominal Vagus of the rat. Anat. Embryol. 181, 101–115. doi: 10.1007/BF00198950

PubMed Abstract | CrossRef Full Text | Google Scholar

Prodanov, D., Nagelkerke, N., and Marani, E. (2007). Spatial clustering analysis in neuroanatomy: Applications of different approaches to motor nerve fiber distribution. J. Neurosci. Methods 160, 93–108. doi: 10.1016/j.jneumeth.2006.08.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Ripley, B. D. (1976). The second-order analysis of stationary point processes. J. Appl. Probabil. 13, 255–266. doi: 10.2307/3212829

CrossRef Full Text | Google Scholar

Ripley, B. D. (1977). Modelling spatial patterns. J. R. Stat. Soc. Ser. B 39, 172–212.

Google Scholar

Ronneberger, O., Fischer, P., and Brox, T. (2015). “U-net: convolutional networks for biomedical image segmentation,” in Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015, Lecture Notes in Computer Science, eds N. Navab, J. Hornegger, W. M. Wells, and A. F. Frangi (Cham: Springer International Publishing), 234–241. doi: 10.1007/978-3-319-24574-4_28

CrossRef Full Text | Google Scholar

Rubner, Y., Tomasi, C., and Guibas, L. J. (2000). The earth mover's distance as a metric for image retrieval. Int. J. Comput. Vis. 40, 99–121. doi: 10.1023/A:1026543900054

CrossRef Full Text | Google Scholar

Safi, S., Ellrich, J., and Neuhuber, W. (2016). Myelinated axons in the auricular branch of the human vagus nerve. Anat. Rec. 299, 1184–1191. doi: 10.1002/ar.23391

PubMed Abstract | CrossRef Full Text | Google Scholar

Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., et al. (2012). Fiji: an open-source platform for biological-image analysis. Nat. Methods 9, 676–682. doi: 10.1038/nmeth.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Sengupta, P., Jovanovic-Talisman, T., and Lippincott-Schwartz, J. (2013). Quantifying spatial organization in point-localization superresolution images using pair correlation analysis. Nat. Protoc. 8, 345–354. doi: 10.1038/nprot.2013.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Sengupta, P., Jovanovic-Talisman, T., Skoko, D., Renz, M., Veatch, S. L., and Lippincott-Schwartz, J. (2011). Probing protein heterogeneity in the plasma membrane using PALM and pair correlation analysis. Nat. Methods 8, 969–975. doi: 10.1038/nmeth.1704

PubMed Abstract | CrossRef Full Text | Google Scholar

Settell, M. L., Skubal, A. C., Chen, R. C. H., Kasole, M., Knudsen, B. E., Nicolai, E. N., et al. (2021). In vivo visualization of pig vagus nerve “vagotopy” using ultrasound. Front. Neurosci. 15:676680. doi: 10.3389/fnins.2021.676680

PubMed Abstract | CrossRef Full Text | Google Scholar

Sinkhorn, R., and Knopp, P. (1967). Concerning nonnegative matrices and doubly stochastic matrices. Pac. J. Math. 21, 343–348. doi: 10.2140/pjm.1967.21.343

CrossRef Full Text | Google Scholar

Soltanpour, N., and Santer, R. M. (1996). Preservation of the cervical vagus nerve in aged rats: morphometric and enzyme histochemical evidence. J. Auton. Nerv. Syst. 60, 93–101. doi: 10.1016/0165-1838(96)00038-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Stoyan, D. (2006). “Fundamentals of point process statistics,” in Case Studies in Spatial Point Process Modeling, eds A. Baddeley, P. Gregori, J. Mateu, R. Stoica, and D. Stoyan (New York, NY: Springer New York), 3–22. doi: 10.1007/0-387-31144-0_1

CrossRef Full Text | Google Scholar

Thompson, N., Mastitskaya, S., and Holder, D. (2019). Avoiding off-target effects in electrical stimulation of the cervical vagus nerve: neuroanatomical tracing techniques to study fascicular anatomy of the vagus nerve. J. Neurosci. Methods 325:108325. doi: 10.1016/j.jneumeth.2019.108325

PubMed Abstract | CrossRef Full Text | Google Scholar

Thompson, N., Ravagli, E., Mastitskaya, S., Iacoviello, F., Stathopoulou, T.-R., Perkins, J., et al. (2023). Organotopic organization of the cervical vagus nerve. bioRxiv [Preprint]. doi: 10.1101/2022.02.24.481810

PubMed Abstract | CrossRef Full Text | Google Scholar

Veatch, S. L., Machta, B. B., Shelby, S. A., Chiang, E. N., Holowka, D. A., and Baird, B. A. (2012). Correlation functions quantify super-resolution images and estimate apparent clustering due to over-counting. PLoS ONE 7:e31457. doi: 10.1371/journal.pone.0031457

PubMed Abstract | CrossRef Full Text | Google Scholar

Villani, C. (2009). “The Wasserstein distances,” in Optimal Transport: Old and New, ed C. Villani (Berlin; Heidelberg: Springer), 93–111. doi: 10.1007/978-3-540-71050-9_6

CrossRef Full Text | Google Scholar

Waller, L. A., Särkkä, A., Olsbo, V., Myllymäki, M., Panoutsopoulou, I. G., Kennedy, W. R., et al. (2011). Second-order spatial analysis of epidermal nerve fibers. Stat. Med. 30, 2827–2841. doi: 10.1002/sim.4315

PubMed Abstract | CrossRef Full Text | Google Scholar

Walter, U., and Tsiberidou, P. (2019). Differential age-, gender-, and side-dependency of vagus, spinal accessory, and phrenic nerve calibers detected with precise ultrasonography measures. Muscle Nerve 59, 486–491. doi: 10.1002/mus.26412

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Ozolek, J. A., Slepčev, D., Lee, A. B., Chen, C., and Rohde, G. K. (2011). An optimal transportation approach for nuclear structure-based pathology. IEEE Trans. Med. Imaging 30, 621–631. doi: 10.1109/TMI.2010.2089693

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Slepčev, D., Basu, S., Ozolek, J. A., and Rohde, G. K. (2013). A linear optimal transportation framework for quantifying and visualizing variations in sets of images. Int. J. Comput. Vis. 101, 254–269. doi: 10.1007/s11263-012-0566-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Watson, G. S. (1964). Smooth regression analysis. Sankhyā 26, 359–372.

Google Scholar

You, K. (2022). T4transport: Tools for Computational Optimal Transport. R package version 0.1.1. Available online at: https://CRAN.R-project.org/package=T4transport