cytoNet: Spatiotemporal network analysis of cell communities.

We introduce cytoNet, a cloud-based tool to characterize cell populations from microscopy images. cytoNet quantifies spatial topology and functional relationships in cell communities using principles of network science. Capturing multicellular dynamics through graph features, cytoNet also evaluates the effect of cell-cell interactions on individual cell phenotypes. We demonstrate cytoNet's capabilities in four case studies: 1) characterizing the temporal dynamics of neural progenitor cell communities during neural differentiation, 2) identifying communities of pain-sensing neurons in vivo, 3) capturing the effect of cell community on endothelial cell morphology, and 4) investigating the effect of laminin α4 on perivascular niches in adipose tissue. The analytical framework introduced here can be used to study the dynamics of complex cell communities in a quantitative manner, leading to a deeper understanding of environmental effects on cellular behavior. The versatile, cloud-based format of cytoNet makes the image analysis framework accessible to researchers across domains.

Introduction

Cell-cell interactions affect the behavior of individual cells through intricate paracrine and juxtracrine signaling. To observe these effects, many research groups have developed in situ profiling techniques to extract highly multiplexed single-cell data while preserving the spatial characteristics of biological samples [1-6]. The increasing prevalence of spatially detailed imaging datasets has led to the proliferation of spatial analysis pipelines for biological research (). While these methods have enabled principled exploration of spatial hypotheses, the majority of the pipelines (with a few exceptions) have been developed for spatial molecular expression data obtained through methods such as mass cytometry, specialized high-resolution imaging, and/or scRNA-seq, with inherent idiosyncrasies. Others have focused on histology and samples obtained for medical applications. As a result, these techniques are not applicable to many standard imaging datasets obtained through routine biological experiments. Further, many pipelines require the user to be familiar with programming and involve the use of customized scripts. All of these limitations mean the most advanced spatial analysis platforms are not commonly employed by biologists. Instead, the spatial analysis platforms are largely used by a subset of labs heavily invested in computational analysis, by core users of specialized microscopy, or by imaging experts themselves. There remains a need for a generalizable, easy-to-use analysis method to test spatial hypotheses applicable to a wide variety of biological imaging data.

To address these limitations, we introduce cytoNet, a user-friendly method to analyze spatial and functional cell community structures from microscope images, using the formalism of network science (). cytoNet enables the time-dependent study of structure-function relationships in cell biology. While many tools automate the analysis of functional read-outs like calcium signaling (), cytoNet’s strengths lies in its ability to integrate spatial and functional topology to provide insight into the nature of structure-function coupling in biological systems. cytoNet’s use of graph theory also enables a single modeling framework to represent multiple descriptors of cell community and provides continuity across spatial and temporal scales. Network models provide a flexible method where objects are nodes and relationships between them are modeled as edges [7,8]. Among existing image-based analyses that employ network science, the cell-graph technique [9] has been employed to great effect in analyzing structure-function relationships in fixed tissue sections. Our early work applying network analysis to fixed samples also enabled rapid classification of cell phenotypes [10,11]. cytoNet advances from these platforms to enable the study of cell community structure and dynamics in multiple local microenvironments.

cytoNet workflow.

(a) The cytoNet pipeline begins with masks and optionally microscope images, which can be static immunofluorescence images or calcium image sequences. (b) Spatial proximity is determined either by measuring shared pixels between cell pairs–type I networks, or by comparing the distance between cell centroids to a threshold distance–type II networks (right panel). Functional networks are estimated from correlations in calcium time series data. (c) Cell community descriptors provide information on local neighborhood characteristics of individual cells, like degree and centrality measures, and global neighborhood characteristics like modularity and path lengths.

cytoNet is available as a web-based interface run on Amazon cloud. Users can choose to analyze image files from their desktops or online servers. Coupled with its ease-of-use, cytoNet’s versatility makes it accessible to researchers across domains. We originally designed the network modeling approach to study populations of developing neurons [12] and characterize how vascular cells respond to neurotrophic factors [10,11]. Here we extend the approach to case studies in a number of other biological systems, and provide cytoNet as a standalone tool for the community. We partnered with labs from across research domains to illustrate applications of the cytoNet platform to stem cell biology, tissue engineering, and neuroscience in both in vitro and in vivo settings. The case studies demonstrate the broad utility of the network modeling approach in studying spatial and functional community structure in complex biological systems.

Results

The cytoNet pipeline enabled us to investigate spatial and functional topology of cell communities in a variety of biological systems. Four case studies are described in the sections below.

dynamic analysis, coupled multicellular spatial and functional topology

immunostain fluorescence images and time-lapse calcium imaging videos (frames in TIF format)

human neural stem cells, human neurons

We designed an in vitro model of neural differentiation to analyze the dynamics of spatial and functional topology during formation of neural circuits from neural progenitor cells (NPCs) [12]. NPCs are known to display structured intercellular communication prior to formation of synapses, which plays an important role in controlling self-renewal and differentiation [23-25]. By leveraging the cytoNet method, we sought to capture the dynamic structure of NPC communities and the effect of such community structure on the phenotypes of individual cells.

In this case study, we describe data obtained using ReNCell VM human neural progenitor cells, in which spontaneous differentiation was triggered through withdrawal of growth factors, leading to rapid cell cycle exit and formation of dense neuronal networks in 5 days [12]. We captured spontaneous calcium activity at days 1, 3, and 5 after withdrawal of growth factors. Following calcium imaging, cells were fixed, and nuclei were stained and reimaged. Nuclei images were then manually aligned by fiducial markers with their corresponding calcium images. The paired image sets allowed the creation of both functional and spatial graphs for the same communities of cells.

