Calciumnetexplorer: an R package for network analysis of calcium imaging data

Background

Analyzing calcium imaging data to understand complex functional networks can be challenging, often requiring multiple tools, custom scripts, and some coding expertise. To address these challenges, we present CalciumNetExploreR (CNER), an R package designed to streamline and standardize the analysis of time-series data from neuronal populations.

Results

CNER integrates essential steps-normalization, binarization, population activity visualization, network construction, degree distribution analysis, principal component analysis, power spectral density evaluation, and event frequency calculations-into a single, cohesive pipeline. This comprehensive approach enables users to efficiently extract and compare network metrics, including clustering coefficients, global efficiency, community structures, and principal component variances. By offering a flexible and customizable framework, CNER simplifies the examination of functional connectivity and network topology, effectively providing the means to characterize a cellular functional network or analogous structures in other modalities.

Conclusion

Designed as a user-friendly package, CNER allows both experimental and computational neuroscientists to incorporate robust statistical and graphical analyses into their workflows without extensive coding knowledge. By unifying key analytical components into one pipeline, CNER reduces barriers associated with large-scale data analyses, ultimately facilitating deeper insights into the functional organization and dynamic properties of neuronal networks across diverse recording techniques.

Calcium live imaging is a pivotal technique in neuroscience and cellular biology for monitoring intracellular calcium dynamics, which are essential to numerous physiological processes such as neurotransmission, muscle contraction, and signal transduction pathways [1, 2]. By utilizing fluorescent calcium indicators, researchers can visualize and quantify calcium fluctuations within living cells and tissues, providing real-time insights into cellular activity and intercellular communication networks.

The analysis of calcium imaging data poses significant challenges due to the large volume and complexity of the datasets. After calcium traces are extracted from imaging data researchers face the task of interpreting these signals to uncover underlying biological phenomena. This requires sophisticated tools for data preprocessing, statistical analysis, network construction, and visualization to extract meaningful insights.

While several software tools focus on the initial extraction of calcium traces from raw imaging data, such as Suite2p [3], CaImAn [4], and Fiji [5] -which can be extended with plugins like TurboReg and Template Matching for motion correction and signal extraction- there is a relative scarcity of tools dedicated to the downstream analysis of these extracted time-series data. Many researchers resort to custom scripts or general-purpose software, which may not be tailored to the specific needs of calcium imaging data analysis. There is a lack of software solutions that enable graph-modeling analyses, such as measuring network clustering, efficiency, and other key topological properties, which are crucial for understanding the dynamic interactions within cellular networks. Existing tools for time-series analysis, such as MATLAB’s Signal Processing Toolbox [6] and Python libraries like NumPy and SciPy [7, 8], offer general-purpose functions but require extensive programming and lack specialized functions. Moreover, they do not provide integrated workflows for network analysis tailored to calcium imaging. Packages like the Brain Connectivity Toolbox [9] and BRAPH [10] offer a wide array of network analysis methods but are not specifically designed for calcium imaging data and often necessitate complex data preparation and manipulation. BRAPH (Brain Analysis using Graph Theory) is a MATLAB-based software package designed for the analysis of brain connectivity using graph theory [10]. While BRAPH offers extensive functionalities for brain connectivity analysis, it is primarily geared towards processing neuroimaging data such as MRI and fMRI. Applying BRAPH to calcium imaging time-series data may not be straightforward due to differences in data structure and may require significant adaptation or preprocessing steps. Moreover, since BRAPH is built on MATLAB, it requires access to proprietary software, which may not be readily available to all researchers due to licensing costs. This can limit accessibility, especially in educational or resource-constrained settings.

To address these gaps, we introduce CalciumNetExploreR (CNER), an R package designed to streamline the analysis of calcium imaging time-series data with a particular focus on understanding cellular networks and dynamic behaviors in biological systems. CNER is developed in R, an open-source programming language and environment widely used for statistical computing and graphics. By being free to use, it lowers the barrier to entry and promotes wider adoption across the research community. Additionally, R has a vast ecosystem of packages and a strong community support, facilitating collaboration and extension of functionalities.

CNER operates on pre-extracted calcium traces and offers a comprehensive suite of tools for data preprocessing, network analysis, and visualization. CNER offers several key advantages:

• Preprocessing: Accepts time-series (calcium fluctuations) and cell/ROI coordinates, with options for normalization and binarization.

• Network analysis: Builds correlation-based networks and extrapolates key features.

