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

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¹¹1at https://github.com/Banus/umf_unet 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 Figure 2A and Figure 2D 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 analysis of SPP aims to study the spatial arrangement of the points in an SPP and establish trends to describe the point pattern. Two fundamental descriptive characteristics of an SPP are intensity and interaction. The intensity ($\rho$ or $\rho(u)$) of a point pattern, a first-order statistics, is the average number of points per unit area, and it could be uniform across the observation window (homogeneous), or it could vary according to an intensity function (inhomogeneous). A point pattern’s intensity is usually denoted by $\lambda$ in literature, but we use $\rho$ to avoid confusion with another notation related to the optimal transport problem. 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 $\mathbf{X}$ be a point pattern and $t(u,r,\mathbf{X})$ be the number of points in $\mathbf{X}$ which lie within distance $r$ of the location $u$. Assuming $\mathbf{X}$ is a homogeneous point pattern with intensity $\rho$, the number of points within distance $r$ of a specific point is represented by $\rho K(r)$ (Ripley, 1976).

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

where $\nu(.)$ 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_{\textrm{inhom}}(r)$) can be evaluated. The estimator for the inhomogeneous $K$-function ($\hat{K}_{\textrm{inhom}}(r)$) is formulated in Equation 3 below.

where $\hat{\rho}(u)$ is an estimator of the intensity function $\rho(u)$, obtained using a kernel-smoothed (described later in the paper) intensity estimator (Baddeley et al., 2015). Figure 3 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. It illustrates spatial inhibition for a small range of interaction distance at the beginning, then spatial attraction for both the samples. The spatial attraction, in other words, the clustering tendency is more apparent in the pelvic sample which agrees with the images shown in Figure 2.

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; Stoyan et al., 2013). Computing the cumulative distribution of the neighbors within a section of the disc of $r$ radius between two directional preferences $\theta_{1}$ and $\theta_{2}$, instead of the entire disc of $r$ radius, gives the sector $K$- and $L$-functions (Baddeley et al., 2015).

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 ($\hat{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 $\hat{K}(r)$ is simply the average of all the $\hat{K}(r,x_{i})$s for $i=1,\dots,n$. Arbitrarily, the centered local $L$-functions are formulated as : $\hat{L}(r,x_{i})=(\hat{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).

We intend to utilize the spatial statistics discussed above, specifically spatial intensity, and local inhomogeneous and anisotropic $L$-functions to describe the spatial arrangement of the point patterns constructed from the unmyelinated axons in the nerve cross-sections.

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\times d$ matrices with row and column summing to $r$ and $c$, respectively. Any matrix $P\in U(r,c)$ describes a transportation plan that transports the mass in $r$ to $c$. Given a $d\times 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:

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 (Villani, 2009; Avis, 1980; Brickell et al., 2008) satisfying the conditions shown in Equation 6.

An OT problem is a convex optimization problem that can be solved using various approaches (Ahuja et al., 1993; Orlin, 1993) and for a general cost matrix the computational cost scales as $\mathcal{O}(d^{3}log(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 3.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 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.

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 $\theta$ values ($\theta=$45°in our experiments). We compute the Sinkhorn distance for all possible values of $\theta$ 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, $\lambda$, as the entropic regularization parameter. As $\lambda\rightarrow 0$, Sinkhorn distance approaches the optimal transport distance (Wasserstein distance, provided that the cost is Euclidean distance). As $\lambda$ 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. We experimented with $\lambda$ = {0.01, 0.05, 0.1, 0.5, 1.0, 2.0, 5.0}, and report the the results for $\lambda$=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_{1}$ and $S_{2}$ be two spatial point patterns, with $n_{1}$ and $n_{2}$ number of points respectively. Considering an optimal transport problem between $S_{1}$ and $S_{2}$, we assume that each point in $S_{1}$ 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_{2}$ 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_{1}$ to $S_{2}$ (balanced). The inputs for the computations are the spatial locations of the points in $S_{1}$ and $S_{2}$. Therefore, we can compute the Euclidean distance matrix between them to be the cost matrix $M$, of dimension $n_{1}\times n_{2}$. Here, the cost matrix $M$ captures the spatial intensity of the point patterns. Further, we compute the transportation plan $P$, of dimension $n_{1}\times n_{2}$, using the formulation described in Section 3.3.1 and obtain the Sinkhorn distance, which provides a measure of similarity of the spatial intensity between $S_{1}$ and $S_{2}$.

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_{1}$ and $S_{2}$ 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_{1}$ and $S_{2}$. 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$\times$ 29, for each of the four spatial features. We use muti-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.

Figure 4B shows an embedding of the spatial intensity of the spatial point patterns (used directly) in the Sinkhorn space for entropic regularization parameter $\lambda$=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 4H 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. Figure 4A and Figure 4C display Image 12 and Image 3 (both vagus) respectively. They are collected from the right cervical trunks. Figure 4I 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. Figure 4(D, 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 Figure 4(F, G) and are embedded closer in the Sinkhorn space.

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 4B, 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. Figure 5(C-E) show three embeddings of the spatial features in the Sinkhorn space for $\lambda$=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 4B. Figure 5A, 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 5B) quite sparse and direction-oriented (anisotropic) in certain regions. The different positioning of Image 18 in the embeddings can reflect these characteristics.

Here, we compute the Sinkhorn distance between every pair of point pattern using the map of the spatial features constructed by kernel-smoothing. When we consider spatial point patterns directly to compute the Sinkhorn distance between nerve cross-sections (in Section 4.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 the points. As we apply kernel-smoothing to the point pattern, the concentrated mass at any point diffuses into its neighborhood. This phenomenon can help capture the essence of locality in the kernel-smoothed maps while computing Sinkhorn distances. In addition, kernel-smoothing can reduce the influence of the differing number of points in the point patterns on the Sinkhorn distance, by constructing maps of the same dimensions. Also, transportation-based metrics are well-suited for quantifying differences between bitmap images in which pixel values can be interpreted as transportable mass without strict geometric constraints (Rubner et al., 2000; Haker et al., 2004; Chefd’Hotel and Bousquet, 2007; Grauman and Darrell, 2004; Wang et al., 2011; wang_linear_2013).

The kernel-smoothed intensity and marks of a point pattern can be represented as bitmaps, where the pixel values represent kernel-smoothed intensity or other derived quantities from the marks, such as the 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 6. The bitmaps for the local inhomogeneous $L$-function, shown in Figure 6(B, F), have higher values of the spatial feature compared to their anisotropic counterparts shown in Figure 6(C, 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 just like quantifying similarities between the spatial features of the original point patterns, using Sinkhorn distance.

Let $I_{1}$ and $I_{2}$ be the centered ($0$-padded as necessary) kernel-smoothed maps (same dimension) of the spatial intensity of the point patterns $S_{1}$ and $S_{2}$, respectively. The pixel values in $I_{1}$ and $I_{2}$ 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]^{2}$ over which the pixel locations of $I_{1}$ and $I_{2}$ are defined. The cost matrix $M$ is the Euclidean distance matrix computed from the $[0,1]^{2}$ grid. The transportation plan $P$ and the Sinkhorn distance between $I_{1}$ and $I_{2}$ 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 4.4.

Figure 7A shows an embedding of the kernel-smoothed maps of the spatial intensity of the point patterns in the Sinkhorn space of $\lambda$=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 4B, when we experimented on point patterns directly. However, vagus sample 6 and 15, and pelvic samples 19 and 26 (see Figure 7(E-H)), that were positioned close to the rest of the samples in Figure 4B, are embedded far apart in the right-most region of the embedding in Figure 7A. Therefore, some behavior in the spatial intensity that were not noticed while dealing with point patterns, have surfaced when maps of the spatial features are used. The embeddings of kernel-smoothed local inhomogeneous and anisotropic $L$-function are illustrated in Figure 7(B-D), portraying a similar trend overall, where the visually and spatially comparable vagus and pelvic samples are positioned in proximity, the rest are far apart, and the vagus samples are more spread out.

We compute the Sinkhorn distances between every pair of nerve cross-sections for the four spatial features using spatial point patterns directly (in Section 4.4) and kernel-smoothed bitmaps representation of the spatial features (in Section 4.5). The resulting Sinkhorn embeddings are shown in Figure 4, Figure 5 and Figure 7. 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 $\Delta=0.91$). However, the sample size is insufficient to determine whether this observed difference reflects biological reality or was the result of random chance. With $n_{\textrm{pelvic}}=11$ and $n_{\textrm{vagus}}=18$, the achieved power ($1-\beta$) is only $0.6$. 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-\beta=0.8$ and $\alpha=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 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 ($\Delta=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 4B). 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 (Figure 7 and 8).

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 9 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. 1.

   within vagus samples

   1. (a)

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

   2. (b)

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

   3. (c)

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

2. 2.

   between pelvic and vagus samples

   1. (a)

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

   2. (b)

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

In Figure 9, 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 (Figure 9A and B) compared to the analysis of the kernel-smoothed bitmaps representing spatial features (Figure 9C and D). This observation applies not just to spatial intensity but also statistical second-order spatial statistics. Figure 9 shows that variability within the abdominal vagus samples is substantially lower than the spatial variability within the cervical vagus cross-sections (p $<0.001$). However, this notion has not been confirmed when looking at Figure 9C, 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 $\Delta=0.59$. The required number of measurements for such effect size should be at least $45$ per class to achieve $1-\beta=0.8$ with $\alpha=0.05$.

The much smaller standardized difference ($\Delta=0.3$) between two sites of vagus nerves sampling and pelvic nerves shown in Figure 9D 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 the dissimilarity between pelvic and cervical vagus nerves (p $=0.0352$). In other words, abdominal vagus samples resemble pelvic cross-sections to a higher degree than cervical vagus cross-sections.