Spatial type II graphs () showed a fall and subsequent rise in network heterogeneity (tendency of the network to contain hub nodes) ( and The low network heterogeneity at day 3 likely reflects the homogeneous distribution of cells at this time point, which gives way to clustered cell bodies with hub-like network structure at day 5. We hypothesize that these trends in spatial properties, independently confirmed in multiple NPC lines [12], reflect a transition from topologies favoring global to hierarchical information flow. We further explored this possibility through calcium imaging. Functional networks constructed from spontaneous calcium activity () revealed network-wide signal correlations, with trends in spontaneous network activity mirroring spatial network parameters (). These results suggest that spatial topology predicts functional communication patterns in differentiating NPCs, with high spatial network homogeneity at intermediate time points facilitating network-wide communication and hub-like network structures at early and late time points mirroring more clustered communication.

Dynamics of spatial and functional topology in developing neural progenitor cells (NPCs).

(a) Spatial NPC networks at day 1, 3 and 5 of differentiation, overlaid on immunofluorescence images of nuclei stained with Hoescht dye; segmented cells are outlined in red, and spatial proximity edges are shown as yellow lines; scale bar = 50 μm. (b) Network heterogeneity of spatial NPC networks is lowest at day 3; red notches show mean and standard deviation; *p < 0.005 from two-sample t-test. (c) Functional networks obtained through calcium imaging with Fluo-4 in developing NPC networks at days 1, 3 and 5. Correlations between calcium traces from individual cells are shown as a network plot overlaid on the maximum intensity image from calcium image sequences. (d) Fraction of active cells in the network; *p < 0.005 from two-sample t-test. Active cells are defined as cells whose normalized fluorescence traces have three or more calcium transients. (e) Frames from time-lapse movies of differentiating NPCs transfected with FUCCI cell cycle reporters. Borders of mCherry+ nuclei (G1) are outlined in magenta, Venus+ nuclei (S/G2/M) are outlined in green, and mCherry-/Venus-nuclei (quiescent) are outlined in blue, spatial edges are overlaid in yellow; scale bar = 50μm. (f) Neighborhood similarity score for low-density culture across time. (g) Neighborhood similarity score across time for medium-density culture. Fig 2A–2D adapted from reference [12].

We next studied the role of cell-cell communication on cell cycle regulation of NPCs. Cell cycle regulation in NPCs is of interest as it has implications for the genetic basis of brain size in different species [26] and aberrant regulation can cause diseases like microcephaly [27]. Studies in the ventricular zone of the embryonic mouse neocortex have shown that clusters of clonally-related neural progenitor cells go through the cell cycle together [28,29]. However, it is unclear whether this community effect is a ubiquitous feature of neural progenitor cells. To this end, we employed the cytoNet workflow to determine whether cell cycle synchronization is a feature of differentiating NPCs cultured in vitro.

For this part of the investigation, ReNCell VM human neural progenitor cells were stably transfected with the FUCCI cell cycle reporters [30] to generate Geminin-Venus/Cdt1-mCherry/H2B-Cerulean (FUCCI-ReN) cells. We captured time-lapse movies of FUCCI-ReN cells after withdrawing growth factors to induce differentiation and built network representations from nucleus images. Adjacency was determined by comparing centroid-centroid distance to a threshold (type II graphs).

In order to evaluate spatiotemporal synchronization in cell cycle, for each individual cell in a frame, we evaluated the average fraction of neighboring cells in a similar phase of the cell cycle (G1 phase–mCherry+ and S/G2/M phases–Venus+), normalized by total fraction of that cell type in the population. We called the average value of this fraction across all cells in an image the neighborhood similarity score, N. Frames from time-lapse movies for low- and medium-density cultures are shown in (see also ). We observed that groups of cells in the low-density culture moved through the cell cycle in unison, which was reflected in periodically high values of the neighborhood similarity score (). In contrast, the composition of cell clusters in the medium density culture was relatively heterogeneous, resulting in relatively low values of the neighborhood similarity score over time (). Neighboring cells in very low-density cultures are likely to be derived from the same clonal lineage, which explains the high level of synchronization in these cultures [19]. This example highlights how cytoNet can be used to test hypotheses on the role of cell-cell interactions in regulating dynamic cell behavior.

dynamic analysis, coupled multicellular spatial and functional topology

time-lapse microscopy videos (frames in TIF format)

mouse brain

In vivo calcium analysis is an avenue for exploring and understanding the role that individual cells of the nervous system play in processing external stimuli including pain. Pain is mainly mediated by a subset of primary sensory neurons known as nociceptors in Dorsal Root Ganglia (DRG) and Trigeminal Ganglia (TG) [46,47]. How DRG neurons function at a population level under physiological and pathological conditions is unknown. Imaging methods developed to record from hundreds to thousands of neurons simultaneously in the brains of live mice are helping elucidate this [47,48]. To investigate population characteristics of pain-sensing neurons, we used cytoNet to evaluate spatial and functional networks from calcium image sequences obtained in a Pirt-GCaMP3 mouse DRG model.

Calcium image sequences, along with single masks identifying individual cells, were inputs to cytoNet (see Methods section “Generation of functional networks,” for details on generation of masks) (). Only cells whose Ca2+ activity changed over the course of the experiment contributed to the functional network analysis displayed. Sensory stimulation experiments produced a single, major signal spike in each segmented cell [47]. Measurement of the magnitude (ΔF/F0) of each spike is sensitive to the quality of segmentation; to mitigate this, we characterized each cell not by its spike magnitude, but by the time a cell took to reach its peak value from 20% of that value (ramp-up) and the time needed for the signal to return to 20% (ramp-down). Inspection of 44 segmented cells revealed 6 unique combinations of ramp-up times and ramp-down times (). Ramp-up times were either 5 or 10 seconds while ramp-down times varied between 5 and 35 seconds. This categorization of cells according to functional similarity was combined with the spatial graph of the segmented cells in order to identify spatial patterns of cells with similar behavior (). In addition, we note that although the vast majority of segmented cells reached their peak intensity at 20 seconds, a small group of cells along the left side of the tissue peaked at 25 seconds suggesting a right-to-left wave of response (). This case study highlights the utility of cytoNet in analyzing spatial patterns of neural populations with unique functional signatures in an in vivo model.

Dynamics of Coupled Functional & Spatial Analysis In Vivo.

cytoNet captures relationships between spatial proximity of neurons and functional features of multicellular modules in vivo. (a) Cells classified according to the time required to first reach their maximum ΔF/F0 values from 20% of that value (ramp-up) and the time required to return to 20% (ramp-down). Edges connect similarly classified cells that are within 10 cell diameters of each other. All cells reached their peak values at 20 seconds except for those circled which reached their peak values at 25 seconds. (b) Calcium time series (ΔF/F0) plotted for 6 categories of cells with unique combinations of ramp-up and ramp-down times. The blue braces indicate a cell’s ramp-up and ramp-down. Each inset image is a spatial pattern of cells with the same ramp-up and ramp-down times. RU = ramp-up time; RD = ramp-down.

time-point analysis, coupled single cell morphology and multicellular spatial topology, compensation for community effects on drug / stimuli

immunostain fluorescence images (TIF format)

human endothelial cells

In a second application to studying human cells in vitro, we used cytoNet to evaluate the relative influence of local neighborhood density and growth factor perturbations on endothelial cell morphology. From a regenerative medicine perspective, studying the morphological response of endothelial cells to neurotrophic stimuli can help assess the cells’ potential angiogenic response following brain injuries that induce the secretion of neurotrophic factors, like ischemic stroke or transient hypoxia [49,50]. Common high-throughput angiogenic assays focus on migration and proliferation as the main cell processes defining angiogenesis, or the growth of new capillaries from existing ones [51]. Distinct morphology and cytoskeletal organization of endothelial cells indicate the cell’s migratory or proliferative nature, and hence their angiogenic contribution within a sprouting capillary [52]. Reproducibly quantifying the morphological response of endothelial cells to neurotrophic factors would enable more targeted approaches to enhancing brain angiogenesis.

We took an image-based approach to this problem, building a library of immunofluorescence images of human umbilical vein endothelial cells (HUVECs) stained for cytoskeletal structural proteins (actin, α-tubulin) and nuclei, in response to various combinations of vascular endothelial growth factor (VEGF) and brain-derived neurotrophic factor (BDNF) treatment. Cell morphology was annotated using 21 metrics described in our previous study [53] (), which included cell shape metrics like circularity and elongation, and texture metrics for cytoskeletal stains such as actin polarity, smoothness etc. Network representations were designated based on shared cell pixels (type I graphs) and local network properties were described using the metrics in

First, we quantified density-dependent effects on endothelial cell morphology in control cultures (without any growth factor perturbation). Our analysis showed correlations between cell morphological features and local network properties (). Some of these relationships were expected, for instance the positive correlation between shared cell border and cell size. Other relationships, such as the negative correlation between cell circularity and closeness centrality, capture intuitive notions of the influence of cell packing on morphology (). The closeness centrality of a cell () describes its relative position in a colony–cells in the middle of a colony will have higher centrality values than cells at the edge of a colony or isolated cells. The negative relationship between circularity and closeness centrality implies that isolated cells and cells located at the edge of colonies are more likely to have a circular morphology, while cells located at the center of colonies tend to be less circular (). Thus, our analysis revealed that local network properties have a quantifiable effect on cell morphology.

Influence of local neighborhood density on primary human endothelial cell (HUVEC) morphology.

(a) Distribution of cell circularity values grouped under different levels of closeness centrality; sample size, n = 786 cells (group 1; cn < 0.025), 741 cells (group 2; 0.025 < cn < 0.05) and 782 cells (group 3; cn > 0.05); Cohen’s d effect size: groups (1, 2) = 0.34, groups (1, 3) = 0.62 (b) Sample immunofluorescence image with graph representation overlaid; scale bar = 50 μm. (c) Heatmap depicting closeness centrality of each cell, with circularity values overlaid in text. (d) Representative cells from cluster analysis, highlighted in magenta. (e) Cell size, closeness centrality and circularity distribution plots for each cluster. (f) Bar plot of variance explained by growth factor treatment and local network metrics. (g) Box plot of cell size as a function of growth factor treatment. (h) Box plot of mean actin intensity as a function of growth factor treatment. Legends and axes in (f-h) contain information on treatment (BDNF, VEGF), concentration (50ng/ml, 100ng/ml) and time of treatment (6 hours and 12 hours). Cohen’s d effect size for (f-h) is shown in .

To determine dominant cell phenotypes, we performed cluster analysis on our dataset consisting of 25,068 cells. This analysis revealed 3 major categories of endothelial cells, with unique morphological and network signatures (). Cluster 1 comprised cells with migratory features, including low circularity and intermediate centrality indicative of their position at the edges of colonies. Cluster 2 contained small, circular cells with low centrality indicative of their isolation. Cells in cluster 3 showed proliferative features with large non-circular shapes, and high centrality indicating their positions in the center of colonies. Through this cluster-based phenotyping, we show how cytoNet can be used to infer the local environment and topological arrangement of distinct cell categories within a culture.

Next, we developed a workflow to analyze the effect of growth factor treatments on cell morphology, while correcting for the effect of local network properties. We did this to infer the independent effects of chemical perturbation and local cell crowding on cell morphology. First, we applied a quantile multidimensional binning approach [14,54] to calculate the variance in morphology metrics that could be individually explained by all local network metrics and growth factor treatments (). We then calculated the values for each morphology metric after correcting for the effect of local network metrics (see Methods). The raw and network-corrected values for two metrics, cell size and mean actin intensity, are shown in The influence of network properties can be clearly seen on cell size, where at 6 hours, large cell sizes are seen in the uncorrected but not corrected plots (). The effect of growth factor treatment can be clearly seen in network-corrected mean actin intensity ( and ), where VEGF and BDNF treatment have dose-dependent effects on mean actin intensity independent of cell crowding effects. Thus, this case study demonstrates the utility of cytoNet in detecting the independent effects of local cell crowding and growth factor perturbations on cell morphology.

coupled single cell morphology and signaling, and complex tissue topology

microscopy color images (TIF format)

mouse adipose tissue

In a second illustration of cytoNet’s utility to analyze intact tissue, we used cytoNet to characterize the pericapillary niche within adipose tissue. Specifically, we sought to understand the role of laminin α4, an extracellular matrix glycoprotein, in adipose tissue. Mice with a null mutation in the laminin α4 gene exhibit resistance to obesity and enhanced insulin sensitivity [55,56]. Understanding how the deletion of laminin α4 affects the spatial distribution of cells present in the adipose tissue can provide insight into the mechanisms underlying the functional change, and guide biomimetic models of the adipose perivascular niche [57,58]. In this Case Study example, the confocal images of adipose tissue and capillaries were segmented by manual tracing on the computer, and provided as input to cytoNet. Because blood vessels have noncircular shapes, the distance between the centroids of vessels and other objects may not give a good sense of proximity. As an alternative graph-generation approach, cytoNet can compute the minimum distance between object perimeters in order to define graph edges. The cell-capillary distances are measured in the plane of each confocal z-stack, which is compressed into its highest intensity projection. The resulting cell-to-cell perimeter distance table and cell area computations were used to determine differences between wild-type and knockout cells (). Differences were observed in the distances among adipocytes stained with the BODIPY lipid dye () and among integrin α7 expressing cells ( between the laminin α4 knockout and wild-type mice models. The observed adipocytes tended to be smaller in knockout tissue compared to wild type, though this was not statistically significant (). These characterizations are consistent with the observation that adipose in knockout mice is more similar to beige adipose tissue. Overall, these observations align with findings that the absence of laminin α4 leads to changes in stromal cell distribution in pericapillary niches within adipose tissue [59]. The resulting data can be used to guide studies into understanding the mechanisms underlying the effect of laminin α4 on adipose tissue function. Thus, this case study demonstrates the utility of cytoNet in detecting regional variations of cell structure within tissues and in addressing testable spatial hypotheses about tissue function.

Spatial Analysis of the Pericapillary Niche in Adipose Tissue.

Example confocal images of wild type (a) and knock out (b) adipose tissue and the corresponding output graph for the wild type image (e). Red = lectin (capillaries). Green = Bodipy (adipocytes). Yellow: integrin α7 positive cells. Violin plots of cell properties comparing wild-type and knockout (c, d, f-h). Distances are measured between the closest border pixels of pairs of objects. Fig 5F is adapted from reference (1). Values represent the results of analyzing confocal images from 17 samples (9 wild-type and 8 knockout). Error bars are mean +/- standard deviation. p-values were computed using the Wilcoxon rank sum test (*: p ≤ 0.05, ***: p ≤ 0.001).

Discussion

Advances in in situ profiling techniques have led to the generation of highly multiplexed imaging datasets describing tissue archicture in great spatial detail [1-6]. Spatially detailed imaging datasets have led to a proliferation of computational pipelines designed to test spatially driven biological hypotheses (). However, many of these analysis pipelines are designed specifically for spatial molecular expression data and are not generalizable to data obtained from other microscopy techniques. Further, due to their reliance on specialized scripts, many pipelines are not readily accessible to biological researchers without programming background.

Here we present cytoNet, a user-friendly pipeline for investigation of spatial hypotheses in cell-and tissue-based biological experiments. cytoNet is available through a web interface, eliminating the need to download and install software. Source code is also provided as MATLAB scripts for more advanced users. cytoNet uses Matlab version 2021a, with the following toolbox dependencies: Bioinformatics Toolbox, Signal Processing Toolbox, Distributed Computing Toolbox, Image Processing Toolbox and Statistics Toolbox. Pre-segmented masks provided as input to cytoNet are used to build network representations of spatial topography. Accompanying epifluorescence or confocal images are used to extract single-cell features and functional relationships. Lastly, network descriptors are combined with single-cell features to explore cell community effects on cell phenotypes.

We demonstrate the utility of cytoNet through four case studies. As shown in detail in our previous study [12], we harness an in vitro model of neuronal network formation from neural progenitor cells (NPCs) to demonstrate a fall and rise in network heterogeneity during neural differentiation. Accompanying functional network analysis through calcium imaging shows that these trends in community structure likely reflect a transition from global to hierarchical communication during the formation of neural circuits. We further use local neighborhood measures to explore the effect of cell community on cell cycle regulation, showing a density-dependent effect on cell cycle synchronization. Further experimental perturbations, for instance addition of chemicals to inhibit cell-cell communication, can be used along with the cytoNet pipeline to investigate more specific hypotheses on the role of cell community in regulating dynamic cell behavior.

Our second case study showed cytoNet’s capability for analyzing time-varying functional image sets. In this case, we characterized spatiotemporal calcium signaling recorded from intact brain tissue. Networks can be constructed based on the similarity of temporal behaviors of cells. The combination of the functional networks and spatial networks reveals local groups of cells with similar behaviors and assists in the development and testing of hypotheses of functional subsystems in neuronal tissue.

We also explored the differential effects of cell density and growth factor stimulation on human endothelial cells using cytoNet. By applying unsupervised clustering approaches on a suite of cytoNet-generated metrics describing cell morphology and local neighborhood, we show the presence of three cell phenotypes. These phenotypes reflect different cytoskeletal states and multicellular interactions indicative of collective behaviors like migration and proliferation. Further, we leverage a quantile multidimensional binning approach to investigate the differential effects of cell density and growth factor perturbations on cell morphology. This workflow can be used to comprehensively characterize the response of cells to chemical perturbations and aid in drug discovery. Case Study 4 illustrated another translational application of cytoNet: this time to study the effect of an extracellular matrix protein on the phenotype of adipose cells within perivascular niches.

Notably, two of the case studies were applied in vitro to human cells, and two were applied to in vivo image sets. Case Study 1 and 2 capitalized on cytoNet’s ability to integrate functional and structural graphs across time in a single mathematical framework. The other two cases illustrated the how cytoNet can be applied to optimize cell phenotyping (Case Study 3 and 4). All of the cases show how cytoNet can help guide hypotheses, inform biomimetic models or tailor therapeutic interventions that reflect a cell’s microenvironment.

The network model utilized by cytoNet is a versatile modeling framework that can be used to incorporate many hypotheses on cell-cell interactions and their role in cellular behavior. In future iterations, this framework can be expanded to incorporate non-binary interactions through weighted networks, shift the focus from individual nodes to motifs through simplicial complexes, and include dynamic reconfiguration of networks over time through multilayer networks. Further, once graphs have been defined, graph theory affords a rich array of metrics that can be used to probe network structure, only some of which were studied here. These include a variety of null graph models that can be used to test specific spatial hypotheses.

In summary, the cytoNet method provides a user-friendly spatial analysis software, leveraging network science to model spatial topography and functional relationships in cell communities. This framework can be used to quantify the structure of multi-cellular communities and to investigate the effect of cell-cell interactions on individual cell phenotypes.

Methods

Software

cytoNet is available as a web-based interface at https://www.QutubLab.org/how and associated scripts are available at https://github.com/qutublab/cytoNet.git. An overview of cytoNet as a resource for the Brain Initiative Alliance community is provided here, along with video tutorials: https://www.braininitiative.org/toolmakers/resources/cytonet/. See for instructions on using cytoNet.

cytoNet web interface: QutubLab.org/how

cytoNet contact: cytoNetProject@gmail.com

cytoNet Github code: https://github.com/qutublab/cytoNet.git

Brain Initiative Alliance: https://www.braininitiative.org/toolmakers/resources/cytonet/

cytoNet image analysis pipeline

The cytoNet pipeline begins with masks and accompanying microscope images. The microscope images may be any color or gray-scale based microscopy images (e.g., immunofluorescence, confocal) or a sequence of calcium images (). The provided mask is used to extract features of cells and to construct spatial and functional graphs (). Spatial graphs are created by having nodes represent mask objects and edges determined by object distance. Edges can be found by one of two methods for spatial graphs: by evaluating the distance between cell boundaries (type I graphs), or by evaluating the proximity of cells in relation to a threshold distance (type II graphs) (). The type I graphs are useful when detailed information of cell boundaries and morphology is available, such as in the case of membrane stains or cells stained for certain cytoskeletal proteins. The type II graphs work well with images of cell nuclei, where detection of exact cell boundaries is not possible. In both approaches, cells deemed adjacent to each other are connected through edges, resulting in a network representation. If calcium imaging sequences are provided as input, a functional graph is created based on correlations among calcium time series of different mask objects ().

Image segmentation

Image segmentation–the identification of salient foreground objects such as cells–is often the first step in image analysis. The quality of image segmentation can have a large effect on the ability of cytoNet to derive meaningful results. Precise detection of cell boundaries is important to derive accurate type I graphs, and precise detection of cell centroids and areas are important to derive accurate type II graphs. For functional networks alone, the dependency on segmentation of cell centroids or cell boundaries is minimal (). The cytoNet pipeline works with pre-segmented masks of images and accompanying microscope images. For users who do not have mask files, cytoNet includes basic image segmentation algorithms including thresholding and watershed operations to generate these masks. The segmentation algorithms included in cytoNet can be parameterized to work well for images with clear delineation of nuclei and cell borders, like the examples provided on the cytoNet website. The cytoNet code also provides frequency detection of cells, where a change in a functional marker (e.g., Ca2+ or FUCCI) delineates cell location. For object detection in more complex image sets, we point the user to programs that focus on cell segmentation [60-62]. Multiple research teams have made significant inroads into designing generalizable image segmentation algorithms, among them classic thresholding and watershed approaches [63], pixel-based classifiers [60] and more recently deep learning approaches [1,61,62]. These programs generate masks as output. Users may wish to implement them prior to analyzing community structure through cytoNet. Image segmentation and graph creation are handled separately by cytoNet, enabling flexibility for the user.

Generation of spatial networks

Type I graphs are generated as follows. Mask boundaries are expanded by 2 pixels and overlap of expanded masks is used to assign edges and build an adjacency matrix. Cells touching the image border are included in calculations of local network properties () for cells not touching the boundary but are excluded for the construction of the adjacency matrix. Type II graphs are generated as follows: for each pair of objects (nuclei), a threshold distance for proximity is defined as the average of the two object diameters, multiplied by a user-defined scaling factor (S). If the Euclidean distance between the object centroids is lower than the threshold distance computed, the pair of objects is connected with an edge. A scaling factor of 3 was consistently used for images of cell nuclei analyzed in this manuscript.

Generation of functional networks

Functional networks are generated from calcium imaging sequences as follows. Users may provide an optional mask to define regions of interest (ROIs). If a mask is not provided, cytoNet calculates a mask based on the maximum intensity projection of the image sequence provided. A time-varying fluorescence trace (F) is calculated for each ROI as the average of the pixel intensities within that ROI. The MATLAB detrend function is used to remove linear trends. Next, a moving baseline fluorescence (F0) for each ROI is calculated as the 8th percentile of intensity values in a window one-hundredth the total number of frames, before and after each frame. These parameters are important variables, which can be modified by the user as needed. The normalized fluorescence trace for the ROI is then calculated as F–F0/F0, followed by smoothing. In order to avoid false positives in the cross-covariance analysis, ROIs with low activity are filtered out by discarding traces with less than three spikes. Spikes are defined as local maxima whose normalized fluorescence is at least 25% of the maximum signal or 0.05, whichever is higher.

Functional networks are created using the method described by Smedler et al, [64] where cross-covariance between signals with a lag time of 0 is used to assign functional connections between pairs of cells (). A randomized dataset is generated by shuffling each signal in the original dataset at a random time point. The 99th percentile of cross-covariance values for the randomized dataset is used as a threshold for determining significant correlations.

Network metric computation

For both spatial and functional graphs, connectivity is denoted mathematically using an adjacency matrix, A, where A = 1 if there exists an edge between cells i and j, and 0 otherwise. This concise representation of hypothesized interactions among cells can be used to generate multiple descriptors at a local level for individual nodes and at a global level for the entire graph (). Extracted metrics are used to visualize and analyze local neighborhood effects on individual cell phenotypes (), as well as global cell community characteristics (). Examples of local metrics are number of connections (degree) or notions of centrality, such as ability to act as a bridge between different cell communities (betweenness centrality). Examples of global metrics include measures of modularity such as the number of connected components, and measures of information flow such as path length. All the network metrics described in Tables and were computed using custom-written code, building upon routines provided in [65].

Generation of random networks

For the purposes of evaluating density-dependent effects on network metrics and hypothesis testing, we provide random null network models. For spatial type II graphs, random spatial networks are created by placing spheres of diameter equal to the average equivalent diameter of all ROIs at random locations in the grid and re-applying connection rules (). For functional networks, random graphs are created by randomly rewiring existing graphs while preserving the degree distribution. Network metrics are computed for 100 iterations of random networks.

Cell culture

Human umbilical vein endothelial cells (HUVEC) were obtained from Lonza and cultured in EBM-2 medium (Lonza) supplemented with penicillin-streptomycin (Fisher Scientific) and EGM-2 SingleQuot bullet kit (Lonza). For imaging experiments, cells were cultured for different periods (6, 12 or 24 hours) in different combinations of vascular endothelial growth factor (VEGF, human recombinant; Millipore) and brain-derived neurotrophic factor (BDNF, human recombinant, Sigma-Aldrich). Concentrations used were in the range 50–100 ng/ml. Controls were the same culture period without growth factor treatments.

Immortalized human neural progenitor cells derived from the ventral midbrain (ReNCell VM) were obtained from Millipore. Cells were expanded on laminin-coated tissue culture flasks, in media containing DMEM/F12 supplemented with B27 (both Life Technologies), 2 μg/ml Heparin (STEMCELL Technologies), 20 ng/ml bFGF (Millipore), 20 ng/ml EGF (Sigma) and penicillin/streptomycin. For differentiation experiments, cells were cultured in medium lacking bFGF and EGF.

Dorsal root Ganglion mouse model

Dorsal laminectomies were performed on anesthetized Pirt-GCaMP3 mice exposing the dorsal root ganglia in the spinal L5 region. The spinal columns were stabilized under a laser-scanning confocal microscope. Stimuli were applied to the hind paw in one of four ways: 1) pressure (rodent pincher analgesia meter), 2) gentle mechanical stroke (brush or von Frey filament), 3) thermal stimuli (immersion in hot or cold water), 4) chemical stimuli (KCl, capsaicin, or TRPV1 agonist applied subcutaneously). Calcium image sequences were acquired at depths of up to 100 μm at 1–3 Hz at intervals of 4–6 seconds.

