Generative modeling of the enteric nervous system employing point pattern analysis and graph construction

We describe a generative network model of the architecture of the enteric nervous system (ENS) in the colon employing data from images of human and mouse tissue samples obtained through confocal microscopy. Our models combine spatial point pattern analysis with graph generation to characterize the spatial and topological properties of the ganglia (clusters of neurons and glial cells), the inter-ganglionic connections, and the neuronal organization within the ganglia. We employ a hybrid hardcore-Strauss process for spatial patterns and a planar random graph generation for constructing the spatially embedded network. We show that our generative model may be helpful in both basic and translational studies, and it is sufficiently expressive to model the ENS architecture of individuals who vary in age and health status. Increased understanding of the ENS connectome will enable the use of neuromodulation strategies in treatment and clarify anatomic diagnostic criteria for people with bowel motility disorders.


Architecture and the connectome of the enteric nervous system
We describe a generative network model of the enteric nervous system (ENS) in the colon, developed with data obtained from confocal images of human and mouse tissue samples. Our models combine spatial point pattern analysis (SPP) with graph generation to characterize the spatial and topological properties of the ganglia (clusters of neurons and glial cells), the interganglionic connections, and the neuronal organization within the ganglia. We show that our approach, initially informed by the mouse colon ENS neuroanatomy (Figure 1), can simulate a variety of human samples differing in age and anatomical pathologies (Graham et al., 2020;Nestor-Kalinoski et al., 2022). We also show that the combination of SPP and graph generation approaches suffices to produce models exhibiting crucial properties of ENS anatomy, morphology, and circuitry. Despite the ongoing focus on mapping the central nervous system (CNS) connectome (Bassett and Sporns, 2017;Vecchio et al., 2017;Swanson and Bota, 2010), the peripheral nervous system (PNS) architecture in general, and enteric nervous system (ENS) in particular, have been largely disregarded. Besides the studies performed using simple biological models (Bezares-Calderón and Jékely, 2016), not much effort has been directed towards describing, quantifying, modeling, and understanding the topology and architecture of the PNS (Irimia and Van Horn, 2021). While the ENS has been described in terms of neuron and glial density or other fundamental characteristics, anatomic descriptors of ENS networks have not been well defined, so it is difficult to recognize when these ENS networks are disorganized (Graham et al., 2020).
Defining ENS network properties is valuable because diverse genetic, toxic, infectious, metabolic, nutritional, and inflammatory factors can disrupt ENS function to cause human diseases collectively called "bowel motility disorders" (Schneider et al., 2019;Furness, 2016). When the ENS is abnormal, bowel dysfunction can be debilitating or deadly since the ENS controls most aspects of bowel function in response to local stimuli. To monitor and regulate bowel function, the ENS is equipped with approximately as many nerve cells as the spinal cord, at least fourteen different neuron types, and nearly four times as many glial cells (Schneider et al., 2019;Graham et al., 2020;Wright et al., 2021). These ENS neurons and glia are distributed along the bowel in interconnected clusters called ganglia, embedded in two ganglionated plexi (intersecting networks): the myenteric plexus and the submucosal plexus. The primary function of the former is the control of muscle contraction and relaxation, while the role of the latter is to regulate epithelial function, blood flow, and the bowel immune system. ENS sensory neurons respond to luminal nutrients and stretch (Neunlist and Schemann, 2014). The ENS influences almost every bowel cell type and ENS cells respond to signals from many other cells (Schneider et al., 2019). These include close interactions of ENS with smooth muscle, interstitial cells of Cajal, PDGFRα cells, enteroendocrine cells, blood vessels, muscularis macrophages, and other bowel immune system cells.
The most dangerous ENS diseases are Hirschsprung's disease (HSCR) (Heuckeroth, 2018;Gosain and Brinkman, 2015;Langer, 2011), where ENS is absent from distal bowel and chronic intestinal pseudoobstruction (CIPO) (Di Nardo et al., 2017), where ENS is present but defective. HSCR and CIPO cause profound bowel dysmotility that can be life-threatening, requiring surgery or intravenous nutrition starting in infancy. Bowel dysmotility also occurs in Parkinson's disease, gastroparesis, achalasia, inflammatory bowel disease, and irritable bowel syndrome. These diseases are challenging to treat, and mechanisms are incompletely understood, limiting targeted therapy. A clear understanding of how these human diseases impact ENS anatomy or how anatomy might predict function would be valuable as we design new diagnostic and therapeutic strategies for people with bowel motility disorders.
Thus, our long-term goal is to describe and analyze the ENS connectome structure using graph-theoretic, network-based, and spatial-point pattern (SPP) process-related analytic approaches to develop robust mathematical models that distinguish normal from abnormal ENS anatomy. This is important because there is currently no definition of what morphology and connectivity patterning constitute normal or disorganized ENS networks. The work on the CNS connectome has demonstrated the value of the graph-analytic perspective for network quantification. Although graph-analytic approaches have been used in biology to model relationships between abstract entities (organisms, individuals, metabolites, regulatory systems, etc.), with a few simple exceptions, graph analysis has rarely been employed to characterize the spatial organization of anatomical structures. On the other hand, SPP techniques have been occasionally used in neuroscience since the seminal work of Ripley (1977). For instance, Diggle et al. (1994) utilized spatial analysis to investigate the distribution of pyramidal neurons in the cingulate cortex of normal subjects and schizophrenic patients.

