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.

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.

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.

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).

An estimator for the empirical K-function 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.

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 K^inhom(r) is formulated in Equation (3) below.

where ρ^(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.

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.

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.

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 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:

The estimator K^(r) is simply the average of all the K^(r,xi)s for i = 1, …, n. The centered local L-functions are formulated as L^(r,xi)=(K^(r,xi)/π)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.

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 ∑i,jPijMij. 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:

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).

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[d3log(d)] (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.

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).

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.

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.

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 1n1 amount of mass (r), therefore the total mass =1n1×n1=1. Similarly, each point in S₂ requires 1n2 amount of mass (c), so the total mass =1n2×n2=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).

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.

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.

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).

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:

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.