Laminin α4 knockout mouse model

Subcutaneous fat was separately collected from laminin α4 knock out mice (8 samples) and wild-type mice (9 samples). The samples were processed and incubated with integrin α7 antibody (1:100, Novus Biologics NBP1-86118) and Griffonia simplicifolia isolectin conjugated with Rhodamine (labels endothelial cells/blood vessels) followed by incubation with a second antibody (Alexa Fluor 647 Donkey Anti-Rabbit IgG, Abcam ab150075) and BODIPY to stain lipid. Images were collected by a Leica TCS SP8 Confocal Microscope.

NPC calcium image acquisition and processing

ReNCell VM neural progenitor cells were plated on LabTek chambered cover glasses for calcium imaging experiments. Cells were loaded with culture medium containing 3 μM of the fluorescent calcium indicator Fluo-4/AM (Life Technologies) and Pluronic F-127 (0.2% w/v, Life Technologies) for 30 min at 37°C. Imaging of spontaneous calcium activity was performed at 37°C using a 20X objective lens (N.A. = 0.75), with 488 nm excitation provided through a SOLA SE Light Engine (Lumencor). 16-bit fluorescence images were acquired at a sampling frequency of 1 Hz for a total duration of 15 min, using a Zyla 5.5 sCMOS camera (Andor). Following calcium imaging, samples were fixed, and nuclei were stained using DAPI. By navigating to the locations where calcium imaging was performed, manual co-registration was done to obtain immunofluorescence images of nuclei for the same fields of view.