Spatial point patterns and spatial networks
We briefly introduce the tools utilized in the process of building our model. The two critical conceptual components are the notions of spatial point patterns and spatial graphs or networks.
A spatial point pattern (SPP) X is a collection of spatial locations assigned to objects or events of interest in 2D or 3D space (Baddeley et al., 2015;Møller and Waagepetersen, 2003;Stoyan, 2006;Jafari Mamaghani et al., 2010). The points in an SPP (i) may be of different types (multitype point pattern), (ii) could be linked to auxiliary information or characteristics (marked point pattern), or (iii) could be linked to the space of interest (covariates). SPP analysis studies the points' spatial arrangement and aims to identify trends that characterize the patterns. It is common practice to compute summary statistics like nearest-neighbor distance, empty-space distance, or pair correlation to examine point patterns. A key component of an SPP is the underlying method or model for generating points. Every SPP can be thought of as an outcome of a spatial point process, a generative statistical model with predefined parameters describing a trend's formation. Two crucial parameters for a point process are intensity and interaction. The intensity (λ) of a point pattern is the average number of points per unit area, and it could be constant across space (stationary point pattern), or it might vary according to an intensity function (non-stationary point pattern). The interaction parameter describes the influence that the points have on each other. If the points are independent, the outcome is referred to as complete spatial randomness (CSR) and could be modeled with the Poisson point process. If the points exhibit positive interaction (spatial attraction), they can be modeled with various cluster and Cox processes. If they display negative interaction (spatial inhibition), they may be modeled with Gibbs processes (Baddeley et al., 2015).
A spatial network is a graph G(V, E) where the set of vertices V may be conceptualized as a spatial point pattern, and the set of edges E represents some context-dependent interaction between the points, constrained by geometry and linked to a distance (Barthelemy, 2018;Huang et al., 2014;Guimera et al., 2005;Marshall et al., 2018). The class of random networks developed by considering the underlying spatial constraints is sometimes referred to as spatially embedded random networks (SERNs) (Hackl and Adey, 2017;Barnett et al., 2007;Parsonage and Roughan, 2017). Spatial networks often share similar traits such as scale-free nature (Albert and Barabási, 2002) and small-world characteristics (Watts and Strogatz, 1998). Typically in these networks, closely located vertices are more likely to be connected than more spatially distant vertices. Vertices also tend to have low degrees. A specific class of spatial networks is a spatial planar network that can be drawn in the plane so that its edges do not intersect (Barthélemy, 2011).
A number of measures have been defined to describe spatial networks. A key measure is the clustering coefficient, which for a vertex i of degree k i is defined as: where E i is the number of edges among the neighbors of vertex i. The average value of the clustering coefficients of all the vertices is the clustering coefficient of the network. This metric is highly relevant for spatial networks with planarity, as this property increases the probability that closer vertices will be connected. Planarity results in a higher clustering coefficient compared to the values in a random Erdös-Rényi graph. Other relevant measures adopted from geographical networks are α (meshedness) (Buhl et al., 2006), g (density of network) (Kansky, 1963) and ψ (compactness) (Courtat et al., 2011), which for planar networks are defined as: where N is the number of vertices, E is the number of edges, A is the total area occupied by the spatially embedded network, and T is the total length of all the edges. Note that a planar graph can have at most 3N − 6 edges, which is the denominator of g; the numerator and denominator of α have the value (N − 1) subtracted from them (it is the number of edges in a tree with N vertices). All three measures have values in the interval [0, 1]. The value of meshedness (α) is 0 for trees and 1 for maximal planar graphs (planar graphs with the maximum number of edges). The density of the network (g) is the ratio of the number of edges present in the graph to the maximum number of edges that could be in the planar graph. Compactness (ψ) represents how much area is filled with edges (Barthelemy, 2018).