• Subset analysis: Isolates user-defined cell groups and measures their interactions with the main network.

• User-Friendly workflow: Features an all-in-one pipeline requiring minimal coding.

• Modular Design: Allows for customization and integration with other R packages.

By providing a cohesive platform specifically designed for the downstream analysis of calcium imaging time-series data, CalciumNetExploreR fills a critical niche in the current landscape of analytical tools. Unlike other software that may require extensive programming, lack specialized functions for calcium imaging data, or do not integrate network analysis capabilities tailored to such data, CNER enables advanced exploratory data analysis with just a few functions. Researchers can isolate a subset of neurons with CNER to compare their activity and interactions against the larger network, offering insights into specialized functions, cell-type-specific behaviors, and how local changes affect global dynamics. This is especially useful when certain cells are labeled differently (e.g., with a distinct fluorescent reporter) or exposed to targeted stimuli, making it easier to study disease, development, or learning processes.

In this paper, we present the design, functionalities, and applications of CNER. We highlight its unique capabilities in subset network analysis and the measurement of interactions between main networks and subsets. Through application on sample datasets, we demonstrate how CNER facilitates advanced analysis and visualization of calcium imaging time-series data, ultimately contributing to a better understanding of complex biological systems.

Data preprocessing

Normalization

To account for potential heterogeneity in baseline levels of the Ca²⁺ indicator (e.g. GCaMP) expression among cells, we normalize the raw calcium signals for each cell individually. This ensures comparability across the cell population by scaling each cell’s signal to a common range and across experiments. However, while this method is straightforward and enhances interpretability, it is important to be aware it might be susceptible to the influence of outliers, which may bias the scaled signal. Nevertheless, it is included as a default in the package as it provides a useful normalization baseline, especially when comparing across cells that may express different levels of calcium indicators such as GCaMP. The normalization is performed using min-max scaling:

where