Regions of interest (ROIs) were obtained by segmenting nucleus images using a local thresholding approach followed by the watershed algorithm. Functional networks were then created as already described.

Generation of FUCCI reporter neural progenitor cell lines

Stable reporter cell lines (FUCCI-ReN) were generated by sequentially nucleofecting ReNcell VM neural progenitor cells with an ePiggyBac [66] construct encoding mCherry-Cdt, Venus-Geminin, or Cerulean-H2B. Each construct introduced to the cells was driven by a CAG promoter containing a blasticidin (ePB-B-CAG-mCherry-Cdt1), puromycin (ePB-P-Venus-Geminin), or neomycin (ePB-N-Cerulean-H2B) resistance gene. Following each round of nucleofection, cells were cultured in the presence of appropriate antibiotics (2 μg/ml blasticidin, 0.1 μg/ml puromycin and 100 μg/ml neomycin).

Acquisition and processing of FUCCI-ReN time lapse videos

FUCCI-ReN cells were plated at different densities on chambered cover glasses (Fisher Scientific) coated with laminin. Cells were imaged after switching to differentiation medium containing phenol red-free DMEM/F12. Time-lapse imaging was performed using a Nikon Ti-E microscope equipped with a motorized stage, a cage incubator for environmental control (Okolab), a 20X objective lens (N.A. = 0.75), SOLA SE Light Engine for LED-based fluorescence excitation (Lumencor), appropriate filters for visualizing mCherry, Venus and Cerulean fluorescent proteins and a Zyla 5.5 sCMOS camera (ANDOR). 16-bit composite fluorescence images were acquired at 10-minute intervals for a total duration of 57.5 hours.