METHODOLOGY
In the simplest conceptual approximation of the ENS network, its structure resembles a planar spatial network with the positions of the ganglia generated from a spatial point process with negative interactions (see Figure 1). The ganglia centroids constitute the vertices of the network, and interganglionic connections are the network edges. Macroscopically, in terms of geometric embedding, the ENS consists of two major interconnected spatial networks (myenteric and submucosal plexi) embedded into concentric cylinders, with the second network internally comprised of highly linked subnetworks (Nagy and Goldstein, 2017;Furness, 2012). However, for the purpose of the modeling, we focus on the myenteric plexus only.

Biological data
For the development of our initial models, we used confocal images of mouse and human ENS in which the tissue was immunolabeled, allowing for easy visualization of different anatomical structures (see Figure 1A and Figure 2A). The interganglionic network in the mouse middle colon resembles the architecture of connections observed in humans. Figure 2A shows the flattened optical z-stacks (combination of multiple images taken at different focal distances) of the myenteric plexus collected from the proximal margin of a region of left colon removed from a five-month-old child with Hirschsprung disease (sample ID: 2475-P). The ENS structures were labeled using antibodies to the human neuron protein HuC/D (RRID: AB 221448) and to neuronal nitric oxide synthase nNOS (RRID: AB 91824) and then visualized using secondary antibodies Alexa Fluor 594 (purple enteric neuron cell bodies, RRID: AB 141633) and Alexa Fluor 488 (green inter-ganglionic connections, RRID: AB 253579). Higher magnification confocal images are in Figures 2B-E. The details of the methods used to capture the images of human ENS samples are described in previously published reports (Graham et al., 2020). The flattened z-stacks of the myenteric plexus were further processed using ImageJ (Fiji) software for manual segmentation and semi-automated extraction of the location of the ganglia centers and the connections of the networks (Schindelin et al., 2012). The semi-automated procedure involves denoising, Fourier filtering, and morphology-based segmentation. Subsequently, the pre-preprocessed network images were processed with Fiji Skeleton Analyzer, which provided the length of the interganglionic connections (network edges), location of the edges, distances between endpoints, and other helpful network features.

ENS ganglia
We use measures commonly employed in the exploratory analysis of spatial point patterns (Baddeley et al., 2015;Møller and Waagepetersen, 2003;Stoyan, 2006;Jafari Mamaghani et al., 2010) to characterize the location of ENS ganglia. A useful summary function, the Ripley K-function, in Besag's L-function format (shown in Equation 3) allows us to demonstrate the non-random spatial arrangement of ganglia by computing the cumulative distribution of the neighbors of the observed data points and comparing it with known alternative arrangements, such as complete spatial randomness (CSR).
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 observation window; e ij denotes weights for edge correction; and r is the distance. Positive values of the L-function depict spatial attraction, and negative values describe spatial inhibition. Figure 2H shows the L-function plot for the spatial point pattern of ganglia centers of the ENS network shown in Figure 2A, where the black curve falls visibly outside of CSR bounds demonstrating that the arrangement of ENS ganglia is non-random. The curve also displays significant spatial interaction (represented by one upward and one significant downward peak) up to a certain distance. Hence the organization of ENS ganglia shows a hybrid form of spatial inhibition occurring up to a certain distance. This behavior suggests that the arrangement of the ganglia might be modeled with Gibbs processes like a hybrid hardcore-Strauss defined by the intensity, interaction parameter, interaction distance, and hardcore distance (Baddeley et al., 2015).

ENS network
The spatial graph in Figure 2G represents the ENS network in Figure 2A with the ganglia centers as vertices and corresponding inter-ganglionic connections as edges. The network is planar, and the vertices have a low degree (with a mode of 3), resembling a small-world model. The network has low meshedness (α=0.19), moderate network density (g=0.46), and high compactness (ψ=0.995). The clustering coefficient of the network is 0.0097. The degree-one vertices are located at the sample boundary because of the dissection artifact. We also characterize the network in terms of the angle and the length of the edges. The distribution of the angle of the edges (see Figure 2I) is multimodal, with modes approximately at 10, 96, and 169 degrees, which gives the entire network a certain dominant orientation. Figure 2J shows the distribution of the length of the edges with an average length of 109 pixels (unit pixel = 1.5743 µm).

Generative model of the ENS network
The characteristics described above provide a mathematical description of the ENS network. Now we develop a well-parameterized generative model that provides realizations closely matching those observed characteristics. Our generative model incorporates a three-step process: (i) generation of starting positions for the ganglia (which we refer to as "ganglia centers" although they may not technically be geometrical centers of the objects), (ii) generation of the edges (inter-ganglionic connections), and (iii) generation of the neurons' centroids (within the body of the ganglia). It is important to emphasize that the sequence of generative steps in the model does not necessarily follow the biological sequence of events during ENS development and maturation. It was chosen for algorithmic convenience. Specifically, we envision a generative model in which the ganglia are created via some pairwise point interaction process (Strauss, 1975). We simulate ganglia centers with the hybrid hardcore-Strauss process (HSP) with intensity parameter β, interaction parameter γ, interaction radius R, and hardcore distance H. In this process, each point contributes a factor β to the probability density of the point pattern, and each pair of points closer than R units contributes a factor γ to the density. 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. The conditional intensity is defined as: The HSP point pattern is stationary when the intensity parameter (β) is a positive numerical constant. The intensity parameter (β) can be a function of certain features of the points (i.e., x-coordinates, y-coordinates) resulting in a non-stationary HSP point pattern. Because ganglia occupy physical space, the incorporated hardcore distance disallows overlap and unrealistically close neighborhoods between structures. Although close pairs of points are allowed above the distance H, they can still be penalized in the model.
For the edge generation process, we choose a combination of a deterministic and a random connection model following the methodology by Hackl and Adey (2017). Let us denote by X an instance of the HSP in R 2 , by C d (X) a deterministic connection model, and by C r (X, c) a random connection model with connection function c(.). We propose the deterministic connection model to be a Delaunay triangulation, which is the dual of Voronoi regions and is based on minimum spanning trees (Hackl and Adey, 2017;Marchette, 2005). We construct the initial template for the inter-ganglionic connections G(V, E), where V = X and E = C d (X) is the Delaunay triangulation on X. This step returns a maximal planar spatial graph (planar graph with most number of edges). Subsequently, we rewire (following the terminology of Watts and Strogatz (1998)) every edge (v i , v j ) ∈ E according to the random connection model C r (X, c). We pick C r (X, c) to be a rejection sampling model so that the vertices are connected with a certain probability (Parsonage and Roughan, 2017). Let G 0 (V 0 , E 0 ) be a real ENS network, and D f (.) the probability distribution of a set of one or more features f of a given network, such as edge length, edge orientation, the degree of vertices, etc. The purpose of this distribution is to compute the probability of the edge (v i , v j ) sampled from the edges of the Delaunay triangulation given the probability distribution of the real ENS network for the selected f . We denote f ij as the value of the selected feature set f for the edge (v i , v j ). The random connection model C r (X, c) keeps the edge (v i , v j ) in E for which the connection function c(v i , v j )=1 and removes the other edges to form the final realization G(V, E ). We formulate the connection function as follows: A generalized outline of the edge generation process is shown in Algorithm 1. After constructing the Delaunay triangulation, we assign a 0-1 normalized weight to each edge as follows: where w ij is the unnormalized weight of the edge (v i , v j ); k v is the degree of the vertex v; w is the set of all unnormalized edge weights; w ij is the 0-1 normalized weight of the edge (v i , v j ); and min(.) and max(.) are the functions returning the minimum value and the maximum value of a given set, respectively. The 0-1 normalized edge weights work during the sampling process as a probability to pick an edge whose end vertices (either or both) have high degrees. For the stopping condition of the sampling process (mentioned as no more edges can be rejected in Algorithm 1) we incorporate two criteria: (i) we stop if any one of the three network measures, meshedness, network density or compactness, for the spatial network under development become lower than the measures of the real ENS network, and (ii) we stop if no edges can be rejected for a significantly large number of iterations. To address the dissection artifact, we only allow vertices at a certain distance ( ) from the boundary of the point pattern to have a degree of one or two. To compare the spatial networks resulting from Algorithm 1 for different choices of the feature set, we use the earth mover's distance (EMD). The EMD, first introduced for image analysis by Rubner et al. (2000;, has an equivalent metric on probability distributions known as Mallows, or Kantorovich-Wasserstein distance (Mallows, 1972;Levina and Bickel, 2001).
Although ganglia centers can play the role of dimensionless vertices in the simplest model, we extend the approach further by simulating the neurons present in the ganglia via another use of spatial point processes. Therefore, in the simulation's last step, our model produces the locations of neurons' centers in the ganglia. Figures 2B-E give close views of the arrangement of the neurons (cell bodies) within the ganglia of a real ENS. The ENS network structure appears to correspond to the neurons' spatial distribution. Therefore, we generate the points simulating the neurons' centers using an inhomogeneous HSP process. The edges of the simulated inter-ganglionic network serve as a basis for an underlying intensity profile, which introduces the inhomogeneity. The objective is to create intensity profiles representing higher densities near ganglia locations, gradually decreasing along the edges of the network. To achieve this, we generate a 2D binary image (pixel matrix) of the network G(V, E ), and then execute a series of morphological image transformations (several dilations followed by Euclidean distance mapping using Fiji to construct a greyscale mask used as the probability profile for the inhomogeneous HSP (Leymarie and Levine, 1992). This image transformation produces an intensity profile that influences the spatial extent of the ganglia. With the profile generated, the inhomogeneity of the process is achieved by applying multiple random thinning operations using the generated thinning surface Jensen and Nielsen, 2001).