• \( x'_t \) is the normalized calcium signal at time point \( t \),

• \( x_t \) is the raw calcium signal at time point \( t \),

• \( \min (x) \) and \( \max (x) \) are the minimum and maximum calcium signals of the cell over the entire recording period.

This transformation scales the calcium signals of each cell to a range between 0 and 1.

Binarization

The binarize() function in CalciumNetExploreR is flexible and allows users to define their own criteria for distinguishing between active and inactive states of cells. Users can specify arbitrary thresholds through custom functions to suit their specific experimental conditions. In the context of this study, we binarized the normalized calcium signals based on a threshold derived from the cell’s signal variability. A cell is considered active at time point \( t \) if its normalized signal exceeds two standard deviations above the mean:

where

• \( x'_{t,\text {bin}} \) is the binarized signal,

• \( x'_t \) is the normalized signal at time \( t \),

• \( \mu _{x'} \) is the mean of the normalized signal \( x' \),

• \( \sigma _{x'} \) is the standard deviation of \( x' \).

This approach identifies significant calcium events relative to the cell’s baseline activity, allowing for the detection of spontaneous or evoked calcium transients. Binarization can be particularly useful in low signal-to-noise conditions, as it suppresses minor fluctuations and emphasizes meaningful activity patterns. As noted in previous studies [11], digitization of signals is often necessary to reduce noise and improve downstream interpretability, particularly when performing correlation-based network analyses.

Population activity analysis

Active cells percentage

We calculate the proportion of active cells over time by summing the binarized signals across all cells at each time point and dividing by the total number of cells:

where:

• \( N \) is the total number of cells,

• \( x'_{t,\text {bin},i} \) is the binarized signal of cell \( i \) at time point \( t \).

This metric provides insight into the overall network activity and can reveal synchronized activity patterns.

Hierarchical clustering

To detect potential clusters of calcium activity patterns among cells, we perform hierarchical clustering directly on the normalized calcium signals \( x' \) using Ward’s method [12]. Nevertheless, by prioritizing compact and homogeneous clusters, the default Ward’s method may bias results toward detecting tightly-knit groups, potentially overlooking broader or less dense patterns of coordination. Thus, users are advised to choose the best clustering method based on their data; can alternatively select from a range of linkage methods (e.g., "complete", "average") to explore different structural assumptions and sensitivities in clustering outcomes.

Network creation and analysis

Functional connectivity

Functional connectivity between cells is assessed using pairwise cross-correlation of the normalized calcium signals. In practice, CNER uses the R function ccf() present in the stats package to compute the cross-correlation coefficients. The ccf() function calculates the cross-correlation between two time series using the formulae as described in [13]: let \( x'_i(t) \) and \( x'_j(t) \) denote the normalized calcium signals of cell \( i \) and cell \( j \) at time \( t \). The cross-correlation coefficient \( \rho _{ij}(\tau ) \) between these two time series at lag \( \tau \) is calculated as:

where

• \( K_{ij}(\tau ) \) is the cross-covariance between \( x'_i(t+\tau ) \) and \( x'_j(t) \);

• \( \sigma _i \) and \( \sigma _j \) are the standard deviations of \( x'_i(t) \) and \( x'_j(t) \), respectively.

The cross-covariance \( K_{ij}(\tau ) \) is defined as:

where

• \( \mu _i \) and \( \mu _j \) are the mean values of \( x'_i(t) \) and \( x'_j(t) \), respectively (the time series are modeled as stochastic processes);

• \( E \) is the expectation value operator. As time in this case is fixed, E is the average value from multiple samples.

To account for short delays in cellular responses, cross-correlation coefficients are computed for user-defined lags \( \tau \). By default, we consider lags \( \tau = -1, 0, +1 \) time points, but the user can specify different lag values. The maximum absolute value of \( \rho _{ij}(\tau ) \) across these lags is selected as the functional connectivity measure between cells:

where

• \( \rho _{ij}(\tau ) \) is the cross-correlation coefficient between cell \( i \) and cell \( j \) at lag \( \tau \),

• \( \text {Lag} \) is the set of lags specified by the user.

Network construction

The functional connectivity matrix \( \rho _{ij,\text {max}} \) is used to construct an undirected weighted graph \( G = (V,E,W) \), where

• \( V \) is the set of nodes representing cells,

• \( E \) is the set of edges representing significant functional connections,

• \( W \) is the set of edge weights given by \( \rho _{ij,\text {max}} \).

A threshold \( \theta \) is applied to \( \rho _{ij,\text {max}} \) to determine significant edges:

Edges with weights below the threshold are excluded to focus on the most significant connections.

Graph metrics

We compute several key topological properties of the network to characterize its structure:

• Degree (\(k_i\)): The degree of node i is the number of edges connected to it:

  $$\begin{aligned} k_i = \sum _{j} a_{ij} \end{aligned}$$

  where \( a_{ij} = 1 \) if there is an edge between nodes \( i \) and \( j \), and \( 0 \) otherwise.

• Clustering Coefficient (\( C(g) \)): Measures the tendency of the entire network to form tightly knit clusters. It is defined as the average clustering coefficient across all nodes in the network:

  $$\begin{aligned} C(g) = \frac{1}{N} \sum _{i=1}^{N} \frac{2 e_i}{k_i (k_i - 1)} \end{aligned}$$

  where \( e_i \) is the number of edges between the neighbors of node \( i \), \( k_i \) is the degree of node \( i \), and \( N \) is the total number of nodes in the network.

• Global Efficiency (\( G(g) \)): Reflects how efficiently information is exchanged over the entire network. It is calculated as:

  $$\begin{aligned} G(g) = \frac{1}{N(N - 1)} \sum _{i \ne j} \frac{1}{d_{ij}} \end{aligned}$$

  where \( d_{ij} \) is the shortest path length between nodes \( i \) and \( j \), and \( N \) is the total number of nodes in the network.

• Mean degree: The average degree across all nodes:

  $$\begin{aligned} \langle k \rangle = \frac{1}{N} \sum _{i=1}^{N} k_i \end{aligned}$$

• Community detection: Communities within the network are identified using the leading eigenvector method [14], which calculates the leading non-negative eigenvector of the modularity matrix of the graph. This method partitions the network into clusters of densely connected nodes by maximizing the modularity.

These metrics provide insights into the network’s topology, revealing patterns of connectivity and potential functional modules within the cellular network.

Subset network analysis

The package allows for the selection and analysis of subsets of cells within the global network. Subsets can be defined based on various criteria, such as anatomical location, functional characteristics, or genetic markers.

For each subset, a functional network is constructed using the same methods described above. The topological properties of the subset network are calculated and compared to those of the main network. Interactions between the subset and the main network are assessed by analyzing the connections between cells in the subset and those outside it.

This analysis provides valuable information about the role of specific cell groups within the overall network and how they influence or are influenced by the broader cellular community.

Principal component analysis (PCA)

CNER performs PCA on the normalized calcium signals using the prcomp function from the stats R package. This process reduces dimensionality and identifies principal components that capture the most variance in the data.

Scree plots are generated to visualize the variance explained by each principal component, aiding in the selection of components for further analysis.

Power spectral density (PSD) analysis

We conduct PSD analysis to examine the frequency content of the calcium signals, identifying dominant frequencies and oscillatory patterns. We use the Welch method [15] to estimate the PSD for each cell’s normalized signal.

The PSD plots display power as a function of frequency, allowing for the identification of periodic behaviors and synchronization phenomena within the cellular population.

Integrated analysis pipeline

An integrated analysis pipeline is provided through the pipeline() function in the CalciumNetExploreR package. This function automates the sequential execution of the analysis steps with minimal user intervention. The pipeline includes:

1. 1.

   Data loading: Reads the time-series data and cell coordinates provided by the user.

2. 2.

   Normalization and binarization: Applies the normalization and binarization procedures as described above.

3. 3.

   Population activity visualization: Generates plots of the proportion of active cells over time, with options for hierarchical clustering to detect activity patterns.

4. 4.

   Network construction and visualization: Builds the functional network using cross-correlation coefficients and visualizes it using graph plotting functions. Node positions are determined by the provided cell coordinates.

5. 5.

   Graph metrics calculation: Computes key topological properties of the network, including degree distribution, clustering coefficient, global efficiency, and identifies communities using the leading eigenvector method.

6. 6.

   Subset analysis: Allows for the selection of cell subsets for specialized analysis, repeating the network construction and metric calculations for the subset.

7. 7.

   PCA and PSD analyses: Performs PCA and PSD analyses on the normalized signals, providing visualizations such as scree plots and PSD graphs.

8. 8.

   Reporting: Compiles the results into summary tables and generates comprehensive visualizations for interpretation (see Figs. 4, 5, 8 and 9 generated using CNER).

The pipeline function enhances reproducibility and efficiency by encapsulating the analysis workflow into a single command, which can be customized with function arguments to suit specific research needs. The function is fully deterministic and does not rely on any random initialization, ensuring stable results across repeated executions without the need to set a random seed.

Software and implementation

The CalciumNetExploreR package is developed using R (version 4.4.0) and can be found, with its documentation, on https://github.com/simo-91/CalNetExploreR.

Data visualization

Images in Figs. 1 and 7 were generated using ImageJ/Fiji, while plots in Figs. 2 and 3 were created using the ggplot2 package in R. Figures 4, 5, 8, and 9 were produced using CNER visualization functions (included in the pipeline() function).

Data requirements

The package requires:

• Time-series data: A matrix or data frame where rows represent cells and columns represent time points, containing the raw or pre-processed calcium signals. Mathematically, this can be represented as \( M_{n \times t} \), where \( n \) is the number of cells (rows) and \( t \) is the number of time points (columns). Each element \( M_{ij} \) represents the calcium signal of cell \( i \) at time point \( j \). It is essential that this matrix be properly formatted so that each row corresponds to a unique cell and each column to a time point in the recording.

• Cell coordinates: A matrix or data frame containing the spatial coordinates of each cell (e.g., X and Y positions), necessary for network visualization and spatial analyses. This is represented as a matrix \( C_{n \times 2} \), where \( n \) is the number of cells and each row contains the \( (x_i, y_i) \) coordinates of cell \( i \). Proper formatting ensures that each cell’s coordinates are aligned with the rows of the time-series matrix, so that spatial information corresponds directly to the calcium data.

• Labels: An additional column in the cell coordinates data frame that specifies the subset classification of each cell. This column, named Label, contains identifiers (a categorical label) indicating whether a cell belongs to a specific subset.

Experimental design

Zebrafish maintenance

Zebrafish (Danio rerio) were maintained at \(28.5^{\circ }\hbox {C}\) on a 14-hour light/10-hour dark cycle in the aquatics facility at BVS Aquatics department based at the Queen Medical Research Institute (University of Edinburgh). Embryos were obtained by natural spawning from adult transgenic Tg(NBT:H2B-GCaMP6s; nacre^-/-) and raised in E3 embryo medium at \(28.5^{\circ }\hbox {C}\). To inhibit pigmentation in pigmented strains, embryos were treated with 200\(\mu \)M 1-phenyl-2-thiourea (PTU) (Sigma) starting at 24 h post-fertilization until the end of the experiment. All animal procedures were approved by the Institutional Animal Care and Use Committee (IACUC) and conducted in accordance with the Animal (Scientific Procedures) Act 1986 and institutional guidelines.

Zebrafish oocyte injections and sample generation

We utilized a transgenic zebrafish line Tg(NBT:H2B-GCaMP6s; nacre^-/-), in which all neurons express the nuclear calcium indicator GCaMP6s, enabling the visualization of neuronal activity in vivo. To achieve mosaic overexpression of the oncogene in the brain, embryos at the single-cell stage were coinjected with a pDEST-NBT:\(\Delta \)lexPR-lexOP-pA driver plasmid and either a Tol2-pDEST-lexOP:AKT1-lexOP:tagRFP plasmid for the AKT1 condition (TREAT) or a Tol2-pDEST-lexOP:tagRFP-pA control plasmid for the control condition (CTRL). Injection solutions contained 40 ng/uL of plasmid DNA, 0.2% Phenol Red dye and Tol2 capped mRNA (20 ng/\(\mu \)L) to enable transient mosaic expression. Briefly, the transcriptional activator \(\Delta \)lexPR is expressed in neurons and activates gene downstream of lexOP sites. Hence, this approach resulted in a subset of neurons overexpressing AKT1 and being labeled with red fluorescent protein (RFP) for the TREAT condition and neurons simply being labeled with tagRFP for the CTRL condition.

To test the pipeline on neurons of the contralateral tectum, we analyzed a sample from a zebrafish line expressing the nuclear-localized calcium indicator H2B-GCaMP6s under the control of the panneuronal promoter elavl3 (Tg(elavl3:H2B-GCaMP6s)). Imaging was performed on a live, awake larva at 4 days post-fertilization (4 dpf).

Calcium live imaging

Six days post-fertilization zebrafish larvae were immobilized by treatment with 0.5 mg/mL Mivacurium Chloride (Abcam, ab143667) to inhibit movement. The larvae were then embedded in 1.5% low melting point agarose (Life Technologies) within glass-bottom dishes (MatTek) filled with E3 medium. Time-lapse calcium imaging was performed using a Zeiss LSM880 confocal microscope equipped with a 20x water-dipping objective lens, focusing on the hindbrain area. Time-lapses lasted for 15 min, capturing frames every 2 s (resulting in a 0.5 Hz acquisition rate). Fluorescence excitation was achieved using 488 nm laser lines for GCaMP6s and 543 nm for tagRFP. To image neurons of the contralateral tectum, a 4dpf Tg(elavl3:H2B-GCaMP6s) larva was embedded in low-melting-point agarose, and the gel was carefully cut to free one of the eyes. Imaging was conducted on the contralateral optic tectum, at a single focal plane located approximately 70 \(\mu \)m deep, for a duration of 5 min at a sampling rate of 2 Hz.

Statistical analysis

Statistical analyses are conducted using built-in R functions and validated methods. When comparing network metrics between groups or conditions, appropriate statistical tests (e.g., t-tests, ANOVA) are applied, and p-values are adjusted for multiple comparisons where necessary.

Explorative and comparative analysis of network features using CNER

To demonstrate the capabilities of CalciumNetExploreR (CNER), we explored Ca²⁺ imaging data from exemplary datasets generated using larval zebrafish. We utilized a transgenic zebrafish line NBT:H2B-GCaMP6s; nacre^-/-, in which all neurons express the nuclear calcium indicator GCaMP6s, enabling the visualization of neuronal activity in vivo. To perform a comparative analysis using CNER, we generated two experimental datasets: 1) A control dataset (CTRL) in which in addition to the GCaMP6s signal a subset of neurons expressed RFP, and a treatment dataset (TREAT) in which a subset of neurons expressed RFP and the human oncogene AKT1 (see Methods for details) [16, 17]. As AKT1 overexpression has been shown to influence Ca²⁺ levels before, we deemed this dataset suitable for a comparative analysis using CNER. Datasets consisted of calcium fluorescence time-series data from hundreds of neurons recorded over T = 450 time points (Fig. 1). Finally, we applied the pipeline() function integrated in CNER to the control (CTRL) and the treatment (TREAT) datasets in order to compare their network features. The features analyzed included:

• Mean event frequency of calcium events per minute for the entire population (Freq.)

• Mean event frequency per minute for the labeled subset of neurons (Freq. labeled), which are the neurons overexpressing AKT1 or carrying the control plasmid.

• Clustering coefficient (\(C_{g}\))

• Global efficiency (\(G_{g}\))

• Percentage of variance explained by the top five principal components (Top5PC Var (%))

• Proportion of connections from labeled-to-unlabeled neurons over total possible labeled-to-unlabeled connections (LtU)

The pipeline() function is a comprehensive analysis tool that automates the processing of calcium imaging data through several key stages, integrating all of the following functions that allow for customization and detailed analysis. It essentially acts as a wrapper for several functions present in the package, namely normalize(), binarize(), population_activity(), make_network(), plot_network(), degrees(), pca(), PSD.plt(), get_top5pc_variance(), subset_connections(), events_per_min(), and also igraph functions transitivity() and global_efficiency().

The results are summarized in Table 1.

Firstly, the datasets were normalized using the normalize() function, integrated in the pipeline, which scales the raw calcium signals across cells to ensure comparability by centering and scaling to unit variance.

Next, we converted the normalized, continuous signals into binary states (active/inactive) using the binarize() function to facilitate the analysis of neuronal activity patterns.

Pipeline() function unravels latent features and differences between treatment and control networks

One of the outputs of the pipeline() function is the mean event frequency per minute for the entire population of neurons (Freq.). According to the analyzed samples, CTRL samples had a Freq. ranging between 2.75 and 3.00 events/min, while TREAT samples showed lower values (around 2.15\(-\)2.25 events/min). This difference was statistically significant (p <0.01, Fig.2A). Another output of the pipeline is the mean event frequency per minute for the labeled subset of neurons (Freq. labeled). The TREAT samples showed a significantly higher mean event frequency (\(4.00 \pm 0.10\)) compared to CTRL samples (\(2.89 \pm 0.09\)), representing an increase of approximately 39% (\(p < 0.001\), independent t-test, Fig. 2B). This suggests that the treatment enhances neuronal excitability or firing rate, which is reflected in the increased frequency of calcium events detected.

The subset_connections() function integrated within the pipeline() quantified the proportion of connections between specified subsets of neurons. We found that the proportion of connections from labeled to unlabeled neurons (LtU) was higher in CTRL samples (\(0.01207 \pm 0.00143\)) compared to TREAT samples (\(0.00600 \pm 0.00040\)), indicating a two-fold decrease in the TREAT condition (\(p < 0.05\), Mann–Whitney U test, Fig. 2C). This result suggests that neurons in the treatment condition are less connected to the surrounding neuronal network, potentially leading to functional isolation or altered network integration.

To construct functional connectivity graphs, the pipeline uses the make_network() function, which performs cross-correlation analysis on the binarized calcium signals to infer connections between neurons. Customization options in this function allowed us to adjust parameters such as lag.max (numeric, default 1) specifying the maximum lag for cross-correlation, and correlation_threshold (numeric or "none", default 0.3) to filter edges by weight; setting it to "none" disables filtering. This flexibility enabled us to tailor the sensitivity of the network construction based on correlation strength and lag. For this use case, we used a correlation threshold of 0.6 and a lag of 1.

Visualization of the networks was done using the integrated plot_network() function, which displays the network structure and highlights community memberships by labeling nodes accordingly. Customization options included coordinates (a data frame with columns "X", "Y", "Cell", and "Label"), cell_ID (character, options "cell" or "communities") to label nodes by cell ID or community, and label (logical, default FALSE) to include labels on the plot.

In the TREAT condition, the network exhibited reduced clustering and global efficiency compared to the CTRL condition. The mean clustering coefficient (\(C_{g}\)) in TREAT samples was \(0.547 \pm 0.011\), significantly lower than in CTRL samples (\(0.772 \pm 0.010\)), representing a 41% decrease (\(p = 0.0009\), independent t-test, Fig. 3A). Global efficiency (\(G_{g}\)), calculated using built-in functions from the igraph package, was reduced by 33% in TREAT samples (\(0.00927 \pm 0.000750\)) compared to CTRL samples (\(0.0138 \pm 0.00125\)) (\(p = 0.004\), independent t-test, Fig. 3B). These metrics indicate that the treatment leads to a less interconnected and less efficient network.

The network graphs generated by the plot_network() function, integrated in the pipeline, visually demonstrate these differences (Fig. 4). The TREAT network (Fig. 4A) shows reduced clustering with only one large community (defined as a community with more than five members), whereas the CTRL network (Fig. 4B) displays increased clustering with 14 large communities. Nodes are color- and number-coded by community membership, illustrating the fragmentation in the TREAT condition.

We visualized population activity using the population_activity.plt() function, which generates raster plots of neuronal activity over time, sorted through hierarchical clustering to reveal patterns of co-activation. Customization options included the dendrogram argument (logical, default FALSE) to include a dendrogram in the plot. In the TREAT condition, the plots showed increased variability and less organized activity patterns compared to CTRL (Fig. 5), suggesting a lack of functional clusters and disrupted synchronization among neurons in the treatment group.

We assessed network complexity using the pca() function also integrated within the pipeline(). This function performs principal component analysis to identify the principal components that capture the most variance in the data, thereby reducing dimensionality and highlighting significant patterns of neuronal activity. The percentage of variance explained by the top five principal components (Top5PC Var (%)) was significantly lower in TREAT samples (\(8.80\% \pm 0.35\%\)) compared to CTRL samples (\(12.38\% \pm 1.02\%\)), indicating a reduction of approximately 29% (\(p < 0.01\), independent t-test). This suggests that neuronal activity in the TREAT condition is more complex and less dominated by a few principal patterns.

The PCA scree plots generated by the pca() function illustrate this difference in network complexity (Fig. 6). In the TREAT condition (Fig. 6A), the variance explained by each successive principal component declines more gradually, indicating that a larger number of components are needed to capture the same amount of variance as in the CTRL condition (Fig. 6B).

By utilizing the pipeline() wrapper function of CalciumNetExploreR, we effectively analyzed and visualized differences in neuronal network features between the TREAT and CTRL conditions. Significant alterations associated with the treatment were identified, demonstrating the utility of CalciumNetExploreR in uncovering functional network changes in calcium imaging data.

Pipeline() function tested on neurons of the contralateral tectum

To demonstrate the generalizability of the pipeline() function, we run it to an independent dataset obtained from an elavl3:H2B-GCaMP6s zebrafish larva at 4 days post-fertilization (4 dpf). A representative frame from the calcium imaging timelapse is shown in Fig. 7. The extracted time-series data were processed using the pipeline() function with default parameters.

The global population activity plot (Fig. 8) reveals a consistent level of neuronal activity throughout the acquisition window, with an average frequency of calcium events estimated at 6.69 events per minute. Network construction based on cross-correlation and subsequent community detection yielded a functional connectivity graph with visible modular organization (Fig. 9).

The topological analysis of the resulting network reported a clustering coefficient of 0.497, suggesting a moderate degree of local interconnectedness among neurons. The global efficiency was measured at 0.0105, indicating relatively sparse long-range communication within the network. Dimensionality reduction via PCA showed that the top five principal components accounted for 22.84% of the total variance, reflecting structured but diverse activity patterns.

These results confirm the capacity of the pipeline() function to extract interpretable and biologically relevant metrics from independent datasets, supporting its applicability across different experimental conditions.

In this study, we introduced CalciumNetExploreR (CNER), a comprehensive R package designed for the analysis of calcium imaging data and the exploration of neuronal network features. By applying CNER to calcium imaging datasets from zebrafish larvae under AKT1 overexpression (TREAT) and control (CTRL) conditions, we demonstrated the package’s capability to extract meaningful network metrics and visualize complex neuronal interactions. The use of the pipeline() function allowed for a simple and standardized analysis workflow, encompassing data normalization, binarization, network construction, feature extraction, and visualization.

By analyzing the TREAT and CTRL datasets, we were able to characterize significant differences in network topology and function using CNER. Specifically, the package revealed that in the TREAT condition there is a substantial decrease in the clustering coefficient (\(C_{g}\)) and global efficiency (\(G_{g}\)). The mean \(C_{g}\) was reduced by approximately 41% in the TREAT condition compared to CTRL, indicating a disruption in local connectivity and network clustering. Similarly, the global efficiency decreased by about 33%, suggesting impaired overall network integration and efficiency in information transfer. The network graphs generated by CNER (Fig. 4) visually illustrated the reduced clustering and altered community structures in the TREAT networks compared to CTRL. Additionally, the population activity plots provided insights into the temporal dynamics of neuronal activity, revealing increased noise and variability in the TREAT condition networks (Fig. 5). The application of principal component analysis (PCA) within the package offered insights into the complexity of neuronal activity. By analyzing the variance explained by principal components, we found that the top five principal components accounted for a lower percentage of variance in the TREAT condition compared to CTRL. This suggests increased complexity and heterogeneity in neuronal firing patterns in the treatment group, which could be further investigated using CNER’s tools. Furthermore, CNER facilitated the analysis of subpopulations within the neuronal network. By examining the connections from labeled (RFP⁺) to unlabeled neurons and calculating the mean event frequencies, we observed that neurons of the treatment condition had elevated activity levels and altered connectivity patterns. Specifically, the mean event frequency per minute for labeled neurons was significantly higher in the TREAT condition, and the proportion of connections from labeled to unlabeled neurons was decreased. These analyses demonstrate how CNER can be used to explore the functional roles of specific neuronal subpopulations within the broader network. Importantly, these findings highlight how CNER can uncover alterations in latent network properties that may be associated with experimental manipulations or disease models.

Graph theory, on which CNER is based, is a well-established and powerful framework for modeling and analyzing cellular population activity in neuroscience. Its continued relevance is underscored by recent advancements in its application to calcium imaging data [18,19,20,21,22]. Several elegant studies have demonstrated the effectiveness of graph theoretical approaches in elucidating network structures and dynamics [11, 23,24,25,26].

Finally, the adaptability and extensibility of CNER allow users to incorporate additional analyses or customize functions to suit their specific research needs. The package’s reliance on widely used R packages such as igraph, ggplot2, and ggraph ensures compatibility and ease of integration with custom analytical workflows.

In summary, our application of CalciumNetExploreR to the TREAT and CTRL datasets demonstrates how the package can be utilized to characterize and compare neuronal networks under different experimental conditions. By providing comprehensive tools for network analysis and visualization, CNER enables researchers to uncover subtle changes in network topology and function, contributing valuable insights into the mechanisms underlying neuronal dynamics.

Limitations and future directions

Despite its capabilities, CalciumNetExploreR has limitations that warrant consideration. The package currently assumes a basic level of proficiency in R programming and familiarity with network analysis concepts, which may present a barrier for some users. A graphical user interface (GUI) is in the works, as it could improve accessibility for a broader range of researchers.

Additionally, while the package offers several network metrics and visualization options, the incorporation of a baseline correction in pre-processing and advanced machine learning tools could further enrich the analyses. Future developments may focus on expanding the package’s functionality to include these features, thereby providing a more comprehensive tool for neuroscientific research.

CalciumNetExploreR is a powerful and versatile R package that facilitates the analysis of calcium imaging data and the exploration of neuronal network properties. By streamlining complex analytical processes and providing intuitive visualization tools, it enables researchers to uncover meaningful insights into neuronal dynamics and network architecture. Our application of CalciumNetExploreR to AKT1 overexpression data illustrates its utility in identifying significant differences in network features between experimental conditions. We believe that CalciumNetExploreR will be a valuable resource for the neuroscience community, promoting advanced network analyses and contributing to a deeper understanding of neuronal function.

The datasets used and analysed during the current study are available from the corresponding author on reasonable request. The source code for CalciumNetExploreR is available on GitHub at: https://github.com/simo-91/CalNetExploreR. The repository includes installation instructions, usage examples, and documentation.

The authors thank the Bioresearch and Veterinary Services (BVS) Aquatics facility (QMRI, University of Edinburgh) for maintenance and care of the zebrafish. The authors are grateful to Dr. Daniel Soong and the UK Zebrafish Imaging and Screening Facility for assistance with microscope imaging. Special thanks to Dr. Lian Hollander-Cohen and Prof. Tim Czopka (Centre for Clinical Brain Sciences, University of Edinburgh)for providing the contralateral tectal neurons timelapse, and Dr. Nicola Romanò for critical reading of this article.

This work was supported by the Medical Research Council with a Precision Medicine PhD fellowship to SM (MR/N013166/1, MR/W006804/1) and by Cancer Research UK with Programme Foundation Award to DS (DRCPFA-May21\(\backslash \)100011).

Authors and Affiliations

1. Centre for Discovery Brain Sciences, University of Edinburgh, 49 Little France Crescent, Edinburgh, EH16 4SB, UK

   Simone Lenci & Dirk Sieger

Contributions

S.L. and D.S. conceptualized the study. S.L.developed the software, conducted the data analysis, and prepared the manuscript. D.S.provided project supervision, manuscript feedback, and editorial support.

Corresponding authors

Correspondence to Simone Lenci or Dirk Sieger.

Ethics approval and consent to participate

Animal experiments were reviewed and approved by the ethical review committee of the University of Edinburgh and conducted under Home Office project license PP8119722, in accordance with the Animals (Scientific Procedures) Act 1986.

Consent for publication

Not applicable.

Competing interests

The authors have no Conflict of interest to declare.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions