Network analysis of timelapse microscopy recordings.  
Jump to Full Text  
MedLine Citation:

PMID: 25278844 Owner: NLM Status: InDataReview 
Abstract/OtherAbstract:

Multicellular organisms rely on intercellular communication to regulate important cellular processes critical to life. To further our understanding of those processes there is a need to scrutinize dynamical signaling events and their functions in both cells and organisms. Here, we report a method and provide MATLAB code that analyzes timelapse microscopy recordings to identify and characterize network structures within large cell populations, such as interconnected neurons. The approach is demonstrated using intracellular calcium (Ca(2+)) recordings in neural progenitors and cardiac myocytes, but could be applied to a wide variety of biosensors employed in diverse cell types and organisms. In this method, network structures are analyzed by applying crosscorrelation signal processing and graph theory to singlecell recordings. The goal of the analysis is to determine if the single cell activity constitutes a network of interconnected cells and to decipher the properties of this network. The method can be applied in many fields of biology in which biosensors are used to monitor signaling events in living cells. Analyzing intercellular communication in cell ensembles can reveal essential network structures that provide important biological insights. 
Authors:

Erik Smedler; Seth Malmersjö; Per Uhlén 
Related Documents
:

6654524  Establishment and characterization of a human colon carcinoma cell line (kms4) from a ... 25433604  Characterization of sequential collagenpoly(ethylene glycol) diacrylate interpenetrati... 17644404  Comparison of different cell lines and incubation times in the isolation by the shell v... 3080224  Chromosome 14 marker appearance in a human b lymphoblastoid cell line of nonmalignant o... 11447194  Argininespecific protease from porphyromonas gingivalis activates proteaseactivated r... 25339614  Reduced myotube diameter, atrophic signalling and elevated oxidative stress in cultured... 
Publication Detail:

Type: Journal Article Date: 20140917 
Journal Detail:

Title: Frontiers in neural circuits Volume: 8 ISSN: 16625110 ISO Abbreviation: Front Neural Circuits Publication Date: 2014 
Date Detail:

Created Date: 20141003 Completed Date:  Revised Date:  
Medline Journal Info:

Nlm Unique ID: 101477940 Medline TA: Front Neural Circuits Country: Switzerland 
Other Details:

Languages: eng Pagination: 111 Citation Subset: IM 
Export Citation:

APA/MLA Format Download EndNote Download BibTex 
MeSH Terms  
Descriptor/Qualifier:

Full Text  
Journal Information Journal ID (nlmta): Front Neural Circuits Journal ID (isoabbrev): Front Neural Circuits Journal ID (publisherid): Front. Neural Circuits ISSN: 16625110 Publisher: Frontiers Media S.A. 
Article Information Download PDF Copyright © 2014 Smedler, Malmersjö and Uhlén. openaccess: Received Day: 08 Month: 7 Year: 2014 Accepted Day: 27 Month: 8 Year: 2014 Electronic publication date: Day: 17 Month: 9 Year: 2014 collection publication date: Year: 2014 Volume: 8Elocation ID: 111 PubMed Id: 25278844 ID: 4166320 DOI: 10.3389/fncir.2014.00111 
Network analysis of timelapse microscopy recordings  
Erik Smedler^{1}^{†}  
Seth Malmersjö^{2}^{†}  
Per Uhlén^{1}^{*}  
^{1}Unit of Molecular Neurobiology, Department of Medical Biochemistry and Biophysics, Karolinska InstitutetStockholm, Sweden 

^{2}Department of Chemical and Systems Biology, School of Medicine, Stanford UniversityStanford, CA, USA 