Grayscale images for each channel (H2B-Cerulean, Geminin-Venus and Cdt1-mCherry) were binarized using locally adaptive thresholding. Seeds for the watershed transform were generated using the regional minima from the distance transform of the grayscale images. Next, the watershed algorithm was applied to detect boundaries between overlapping cell nuclei. Finally, information from different channels were used to correct undersegmented nuclei ().

Acquisition and processing of HUVEC immunocytochemistry images

For imaging experiments, HUVECs were cultured on glass dishes coated with fibronectin (Sigma-Aldrich). After appropriate growth factor treatments, cultures were fixed with 4% paraformaldehyde, free aldehyde groups were quenched using 1 mg/mL sodium borohydride, and membranes were permeabilized with 0.2% Triton-X-100 solution in PBS. Actin fibers were visualized using an Alexa Fluor 488-phalloidin antibody (1:40, Molecular Probes) and microtubules were visualized using a mouse monoclonal anti-α-Tubulin antibody (1:250, Sigma-Aldrich) followed by a goat anti-mouse Alexa Fluor 647 secondary antibody. Nuclei were stained using Hoescht (Molecular Probes). 16-bit composite immunofluorescence images were acquired through a 20X objective (N.A. = 0.75) on a Nikon Ti-E epifluorescence microscope. Physical pixel size was 0.32 μm.