Simulation of the generative model
In our work, we use existing comprehensive tools for spatial point process statistics and graph modeling, including R libraries spatstat (Baddeley et al., 2015(Baddeley et al., , 2020 and igraph (Csardi and Nepusz, 2006). We fit the hybrid hardcore-Strauss process (HSP) to the spatial point patterns formed by the centers of the ENS ganglia in such a way that the interaction parameter confirms spatial inhibition (γ ≤ 1). We use the function spatstat::ppm which includes a computationally efficient technique, based on logistic regression, for fitting Gibbs point process models to SPP data (Baddeley et al., 2014;Diggle et al., 1994;Jensen and Künsch, 1994). Although we experimented with both stationary and non-stationary HSP models, we report the non-stationary HSP results only as they are more appropriate for the generation of ENS ganglia centers. Equation 7 shows the intensity function of the non-stationary HSP model.
where (u x , u y ) are the spatial coordinates of the location u, and θ 0 , θ x and θ y are the model parameters to be estimated. To illustrate the process, we compute the parameters for the ENS image shown in Figure 2F. We simulate the ENS network's realization by constructing the planar spatial network on top of the simulated ganglia centers. The vertices of this initial Delaunay triangulation have higher degrees than the real ENS networks (see Figure 3C). Then, we sample edges from the triangulation using Algorithm 1. We experimented with three constraints imposed on sampling: the limitation in angles of the edges, the restriction in the length of the edges, use of both of them simultaneously. The evaluation showed that the simultaneous use of angle and length leads to realizations that are most similar (in the EMD metric) to those observed in real networks. Figure 3D shows one example of the simulated ENS realization. Though there are a few high-degree vertices (with degrees 6 and 7), the vertices have an overall low degree (a mode of 3), similar to the real ENS network. All the degree-one vertices are located at the sample window boundary, recreating the sample dissection artifact. Figures 3E-F show two more simulated realizations of the ENS network built on the same simulated ganglia centers. The distribution of the angle and the length of the edges of the simulated spatial networks are shown in Figures 3G-H, respectively, along with the distributions of the real ENS network and the Delaunay triangulation-based network. The distribution of the edge angles and the edge lengths of the three simulated spatial networks are similar, and they are very close to the corresponding distributions of the real ENS network. The plots also illustrate that the starting planar network (Delaunay triangulation) is significantly different from the real one. As the rejection sampling process terminates, the angle and length distributions get closer to the real ENS network.

B C A
To generate realizations of the neurons' centers, we construct the greyscale intensity profile using the simulated ENS network and generate a stationary HSP point pattern over the simulated ENS network's window. The parameter values are tuned by grid search. For the purpose of practical demonstration, we used the intensity parameter β=0.003 points per unit squared, interaction parameter γ=0.78, interaction radius R = 0.0285 and hardcore distance H=0.01. Choosing different values for the HSP model parameters results in different spatial arrangements of the neurons, hence in different appearances of the ganglia. Finally, we introduce inhomogeneity by thinning the points with the  Table 2 for the pixel size of the sample ENS images). *γ=N/A and R=N/A indicates hardcore process with hardcore distance H. greyscale intensity profile. Figures 4A-C show a simulated realization of the ENS network, the generated intensity profile, and a complete ENS network with neurons, respectively. The correlation between the network structure and the arrangement of the simulated neurons is better visualized in Figure 5.  Table 1 shows the physiological information about the used set of human ENS samples. These samples are highly diverse in terms of age and health status. The table shows also the values of the selected parameters of our generative model: θ intensity coefficients, interaction parameter γ, interaction radius R, hardcore distance H (see Equation 4 and Equation 7) and average edge length of the network. The parameters H and R determine a disk-like neighborhood for each ganglia center. A ganglion mostly interacts (has inter-ganglionic connections) with other ganglia in its neighborhood, which means the average length of the inter-ganglionic connections should be close to R. This interpretation is reflected in the parameter values recorded in Table 1. We note some correspondence between the parameter values and the physiology. For instance, the interaction radius (R) is higher in adults than in the pediatric sample 2475-P. However, at this stage of research, it is difficult to discern whether the different ENS architectural parameters in different samples reflect distinct disease processes or age differences, or both. The generative model developed with sufficient ENS samples could give a range for interaction distance (R) and intensity parameters (β or θ) that are separable in healthy adults and healthy children or healthy and diseased individuals. Due to the limited amount of biological data and difference in sample size, we cannot yet conclude that the recorded parameter values are the best-fitted ones for all ENS samples belonging to the same age group or having similar pathology. However, we argue that our generative model is expressive enough to capture the diversities in such ENS samples.

CONCLUSIONS AND FUTURE WORK
We describe the initial development of a generative model of the enteric nervous system network. The model combines concepts from spatial point pattern analysis and graph analysis, two fields that have rarely been used together to solve practical problems from biological domains. In our model, we characterize the spatial location of ganglia and the neurons present in them with hybrid hardcore-Strauss processes (HSP). We describe the inter-ganglionic connections as spatial planar graphs and use measures such as the degrees of the vertices and angles and lengths of the edges to provide the spatial network characteristics. The model operates by generating realizations of ENS ganglia from the HSP fitted to actual biological samples. The edges are created by combining deterministic and random connection models. Finally, the model operation concludes by generating the location of ganglionic neurons using a second HSP. The computed features of the simulated spatial networks are similar to that of the real ENS networks. Our research demonstrated that the combined graph/SPP modeling technique presented, utilizing existing statistical approaches, offers a parametric space suitable for describing a wide variety of real ENS organizations. Our models were appropriate for ENS architectures present in samples originating from adults, children, newborns, and healthy and sick individuals. The ability to represent ENS with a generative model dependent on a relatively low number of ante-hoc explainable parameters (some related to the SPP aspect and others from the connection model) will allow us to create compact and easy-to-process ENS feature vectors. We envision using our pipeline as a tool to rapidly and robustly characterize ENS architecture to define normal parameters and to identify disorganized ENS such as may occur in people with Hirschsprung disease, chronic intestinal pseudo-obstruction, or other severe bowel motility disorders. The proposed generative modeling approach and feature extraction from the model could complement traditional techniques used in pathology, which are based on qualitative descriptions, basic neuron counting, and other simple image analysis methods. These latter techniques might be operator-dependent and are often challenging to reproduce (Swaminathan and Kapur, 2010). We see our modeling approach as a convenient alternative, which could be easily paired with computer-aided diagnostic (CAD) tools in downstream analysis.
We are also aware of the model's limitations. Firstly, the number of accessible clinical samples utilized in the model development and validation is relatively low. This is due to the difficulty in obtaining the biological samples and technological challenges associated with the pre-processing required for sophisticated laser-scanning microscopy imaging. Consequently, we face a tough choice of whether the model development should be focused on a narrow range of conditions and patients or should be aimed at a more generalized ENS description. The first route would allow a demonstration of applicability in computer-aided diagnostics, but it would also be very self-limiting. Secondly, our model does not capture the diversity of neurons and glial cells present in ENS in the current version. It focuses on the architectural rather than functional properties of the network. Therefore, any biological conclusions based on the ENS descriptions recovered by the model must be paired with cell-functional descriptions obtained using immunostaining and manual examination of the samples conducted by neuroscientists and pathologists. Naturally, we plan to expand our model to multiple cell types.
Despite these limitations, we believe that our approach opens the door to a significant breakthrough in the use of ENS imaging for diagnostics. By introducing a well-defined procedure for characterizing ENS networks to describe overall network properties, including the neuronal clusters of ganglia and the inter-ganglionic connections, we are moving ahead from observer-dependent scoring and qualitative descriptors of ENS complexity, which do not readily capture network architectural features defined in our models.