^{1}

^{2}

^{2}

^{1}

^{3}

^{2}

^{2}

^{1}

^{1}

^{2}

^{3}

Multispectral three-dimensional (3D) imaging provides spatial information for biological structures that cannot be measured by traditional methods. This work presents a method of tracking 3D biological structures to quantify changes over time using graph theory. Cell-graphs were generated based on the pairwise distances, in 3D-Euclidean space, between nuclei during collagen I gel compaction. From these graphs quantitative features are extracted that measure both the global topography and the frequently occurring local structures of the “tissue constructs.” The feature trends can be controlled by manipulating compaction through cell density and are significant when compared to random graphs. This work presents a novel methodology to track a simple 3D biological event and quantitatively analyze the underlying structural change. Further application of this method will allow for the study of complex biological problems that require the quantification of temporal-spatial information in 3D and establish a new paradigm in understanding structure-function relationships.

Traditional reductionist biology has developed a wealth of knowledge concerning individual molecules and their function. This data set, however, has exceeded our capabilities for integration and analysis. Attempts to step outside these subcellular networks, into the equally complex world of three dimensional (3D) tissue development, are limited by our ability to deal with the spatial and temporal information that can be gathered in such studies [

Graph theory has emerged as a method to characterize the structure of large complex networks leading to a better understanding of the dynamic interactions that exist between their components. A graph is a mathematical structure that represents the relationships between members of a given set and is depicted as discrete points (nodes) connected by interactive links (edges) [

Graph theory has been applied to a variety of real networks from biology to social networking the internet and literature citations. Biological networks have common designs that are governed by simple and quantifiable organizing principles [

Natural and engineered tissues are a collection of cells arranged within a structural scaffold of extracellular matrix (ECM) proteins that provide biochemical and mechanical cues to direct the function of those cells. Tissue function is therefore dependent upon the spatiotemporal resolution of matrix proteins, cells, and signaling molecules. Understanding the interactions of this complex set of components over time requires novel techniques of extracting data from the intact tissue. In fact, the three dimensionality of tissue itself is critical for the maintenance of normal cellular function [

This paper presents a novel application of graph theory to extract and quantify the tissue structure/function relationship. We construct cell-graphs that place cell nuclei as the nodes of the graphs and link those nuclei as a function of their pairwise distance based on the assumption that cells that are in physical contact (i.e., membrane proximity) interact. Previous work demonstrates that graph theory-based analysis of the locations of cell nuclei provides important information concerning the (dys)functional states of the tissue [

Cryopreserved hMSCs (Lonza) were grown according to manufacturer’s instructions. hMSCs were cultured in Dulbecco’s Modification of Eagle’s Medium 1x (DMEM) supplemented with 10% fetal bovine serum (FBS) and fungizone/penicillin/streptomycin (FPS) (10,000 units/mL). Medium was changed every three days and cultures were incubated at _{2}. Cells were detached using trypsin-EDTA and passaged into fresh culture flasks upon reaching confluence. hMSCs were used between passages 6 and 8. In preparation for incorporation into 3D constructs, cells were washed with PBS, detached with trypsin-EDTA, collected, and counted. All reagents were purchased from Fisher Scientific unless otherwise noted.

Three-dimensional collagen I gels were prepared by mixing cells with the following reagents: DMEM (14%), FBS (10%), 5^{6} cells per mL ECM and 3.0 ^{5} cells per mL, for their respective experimental sets. The constructs were incubated at

Collagen I constructs were prepared for fluorescence staining at hours 0, 1, 2, 4, 6, 10, 16, and 24. Constructs were washed with PBS followed by a 30-minute incubation at

The 3D confocal images were segmented using the Otsu Thresholding algorithm [

A graph is given by

For each pair of vertices the pairwise Euclidean distance is calculated and a link is set when that distance is smaller than the linking (edge) threshold. This threshold value is determined by considering the nucleus-membrane ratio and estimating cell diameter. We have identified thresholds, corresponding to the approximate radii of spread mammalian cells, of 65, 70, and 75 microns (from a distortion free sphere representation of the cell membrane to observable distortion) and executed our algorithms for each.

Topographical properties for each cell graph were extracted and averaged over the entire graph to represent global measures of “tissue organization”; see Table

Description of graph metrics.

Number of nodes | The number of vertices in the graph, in this case cell nuclei as defined by image segmentation. |

Number of edges | The number of links between the vertices of the graph, in this case the probability of cellular interaction based on distance. |

Average degree | The number of neighbors/connections that a given node has averaged over all the nodes in the graph. |

Clustering coefficient (C) | Is defined as |

Path length or hop count: | The shortest distance between two nodes reflective of the weight of each link. |

Number of nodes | The number of vertices in the graph, in this case cell nuclei as defined by image segmentation. |

Eccentricity and closeness | The maximum and average of the shortest lengths, respectively. The average of these two features provides global measure of average eccentricity and average closeness. |

Diameter | Defined as the maximum eccentricity, greatest distance between two nodes in the graph. |

Central points | Defined as nodes that have an eccentricity equal to the average eccentricity. |

Hop plot value | Reflects the size of a neighborhood between any two nodes with a “hop”, a single edge. |

Hop plot exponent | Slope of the hop plot values as a function of |

Effective diameter | Defined as |

Giant connected component | Defined as the largest set of nodes that are reachable from each other. |

Giant connected ratio | The ratio of the size of the giant connected component to the number of nodes in the graph. |

Isolated node | A node with no edges and therefore a degree of 0. |

End node | A node with one edge. |

The Cumulative Nearest Neighborhood Distribution (CNND) function [

This function value represents the fraction of the cell counts for which their nearest neighbor falls within the threshold distance

The randomness of cell-cell relationships was tested by defining an edge function according to the Erdos-Renyi model, calculating the same set of features on our random graphs, and comparing the random graph metrics to that of the cell-graphs:

We compared the cell-graph technique to Voronoi diagrams and the Delaunay triangulation method [

Classical graph mining [

For the cell graphs we sought subgraphs that are frequent with respect to their embedding count in one large graph. Support of a subgraph,

For each subgraph we detected the frequency of its occurrence at each time point and for each cell density using edge-disjoint embeddings. Subgraph pattern size was restricted to an upper limit of four edges, and each vertex is labeled by an integer number corresponding to the proportional degree of that vertex to capture local network information. The labeling scheme is shown in Table

Labeling scheme to capture local degree of subgraph networks.

Degree range | Label |
---|---|

0 | 0 |

1-2 | 1 |

3-4 | 2 |

5-6 | 3 |

7-8 | 4 |

9-10 | 5 |

13-14 | 7 |

15-16 | 8 |

17–19 | 9 |

20–22 | 10 |

23–26 | 11 |

More than 26 | 12 |

To provide proof of concept for our cell-graph method we encapsulated hMSC in 3D collagen I hydrogels at a concentration of

Macro- and micro scale compaction occurs over 24 hours of 3.0 ^{5} cell/mL gels. hMSCs at a cell density of 3.0 ^{5} cells/mL were embedded in 2 mg/mL collagen I gels. Collagen I gels were fixed and stained for nuclei with SYTOX Green (a) Macroscale images of hydrogel compaction (Scale bar

Cell-graphs, Figure

Graph theory applied to fluorescent confocal z-stack images to generate distance dependent cell-graphs for 3.0 ^{5} cell/mL gels. (a) 3D representation of fluorescent image z-stack. (b) 3D cell-graph generated from fluorescent image z-stack.

The utility of these cell-graphs comes from the ability to extract meaningful and quantitative metrics to describe the relationships between nodes. The metrics (Table

Quantitative metrics for 3.0 ^{5} cells/mL collagen I hydrogel 24-hour timecourse.

Time (hours) | |||||||

0 | 1 | 2 | 6 | 10 | 16 | 24 | |

Average degree | 1.18 | 1.53 | 2.45 | 3.62 | 5.56 | 9.40 | 14.40 |

Clustering Coefficient | 0.26 | 0.37 | 0.47 | 0.56 | 0.63 | 0.59 | 0.61 |

Average eccentricity | 1.11 | 2.08 | 2.16 | 4.27 | 13.52 | 14.38 | 14.24 |

Maximum eccentricity (diameter) | 5.00 | 8.00 | 7.00 | 9.00 | 23.00 | 20.00 | 19.00 |

Minimum eccentricity (radius) | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | 10.00 | 10.00 |

Average eccentricity 90 | 1.11 | 1.89 | 2.16 | 3.77 | 11.28 | 10.97 | 11.28 |

Maximum eccentricity 90 | 5.00 | 7.00 | 7.00 | 8.00 | 20.00 | 16.00 | 17.00 |

Minimum eccentricity 90 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | 8.00 | 8.00 |

Average path length (closeness) | 0.84 | 1.31 | 1.41 | 2.39 | 6.97 | 7.07 | 7.36 |

Hop plot exponent | 0.20 | 0.49 | 0.41 | 0.70 | 1.08 | 1.51 | 1.46 |

Effective hop diameter | 14410546.75 | 961.44 | 3191.99 | 164.88 | 26.67 | 10.34 | 11.29 |

Giant connected ratio | 0.13 | 0.22 | 0.20 | 0.20 | 0.68 | 0.99 | 1.00 |

Number of connected components | 35.00 | 34.00 | 35.00 | 26.00 | 12.00 | 3.00 | 1.00 |

Percentage of isolated points | 0.38 | 0.32 | 0.22 | 0.05 | 0.02 | 0.01 | 0.00 |

Percentage of end points | 0.28 | 0.19 | 0.17 | 0.12 | 0.05 | 0.02 | 0.01 |

Number of central points | 23.00 | 23.00 | 21.00 | 9.00 | 11.00 | 1.00 | 14.00 |

Percentage of central points | 0.38 | 0.32 | 0.22 | 0.05 | 0.05 | 0.00 | 0.03 |

Quantitative metric statistics for different link thresholds for 3.0 ^{5} cells/mL collagen I hydrogel 24-hour timecourse.

Time (hours) | |||||||

0 | 1 | 2 | 6 | 10 | 16 | 24 | |

Average degree | |||||||

Clustering coefficient | |||||||

Average eccentricity | |||||||

Maximum eccentricity (diameter) | |||||||

Minimum eccentricity (radius) | |||||||

Average eccentricity 90 | |||||||

Maximum eccentricity 90 | |||||||

Minimum eccentricity 90 | |||||||

Average path length (closeness) | |||||||

Hop plot exponent | |||||||

Effective hop diameter | |||||||

Giant connected ratio | |||||||

Number of Connected Components | |||||||

Percentage of isolated points | |||||||

Percentage of end points | |||||||

Number of central points | |||||||

Percentage of central points | |||||||

Number of nodes | |||||||

Number of edges |

To demonstrate the metrics’ sensitivity to perturbations in the biological system we repeated the experiment using a higher cell concentration (

Changes in cell density alter compaction and nuclear organization. hMSC at a cell density of 1 ^{6} cells/mL embedded in 2 mg/mL collagen I gel. Collagen I gels were fixed and then stained for nuclei with SYTOX Green. (a) Macroscale images of hydrogel compaction (Scale bar

Cell-graphs, Figure

Quantitative metrics for 1 ^{6} cells/mL collagen I hydrogel 24-hour timecourse.

Time (hours) | |||||||

0 | 1 | 2 | 6 | 10 | 16 | 24 | |

Average degree | 12.47 | 6.30 | 7.05 | 11.47 | 14.70 | 22.68 | 50.07 |

Clustering Coefficient | 0.71 | 0.64 | 0.67 | 0.67 | 0.65 | 0.62 | 0.66 |

Average eccentricity | 15.93 | 10.66 | 12.90 | 15.59 | 13.41 | 11.68 | 11.29 |

Maximum eccentricity (diameter) | 21.67 | 18.33 | 21.00 | 21.67 | 17.33 | 15.00 | 14.67 |

Minimum eccentricity (radius) | 1.00 | 0.67 | 1.00 | 8.00 | 6.33 | 8.00 | 7.67 |

Average eccentricity 90 | 12.28 | 8.56 | 11.08 | 12.51 | 10.82 | 9.35 | 9.02 |

Maximum eccentricity 90 | 19.33 | 16.33 | 18.67 | 18.67 | 14.67 | 12.67 | 12.33 |

Minimum eccentricity 90 | 1.00 | 0.67 | 1.00 | 6.33 | 5.00 | 7.00 | 6.00 |

Average path length (closeness) | 7.84 | 5.23 | 6.65 | 7.58 | 7.00 | 6.11 | 5.85 |

Hop plot exponent | 1.33 | 1.01 | 1.07 | 1.24 | 1.39 | 1.54 | 1.52 |

Effective hop diameter | 12.47 | 35.15 | 22.77 | 16.26 | 10.99 | 9.44 | 9.12 |

Giant connected ratio | 0.98 | 0.65 | 0.73 | 0.88 | 0.99 | 1.00 | 1.00 |

Number of connected components | 4.00 | 14.67 | 11.33 | 3.33 | 3.00 | 1.00 | 1.00 |

Percentage of isolated points | 0.00 | 0.02 | 0.02 | 0.00 | 0.01 | 0.00 | 0.00 |

Percentage of end points | 0.01 | 0.09 | 0.04 | 0.01 | 0.01 | 0.00 | 0.00 |

Number of central points | 4.67 | 10.33 | 10.33 | 5.33 | 9.67 | 13.33 | 34.67 |

Percentage of central points | 0.01 | 0.06 | 0.05 | 0.02 | 0.02 | 0.02 | 0.02 |

Number of nodes | 380.33 | 190.00 | 218.33 | 350.33 | 440.00 | 743.33 | 1458.67 |

Number of edges | 2378.67 | 604.33 | 769.00 | 2049.33 | 3339.00 | 8454.33 | 36770.33 |

Graph theory applied to fluorescent confocal z-stacks from the 1 ^{6} cells/mL time course to generate distance dependent cell-graphs. (a) 3D representations of fluorescent image z-stacks of 1 ^{6} cell densities. (b) 3D cell-graph generated from fluorescent image z-stacks of 1 ^{6} cell density.

The sensitivity of our metrics to gross changes in biological function is a prerequisite to the extension of this method into more complex problems. Quantification of gel compaction revealed a relationship between collagen I compaction and cell density; see Figure

Collagen I hydrogels compact as a function of cellular density. The diameter of the collagen I gels at each time point was measured using Image J and plotted as a percentage of the original volume. Changes in cellular density resulted in significant changes in compaction over the observed 24 hours time course, where an increase in cell density correlated with an increase in gel compaction.

Figure

Cell-graph mining extracts metrics which represent the topological changes in 3D cell culture as a function of cell density. The effect of change in cell density from 3.0 ^{5} to 1.0 ^{6} on our graph features was plotted. (a) Average degree. (b) Clustering coefficient. (c) Diameter. (d) Number of edges. (e) Giant connected ratio. (f) Number of connected components over the 24-hour time course. The extracted metrics have been normalized from

To show that tissue has a structure which is different from the random organization of the same number of cells, we established two properties: (i) cells are not randomly scattered in 3D space, and (ii) cells interact in a nonrandom fashion. For the first property, we randomly distributed points representing the cells within 3D space with the same density and geometry seen in the images. The CNND function was implemented explained for both random scattering and for actual cell locations obtained from image segmentation. The results, Figure

Point Pattern Analysis indicates that cells do not pick random locations. In order to show that cell locations obtained from segmented images are different from random scattering, we implemented CNND function

Cell-graph features are distinct from random graphs. Random graphs were generated keeping the number of nodes and links constant. Features were extracted from the random graphs and compared to features from cell-graphs to determine significance. (a) Clustering Coefficient. (b) Giant Connected Component. (c) Average Eccentricity.

To show that cell-cell interactions cannot be modeled by randomly picking pairwise interactions between cells, we defined an edge function according to Erdos-Renyi (ER) model and generated random graphs to distinguish their features from those extracted from our cell-graphs. Consider three cases for generating random graphs: (i) cell locations are chosen randomly and the ER edge function is used, (ii) cell locations, obtained from image segmentation, are taken as node locations and ER edge function is used, and (iii) node locations are randomly chosen but the cell-graph edge function is used. Since case (ii) is subsumed in (i) (recall that ER model takes number of nodes

For case (i) random graphs were generated with equal numbers of nodes and edges as found in the cell-graphs at each given time point. We extracted the same features calculated for our cell-graphs from these random graphs. Some of the metrics are defined by the number of nodes and links present in the tissue (i.e., average degree) and therefore must be consistent between the cell and random graphs. The remaining values in our cell-graphs were significantly different from those found in the random graphs, indicating that the cell-graph metrics and random metrics did not come from the same distribution; see Figure ^{5 }cell/mL samples fluctuated between 0.02 and 0.03, but gradually increased from 0.35 to 0.62 for the cell-graphs over the 24-hour time period. This difference captures the compaction of the tissue samples as the value gradually increases. Also, the giant connected component of the high density random graphs maintained a value of 1 throughout the 24 hour period, whereas the values for the cell-graphs steadily increased to 1, matching the compaction of the tissue samples. The average eccentricity increased throughout compaction in cell-graphs whereas in random graphs this metric was almost constant.

For case (iii) we also randomly distributed cells in a 3D space with the same density and geometry of the corresponding samples and used the cell-graph edge function to construct the cell-graphs. We compared the graph metrics on these random-cell-graphs with graphs from the 3D image segmentation. The cell-graphs defined a relationship (edge) between the nearest random nodes and the graph features correlated with many of the tissue cell-graph features. However, certain graph metrics such as the graph radius were different. In Figure

Kolmogorov-Smirnov test for graph radius. The plot (although the sample size is small) shows the difference between a cell graph with randomly scattered nodes and one with segmented nodes.

In order to determine the relevance of our approach we assessed the performance of our cell graph method to the more traditional Delaunay triangulation model. The tissue was partitioned into Voronoi cells and the corresponding Delaunay graph was generated for each set of images. Graph features were extracted from these Delaunay graphs as shown in Table

Quantitative delaunay graph metrics for 1 ^{6} cells/mL collagen I hydrogel 24-hour timecourse.

Time (hours) | |||||||

0 | 1 | 2 | 6 | 10 | 16 | 24 | |

Average degree | 13.55 | 13.24 | 13.17 | 13.9 | 13.93 | 14.3 | 14.33 |

Clustering Coefficient | 0.48 | 0.49 | 0.49 | 0.47 | 0.47 | 0.45 | 0.45 |

Average eccentricity | 4.7 | 4.4 | 4.45 | 5.23 | 5.44 | 5.92 | 5.85 |

Maximum eccentricity (diameter) | 5.8 | 5.4 | 5.2 | 6.4 | 6.8 | 7 | 7 |

Minimum eccentricity (radius) | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

Average eccentricity 90 | 4.04 | 3.85 | 3.88 | 4.45 | 4.71 | 5.03 | 5.01 |

Maximum eccentricity 90 | 5 | 4.6 | 4.6 | 5.4 | 5.8 | 6 | 6 |

Minimum eccentricity 90 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

Average path length (closeness) | 3.01 | 2.85 | 2.89 | 3.38 | 3.58 | 3.9 | 3.93 |

Hop plot exponent | 2.03 | 1.83 | 1.91 | 2.16 | 2.26 | 2.46 | 2.54 |

Effective hop diameter | 3.94 | 4.02 | 3.91 | 4.43 | 4.71 | 4.84 | 4.83 |

Giant connected ratio | 1 | 1 | 1 | 1 | 1 | 1 | 1 |

Number of connected components | 2 | 2 | 2 | 2 | 2 | 2 | 2 |

Percentage of isolated points | 0 | 0.01 | 0.01 | 0 | 0 | 0 | 0 |

Percentage of end points | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

Number of central points | 1 | 1 | 1 | 1 | 1 | 1 | 1 |

Percentage of central points | 0 | 0.01 | 0.01 | 0 | 0 | 0 | 0 |

Number of nodes | 239.8 | 204 | 205 | 396 | 493.6 | 748.2 | 849.6 |

Number of edges | 1628.6 | 1375.6 | 1371.2 | 2770 | 3438.4 | 5352.2 | 6100.8 |

The features extracted from the tissue cell-graphs were representative of the global (average) structure of our “tissues.” To extract meaningful local structures and motifs we used subgraph mining techniques that identify patterns of node interaction and measure their frequency of occurrence over time. Thousands of subgraph patterns were generated that represent repeated local motifs. The frequency of common subgraphs was also measured in random graphs, to eliminate the possibility that the local structures extracted were simply structures present in any random organization of nodes and links. Distinguishable patterns were extracted that existed only in the cell-graphs, Figure

Common cell subgraph patterns are distinct from patterns extracted from random graphs. Subgraph patterns were identified that represented local graph topologies and compared to their frequency in random graphs. (a) Significant cell-graph patterns not found in random graphs. (b) Patterns found equally in random and cell-graphs and are therefore not significant. (c) Histogram of all patterns common among cell-graphs plotted as the absolute difference between frequencies in random versus cell-graphs. Patterns at the extremes indicate significant patterns.

The appearance of specific subgraphs was observed at distinct time points and cellular densities. To reveal the time evolution of the significant subgraphs we calculated the number of patterns that first appeared in the 3.0 ^{5} cell/mL graphs at each time point as a function of their frequency at a given time point in the 1.0 ^{6} cell/mL graphs; see Table ^{5} cell/mL time series correlated with patterns frequently found at early time points in the 1.0 ^{6} cell/mL series. Some time points showed a very strong correlation, with a high number of similar patterns, specifically time 10 hours in the 1.0 ^{6} series and time 16 hours in the 3.0 ^{5} series (bold value in Table ^{5} subgraphs that appeared before, after and during respective 1.0 ^{6} subgraphs.

Time evolution of subgraphs.

time (hrs) | 0 | 1 | 2 | 6 | 10 | 16 | 24 | |

0 | 7 | 1 | 0 | 0 | 0 | 0 | 0 | |

1 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | |

2 | 9 | 0 | 0 | 18 | 3 | 0 | 0 | |

6 | 26 | 0 | 10 | 5 | 42 | 0 | 0 | |

10 | 122 | 3 | 0 | 19 | 76 | 1 | 36 | |

16 | 143 | 2 | 9 | 46 | 1303 | 1274 | ||

24 | 8 | 0 | 0 | 0 | 457 | 842 | 515 |

The extraction of nonrandom repeating motifs in biological structures presents a significant opportunity to define organizing design principles and use those principles to inform and accelerate natural and synthetic processes of tissue development. Here we present the spatial and temporal analysis of a simple biological observation, the compaction of collagen I hydrogels by hMSC. The progression of this event over time was dependent upon the cellular density of the starting material and represents a direct, cell-mediated organizational event that can be observed through time and space. We observed that an increase in cell density caused an accelerated compaction event that resulted in a decrease in the final volume as compared to the original samples. We tracked this organization event using confocal microscopy and employed graph theory as a method to extract the 3D topography of our “tissue constructs” over the time course of the compaction event. As the hMSCs interact with and organize the naïve collagen matrix, we hypothesized that structural motifs would appear frequently and be evident at earlier time points in the dense samples. This predicts that cells interact with one another and their extracellular environment in nonrandom, defined patterns that we can quantify, and that these patterns are consistent in terms of their temporal organization.

To address this hypothesis we constructed cell-graphs from our 3D micrographs and extracted topographical features that assess the global 3D structure of the “tissue constructs” as they organize over time. Such techniques have previously been reported to extract meaningful information concerning the (dys)functional state of 2D tissue and can play a significant role in diagnosis of cancer pathologies [

The use of random graphs as a negative control was critical to ensure the extraction of meaningful and nonrandom metrics. The divergence of features extracted from our cell-graphs from those found in random graphs indicates that there exists structure within our “tissue” that is different from a random organization of the same number of cells and links. This is consistent with the theory that random networks do not accurately describe real networks, that networks are highly ordered structures, and that keys to their function lie within that inherent structure [^{6} cell/mL the graph is connected, and as the density increases, the diameter decreases as more links are added to the graph. Within our cell-graphs, however, at this same cell density, we see a steadily increasing diameter. This increase in diameter is due to the fact that the connectivity at time zero is sufficiently low to observe increases as the gels condense over time, unlike randomly generated graphs, and therefore the diameter of the entire graph responds by growing as the giant connected component size grows as well. The random graph is connected as a function of the number of nodes and links incorporated into its generation, but there is no true structural relationship evident, as there is no change observed in eccentricity or clustering coefficient. Instead what we see in our cell-graphs is a steadily emerging structural relationship between nodes that is quantified by changes in eccentricity and the clustering of the nodes.

We also used subgraph mining to extract a finer picture of the local structure present within our cell-graphs. Through a multilevel analysis of the cell-graphs we can project a quantitative, sensitive, and rigorous fine structure on the global features extracted above. Searching the micro level space of the graph we find local structures that are not represented equally over the whole of the graph. In this way, functional components of the tissue, or modules, are identified. For example, a sample in which only a portion represents cancerous tissue might be averaged out and classified as inflamed or normal. However, through the use of subgraph mining, and local feature analysis, an area of densely connected cells would be identified as a specific tissue compartment with unique properties and lead to reclassification as cancer [

Over time these subgraphs became increasingly complex and incorporated nodes with higher degrees and therefore more interaction with neighboring cells, consistent with the conclusion from the global features. In addition, the frequent subgraphs appeared at specific points in time, and not all subgraphs were present at all time points for all densities examined in this study. We hypothesize that the formation of distinct local structures occurs in a specific temporal order. The time evolution of our subgraphs indicates that there is a “route” that a population of cells takes while interacting with and organizing a naïve matrix, and if we change a fundamental property of that population, then we change the rate at which the “tissue” is able to complete that “route.” This hypothesis was confirmed by quantifying the frequency of subgraph occurrence as a function of time. When compared across the two densities, we observed that the appearance of new subgraphs in the less dense sample at late time points consistently correlated with high frequencies of those same patterns at earlier time points in the dense sample. We conclude that a population of cells follows a consistent organizing pattern over time, and while the rate of that organization can be altered by changes in the cell population, the route or order in which that organization takes place remains constant.

Point pattern analysis and its variants [

Therefore, there are several advantages of our cell-graph approach to tissue modeling: First, cell-graphs enable us to benefit from well established principles of graph theory and provide a rich set of features defined precisely by these principles to be used as quantitative descriptor features. These features could be defined and computed locally from a single node’s point of view (e.g., number of its neighbors) or globally for the entire tissue sample (i.e., the shortest or longest distance in the cell-graph between any two nodes). Second, cell-graphs benefit from sophisticated segmentation algorithms that can provide cell level attributes such as convexity, size, physical contact, and shape. to establish links between a pair of nodes. Furthermore, molecular details of a particular cell type (e.g., is it diseased?) do not need to be resolved, and a specific textural change in the image is not required. Third, cell-graphs provide a precise mathematical representation of cellular organization and the ECM that surrounds cells. If the images carry multichannel information by applying more sophisticated staining techniques (e.g., multispectral fluorescence imaging), it is possible to build cell-graphs that have different types of nodes corresponding to different types of cells that coexist (e.g., epithelial versus fibroblast) and other ECM entities (e.g., basement membrane underlying epithelial cell layers and blood vessels). With 3D images and 3D cell-graphs, such representation becomes more accurate and powerful.

The demonstration of higher order in biological structures is evidence of recurring principles that govern the design and therefore development of metabolic networks, protein-interactions, neural circuitry, and tissue. Cells interact with one another and their extracellular environment in non-random, defined patterns, and these patterns are consistent in terms of their spatial and temporal organization. Understanding the temporal order of this patterning provides organizing principles that can be used to help biologists and tissue engineers understand or optimize processes of natural and synthetic tissue development.

This work was supported in part by the National Institutes of Health, Grant no. RO1 EB008016 and no. RO1 AR053231. The first and the second authors contributed equally to this work.