Fluorescence images were processed as described previously [67] (). Briefly, the following steps were used.

Contrast was enhanced using histogram equalization.

Images were smoothed using a 2D Gaussian lowpass filter.

Initial binarization was performed using Otsu’s method.

The binary image was dilated to fill in individual cell areas.

All objects <1% of the total image area were removed. This was called the final binary image.

A binary representation of the nuclear and microtubule image layers was generated using a high input threshold value. This was called the marker image.

Another binary image was created with values of 0 where either the final binary image (step 5) or the marker image (step 6) had a value of 1.

Watershed markers were generated by imposing the minimum of the complement of images obtained in steps 2 and 7. This image had black markers contained within cells to serve as basins for flooding, while cell areas themselves were represented by lighter pixels that served as the rising contours of the basins.

The watershed algorithm was implemented using Matlab’s built-in function to generate cell boundaries.

Masks generated in step 9 were refined by using composite images of microtubules and actin as the marker image (step 6).

In order to automate the threshold generation, the area of cell masks obtained from segmentation were compared to those obtained through thresholding with a high threshold. The entire process was then iterated until an acceptable area ratio was achieved.

Processing of in vivo calcium image sequences

Calcium image sequences from dorsal root ganglion models were processed as follows. To generate a mask, the calcium image sequence was first decomposed into individual grayscale frames. Next, for each pixel location, the maximum and minimum intensities were found across all frames. The differences between the maximum and minimum intensities were stored in an array (of delta values) and normalized. An initial segmentation of the delta values was done by thresholding using Otsu’s method, resulting in an initial binary mask. The initial mask was refined by computing a new threshold by applying Otsu’s method to only those delta values that were identified as foreground objects in the initial segmentation. The resulting binary image underwent a morphological closing with a disk of radius 3, and objects of fewer than 10 pixels were removed to generate the final mask. By this method, only active cells are segmented in the Case Study 2 functional analysis (). Cells whose Ca2+ activity did not change over time were not analyzed further.