[editedby] Edited by: Mark A. Frye, University of California Los Angeles, USA [editedby] Reviewed by: Emre Yaksi, Vlaams Instituut voor Biotechnologie, Belgium; Jacob Aptekar, Howard Hughes Medical Institute, USA Correspondence: *Correspondence: Per Uhlén, Department of Medical Biochemistry and Biophysics, Karolinska Institutet, Scheeles väg 2, SE171 77 Stockholm, Sweden email: per.uhlen@ki.se [other] This article was submitted to the journal Frontiers in Neural Circuits. [presentaddress] †These authors have contributed equally to this work. 
Intercellular communication has been extensively investigated in the adult brain, in which cells form neural network circuits whose activities underlie the basic functions of the brain (Buzsaki, ^{2004}). During the postnatal period, developing networks exhibit spontaneous correlated neuronal activity that plays a central role in their establishment (Katz and Shatz, ^{1996}; Khazipov et al., ^{2004}; Cang et al., ^{2005}; Kandler and Gillespie, ^{2005}; Nicol et al., ^{2007}). In developing neocortical structures, synchronized activity has been observed at scales ranging from correlated pairs of neural progenitor cells (Owens and Kriegstein, ^{1998}) to gap junctionsynchronized cortical columns (Yuste et al., ^{1992}; Kandler and Katz, ^{1998}; Dupont et al., ^{2006}). To understand the nature and physiological roles of these networks, it is essential to identify them and characterize their structures.
Graph theory is the study of graphs (mathematical term for networks) consisting of nodes (also called vertices) and links connecting the nodes (also called edges) (Feldt et al., ^{2011}). Many relations and dynamic processes can be modeled by graphs, in such diverse fields as biology, social sciences, information systems, and transportation systems (Newman, ^{2010}). When studying biological networks, a typical node can represent a cell (Malmersjö et al., ^{2013a}), specific brain area (Biswal et al., ^{2010}), or protein (Jeong et al., ^{2001}). A link between two nodes can represent a functional interaction, e.g., via ion fluxes (Malmersjö et al., ^{2013a}), synapses (Bonifazi et al., ^{2009}), physical connections by axon bundles (Sporns et al., ^{2005}), or protein–protein interaction (Jeong et al., ^{2001}). The method presented here uses crosscorrelation analysis and graph theory to evaluate network topologies among living cells (Figure 1A) exhibiting spontaneous Ca^{2+} signaling (Figure 1B). The analysis identifies cell pairs with highly correlated Ca^{2+} signals (Figure 1C) or with uncorrelated Ca^{2+} signals (Figure 1D). Furthermore, the method can also reveal highly connected “hub cells” (scalefreeness) and high clustering accompanied by short internodal distances (smallworldness). Such functional network designs are effective in many types of biological systems (Barkai and Leibler, ^{1997}; Barabasi and Oltvai, ^{2004}). Furthermore, network analysis has also been used in fMRI based studies of cognitive neuroscience (Sporns, ^{2014}). By comparing the activity of healthy and diseased subjects, changes in network structure and activity can be linked to physiology and pathophysiology (Honey and Sporns, ^{2008}).
We developed a method to identify and characterize network structures from timelapse microscopy recordings, and implemented it in MATLAB. The resultant software tool is easy to extend and use. The method was originally designed to study network formations among large populations of neural progenitor cells exhibiting spontaneous intercellular Ca^{2+} activity (Malmersjö et al., ^{2013a},^{b}). However, this method could potentially be used to analyze any kind of cellular activity that is detectable in living single cells by timelapse microscopy, e.g., contractile activity measured by brightfield microscopy (Kitambi et al., ^{2012}) or activity measured using various biosensors, such as NFκ B, cAMP, or ERK (Nelson et al., ^{2004}; Malmersjö et al., ^{2010}; Aoki et al., ^{2013}; Everett and Cooper, ^{2013}). Thus, this method could be used to identify and investigate the details of uncharacterized network structures.
Cells were loaded with the Ca^{2+}sensitive fluorescence indicator Fluo3/AM (5 μM, Invitrogen) at 37°C for 30 min in KrebsRinger's buffer containing 119.0 mM NaCl, 2.5 mM KCl, 2.5 mM CaCl_{2}, 1.3 mM MgCl_{2}, 1.0 mM NaH_{2}PO_{4}, 20.0 mM HEPES (pH 7.4), and 11.0 mM dextrose. Measurements of cytosolic Ca^{2+} were carried out in KrebsRinger's buffer at 37°C using a heatcontrolled chamber (QE1; Warner Instruments) and a cooled electronmultiplying chargedcoupled camera (QuantEM 512SC, Photometrics) mounted on an upright fixed stage microscope (Axio Examiner.A1, Carl Zeiss) equipped with a 20× 1.0 N.A. lens (Carl Zeiss). Excitation at 480 nm was assessed with an illumination system (DG4, Sutter Instrument) at sampling frequency of 0.2–1 Hz (T = 1–5 s). MetaFluor (Molecular Devices) was used to control all devices and to analyze acquired images. The cellfree area was created by making a cut with a fine syringe (BD Microlance™ 3, 0.4 × 19 mm) in confluent HL1 cells. Dishes were placed in an incubator for 5 h before imaging.
Neural progenitors were derived from mouse embryonic stem cells as described before (Malmersjö et al., ^{2013a}). HL1 cells were cultured as previously described (Claycomb et al., ^{1998}).
Crosscorrelation was used to determine whether two cells were functionally interconnected. Crosscorrelation analysis is a mathematical method for quantifying the linear similarity between two waves as one of them is shifted in time (Brockwell and Davis, ^{1998}). When crosscorrelation analysis is applied in signal processing, the waves are typically time series consisting of discrete sets of data points [X_{t}, t ∈ T], e.g., images acquired by timelapse microscopy. The normalized version of the crosscorrelation function, i.e., the crosscovariance, is commonly used for imageprocessing applications in which the brightness of the image is the quantitative measure. In MATLAB, the crosscovariance is implemented as xcov:
(1)
cxy(m)={∑n = 0N − m −1(x(n+m)−1N∑i = 0N − 1xi) (yn*−1N∑i = 0N − 1yn*) if m≥0cyx(−m)* if m<0 
Here, m is the lag, N is the number of time points, n is the summation index, and x and y are the two time series. Because N is a finite number, the above function (Equation 1) is just an estimation of the real crosscovariance function:
(2)
φxy(μ)=E[(xn + m−μx)(yn−μy)] 
Here, μ_{x} and μ_{y} are the mean values of the stochastic processes (the time series are modeled as stochastic) and E is the expectation value operator (the average value from multiple samples). If time is fixed in the crosscovariance function (Equation 2), it will result in the wellknown correlation coefficient, also known as the Pearson correlation, a real number between 1 and 1. A correlation coefficient equal to 0 indicates no linear relation between the waves, whereas a coefficient equal to 1 or 1 demonstrates a perfect linear relation.
Two time series might be highly correlated even if one of them is shifted in time. Calculating the correlation as a function of lag enables determination of the maximum correlation despite lag. Figure 2A shows two sine waves with identical frequency, but different amplitudes and phases. Figure 2B shows correlation as a function of lag for the two sine waves. The phase shift is 2.5 s. Note that the correlation function is amplitudeindependent and only considers the relative amplitude. In some cases, for example neurons interconnected with synapses, the identified lag could be related to the pausing time between two neurons. However, most often this effect is interpreted as an effective phase shift.
Before calculating the correlation between two signals, they can be filtered by subtracting underlying trends; this process is called trend correction. For instance, bleaching or focus shifts might lead to a gradual decay, superimposed on the actual signal. By fitting the signals to a polynomial function with a certain degree (for example a linear function for linear trends), this effect can be reduced.
It is important to decide a cutoff that filters out insignificant correlations. We have developed a method for determining such cutoff values using a scrambled data set. A scrambled data set f_{scrambled} is created by shuffling the individual time series f to random starting points t (Equation 3). Thus, each original time series is divided into two parts at a random position and then put together again in the opposite order. Figure 3 illustrates a time series between t_{0} and t_{n} (Figure 3A) that is shuffled to t_{x} (Figure 3B). This procedure is then repeated for alltime series in the data set in a cyclic fashion:
(3)
fscrambled[1:N]=(f[t:N] f[1:t−1]) 
The total activity in the original data set and the scrambled data set are thereby conserved. The mean or the 99th percentile of the crosscovariance values of the scrambled data set can then be applied as the cutoff value.
Graph theory can be used to characterize the topology of a network (Feldt et al., ^{2011}). In graph theory, several network parameters are calculated for each network to determine the network's type (Newman, ^{2010}), shown in Figures 4A,B. Connectivity is here defined as the number of cells with a correlation coefficient larger than the cutoff, divided by the total number of cells. This is a measure of the degree of connections for the network, because a high cutoff results in a low connectivity and vice versa. The edge density, also referred to as connectance (Newman, ^{2010}), is defined as the number of edges divided by the maximum number of edges. The neighbors of a node are all the nodes connected to it in one step. The degree of a node is its number of neighbors; hence the degree distribution P(k) of a network is the distribution of nodes with a degree equal to k. P(k) is obtained by counting the number of nodes N(k) with k = 1, 2, 3, … connections and dividing by the total number of nodes.
Classical models of networks occurring in nature are either regular or random, as shown in Figure 4A. In a regular network, each node is connected to k other nodes (Feldt et al., ^{2011}). Figure 4A shows a regular network with k = 4. In a random network, nodes are connected with links stochastically, resulting in a Poissonshaped degree distribution around pN, where p is the probability and N is the total number of nodes, plotted in panel Figure 4C. In this example, 1000 nodes are randomly connected to each other with 50% probability.
Smallworld networks combine features of both regular and random networks, with high clustering as in regular networks but short internodal distances as in random networks (Watts and Strogatz, ^{1998}). According to the WattsStrogatz model, a smallworld network can be generated by randomly rewiring links in a regular network (Watts and Strogatz, ^{1998}). The properties of smallworld networks are assessed by calculating the mean clustering coefficient C and the mean shortest path length L of the network using the MatlabBGL library (http://www.mathworks.se/matlabcentral/fileexchange/10922matlabbgl). Hence, a smallworld network is characterized by the following relations:
(4)
σ=CCrand≫1 and λ=LLrand≈1 
The clustering coefficient is the number of neighbors of a node that are also neighbors of each other divided by the total number of possible links between the neighbors, as shown in Figure 1A. Thus, it reflects the number of groups in a network. The shortest path length is the minimum number of nodes that must be passed to travel from one node to another, as depicted in Figure 1A. The values of C and L are then compared with the corresponding values of C_{rand} and L_{rand} for a randomized version of the network. Figure 4A shows a simplified scheme for classification of smallworld (red), regular (blue), and random (green) networks according to their mean clustering coefficient σ and mean shortest path length λ. The interrelationship between the clustering coefficient (σ), shortest path length (λ), and smallworld parameter (σ/λ) as a function of rewiring probability of a regular network is shown in Figure 4B. In this simulation, 1000 nodes were connected in a circle to their two closest neighbors, and random links were rewired with increasing probability.
A network is defined as possessing smallworld characteristics if the mean path length is as short as in the corresponding random network, whereas the mean clustering coefficient is higher. At the level of the central nervous system, smallworld networks are thought to promote efficient information flow at a low wiring cost (Achard and Bullmore, ^{2007}).
The Barabási–Albert model of preferential attachment states that a scalefree network can be generated by allowing a random network to grow according to preferential attachment (Barabasi and Albert, ^{1999}). If the degree distribution approximately follows a power law (a heavytailed function without a clear mean value or scale), the network is defined as scalefree.
(5)
P(k)∝k−γ↔log(P(k))∝−γlog(k) 
In Figure 4D the degree distribution is plotted in a log–log scale that shows a linear function if possessing scalefree characteristics. The scale factor γ is inferred using the leastsquares method; however, some authors have argued that the maximum likelihood method is superior (Clauset et al., ^{2009}). Scalefree networks have some nodes with many neighbors that can act as hubs (Barabasi and Albert, ^{1999}). On theoretical grounds, these networks are thought to be robust (Albert et al., ^{2000}): as depicted in Figure 4E. In Figure 4D, a network with 100 start nodes with random connections was grown by connecting 900 new nodes in a “richgetricher” fashion.
MATLAB version 7.12.0 (http://www.mathworks.com), or later, including the Signal Processing toolbox, is required to implement the analytical method we describe here. The current versions of MATLAB run on multiple operating systems. NetworkIdentification.m and NetworkAnalysis.m are the files containing the main code, supported by the function files mic2net.m, crosscorrelation.m and pickcells.m. The program is initiated by mic2net.m. All files are included as supplementary material for this article (Supplementary Materials). The MatlabBGL library for MATLAB is required for the topology analysis that calculates the shortest path length and clustering coefficients, and is provided as supplementary material online. The original source for the MatlabBGL can be found elsewhere: (http://www.mathworks.com/matlabcentral/fileexchange/10922).
Fiji (http://www.fiji.sc/Fiji) or any other imageprocessing software that can extract the mean intensities, xcoordinates, and ycoordinates from multiple regions of interest (ROIs) in timelapse microscopy recordings is required to analyze timelapse microscopy recordings.
The first step of the analysis is to extract data from the timelapse microscopy recordings using a fluorescent biosensor to monitor e.g., intracellular Ca^{2+} signaling (Figures 5A,B). This is carried out by calculating the mean fluorescence value of each cell marked with an ROI. Because network structures are spatial structures, we also need the x and ycoordinates of each cell. Each ROI is therefore assigned its x and ycoordinates. The values are saved in a file to be loaded into MATLAB. Next, all possible combinations of cell pairs are examined using crosscorrelation analysis (Figures 5C,D). A correlation value between 0 and 1 is calculated for each cell pair that reflects the strength of interaction. Note that only the absolute value of the crosscorrelation function is saved in the correlation matrix, thus neglecting the difference between correlation and anticorrelation. To discriminate between real and false correlations, a cutoff value should be established (Figures 5E,F). Our method provides an unbiased way to calculate such cutoff values using signals that have been scrambled by randomly shuffling each cell signal in the time domain (see Section Materials and Methods); thus the total activity within the data set is preserved. Thereafter, the cutoff value can be set as the 99th percentile of all correlation coefficients within the scrambled data set (Malmersjö et al., ^{2013a},^{b}). Thus, two cells are defined as interconnected if their correlation coefficient exceeds the cutoff value.
We applied this method to the investigation of functional networks in neural progenitor cells exhibiting spontaneous Ca^{2+} activity. Before imaging, neural progenitor cells derived from embryonic stem cells were loaded with the Ca^{2+} indicator Fluo3/AM for 30 min. Timelapse recording was performed by fluorescence microscopy, using a sampling frequency of 0.2 Hz. After approximately 30 min of recording, the timelapse experiment was terminated, and the network analysis could start. Using the Fiji imageprocessing software, each cell was marked with an ROI and the mean fluorescence value was calculated (see Supplementary Movie 1 for a stepbystep guide). The values were saved in a tabdelimited text file that can be loaded into MATLAB (see Supplementary Movie 2 for a stepbystep guide). The sampling frequency and the physical pixel size were stated explicitly in MATLAB. Signal artifacts, caused, e.g., by focus change and/or drug application, were filtered out by choosing a “clean” time window for the analysis. Because a data set containing cells with no activity will result in a correlation analysis with false positives, only active cells were chosen for further analysis (Figures 6A,B). Cells with very few peaks (Figure 6C) or dying cells that create a sustained plateau increase in the signal (Figure 6D) can also result in false positives. Active cells could be defined, for example, as cells exhibiting three transient intensity increases or more that exceed 15% over the baseline (Malmersjö et al., ^{2013a}) or more than two standard deviations over the noise. Here noise is defined as the baseline fluctuations of a nonactive cell. If the data set has a low signaltonoise ratio, digitization is required. The process of digitization replaces numbers that do not exceed the threshold with 0 and all others with 1. To further optimize the data set before analysis, singlecell signals can be trendcorrected and/or maximal lag can be selected (see Section Materials and Methods). Thereafter, it is advisable to manually check cell pairs with high correlation coefficients in order to filter out false positives. Another source of false positive correlations could arise from defining two separate ROIs that inadvertently overlaps covering parts of the same cell. To avoid this, an alternative approach would be to automatically segment cells with a nuclear marker used as a mask.
The network analysis generates a network structure that can be plotted on top of a microscope image. Figure 7A shows the network plot of correlation coefficients for neural progenitor cells exhibiting spontaneous Ca^{2+} activity. Only connections exceeding a cutoff (see below) were plotted. The network plot suggests that neighboring cells are strongly correlated (Figure 7A). By studying the distance distribution of the correlation coefficients in detail, we indeed observed that cells close to each other were strongly correlated (Figures 7B,C). We previously tested this observation experimentally by pharmacological inhibition of gap junctions with Octanol (1 mM) and shRNA knockdown of the Connexin 43 gene, both of which disrupt the network activity (Malmersjö et al., ^{2013a}). As described above, the cutoff was determined by calculating the mean of the 99th percentile of the scrambled data (see Section Materials and Methods). Plotting the distance distribution of the scrambled data revealed no cell pairs with high correlations and short internodal distances (Figures 7D,E). Thus, the cutoff was set to 0.39 to filter out the bulk of cells with low correlation values (Figures 7B,C,E; red lines). The degree distribution revealed a typical scalefree network structure (Figure 7F). Calculation of the network parameters σ = 7.7 ± 0.92, λ = 0.97 ± 0.056, and γ = −1.2 ± 0.049 (N = 6) yielded a characterization of this network as a smallworld network with scalefree topology.
Network analysis was also performed on cardiac HL1 cells interconnected by gap junctions (Claycomb et al., ^{1998}). These cells exhibited spontaneous Ca^{2+} activity with network properties (Figure 8A). The analysis revealed that these networks did not have smallworld characteristics (σ = 1.06 ± 0.020 and λ = 1.0 ± 0.0083, N = 4). Instead, the HL1 networks showed similarities to a random network. To test the strength of the method, the network was divided into two subnetworks by creating a cut with a syringe (Figure 8A). The distance distribution of HL1 cells revealed a majority of cell pairs with high correlations and short internodal distances (Figure 8B), although not as clearly as in the neural progenitor cells described above (Figure 7B). Studying the maximal correlation as a function of lag revealed a highly synchronized cell population (Figure 8C). Next, we tested the dependency of correlations on distance by randomizing the cell positions of the cells. After randomizing the cells' coordinates, the cut was no longer visible (Figure 8D). Plotting the distance distribution of the data with randomized positions resulted in more evenly spread cell pairs (Figure 8E), revealing that the randomization abolished the dependency of correlation on distance (Figures 8D,E). In this analysis, the cutoff was set to 0.23, using the same method as described above for neural progenitor cells, to filter out the bulk of cells with low correlation values. Overall, the HL1 network exhibited similarities with a random network, as indicated by the degree distribution (Figure 8F).
Usage of our method will generate a number of different parameter values describing the corresponding functional network: mean value of all correlations, mean value of correlations above cutoff, 99th percentile of all correlations, connectivity, edgedensity, λ, σ, smallworld parameter σ/λ, and γ. How to statistically compare and evaluate networks is nontrivial and still being studied. As a first approximation, different networks of different type of cells or treatments can be statistically compared using for instance the Student's ttest, thus assuming normally distributed values (Malmersjö et al., ^{2013a}). If, for example, one treatment is hypothesized to decrease the overall connectivity within a network, mean values of a number of samples of edgedensity and connectivity can be calculated and compared between groups. However, it is more accurate to compare several parameters simultaneously and perform a Bonferroni correction to counteract the problem with multiple comparisons. Using a Bonferroni correction sets a lower bound on the significance threshold. Regarding degree distributions, for example the scale factor for powerlaw distributions can be inferred using standard leastsquare or the maximum likelihood method (Clauset et al., ^{2009}).
The method described here identifies network formations in data sets acquired by livecell imaging and enables objective analyses of these networks by characterizing their network topology. The information provided here, together with the accompanying software tools, should enable the experimental biologists without programming skills to identify and quantitatively analyze networks of timelapse microscopy recordings. The crucial steps of preprocessing the data and determining the cutoff are discussed further below. Comparisons with other network methods and similarities between cell networks and social networks are also discussed.
The outcome of the analysis depends critically on the quality of the input data. It is therefore essential to preprocess the data; for example, by selecting a time window that includes no artifacts or performing trend correction to remove the influences of fluorescence/dye bleaching or dye leakage. It is also crucial to preselect cells that are active: the crosscovariance function only considers the relative amplitude; consequently, inactive cells with normal background noise appear correlated in the program, resulting in falsepositive network links (Figures 6A–C). Therefore, inactive cells should be manually removed from the data set. Cells with sustained plateau increases may also result in falsepositive network links (Figure 6D).
In the two experiments above, different cutoff values were used to filter out only strongly connected cells. The cutoffs were chosen as the mean of the 99th percentile of the correlation coefficients for the scrambled data. Changing the cutoff affected the edge density and the connectivity, two different definitions of the number of connected cells (Figure 9A). In addition, the smallworldness of a network is strongly related to its edge density (Humphries and Gurney, ^{2008}). As illustrated by a functional network of neural progenitors, increasing the cutoff (the same as decreasing the edge density) increases the smallworld parameter (Figures 9B,C). Thus, it is very important to choose an adequate cutoff; for example, by studying the correlation distribution of a scrambled data set. It is worth pointing out that a network of strongly correlated cells does not automatically have the properties of smallworldness or scalefreeness, as is clear from a theoretical point of view (Watts and Strogatz, ^{1998}; Barabasi and Albert, ^{1999}) and from our two examples of highly correlated cells with different network structures (Figures 7A, 8A).
Our method allows identification and analysis of network structures from cellimaging data sets acquired by timelapse microscopy using code written in MATLAB. MATLAB is a widely used numerical computing environment and programming language developed by MathWorks (http://www.mathworks.com). A number of software tools for analyzing networks exist, but these are most often suited for genetic and molecular data, rather than timelapse microscopy recordings (Brohee et al., ^{2008}; Doncheva et al., ^{2012}). LEDA (Mehlhorn and Naher, ^{1995}) (C++) and QuACN (Mueller et al., ^{2011}) (R) are tools for users with programming skills, whereas Pajek (Batagelj and Mrvar, ^{2004}), Cytoscape, and yED are more graphically oriented solutions. A common feature of all of these software tools is that they are primarily used to analyze already identified networks, and are therefore not applicable to data sets consisting of timelapse recordings of unknown network activity.
We hypothesized that the scalefree network structure exhibited by neural progenitor cells (Figure 7A) was a consequence of the fact that dividing cells are more prone to connect to already existing highly connected nodes, in accordance with the Barabási–Albert model of preferential attachment (Barabasi and Albert, ^{1999}). Highly active neural progenitors with many neighbors divide more frequently and tend to connect to their daughter cells. The same phenomenon exists in social networks, in which people with many friends can acquire new friends more easily than those with few. By analogy to the WattsStrogatz model of long shortcuts (Watts and Strogatz, ^{1998}), the smallworldness in the central nervous system could be generated by extensions from neurons that connect to other neurons at relatively long distances, thereby decreasing the mean internodal distance. Examples of smallworld networks include social networks of actors in Hollywood, in which the distance between two random actors is shorter than expected and cliques are present (Watts and Strogatz, ^{1998}). An example of a scalefree smallworld network is how airlines connect the world through nodes of airports (Guimera et al., ^{2005}). A random disruption to one of the thousands of airport around the world would usually not disturb the flow of travelers, but a shutdown of a hub, such as London Heathrow Airport, could severely harm the network. Hence a scalefree smallworld network has a good tolerance for random deletion of nodes, but low tolerance for a directed attack to a hub. Graph theory predicts that such network designs are effective for biological systems, since they enable efficient information transfer and robustness against failure of single cells (Barabasi and Oltvai, ^{2004}).
Studying network structures is becoming increasingly popular in many subfields of biology, in part because network designs of biological systems resemble the Internet, social networks, and transportation systems. Studies of networking in cell biology and social networks will benefit from each other and further our knowledge of how cells interact to build a functioning organism. Obviously, this is even more important in the field of neuroscience where neural circuits perform computations dependent on their structure.
Currently, most analyses are twodimensional. However, biological systems are 3dimensional. Therefore, future analyses should be performed in three dimensions. Modern twophoton laser scanning microscopes and lightsheet microscopy systems should make it possible to image intact organs or whole animals in real time. Experiments using these image techniques generate huge data sets that require fast computers with large storage capacities to perform network analyses.
In the future, network analyses may benefit from employing more sophisticated methods to identify network links. This may include Granger causality test, which is a statistical hypothesis test for determining whether one time series is useful in forecasting another (Seth, ^{2010}). Furthermore, both crosscorrelation analysis and Granger causality detects linear relationships, opening up for further improvement for detection of nonlinear connections such as transfer entropy (Vicente et al., ^{2011}).
Erik Smedler, Seth Malmersjö, and Per Uhlén designed the algorithm. Erik Smedler and Seth Malmersjö wrote the code. Erik Smedler and Seth Malmersjö analyzed the data. Erik Smedler, Seth Malmersjö, and Per Uhlén wrote the manuscript.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank Dr. Zachi Horn, Dr. Markus Kruusmägi, David Larsson, and Hampus Sunner for useful comments on this work and Dr. William Claycomb (Louisiana State University) for kindly providing immortalized HL1 cardiomyocytes. This study was supported by the Swedish Research Council Grants 20056682, 20057204, 20086577, 20093364, 20104392, and 20133189 (Per Uhlén), Swedish Foundation for Strategic Research (Per Uhlén), Linnaeus Center in Developmental Biology for Regenerative Medicine (DBRM) (Per Uhlén), the Knut and Alice Wallenberg Foundation Grant to the Center of Live Imaging of Cells at Karolinska (CLICK) Institutet (Per Uhlén), the Royal Swedish Academy of Sciences (Per Uhlén), and Hjärnfonden Grant FO20110300 (Per Uhlén), and the Swedish Society for Medical Research (Seth Malmersjö and Per Uhlén).
The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fncir.2014.00111/abstract
Click here for additional data file (Video1.WMV)Click here for additional data file (Video2.WMV)
Click here for additional data file (Presentation1.ZIP)
References
Achard S.,Bullmore E.. (Year: 2007). Efficiency and cost of economical brain functional networks. PLoS Comput. Biol. 3:e1710.1371/journal.pcbi.003001717274684  
Albert R.,Jeong H.,Barabasi A. L.. (Year: 2000). Error and attack tolerance of complex networks. Nature406, 378–38210.1038/3501901910935628  
Aoki K.,Kumagai Y.,Sakurai A.,Komatsu N.,Fujita Y.,Shionyu C.,et al. (Year: 2013). Stochastic ERK activation induced by noise and celltocell propagation regulates cell densitydependent proliferation. Mol. Cell52, 529–54010.1016/j.molcel.2013.09.01524140422  
Barabasi A. L.,Albert R.. (Year: 1999). Emergence of scaling in random networks. Science286, 509–51210521342  
Barabasi A. L.,Oltvai Z. N.. (Year: 2004). Network biology: understanding the cell's functional organization. Nat. Rev. Genet. 5, 101–11310.1038/nrg127214735121  
Barkai N.,Leibler S.. (Year: 1997). Robustness in simple biochemical networks. Nature387, 913–91710.1038/431999202124  
Batagelj V.,Mrvar A.. (Year: 2004). Pajek  analysis and visualization of large networks, in Graph Drawing Software, eds Jünge M.,Mutzel P. (Berlin; Heidelberg: Springer), 77–10310.1007/9783642186387_4  
Biswal B. B.,Mennes M.,Zuo X. N.,Gohel S.,Kelly C.,Smith S. M.,et al. (Year: 2010). Toward discovery science of human brain function. Proc. Natl. Acad. Sci. U.S.A. 107, 4734–473910.1073/pnas.091185510720176931  
Bonifazi P.,Goldin M.,Picardo M. A.,Jorquera I.,Cattani A.,Bianconi G.,et al. (Year: 2009). GABAergic hub neurons orchestrate synchrony in developing hippocampal networks. Science326, 1419–142410.1126/science.117550919965761  
Brohee S.,Faust K.,LimaMendez G.,Vanderstocken G.,Van Helden J.. (Year: 2008). Network analysis tools: from biological networks to clusters and pathways. Nat. Protoc. 3, 1616–162910.1038/nprot.2008.10018802442  
Buzsaki G.. (Year: 2004). Largescale recording of neuronal ensembles. Nat. Neurosci. 7, 446–45110.1038/nn123315114356  
Cang J.,Renteria R. C.,Kaneko M.,Liu X.,Copenhagen D. R.,Stryker M. P.. (Year: 2005). Development of precise maps in visual cortex requires patterned spontaneous activity in the retina. Neuron48, 797–80910.1016/j.neuron.2005.09.01516337917  
Clauset A.,Shalizi C. R.,Newman M. E. J.. (Year: 2009). PowerLaw distributions in empirical data. Siam Rev. 51, 661–70310.1137/070710111  
Claycomb W. C.,Lanson N. A.,Stallworth B. S.,Egeland D. B.,Delcarpio J. B.,Bahinski A.,et al. (Year: 1998). HL1 cells: a cardiac muscle cell line that contracts and retains phenotypic characteristics of the adult cardiomyocyte. Proc. Natl. Acad. Sci. U.S.A. 95, 2979–29849501201  
Doncheva N. T.,Assenov Y.,Domingues F. S.,Albrecht M.. (Year: 2012). Topological analysis and interactive visualization of biological networks and protein structures. Nat. Protoc. 7, 670–68510.1038/nprot.2012.00422422314  
Dupont E.,Hanganu I. L.,Kilb W.,Hirsch S.,Luhmann H. J.. (Year: 2006). Rapid developmental switch in the mechanisms driving early cortical columnar networks. Nature439, 79–8310.1038/nature0426416327778  
Everett K. L.,Cooper D. M.. (Year: 2013). An improved targeted cAMP sensor to study the regulation of Adenylyl Cyclase 8 by Ca(2+) entry through voltagegated channels. PLoS ONE8:e7594210.1371/journal.pone.007594224086669  
Feldt S.,Bonifazi P.,Cossart R.. (Year: 2011). Dissecting functional connectivity of neuronal microcircuits: experimental and theoretical insights. Trends Neurosci. 34, 225–23610.1016/j.tins.2011.02.00721459463  
Guimera R.,Mossa S.,Turtschi A.,Amaral L. A.. (Year: 2005). The worldwide air transportation network: anomalous centrality, community structure, and cities' global roles. Proc. Natl. Acad. Sci. U.S.A. 102, 7794–779910.1073/pnas.040799410215911778  
Honey C. J.,Sporns O.. (Year: 2008). Dynamical consequences of lesions in cortical networks. Hum. Brain Mapp. 29, 802–80910.1002/hbm.2057918438885  
Humphries M. D.,Gurney K.. (Year: 2008). Network ‘smallworldness’: a quantitative method for determining canonical network equivalence. PLoS ONE3:e000205110.1371/journal.pone.000205118446219  
Jeong H.,Mason S. P.,Barabasi A. L.,Oltvai Z. N.. (Year: 2001). Lethality and centrality in protein networks. Nature411, 41–4210.1038/3507513811333967  
Kandler K.,Gillespie D. C.. (Year: 2005). Developmental refinement of inhibitory soundlocalization circuits. Trends Neurosci. 28, 290–29610.1016/j.tins.2005.04.00715927684  
Kandler K.,Katz L. C.. (Year: 1998). Coordination of neuronal activity in developing visual cortex by gap junctionmediated biochemical communication. J. Neurosci. 18, 1419–14279454851  
Katz L. C.,Shatz C. J.. (Year: 1996). Synaptic activity and the construction of cortical circuits. Science274, 1133–11388895456  
Khazipov R.,Sirota A.,Leinekugel X.,Holmes G. L.,BenAri Y.,Buzsaki G.. (Year: 2004). Early motor activity drives spindle bursts in the developing somatosensory cortex. Nature432, 758–76110.1038/nature0313215592414  
Kitambi S. S.,Nilsson E. S.,Sekyrova P.,Ibarra C.,Tekeoh G. N.,Andang M.,et al. (Year: 2012). Small molecule screening platform for assessment of cardiovascular toxicity on adult zebrafish heart. BMC Physiol. 12:310.1186/1472679312322449203  
Malmersjö S.,Liste I.,Dyachok O.,Tengholm A.,Arenas E.,Uhlén P.. (Year: 2010). Ca2+ and cAMP signaling in human embryonic stem cellderived dopamine neurons. Stem Cells Dev. 19, 1355–136410.1089/scd.2009.043620043754  
Malmersjö S.,Rebellato P.,Smedler E.,Planert H.,Kanatani S.,Liste I.,et al. (Year: 2013a). Neural progenitors organize in smallworld networks to promote cell proliferation. Proc. Natl. Acad. Sci. U.S.A. 110, E1524–E153210.1073/pnas.122017911023576737  
Malmersjö S.,Rebellato P.,Smedler E.,Uhlén P.. (Year: 2013b). Smallworld networks of spontaneous Ca(2+) activity. Commun. Integr. Biol. 6:e2478810.4161/cib.2478823986813  
Mehlhorn K.,Naher S.. (Year: 1995). Leda  a platform for combinatorial and geometric computing. Commun. ACM38, 96–102  
Mueller L. A. J.,Kugler K. G.,Dander A.,Graber A.,Dehmer M.. (Year: 2011). QuACN: an R package for analyzing complex biological networks quantitatively. Bioinformatics27, 140–14110.1093/bioinformatics/btq60621075747  
Nelson D. E.,Ihekwaba A. E.,Elliott M.,Johnson J. R.,Gibney C. A.,Foreman B. E.,et al. (Year: 2004). Oscillations in NFkappaB signaling control the dynamics of gene expression. Science306, 704–70810.1126/science.109996215499023  
Newman M.. (Year: 2010). Networks: An Introduction. Oxford: Oxford University Press  
Nicol X.,Voyatzis S.,Muzerelle A.,NarbouxNeme N.,Sudhof T. C.,Miles R.,et al. (Year: 2007). cAMP oscillations and retinal activity are permissive for ephrin signaling during the establishment of the retinotopic map. Nat. Neurosci. 10, 340–34710.1038/nn184217259982  
Owens D. F.,Kriegstein A. R.. (Year: 1998). Patterns of intracellular calcium fluctuation in precursor cells of the neocortical ventricular zone. J. Neurosci. 18, 5374–53889651220  
Brockwell P. J.,Davis R. A.. (Year: 1998). Time Series: Theory and Methods. New York, NY: Springer  
Seth A. K.. (Year: 2010). A MATLAB toolbox for Granger causal connectivity analysis. J. Neurosci. Methods186, 262–27310.1016/j.jneumeth.2009.11.02019961876  
Sporns O.. (Year: 2014). Contributions and challenges for network models in cognitive neuroscience. Nat. Neurosci. 17, 652–66010.1038/nn.369024686784  
Sporns O.,Tononi G.,Kotter R.. (Year: 2005). The human connectome: a structural description of the human brain. PLoS Comput. Biol. 1:e4210.1371/journal.pcbi.001004216201007  
Vicente R.,Wibral M.,Lindner M.,Pipa G.. (Year: 2011). Transfer entropy–a modelfree measure of effective connectivity for the neurosciences. J. Comput. Neurosci. 30, 45–6710.1007/s108270100262320706781  
Watts D. J.,Strogatz S. H.. (Year: 1998). Collective dynamics of ‘smallworld’ networks. Nature393, 440–4429623998  
Yuste R.,Peinado A.,Katz L. C.. (Year: 1992). Neuronal domains in developing neocortex. Science257, 665–6691496379 
Figures
[Figure ID: F1] 
Figure 1
Principle of crosscorrelation analysis of a cell signaling network. (A) Illustration of a group of ten cells, numbered with roman numerals, with links connecting strongly correlated cells. The shortest path length between cell i and cell x, as well as the clustering coefficient for cell v, is stated in the figure. The shortest path between cell i and cell x is 2. Cell v has three connected neighbors, and they are all neighbors of each other; thus, the clustering coefficient is 3/3 = 1. (B) Pseudo Ca^{2+} concentration vs. time traces for all ten cells. Two time series with high (C) or low (D) correlation coefficients. 
[Figure ID: F2] 
Figure 2
Correlation as a function of lag. (A) Two sine functions with the same frequency but different amplitudes and phases, plotted in the same graph. (B) The correlation as a function of lag of the two sine waves in (A). 
[Figure ID: F3] 
Figure 3
Producing a scrambled signal. (A) A pseudosignal with start and end time points t_{0} and t_{n}. (B) A scrambled signal was created by shuffling the signal in (A) to a random time point t_{x}. 
[Figure ID: F4] 
Figure 4
Characterizing network structures. (A) The topologies of a smallworld (red), regular (blue), and random (green) network depicted with their internal relations to the network parameters σ and λ. (B) The network parameters σ and λ for a smallworld (red), regular (blue), and random (green) network, plotted as a function of rewiring. The degree distributions for a random (C) and scalefree (D) network plotted in a log–log scale. (E) The topology for a scalefree network. 
[Figure ID: F5] 
Figure 5
Cartoon illustrating the basic steps of the method. (A) A timelapse cell recording is imported into imageprocessing software that calculates the mean intensity and x,ycoordinates for each ROI and saves all the values in a tabdelimited text file. Time points are indicated by t_{1}, t_{2}, … t_{n}. (B) The file is loaded into MATLAB, which generates a matrix in which the individual cells' time series are divided into columns and the time points are divided into rows. The ten cells are indicated by roman numerals. The entire network structure is plotted (C) from the correlation matrix (D). Colored lines and matrix elements indicate the strength of correlation between cell pairs. Significant correlations were plotted (E) by applying a cutoff to the correlation matrix (F). 
[Figure ID: F6] 
Figure 6
Time series can produce false correlations. (A) Two separate time series consisting of pseudorandom numbers normally distributed. (B) The crosscorrelation calculation for the time series in (A). The mean correlation is −0.0011 and the Pearson correlation is 0.2863. Illustration of falsepositive cell pairs with high correlations due to small numbers of peaks (C) or sustained plateau increases (D). 
[Figure ID: F7] 
Figure 7
Network analyses performed on timelapse Ca^{2+} recordings of neural progenitor cells. (A) Network with a cutoff of 0.39 for neural progenitor cells plotted on top of a microscopic image. (B) Correlation as a function of distance for the experiment shown in (A). Red line indicates the cutoff level. Arrows indicate cell pairs with high correlations and short intermodal distances. (C) Distance distribution of data in (B) with correlations above cutoff (>Cutoff), shown in red, and below cutoff (<Cutoff), shown in black. Note the higher frequency of shorter distances for correlations above the cutoff. (D) Network plot of data from (A) scrambled in time domain, as described in Figure 3. (E) Correlation as a function of distance for the data shown in (D). (F) Degree distribution of the analysis on neural progenitor cells showed in (A). Scale bars, 50 μm. 
[Figure ID: F8] 
Figure 8
Network analyses performed on two separated cell ensembles in cardiac HL1 cells. (A) Network of atrial HL1 cells with a cutoff of 0.23, plotted on top of a microscope image. A cut divided the cells into two populations, as indicated in the figure. (B) Correlation as a function of distance for the experiment shown in (A). (C) Lag distribution for the maximum correlations above cutoff (>Cutoff), shown in red, and below cutoff (<Cutoff), shown in black. Note the small peaks for values below cutoff. (D) Network plot of data from (A) with xycoordinates of all cells randomized. (E) Correlation as a function of distance for the data shown in (D). (F) Degree distribution of the analysis of the HL1 cells shown in (A). Scale bars, 50 μm. 
[Figure ID: F9] 
Figure 9
Network properties are dependent on the cutoff value. (A) Increasing the cutoff for an experiment on neural progenitor cells decreased both the edge density and connectivity. (B) Increasing the cutoff decreased the shortest path length and increased clustering and the smallworld parameter. (C) The smallworld parameter does not exhibit a linear dependency on the edge density. Data were fit to a polynomial function of the first order inverse. 
Article Categories:
Keywords: networks, graph theory, neural networks, cross correlation, smallworld, scalefree, calcium signaling. 
Previous Document: Simultaneous visualization of extrinsic and intrinsic axon collaterals in Golgilike detail for mous...
Next Document: Optical suppression of drugevoked phasic dopamine release.