To generate functional networks, edges were placed between two cells whenever: a) the two cells had the same ramp-up and ramp-down times, and b) the Euclidean distance between the centroids of the two cells was less than or equal to 10 times the mean of the diameter of each of the two cells.

Cluster analysis

We performed cluster analysis on the HUVEC imaging dataset using Shrinkage Clustering [68], a two-in-one clustering and cluster optimization algorithm based on matrix factorization that simultaneously finds the optimal number of clusters while partitioning the data. Cells whose features had the smallest sum of squares distance to the median values for each cluster were identified as representative cells for each cluster.

Correction of morphology metrics for effects of local network properties and treatment conditions

We performed quantile multidimensional binning [17] of cells for all 7 network metrics (5 bins per metric). The mean of each morphology metric was calculated for each multidimensional bin, and this mean was subtracted from the raw measurements to generate the network-corrected measurements for each cell. Treatment-corrected measurements were generated similarly by calculating the mean of each morphology metric under each treatment condition and then subtracting it from the raw measurements.

Variance explained by local network properties and treatment conditions

The variance explained by each factor was calculated using the following formula [56] V is the variance of the corrected measurements, and V is the variance of the uncorrected measurements.

Random graph creation for type II spatial graphs.

(a) Original grayscale image. (b) Binary mask with overlaid spatial graph. (c) Random graph created by placing spheres of diameter equal to the average equivalent diameter of all ROIs in (b) at random locations and reapplying connection rules. (d) Second iteration of random graph creation.

(PDF)

Click here for additional data file.

Image segmentation of HUVEC immunofluorescence images.

(a) Original grayscale image. (b) Image after adaptive histogram equalization and Gaussian filtering. (c) Binary image obtained using Otsu’s threshold, with small objects removed. (d) Complement of filtered image in (b). (e) Watershed basins obtained through imposing minimum of images in (d) and the marker image (obtained by combining the binary image in (c) and the image obtained through binarization of microtubules and nuclei). (f) Final cell borders.

(PDF)

Click here for additional data file.

Image processing steps for FUCCI-ReN nucleus images.

(a) Fluorescence image from H2B-Cerulean channel marking all nuclei. (b) Binary mask obtained through adaptive thresholding. (c) Image obtained through imposing minimum of distance transform of binary image in (b) and local minima. This image serves as seeds for the watershed algorithm. (d) Final mask obtained after watershed transform.

(PDF)

Click here for additional data file.

Correlation heatmap of local network metrics and morphology metrics for immunofluorescence HUVEC images.

All morphology and local network metrics (S1 and S2 Tables) were combined into a single matrix. The cluster dendrogram was obtained through hierarchical clustering of the covariance matrix using Pearson’s correlation as the similarity metric.

(PDF)

Click here for additional data file.

Metrics used to define endothelial cell morphology.

(PDF)

Click here for additional data file.

Cohen’s d effect size for treatment conditions on morphology metrics shown in Fig 4G and 4H) in the main text.

(PDF)

Click here for additional data file.

Time-lapse movie of sparse culture of FUCCI-ReN cells.

Magenta: Cdt1-mCherry, Green: Geminin-Venus. Time stamp is shown on top right corner.

(M4V)

Click here for additional data file.

Time-lapse movie of sparse culture of FUCCI-ReN cells with graph overlay.

Movie displays phase contrast frames from movie in S1 Video, with Cdt1(-)/mCherry(+) nuclei circled in magenta, Geminin(-)/Venus(+) nuclei circled in green and mCherry(-)/Venus(-) nuclei circled in blue. Yellow lines represents proximity edges.

(MP4)

Click here for additional data file.

Time-lapse movie of dense culture of FUCCI-ReN cells.

Magenta: Cdt1-mCherry, Green: Geminin-Venus. Time stamp is shown on top right corner.

(M4V)

Click here for additional data file.

Time-lapse movie of dense culture of FUCCI-ReN cells with graph overlay.

Movie displays phase contrast frames from movie in S3 Video, with Cdt1(-)/mCherry(+) nuclei circled in magenta, Geminin(-)/Venus(+) nuclei circled in green and mCherry(-)/Venus(-) nuclei circled in blue. Yellow lines represents proximity edges.

(M4V)

Click here for additional data file.

Instructions for using the cytoNet web-based user interface.

(PDF)

Click here for additional data file.

Table 1

Software tools for spatial analysis.

Software	Platform	Input	Output	Reference
histoCAT	MATLAB, standalone program	Imaging mass cytometry	User-guided cell neighborhood for selected cells, enrichments/depletion of cell-cell interactions based on comparison to spatially randomized data	[13]
Pelkmans lab	Module compatible with CellProfiler	Cell cultures	Local cell density, population size, cell islet edges	[14–17]
Cell-graph	Standalone tool	H&E stained tissue samples	Multiple graph metrics, e.g. clustering coefficient, network diameter	[9]
PySpacell	Python	Cell cultures	Statistical tests of magnitude and scale of spatial effects	[18]
SpatialDE	Python	Spatial transcriptomics datasets	Statistical tests of genes with spatial variation, spatial gene-clustering	[19]
trendsceek	R	Spatial transcriptomics datasets	Statistical tests of genes with spatial variation	[20]
cytoMAP	MATLAB	Histo-cytometry data	Multi-scale characterization of tissue structure	[21]
MuSIC	Cytoscape	Immunofluorescence and affinity purification mass spectrometry data	Intracellular protein positions and distances	[22]

Table 2

Software tools for calcium signal analysis.

Software	Platform	Input	Output	Reference
unnamed	MATLAB	Images	Segmentation, signal extraction, stimulus response analysis, assembly detection, network dynamics analysis	[31]
CaImAn	Python	Images	Motion correction, source extraction, deconvolution, registration	[32]
EZcalcium	MATLAB	Images	Motion correction, segmentation, signal extraction, deconvolution	[33]
NA 3	ImageJ, R	Images	Total activity value, variance area	[34]
CAVE	MATLAB	Images	Motion correction, ΔF/F calculation, cell detection, calcium trace analysis	[35]
CaSiAn	Java	Signal data	Peak and nadir detection, interspike interval and average period regression, signal correlation	[36]
SIMA	Python	Images	Motion correction, segmentation, signal extraction, ROI registration	[37]
Suite2p	MATLAB, Python	Images	Image registration, ROI detection, cell determination, activity and neuropil extraction, spike deconvolution	[38]
CNMF-E	MATLAB	Images	Contour detection, signal extraction	[39]
ABLE	MATLAB	Images	Contour detection, neuropil correction, signal extraction	[40]
SCALPEL	R	Images	Segmentation, signal extraction	[41]
MIN1PIPE	MATLAB	Images	Motion correction, segmentation, signal extraction, deconvolution	[42]
SamuROI	Python	Images	Image stabilization, event detection	[43]

Table 3

Global graph metrics and their normalization to account for network size.

n = number of nodes, m = number of edges.

Graph Metrics	Symbol	Definition
Node Count	n	Number of nodes
Edge Count	m	Number of edges
Fraction Area Cells	A	Fraction of total surface area in field of view covered by cells
Average Degree	avgeK	Average number of connections for a node in the network
Variance in Degree	varK	Variance of node degree sequence
Network Heterogeneity	NetworkHeterogeneity	Standard deviation of node degree sequence divided by mean of degree sequence–reflects tendency of network to contain hub nodes
Average Neighbor Degree	avgeNeighborK	Average degree of local neighborhood, averaged across all nodes
Variance in Neighbor Degree	varNeighborK	Variance of the average neighbor degree sequence
Network Efficiency	E	The average reciprocal of shortest path length across all pairs of nodes, E
Average Clustering Coefficient	C	Fraction of total possible links among the neighbors of a node that are actually present, averaged across all nodes, C
Number of connected components	nConnectedComponents	Number of disconnected sub-graphs in main graph
Average Size of Connected Components	avgeComponentSize	Average number of nodes in each connected component
Variance in size of connected components	varComponentSize	Variance in component size sequence
Network Diameter	networkDiameter	Longest shortest path length of network
Isolated Node Count	nIsolatedNodes	Number of nodes with no neighbors
Pair Node Count	nPairNodes	Number of independent pairs of nodes
Triangular loop count	nLoops3	Number of loops of 3 nodes
4-star motif Count	nStar4	Number of star motifs with one hub and three spokes
5-star motif count	nStar5	Number of star motifs with one hub and four spokes
6-star motif count	nStar6	Number of star motifs with one hub and five spokes
Rich-Club Metric Average	avgeRichClubMetric	Measure of the tendency of nodes with high number of links to be well connected among each other [44]; Computed for threshold degrees between 1 and (n-1)
Rich-Club Metric Variance	varRichClubMetric	Variance in rich-club metric for thresholds from 1 to (n-1)
Assortativity	Assortativity	Pearson correlation coefficient of degrees between pairs of linked nodes [45].

Table 4

Local neighborhood metrics calculated at the individual cell level.

Graph Metrics	Symbol	Definition
Degree	k	Number of neighbors one link away from cell of interest
Average Neighbor Degree	k n	Average degree of all neighboring cells
Clustering Coefficient	C	Number of edges in local neighborhood of a cell, divided by total possible connections
Local Efficiency	E l	Average shortest path length in local neighborhood
Node Closeness Centrality	c n	Sum of reciprocal distances in number of links to all other nodes
Node Betweenness Centrality	w n	Number of shortest paths that pass through a node
Shared Cell Border 1	S b	Total number of pixels shared with neighbors

1Relevant only for type I graphs