Unraveling the neural code: analysis of large-scale two-photon microscopy data

Abstract

The brain is an intricate neuronal network that orchestrates our thoughts, emotions and actions through dynamic interactions between neurons. If we could record the activity of all neurons simultaneously in detail, it could revolutionize our understanding of brain function and lead to breakthroughs in treating neurological diseases. Recent technological innovations, particularly in large field-of-view two-photon microscopes, have made it possible to record the activity of tens of thousands of neurons simultaneously. However, the size and complexity of the datasets present significant challenges in extracting interpretable information. Conventional analysis methods are often insufficient, necessitating the development of new theoretical frameworks and computational efficiencies. In this review, we describe the characteristics of the data obtained from advanced imaging techniques and discuss analytical methods to facilitate mutual understanding between experimentalists and theorists. This interdisciplinary approach is crucial for effectively managing and interpreting large-scale neural activity datasets, ultimately advancing our understanding of brain function.

large field-of-view two-photon microscopy, neural activity, data processing, single neuron, neural population, neuronal network

Introduction

Neuroscience aims to elucidate the underlying mechanisms of brain functions such as perception, action, decision-making, memory and emotion, among others. Unveiling the intricate neural underpinnings of these functions requires high-fidelity measurements of neural activity. Although the term neural activity is singular, it encompasses a wide range of spatiotemporal and qualitative resolutions, each of which is determined by the recording method. Therefore, the choice of recording method is a critical first step in neuroscience.

Neuroscientists primarily employ two broad categories of recording techniques: electrophysiology and optical imaging [1]. Recent technological advances in both electrophysiology and imaging have enabled large-scale neural recordings across different spatiotemporal scales, offering unprecedented insights into the neural basis of brain function. Herein, we aim to describe the principles of these two recording methods and discuss their respective strengths and limitations. We then aim to examine how the datasets generated by optical imaging differ from those generated by electrophysiology and how these differences affect data analysis and interpretation.

Notably, two-photon microscopy with a large field-of-view (FOV) has emerged as a powerful tool, allowing the simultaneous observation of >10 000 neurons in multiple cortical areas [2–7]. It has opened new possibilities to elucidate the spatial distribution of information coding [8,9], the dynamics of neural populations [10–12], interregional interaction [13,14] and the organization of cortical-wide neuronal networks [5]. However, the vast amount of data generated by these advanced imaging techniques presents significant challenges. Effective data preprocessing is essential to manage noise, correct motion artifacts and extract meaningful signals from raw data. Proper preprocessing steps ensure the accuracy and reliability of subsequent analyses, enabling neuroscientists to draw valid conclusions from complex data sets. In this review, we aim to address the critical role of preprocessing and discuss best practices for optimizing data quality for robust analysis of optical imaging. Furthermore, we aim to explore methods for analyzing neural metrics at various levels: single-neuron metrics, population metrics and network metrics. Comprehensively covering these topics, this review aims to provide an overview of the current landscape in neural activity recording and analysis, supporting future research and innovation in neuroscience.

Conventional in vivo neural recording methods

The most prominent physiological property of neurons is their sub-millisecond electrical activity called action potentials, also referred to as spikes. Therefore, the straightforward approach to recording neural activity is to monitor the electrical activity of neurons directly. In electrophysiology, researchers use electrodes to record intracellular or extracellular electrical signals at various resolutions. Extracellular recordings are widely used to monitor neural activity from many neurons and/or regions [15,16]. For example, electrocorticography (ECoG) is performed using low-impedance electrodes placed on the cortical surface (Fig. 1a). In general, ECoG records field potentials that reflect the population activity of multiple neurons [17]. The array of ECoG electrodes allows the simultaneous monitoring of different cortical regions [18,19]. To probe the activity of individual neurons within these populations, high-impedance multichannel electrodes such as tetrodes or high-density silicon probes are used (Fig. 1b) [20–22]. These electrodes can simultaneously record field potentials and spiking activity from multiple neurons. The activity of individual neurons can be isolated through a spike-sorting process, which involves detecting spikes, extracting their waveforms and clustering similar waveforms [23,24]. Recent advances in high-density probes, such as Neuropixels, enable the recording of ∼100 neurons from a single brain region [25–27].

Fig. 1.

Basic neural recording methods. (a) ECoG grid array records field potentials across multiple cortical regions. (b) High-density silicon probe recording enables isolating the spiking activity of individual neurons through spike sorting. (c) Widefield one-photon imaging records the population of Ca²⁺ dynamics from the entire dorsal cortex. (d) Conventional two-photon imaging records the Ca²⁺ dynamics of individual neurons in a specific brain region.

While electrophysiology provides direct measurements of electrical activity with high temporal resolution, optical imaging methods are also powerful tools for recording neural activity as fluorescent changes [28,29]. In these methods, genetically encoded Ca²⁺ indicators (GECIs) such as GCaMPs are used to monitor neuronal Ca²⁺ dynamics associated with spiking activity [30–35]. Optical imaging, similar to electrophysiology, can be performed at different spatiotemporal scales depending on the combination of the objective and optical system of the microscope. The magnification of the lenses determines the size of the FOV, and the numerical aperture (NA) of the lenses primarily determines the optical resolution (higher NA has higher resolution). The magnification of a lens system is determined by the ratio of the focal length of the objective lens to the focal length of the tube lens. Low-magnification lenses are designed to provide a wide FOV, which requires a longer focal length of objective lenses. However, this involves a trade-off, as the longer focal length limits the light collection angle, resulting in a lower NA and reduced optical resolution. For instance, a one-photon microscope equipped with a low magnification and low NA objective (e.g. 1×, 0.25 NA) is commonly used to record cortical-wide Ca²⁺ dynamics from neuron population (Fig. 1c) [36,37]. In this recording method, a complementary metal-oxide semiconductor sensor has recently been used as a photodetector to collect fluorescence signals, allowing high-speed recording of cortical-wide activity. However, one-photon excitation illuminates the entire tissue, even outside the focal plane. Thus, a photodetector collects contaminated fluorescence signals, making it difficult to detect individual neuron activity, even if the optical resolution is sufficient to identify the neuron. While this scattering problem is solved by using confocal microscopes, it is difficult to illuminate deeper into the tissue. To overcome these limitations, a two-photon microscope equipped with higher magnification and higher NA objective lens (e.g. 16×, 0.8 NA) is employed for high-resolution imaging of individual neuron activity (Fig. 1d) [38,39]. This technique uses pulsed infrared lasers to excite fluorescent indicators in a highly localized manner, allowing for imaging deep (∼1 mm) within the brain tissue. Conventional two-photon microscopy can simultaneously record the activity of tens to hundreds of neurons from ∼500 μm × 500 μm at 10–30 Hz, providing detailed insights into the dynamics of individual neurons and their interactions within a given cortical region [40].

Sampling bias and data properties in electrophysiology and Ca²⁺ imaging

Extracellular electrophysiology provides high temporal resolution, capturing spikes with sub-millisecond precision from both cortical and subcortical regions. However, the data tend to overrepresent neurons with higher spiking rates because spike sorting algorithms require a sufficient number of spikes for accurate clustering [41]. As a result, neurons with lower spiking rates may not be reliably detected or classified, leading to their underrepresentation in the dataset.

Meanwhile, Ca²⁺ imaging provides high spatial resolution, allowing visualization of larger neural populations compared to electrophysiology. This method captures an indirect measure of neural activity, as Ca²⁺ signals reflect the cumulative effect of multiple spikes over time. The temporal resolution of Ca²⁺ imaging is limited by the kinetics of the Ca²⁺ indicators and the frame rate, which may make the activity of neurons with higher spiking rates less accurately represented in the data. Moreover, the Ca²⁺ events are mostly associated with multiple spikes in a short period rather than single spikes [42,43]. Despite these limitations, Ca²⁺ imaging has several advantages, including the ability to (i) track the same neurons over several days, (ii) record neural activity from subcellular compartments such as dendrites, axons and spines, (iii) complement these measurements by using sensors specifically designed to detect other modalities of neural activity, such as neurotransmitter or neuromodulator dynamics, (iv) record from genetically defined specific neural subpopulations by combining transgenic lines that specifically express GECIs and (v) combine optogenetic manipulations with single-cell resolution.

Moreover, recent developments in GECIs, such as the jGCaMP8 series, which have faster kinetics and high sensitivity to detect single spikes, offer the potential to infer the underlying true spiking activity more accurately [32]. Advances in microscopy and biological sensors are addressing current limitations and moving toward the future of optical imaging with a sufficient temporal resolution to capture single-spike dynamics.

Development of a large FOV two-photon microscope

Ca²⁺ imaging with a two-photon microscope is a powerful tool for elucidating the neural dynamics underlying brain function. However, the FOV of conventional two-photon microscopes is limited to recording from a single cortical region. In order to gain a comprehensive understanding of the interactions between cortical regions, it is essential to record neural activity across multiple regions.

To develop a large FOV two-photon microscope, it is necessary to overcome the following trade-offs [7]: (i) Trade-off between FOV and optical resolution – while high-resolution imaging requires a higher NA, the NA of a conventional low-magnification objective for larger FOV imaging is not sufficient to resolve individual neurons. Low magnification and high NA can be achieved by increasing the size of the lens system, but designing a lens that has both low magnification and high NA while suppressing optical aberrations has been a practical challenge. (ii) Trade-off between FOV and frame rate – this trade-off arises from the laser scanning mirror system commonly used in two-photon microscopes. Mirrors scan an excitation laser spot across the entire FOV to capture fluorescence. Therefore, a larger FOV requires the laser spot to travel a greater distance, increasing the scan time and subsequently decreasing the frame rate. (iii) Trade-off between FOV and signal-to-noise ratio (SNR) – as the FOV increases, the duration of laser excitation at each pixel (called dwell time) decreases when maintaining the frame rate. This reduction in dwell time decreases the probability of two-photon excitation and the number of fluorescence photons per pixel, resulting in lower SNR in large FOV imaging.

To address these trade-offs, we developed a customized two-photon microscope called fast scanning high optical invariant two-photon microscopy (FASHIO-2PM) [5]. Our design incorporates several key innovations to overcome these challenges: (i) We designed a large objective lens along with scanning and tube lenses to increase NA. While increasing the pupil diameter of the objective lens increases NA, minimizing optical aberrations throughout the system is critical because optical aberrations reduce resolution and two-photon excitation efficiency, resulting in a decrease in SNR. The large objective lens we developed has an excitation NA of 0.4, sufficient for resolving single neurons, with near diffraction-limited (theoretical limit) performance across the entire 3 mm × 3 mm FOV. (ii) As discussed earlier, the frame rate depends on the laser scanning mirror system. In particular, resonant mirrors are widely used due to their ability to achieve higher scanning speeds compared to traditional galvanometric mirrors. Basically, there are three choices of resonant mirrors for two-photon microscopes (resonant frequency: 4k, 8k and 12k). We adopted the 8k resonant mirror in our system because it is well balanced in terms of scanning angle, scanning speed and SNR. (iii) We developed a pre-chirper and a large aperture GaAsP photomultiplier tube (PMT) sensor to improve the SNR. The pre-chirp corrects the stretching degradation of the laser pulse to improve the two-photon excitation efficiency over the entire FOV, and the GaAsP PMT sensor improves the sensitivity of signal detection. In addition, the NA of the objective for collection is 0.8 enabling efficient collection of fluorescence signals. By combining these technologies, FASHIO-2PM can scan the 3 mm × 3 mm FOV at ∼7.5 Hz for 2048 × 2048 pixels or ∼15 Hz for 1024 × 1024 pixels while maintaining single-neuron resolution with high SNR (Fig. 2a). Other research groups have also achieved large FOV two-photon imaging using different techniques, resulting in an exponential increase in the number of neurons that can be imaged simultaneously from multiple cortical regions [3,4,6].

Fig. 2.

Large FOV two-photon Ca²⁺ imaging and data processing. (a) Schematics of FASIO-2PM recording. This method allows the simultaneous recording of Ca²⁺ dynamics of individual neurons from multiple cortical regions. (b) Data processing workflow. (c) Motion correction. Standard deviation projection image of first 500 frames before motion correction (top). Standard deviation projection of the first 500 frames after motion correction (bottom). (d) ROIs detected with LCCD. (e) Signal processing procedures. Neuropil signal is subtracted from soma raw signal (top). After correcting for baseline variation, ΔF/F₀ is calculated and spikes were estimated by deconvolution (bottom). (f) Examples of subsequent analysis methods.

Preprocessing two-photon microscopy data for accurate analysis

These large-scale imaging techniques allow us to monitor cortical neural dynamics on an unprecedented scale. However, biological data contain a variety of observational noises, such as shot noise, brain motion and sensor expression heterogeneity, and improper preprocessing can greatly affect subsequent analysis and lead to fatally erroneous conclusions. In this section, we describe how to perform preprocessing such as motion correction, neuron detection, signal extraction and spike inference on data acquired by two-photon microscopy (Fig. 2b).

Motion correction

Managing motion artifacts is a significant challenge in two-photon microscopy, especially for the imaging of awake-behaving animals. These artifacts are caused by involuntary movements of the animal, such as breathing or locomotion. Even small movements can cause significant misalignment of fluorescence signals across multiple frames, resulting in correlated noise in the data [44]. This noise can obscure true neural activity and prevent accurate analysis. Several motion correction algorithms have been developed to address this problem [44–48]. These algorithms typically work by creating a reference image, often derived from the initial frames, and then registering each subsequent frame against this reference through rigid or nonrigid transformations. In the context of large FOV imaging, motion correction becomes even more critical due to the increased likelihood of regional variations in motion. These regional variations may be due to slight differences in tissue properties or the subject’s specific motion patterns. In this case, an approach that corrects motion in a segmented manner shows exceptional suitability for large FOV imaging [47,48]. By dividing the large image into smaller, overlapping patches and aligning the smaller patches independently, motion correction is more accurate than algorithms that treat the entire image as a single unit (Fig. 2c).

Region of interest detection

The next step is to identify regions of interest (ROIs) for each neuron. Manually setting ROIs for many neurons in a large FOV two-photon imaging dataset is a time-consuming and impractical task. We developed the low computational-cost cell detection (LCCD) method to address this challenge [49,50]. This algorithm provides an efficient and accurate solution for automatic cell detection from large-scale imaging data, enabling ROI identification in tens of minutes. LCCD works by first dividing the entire Ca²⁺ imaging data into smaller blocks along the time axis. Within each block, LCCD captures the morphology of active neurons from projection images. This approach allows separate detection of spatially adjacent neurons. For example, even if two neurons, A and B, are close to each other, LCCD can detect only A in a specific block when A is active but not B. This block-by-block detection more effectively separates surrounding neurons than using a single projection image that merges all activity over time. Finally, LCCD integrates the information from each block to create the final set of ROIs over the entire data (Fig. 2d). Other groups have also developed other ROI detection pipelines, such as constrained nonnegative matrix factorization (CNMF) and Suite2p [51,52]. These approaches decompose the imaging data matrix into spatial components, background, noise and temporal components. In the CNMF, nonnegative matrix factorization is used to decompose the data, and the Ca²⁺ dynamics are incorporated into an autoregressive model as a temporal constraint. Suite2p incorporates the model of the neuropil (a structure of a mixture of dendrites and axons around the soma) into the background component. This is because the neuropil signal is spatially diffuse and contaminates the background noise as large correlated fluctuations [53,54]. Importantly, it is critical for researchers to compare the ROI detection results of different methods on their own datasets to find the most suitable approach. The effectiveness of these methods depends heavily on the specific data characteristics and parameters chosen.

Signal extraction and normalization

Once neurons are detected, the fluorescence changes of the identified cells are extracted (Fig. 2e, top). To prevent the signal from being contaminated by the neuropil, the signal for a few pixels around the cell body or the neuropil component extracted by matrix factorization is subtracted from the soma signal (Fig. 2e, top). Although the Ca²⁺ dynamics are extracted in this way, they cannot be directly compared because of variations in brightness across neurons. This variation arises from differences in the baseline expression level of GECIs within each neuron. Therefore, to enable comparisons, the change in fluorescence value (ΔF) is normalized to the baseline (F₀) to calculate ΔF/F₀ (Fig. 2e, bottom) [55]. It is better to compute the baseline successively in a moving window rather than averaging over the entire period to compensate for baseline fluctuations due to slow baseline changes such as fading. All of this processing removes various noise components that do not depend on neural activity. The resulting ΔF/F₀ dynamics can then be used for analysis (Fig. 2f).

Spike inference

While the ΔF/F₀ dynamics provides a valuable measure of neural activity, it is an indirect reflection of the underlying spiking events. Since the Ca²⁺ signal is assumed to be a convolution of the actual spiking activity with the GECI kernel function, it should be possible to deconvolve the Ca²⁺ signal and recover the original spiking events (Fig. 2e, bottom). Deconvolution allows for higher temporal precision in neural activity analysis and more accurate estimation of correlations between neurons. Several studies have developed algorithms for deconvolution, such as MLspike, Online Active Set method to Infer Spikes and FastLZeroSpikeInference [56–58]. Although these algorithms estimate the underlying spiking activity under ideal conditions, the accuracy of spike inference depends on the datasets. The slow kinetics of GECIs and the limitations of frame rate restrict the temporal resolution of spike inference, leading to missed spikes for high-spiking rate neurons. Additionally, low SNR recordings can reduce the accuracy of spike detection, resulting in overestimation of spikes. Since experimental conditions such as microscope settings (laser power, sensor gain and frame rate), family of GECIs, expression method of GECIs (virus vector or transgenic line), recording region and layer vary from experiment to experiment, finding a universally optimal set of parameters for these algorithms is challenging. Therefore, it is recommended to use generative models to simulate Ca²⁺ dynamics from spike events that reproduce similar dynamics to real datasets and to explore optimal parameters for spike inference on your datasets practically [43,56,59].

Recently, deep learning tools have shown promise in offering more accurate and noise-resistant inference of spiking activity. One approach is to use a self-supervised deep learning method, such as DeepCAD, to improve the SNR of the acquired images prior to ROI detection [60,61]. DeepCAD can significantly reduce acquisition noise, resulting in high-SNR Ca2+ signals and improved spike inference. Other approaches, like CASCADE and ENS² (effective and efficient neural networks for spike inference from calcium signals), utilize deep neural networks trained on datasets where spike data are recorded simultaneously with Ca²⁺ imaging as ground truth [62,63]. These deep neural networks learn the relationship between ground truth spikes and Ca²⁺ signals, leading to significantly improved inference compared to conventional deconvolution methods. Alongside advancements in biosensors and microscopy, the development of these computational methods is paving the way for enhanced accuracy and reliability in Ca²⁺ imaging.

Quantification of single neuron encoding properties

Understanding how the brain processes information and encodes external and internal states heavily necessitates a characterization of the individual and interactive effects of external and internal factors on single neurons. Previous studies have demonstrated that even a single neuron or a small number of neurons are sufficient to alter the brain states [64] and drive sensory-guided decision-making [65,66]. These findings indicate that individual neurons can have a prominent impact on the entire network of the brain and influence behavior. Hence, investigating the activity patterns and their correlation with behavior or task-related information in individual neurons is still an important and intriguing question [67].

One of the key methodologies for analyzing single-neuron activity is event-triggered analysis, which examines neural responses triggered by specific behavioral/experimental events such as stimulus or movement. Researchers can analyze how neurons respond to events dynamically and compare these responses across different conditions by aligning neural activity to these events. This alignment of the neural activity to the event is often represented by peri-event time histograms (PETHs), which visualize the average neural responses around the event, providing insights into how these events influence neural activity (Fig. 3a). Additionally, spike-triggered averaging can reveal the influence of specific information on the timing of neural responses. The event-triggered analysis can also help determine the selectivity of each neuron. Quantifying neural selectivity is fundamental in neuroscience research, as it provides insights into how neurons encode information under varying conditions over time. To determine the selectivity of the specific information, analytical methods such as tuning curves [68,69], information theory metrics (e.g. mutual information) [70], D-prime (d′) [71], receiver operating characteristic analysis [71,72] and entropy-based measures [73] are commonly used. These methods collectively help elucidate how neurons selectively encode and discriminate information within neural circuits.

Fig. 3.

Decompose single-neuron information encoding. (a) PETHs of neural activity aligned to stimulus onset (left) and movement onset (right), showing average responses around these events. (b) Example neural activity responding to stimulus and movement information. The observed neural response can be split into stimulus and movement-related components. (c) Schematic of a GLM fitted to the neural response to the task and behavioral variables. (d) Example neural activity (top). Structure of predictors design matrix (bottom). The rows in the matrix represent stimulus and movement variables (predictors), which can have a 0 or 1 value for time points relative to the appropriate time onset from the stimulus and action, respectively.

A major challenge in analyzing single-neuron responses is statistically quantifying and separating the contributions of multiple variables encoded within individual neurons. Many neurons often exhibit multivariate response characteristics to multiple variables [74], making it difficult to precisely estimate the contributions of individual information based on changes in average neural activity. For example, if a single neuron encodes both sensory and motor information, its response during stimulus presentation will reflect both types of information if the subject performs a movement simultaneously (Fig. 3b). To resolve such confounding, it is necessary to apply models that mathematically characterize the relationship between neural activity patterns and sensory or motor information, thereby disentangling their respective contributions. This mathematical model can quantitatively describe the temporal dynamics of information encoding and the contributions of multiple information sources. Furthermore, explicitly hypothesizing the dynamics of neural activity through parameterized models can help extract important neural patterns and encoding schemes that might be difficult to characterize using simple statistical methods. The generalized linear model (GLM) is used for modeling statistical relationships between neural activity and experimental and behavioral parameters such as stimulus and movement (Fig. 3c) [75]. This versatile technique facilitates the understanding of how different types of information are encoded in neural activity and can reveal underlying patterns that contribute to our comprehension of brain function. The GLM assumes that the outputs of behavioral/experimental variables are linearly summed and transformed through nonlinear function (link function) to determine the predicted neuron’s response |$r$| as follows:

$$r = f\,\left( {{\beta _0} + {\beta _1}{X_1} + \ldots + {\beta _k}{X_k}\,} \right)$$

(1)

where |$f$| is the nonlinear activation function. In the neuroscience field, exponential and log functions are commonly used because neurons typically exhibit sparse firing patterns and their spike densities show strong nonlinearity, which can be effectively modeled by these functions. If a neuron’s response adheres to a specific probability distribution, such as a Poisson distribution, the predicted response can be transformed using the corresponding distribution function, denoted as |$r^{\prime}\sim Poisson\left( r \right)\,$|⁠.

The GLMs are generally used for modeling neural activity by the parameters within the epoch, but the models can also be expanded to predict time-varying responses. In this case, the neural response is described as a linear sum of temporal filters aligned to behavioral or experimental events. Suppose the subject is presented with a stimulus and reacts to the stimulus by licking the spout to get a reward. Stimulus and movement-related parameters can be aligned to stimulus presentation onset and movement onset (Fig. 3d, top). In the model, the predicted neural response |$r\left( t \right)$| is given as follows:

$$r\left( t \right) = f\,\left( {\mathop \sum \limits_{{t_s} \in \,{T_s}} {\beta _s}\left( {t - {t_s}} \right) + \mathop \sum \limits_{{t_m} \in \,{T_m}} {\beta _m}\left( {t - {t_m}} \right) + {\beta _0}} \right)$$

(2)

where |${\beta _s}\left( t \right)$| and |${\beta _m}\left( t \right)$| represent the weights of stimulus and movement variables at time t, |${T_s}$| and |${T_m}$| represent the set of times for which those variables appeared and |${\beta _0}$| is the intercept. The stimulus and movement weights |$\beta_s\left(t\right)$| and |$\beta_m\left(t\right)$| support the experimenter-defined window relative to the stimulus onset and movement onset, respectively.

Building the model helps statistically assess the time-varying response properties of a single neuron. When the models are parameterized with time-varying response properties as those described earlier, the success of the model depends strongly on the shape of the design matrix (Fig. 3d, bottom), which is the matrix that concatenates all parameters (variables) including the bias term. However, including too many parameters in a model may result in overfitting, especially for neurons with a low SNR and limited number of time samples. The number of parameters that can be safely incorporated is constrained by the sample size of the dataset. As the number of parameters increases, model complexity rises, which may lead to overestimation by incorporating noise as a signal in datasets with low SNR of neuronal response [76]. Various machine learning and statistical techniques are applied to overcome this problem, such as the application of a basis function [77,78], a reduced-rank approach [79,80], regularization [81–83] and cross-validation [84].

Overall, given the multivariate nature of neural activity, using encoding models like the GLMs helps distinguish among the different effects. Combining large-scale imaging with GLMs allows for precise analysis of how widespread neural populations represent different types of information and enables the examination of its spatial distribution. For example, Tseng et al. showed that task- and behavior-relevant information is represented with a specific gradient in widespread regions of the posterior cortex [9]. Such integration of large-scale imaging with encoding models unravels a comprehensive view of how neurons encode and process information, paving the way for deeper insights into neural function.

Uncovering the latent dynamics underlying neural population activity

Neurophysiology has shown that individual neurons have at least two major limitations in encoding information. First, a neuron is noisy and spikes probabilistically. Even if physically identical stimuli are presented, the activity of single neurons exhibits substantial variability in response and timing across repeated presentations. This variability can be due to synaptic noise [85], making consistent and reliable information encoding challenging. Second, the amount of information a single neuron can represent is limited by physiological constraints on spiking rate and temporal response patterns. These constraints limit the encoding capacity of individual neurons, creating a need for population-level averaging of information encoding. For nearly a century, researchers have studied how neurons behave collectively and how information is encoded with the population of neurons. This macroscopic concept of encoding information by the population of neurons has been shaped by numerous theoretical frameworks and experimental evidence, such as Hebb’s assembly hypothesis from the 1940s [86] and Georgopoulos’s study about population vectors in the 1980s [87]. The key premise is that a population, not a single neuron, is a fundamental unit of encoding information in the brain, which is known as the population doctrine [12,88]. In this framework, single-neuron activities are treated as reflecting latent factors that are information encoded by time-varying patterns that are not observed directly at the single-neuron level but appear by the collective and correlated activity underlying neural population activity. This dynamic evolution of population activity can be mathematically modeled by various approaches including a linear dynamical system (LDS) [89,90], Gaussian process models [91] and recurrent neural network (RNN) models [92]. While these approaches have shared and distinct characteristics, here we discuss neural population dynamics based on LDS, where the model infers the current state based on past observations. The LDS is the simplest framework to model dynamical systems due to the linear projection of low-dimensional latent variables with linear dynamics. An LDS is described by a state equation (latent dynamics)

$$z\left( {t + 1} \right)\, = Az\left( t \right) + Bu\left( t \right) + \varepsilon \left( t \right)$$

and observation equation (neural activity)

$$y\left( t \right)\, = Cz\left( t \right) + Du\left( t \right) + \xi \left( t \right)$$

where |${z_t} \in {\mathbb{R}^d}$| corresponds to the latent state vector at time |$t$| over d latent factors, |$A \in {\mathbb{R}^{d \times d}}{\rm{ }}$| is the state transition matrix, describing the dynamics of the latent state, |$u\left( t \right) \in {\mathbb{R}^m}$| is the |$m$| external inputs (e.g. stimulus) at time |$t$|⁠, |$B \in {\mathbb{R}^{d \times m}}$| is the input matrix, describing how the external input or other brain areas influence the latent state, |$y\left( t \right) \in {\mathbb{R}^n}$| is the observed neural activity of |$n$| neurons at time |$t$|⁠, |$C \in {\mathbb{R}^{n \times d}}$| is the observation matrix, mapping the latent state to the observed neural activity (Fig. 4a), |$D \in {\mathbb{R}^{n \times m}}$| is the feedforward matrix, describing the direct influence of external input on the observed neural activity and |$\varepsilon \left( t \right) \sim N\left( {0,Q} \right),\,\xi \left( t \right) \sim N\left( {0,R} \right)$| is the noise assumed to be Gaussian with zero mean and covariance |$Q$| and |$R$|⁠. Typically, |$z\left( t \right)$| is an abstract representation in a low-dimensional subspace (or manifold), which is a lower number than the dimensionality of observed neurons (Fig. 4a, latent dynamics). This low-dimensional subspace is found via many types of linear/nonlinear dimensionality reduction methods or machine learning-based techniques that we introduce later. These equations allow for the decomposition of neural population activity into latent factors that evolve according to linear dynamics through neural state space, constituting a neural population trajectory. A characteristic of the LDS-based approach is that it assumes the latent state at a given time point is a linear function of the previous state (i.e. |$Az\left( t \right)$|⁠) and inputs (i.e. |$Bu\left( t \right)$|⁠). Latent factors, in turn, influence the observed activity of individual neurons through linear mapping. Researchers can infer the underlying dynamics driving neural population responses, revealing how these responses encode information about external stimuli and internal states.

Fig. 4.

Neural population dynamics and observed neural activities. (a) The relationship between latent dynamics and observed neural activities. The latent dynamics |$z\left( t \right)$| are described by the LDS. The observed neural activities |$y\left( t \right)$| are mapped by the latent dynamics. (b) schematic depicting interaction between two brain regions through the communication subspace. In both Regions 1 and 2, observed activities of |$N$| neurons are reduced down to a low-dimensional latent variables |$z$|⁠, such that each latent variable is inferred via linear combination of |$N$| observed neurons. The population activity in Region 1 at each time point corresponds to the inferred latent variables, which can be shown in three (or more) dimensions (low-dimensional dynamics in Region 1). Region 2 is assumed to influence population activity in Region 1 via the communication subspace defined by linear combination of |${B_{2 \to 1}}z\left( t \right)$| (communication subspace in Region 2).

An important analytical aspect under the population doctrine is how to mathematically and statistically infer the latent state |$z\left( t \right)$| in the neural population. Dimensionality reduction methods are commonly used and considered suitable for analyzing neural population activity. These methods facilitate the identification of low-dimensional latent states or variables that capture features of interest in the data. Dimensionality reduction has been applied in various previous studies, revealing the selection and integration of sensory information in the prefrontal cortex [93], population dynamics during reaching behavior in the motor cortex [94], choice-specific dynamics in the parietal cortex [95], visually guided decision-making [96], representation of latent variables in the hippocampus [97], representation of cognitive maps in the hippocampus [98] and more. However, there are many dimensionality reduction techniques such as principal component analysis (PCA), factor analysis (FA), nonnegative matrix factorization (NMF), locally linear embedding (LLE) and isometric mapping (Isomap). Each method captures the statistical structure of the data differently, resulting in a potentially substantial impact on the interpretations. Here, we briefly introduce basic dimensionality reduction methods, but further discussion can be found in [10].

Principal component analysis (PCA)

PCA is the most commonly used method for reducing dimensionality because it is typically simple to apply and relatively easy to interpret. PCA captures the low-dimensional elements of data while preserving data covariance. It aims to identify an ordered set of orthogonal (independent) axes that capture the greatest variances of the data: the first axis (principal component (PC)) summarizes the major axis variation and the second PC captures the next largest. The data can then be projected onto each axis and traced as neural trajectories (⁠|$z\left( t \right)$|⁠) in the low-dimensional space. Since PCA captures greater data variance, data must be well preprocessed using methods such as trial-averaging and denoising. PCA is particularly useful for analyzing population-level activity, such as identifying dominant firing patterns across a large number of neurons during a specific task or behavior.

Factor analysis (FA)

FA is also commonly used in neuroscience to infer the low-dimensional space of the data. While the goals of applying PCA and FA are similar, FA seeks latent variables by separating the shared variance of neurons from the independent variance, which is considered spiking variability as mentioned in the previous section. Thus, FA is essentially PCA with an explicit noise model and is suitable for single-trial analysis [91]. FA is especially beneficial when analyzing neural data from individual trials, where trial-to-trial variability is of interest. It can help distinguish between coordinated neural activity that might represent task-related information and the inherent variability in individual neuron firing rates.

Nonnegative matrix factorization (NMF)

NMF is also a useful and tractable linear dimensionality reduction method for inferring the low-dimensional representation of the data. NMF operates under the constraint that the data are all nonnegative, aligning with the nature of neural activity. This nonnegative constraint results in additive, nonsubtractive combinations, ensuring that each latent factor contributes positively to the whole. Therefore, NMF is considered a method to seek a parts-based representation [99], making it well suited for finding cell assemblies.

Locally linear embedding (LLM)

PCA, FA and NMF assume that the latent state space or manifold is projected via linear mapping from the neural data. However, some studies revealed that latent variables encoded in the low-dimensional space (or manifold) are represented nonlinearly such as curved shapes [100]. Due to this nonlinearity, linear dimensionality reduction may overestimate the number of dimensions to encode a latent variable or fail to capture the complex structure such as a Swiss roll shape. LLE is a commonly used method to extract latent variables by mapping high-dimensional data points into a lower-dimensional space while preserving neighborhood embeddings [101]. By locally and linearly reconstructing the data, LLE captures the nonlinear global structure in the high-dimensional data while preserving local structure. Therefore, LLE is particularly useful for analyzing low-dimensional dynamics or the structure of population activity that follows nonlinear trajectories in state space. It can uncover the underlying structure of neural dynamics that might be missed by linear dimensionality reduction methods.

Isometric mapping (Isomap)

Isomap is another nonlinear dimensionality reduction method widely used in neuroscience [97]. Like the LLE, Isomap preserves the local structure of the high-dimensional data and maps it into low-dimensional space. However, Isomap goes further, preserving the global structure of the data more effectively than LLE, making it suitable for data with a global geometrical shape. Therefore, Isomap is particularly effective when dealing with neural data that might form continuous, curved manifolds in the high-dimensional neural space.

Other approaches

Many previous studies have adopted the LDS framework, extending critical insight into characteristics of the computational framework in the brain [93,102–105]. Recently, various methods have been developed to expand LDS. These include Latent Factor Analysis via Dynamical Systems (LFADS) [106], which combines RNN and LDS concepts; low-rank RNN [107], which learns low-dimensional dynamics directly; and Gaussian process method [91], which smooths temporal dynamics more effectively compared to LDS. Notably, due to the constraints of Ca²⁺ imaging data such as the slow kinematics of GECIs, frame rate and insufficient spike inference, some statistical methods that specifically adapt to the data have developed. For example, the Calcium Imaging Linear Dynamical System addresses these issues by performing deconvolution of the data and dimensionality reduction simultaneously [108]. Another example is the Recurrent Autoencoder for Discovering Imaged Calcium Latents (RADICaL) [109]. RADICaL extended LFADS to infer latent dynamics from Ca²⁺ imaging data. It uses a zero-inflated gamma model to handle the statistical properties of the deconvolved Ca²⁺ signals, and a novel training strategy that uses precise neuron sampling times within each image frame. These methods enhance the recovery of high-frequency features and accurately infer latent dynamics from Ca²⁺ imaging data. As such, the development of these statistical methods continues as imaging techniques advance.

Revealing communication subspace underlying interregional interaction

The advancement of large FOV two-photon imaging has the potential to reveal inter-regional interaction dynamics. Recent findings have shown that widespread cortical areas encode external information and engage in the task, which is often called distributed coding [80,110]. These findings suggest a shift from the traditional view of functional specialization of each brain region to a more holistic perspective of the brain as an integrated system. Additionally, Pinto et al. showed that while widespread cortical regions are causally engaged in the task, these regions become decorrelated if the task demand increases [111], suggesting that different dynamics across regions complement each other by processing information in distinct ways. Since LDS models do not account for such inter-regional influences, incorporating these interactions into the model is crucial to understanding the brain dynamics better (Fig. 4b). Shenoy and Kao proposed simple model of the two regions with interaction terms [90]:

$${z_1}\left( {t + 1} \right)\, = {A_1}{z_1}\left( t \right) + {B_1}{u_1}\left( t \right) + {B_{2 \to 1}}{z_2}\left( t \right)$$

$${z_2}\left( {t + 1} \right)\, = {A_2}{z_2}\left( t \right) + {B_2}{u_2}\left( t \right) + {B_{1 \to 2}}{z_1}\left( t \right)$$

where |${A_1}\,and\,{A_2}$| represent transition matrices for regions 1 and 2, respectively, |${z_1} \in {\mathbb{R}^{{d_1}}}$| and |${z_2} \in {\mathbb{R}^{{d_2}}}$| correspond to the latent state vectors at time |$t$| over |${d_1}$| and |${d_2}$| latent factors in regions 1 and 2 and |${B_{2 \to 1}} \in {\mathbb{R}^{{d_2} \times {d_1}}}$| and |${B_{1 \to 2}} \in {\mathbb{R}^{{d_1} \times {d_2}}}$| represent the inter-regional input matrix, describing how the other brain areas influence the latent state. In particular, |${B_{2 \to 1}}{z_2}$| and |${B_{1 \to 2}}{z_1}$| can be thought of as a population activity of low-dimensional communication subspace that defines which activity patterns are effectively relayed between regions. This subspace is the low-dimensional space to be calculated by the fluctuation between brain regions, and there are several statistical methods to infer the subspace such as reduced rank regression (RRR), canonical correlation analysis (CCA) and partial least squares (PLS). Further discussion of the statistical methods for communication between regions can be found in Semedo et al. [112].

Reduced rank regression (RRR)

RRR is an extended multivariate linear regression model with the dimensionality reduction approach. This model can be extended to infer the communication subspace by using neural activities in one region to predict activities in another region. Assuming that |${y_1} \in {\mathbb{R}^{p \times {n_1}}}$| and |${y_2} \in {\mathbb{R}^{p \times {n_2}}}$| are the neural activities in regions 1 and 2, where |${n_1}$| and |${n_2}$| are the number of neurons and |$p$| is the number of observations (e.g. trial-by-trial activities). RRR model is defined as |${y_1}\sim{B_{2 \to 1}}{y_2}$| and has constraints |${\rm{rank}}\left( {{B_{2 \to 1}}} \right) \lt {\rm{min}}\left( {{n_1},{n_2}} \right)$|⁠. RRR can be solved using singular value decomposition:

$${B_{2 \to 1}} = BV{V^T}$$

where |${B_{2 \to 1}} = {({y_2}^T{y_2})^{ - 1}}{y_2}^T{y_1}$| is the least square solution of the coefficients. The columns of the |${n_2} \times m$| matrix |$V$| contain the top |$m$| principal components. First |$m$| components of |${B_{2 \to 1}}{y_2}$| thus correspond to the communication subspace from Region 2 to Region 1 [13].

Canonical correlation analysis (CCA)

CCA is a powerful and versatile multivariate tool used to investigate relationships among multiple brain regions jointly. CCA facilitates the understanding of inter-regional interactions and their impact on population dynamics by focusing on shared variance and providing a low-dimensional subspace. CCA is designed to maximize the correlation between neural activities between two brain regions |${y_1} \in {\mathbb{R}^{p \times {n_1}}}$| and |${y_2} \in {\mathbb{R}^{p \times {n_2}}}$|⁠, where |${n_1}$| and |${n_2}$| are the number of neurons and |${p_{}}$| is the number of observations (e.g. trial-by-trial activities). CCA aims to find linear combinations of |${y_1}$| and |${y_2}$| that maximize the correlation between these combinations. Mathematically, this is equivalent to solving the following equation:

$\begin{equation*}\max_{a,b}\;\text{corr}(y_1a,\;y_2b)\end{equation*}$

where |$a \in {\mathbb{R}^{{\rm{\,}}{n_1}}}$| and |$b \in {\mathbb{R}^{{\rm{\,}}{n_2}}}{\rm{\,}}$| are the vectors that map neural activity into the low-dimensional subspace. This is solved via eigenvalue decomposition:

$\begin{equation*}\textstyle\sum\nolimits_{y_1y_2}^{}\sum\nolimits_{y_2y_2}^{-1}\sum\nolimits_{y_2y_1}^{}a=\lambda\sum\nolimits_{y_1y_1}^{}a\end{equation*}$

$\begin{equation*}\textstyle\sum_{y_2y_1}^{}\sum_{y_1y_1}^{-1}\sum_{y_1y_2}^{}b=\lambda\sum_{y_2y_2}^{}b\end{equation*}$

where |$\mathop \sum$| is the covariance and the vectors |$a$| and |$b$| corresponding to the largest eigenvalue |$\lambda $| provide the linear combinations that maximize the correlation. Thus, |${B_{1 \to 2}} = {y_1}{a^T}$| and |${B_{2 \to 1}} = {y_2}{a^T}$| are the canonical variables that define the communication subspaces [113,114].

Partial least squares (PLS)

PLS is also a statistical method developed as an extension of the multivariate linear regression model [115,116]. Similar to CCA, PLS also identifies subspaces that are shared between brain regions. While CCA seeks subspaces that maximize the correlation between two sets of data, PLS is specifically designed to maximize the covariance between them. This characteristic makes PLS generally more robust to noise than CCA, as the covariance of the noise is expected to be close to zero. Let |${y_1} \in {\mathbb{R}^{p \times {n_1}}}$| and |${y_2} \in {\mathbb{R}^{p \times {n_2}}}$| represent two datasets, where |$p$| represents the number of observations and |${n_1}$|and |${n_2}$| represent the number of neurons. The general model for PLS can be as expressed as follows:

$\begin{equation*}y_1=TP^\top+E_1,\;y_2=UQ^\top+E_2\end{equation*}$

Here, |$T$| and |$U$| are the matrices of latent components (or scores), |$P$| and |$Q$| are the loading matrices and |${E_1}$| and |${E_2}$| are the residual matrices.

PLS iteratively finds the latent components |${T_i}$| that maximize the covariance between the linear combinations of |${y_1}$| and |${y_2}$|⁠, which are defined as follows:

$\begin{equation*}t_i=y_1w_i\end{equation*}$

and

$\begin{equation*}u_i=y_2c_i\end{equation*}$

where |${w_i}$| and |${c_i}$| are the weight vectors that define the direction of the |$i$|th latent components. Mathematically, this is equivalent to solving the following optimization problem:

$$ma{x_{{w_i},{c_i}}}Cov\left( {{t_i},{u_i}} \right),\,subject\,to\,\left| {\left| {{w_i}} \right|} \right| = \left| {\left| {{c_i}} \right|} \right| = 1$$

After extracting components, PLS provides a low-dimensional subspace that captures the most important shared covariance between |${y_1}$| and |${y_2}$|⁠. Therefore, |${t_i}$| and |${u_i}$| are the variables that define the communication subspaces.

The investigation of communication subspace is still in its early stages. This is primarily because the techniques for recording large-scale neural activity across multiple brain regions have only recently become available. Sophisticated analytical methods are essential to interpret the complex data obtained from these large-scale recordings effectively. The population doctrine is currently the most reliable and tractable framework for bridging neural activity in a high-dimensional space and information encoding and behavior and is made possible by dimensionality reduction and modeling interaction between brain regions. As interest in investigating information encoding at a population level continues to grow, and as large-scale recording technologies advance and become more widely adopted, statistical methodology to model population doctrine is increasingly critical to neuroscience.

Elucidating the topology of cortical functional network

So far, we have discussed information representation at the level of single neurons and populations of neurons. We now delve into the domain of network neuroscience [117,118]. Network neuroscience centers around the concept of functional connectivity, which refers to the statistical dependencies observed between individual neurons during specific tasks or brain states Fig. 5a-c). These dependencies can arise from unidirectional or bidirectional local synaptic connections, shared inputs and intrinsic neural states. Measuring structural connectivity, which focuses on anatomical connections between neurons, is challenging in living subjects. Therefore, researchers often rely on functional connectivity to infer the organization of neuronal networks that underlie brain functions. Furthermore, unlike structural connectivity, functional connectivity captures temporal co-activation and statistical dependencies between neurons, revealing dynamic communication patterns rather than static physical connectivity. This approach allows researchers to infer networks that reflect the interactions underlying specific brain functions or states and provides insights into the underlying mechanisms.

Fig. 5.

Analysis of cortical functional network. (a) Inferred spiking activity of individual neurons recorded with FASHIO-2PM. (b) A functional network is estimated with various methods. Here is an adjacency matrix constructed by a partial correlation of pairs of neurons conditioned by a movement signal. (c) Examples of graph theory-based measurements of network topology.

Correlation-based methods

Within network neuroscience, statistical methods are needed to estimate functional connectivity. However, it is still an open question how to reconstruct the functional connectivity from the observed neural activity without observing the connections themselves. The simplest way to estimate functional connectivity is to compute Pearson’s correlation coefficients between neurons [119]. This method quantifies the linear relationship between the activity patterns of two neurons. However, correlation not only reflects the direct connection between two neurons but also largely reflects the response to shared inputs. Partial correlation, a technique to condition the correlation by accounting for the influence of confounding factors such as common inputs, is used to address this issue [5]. This approach potentially provides a more accurate estimate of the linear relationship between two neurons. However, as the number of recorded neurons increases, it becomes difficult to accurately estimate the connectivity between neurons because the number of estimated coefficients increases quadratically, leading to an increase in spurious correlations, overestimation of shared activity and insufficient conditioning of partial correlations. To improve the estimation of functional connectivity in large-scale neural recordings, a regularization technique is used that reduces the degrees of freedom by imposing constraints on the data. Yatsenko et al. proposed that modeling the correlation matrix by assuming sparse connections between neurons using regularized partial correlations and shared inputs using a low-rank component improves the estimation and interpretation of functional connectivity [120].

Causal inference methods

In addition to these correlation-based methods, causal inference methods are used to estimate directed connections between neurons. For example, Granger causality is a statistical test that determines whether one time series can predict another [121]. Granger causality is based on linear autoregressive models and tests whether the activity history of one neuron can improve the prediction of the activity of another neuron [122,123]. Another method is transfer entropy [124]. Transfer entropy is based on information theory and measures the reduction in uncertainty of the future activity of one neuron given the past activity of another neuron. Unlike Granger causality, transfer entropy does not assume linear relationships and can capture both linear and nonlinear dependencies. Similar to correlation-based methods, conditioning Granger causality and conditioning transfer entropy are proposed to improve estimations by removing confounding factors [125,126].

Statistical inference of generative models

One of the other approaches is to reconstruct a network by statistical inference from generative models. This framework explicitly assumes that observed data are generated by a generative model with a weighted network. It aims to reconstruct the network from observed data by solving it as a statistical inference problem. One approach utilizes a Bayesian framework [127]. This method incorporates prior knowledge about possible network structures and updates it with the observed data to arrive at posterior distribution [128]. This distribution represents the most likely network structures based on both prior knowledge and the observed data. Recently, researchers proposed an algorithm to estimate parameters based on the minimum description length principle [129]. This is a model selection principle where the best model is the one that captures the data with the shortest description. The proposed algorithm offers several advantages, including faster computation and reduced overfitting, which are particularly beneficial for reconstructing large networks [129]. In addition to these methods, researchers have developed data-constrained RNN (dRNN) models that reproduce the obtained neural dynamics by training the RNN [92,130]. The dRNN consists of the same number of units as the actual recorded neurons, with a randomly assigned initial weight between units. The dRNN is trained to reduce the error between the observed neural dynamics and the dRNN output. This method produces the interpretation of how the neural dynamics is generated in the neuronal network and has a possibility to decompose the dynamic neuronal interaction in the brain [92]. However, estimating the ground truth connections in large-scale recordings is challenging. The model’s degrees of freedom increase with the number of units, allowing for different connectivity patterns to reproduce the same dynamics. Thus, it is essential to integrate physiological constraints into the modeling framework to reduce these degrees of freedom [131–133].

Network analysis using graph theory

After reconstructing the functional network, the next step is to analyze the network using various methods [117] (Fig. 5d). Network analysis using graph theory examines the topology of the reconstructed network. This analysis helps identify key properties and principles that govern neural interactions and brain function. The methods introduced earlier capture the pairwise relationships between neurons, resulting in an adjacency matrix that represents them as a graph, where nodes represent neurons and edges represent functional connections. If the matrix is binarized at a certain threshold value, the resulting graph is unweighted. If the values are retained, the graph is weighted. Additionally, a symmetric matrix represents an undirected graph, while an asymmetric matrix represents a directed graph. From the estimated graph, several network properties can be measured as follows: (i) Node degree measures the number of connections each neuron has, indicating its centrality and potential influence within the network. (ii) The clustering coefficient indicates the extent to which neurons form clusters or groups, reflecting connectivity patterns. (iii) Path length represents the average number of steps required to connect any two neurons, providing insights into the efficiency of information transfer across the network. (iv) Centrality measures identify key neurons that play critical roles in network communication. Betweenness centrality indicates how often a node acts as a bridge along the shortest path between two other nodes. Eigenvector centrality considers a node’s influence based on the number and quality of its connections. These measurements produce topological insights into the neuronal network. In our study, for example, we examined the network topology that underlies spontaneous activity of >10 000 neurons recorded from multiple cortical regions [5]. In this study, we obtained a functional connectivity matrix by calculating partial correlation, which conditioned movement information as a confounding factor, and an unweighted undirected graph was constructed by thresholding. We then found that the degree distribution is heavy-tailed, but not a power-law distribution. This finding indicates that the functional network of the cortex is neither a random network nor a scale-free network. Furthermore, we found that there are a small number of hub neurons with a very large number of functional connections. By calculating network metrics such as shortest path length, cluster coefficient and small-world propensity [134], we showed that the cortical functional network has small-world properties. In general, small-world systems are known to process information at low energy cost. This network structure of the brain may be the basis for our brain’s ultra-low energy consumption.

Outlook

Bridging between single neuron, population dynamics and brain function

Recent advances in large-scale Ca²⁺ imaging allow researchers to capture detailed cortical dynamics across various levels, from single neuron level to population level. Future research will focus on bridging these levels to understand how they contribute to brain functions. A powerful tool for bridging them is to use two-photon optogenetics [65,66,135,136]. This technique allows researchers to activate or inhibit specific neurons while simultaneously observing the effects on behavior and neural circuits. Previous studies with conventional two-photon microscopes have revealed that the activation of a small subpopulation of neurons in the sensory cortex is sufficient to drive behavior [66,137]. However, it remains an open question how such manipulation influences cortical-wide population dynamics and how the population dynamics drives behavior. Two-photon optogenetics combined with large FOV imaging offers a unique opportunity to investigate the precise communication mechanisms between cortical regions that orchestrate behavior.

Revealing organizational principles of cortical network

Understanding the organizational principles of cortical networks is crucial for deciphering the complex mechanisms of brain function. Two-photon optogenetics can also directly probe the influence of specific neural populations within these networks. For instance, by manipulating hub neurons – highly connected neurons within the network – researchers can assess their impact on network dynamics and gain insight into the hierarchical organization of neural circuits.

Furthermore, an intriguing question is the relationship between functional connectivity and structural (synaptic) connectivity [138–140]. Recent advances in electron microscopy enable the integration of large-scale functional two-photon imaging with dense reconstruction of neural structures, allowing for a direct comparison of functional and structural connectivity patterns [141,142].

Bridging the gap between data and models

The increasing data volume and complexity generated by advanced imaging techniques necessitate the integration of computational modeling to extract meaningful insights. Utilizing simulation frameworks opens new possibilities for bridging computational models and empirical data [92,129]. Methods for reconstructing neural networks through statistical inference of generative models allow researchers to set various parameters and simulate their effects within the model. This comprehensive view of how different parameters affect network dynamics leads to new hypotheses about neural dynamics and their underlying networks. Researchers can make targeted perturbations based on simulation predictions and observe the effects in vivo. This approach holds significant potential for advancing our understanding of brain function and provides a robust framework for validating and refining computational theories.

Application to disease research

Large-scale Ca²⁺ imaging technologies will revolutionize our ability to study not only healthy brain function but also dysfunction. Ca²⁺ imaging offers several key advantages for studying neural activity in the context of disease. First, its high spatial resolution allows researchers to track the activity of individual neurons over extended periods [143]. This longitudinal tracking capability is invaluable for elucidating the dynamic changes in neural function associated with disease progression. Second, advancements in technology have enabled the integration of Ca²⁺ imaging with spatial transcriptomics [144,145]. This combination facilitates the examination of gene expression patterns within the same neurons whose activity has been recorded and establishes connections between single-neuron genetics and disease-related alterations in neural dynamics. Third, large-scale Ca²⁺ imaging can be combined with other imaging modalities, allowing researchers to simultaneously monitor multiple aspects of neural function, including neuromodulator dynamics, blood flow, and glial cell activity and morphology [146–150]. This comprehensive data collection provides a richer understanding of the complex interplay between different factors that contribute to disease pathology.

Conclusion

In conclusion, future directions in large-scale Ca²⁺ imaging research focus on integrating information across multiple levels of analysis, from single neurons to entire networks. This integration will lead to a more comprehensive understanding of the complex mechanisms underlying brain function and behavior. Advances in microscopy technology and the integration of various new techniques (e.g. two-photon optogenetics, spatial transcriptomics and theoretical approach) are opening exciting doors for the development of more targeted and effective therapeutic interventions for neurological disorders. Achieving these goals requires collaboration between experimentalists and theorists. We hope this review promotes mutual understanding among researchers in different fields and fosters future collaborations.

Acknowledgements

We thank Hiroyuki Uwamori, Kensuke Yoshida and Takeru Matsuda for their constructive comments on this manuscript.

Funding

AMED-Brain/Minds 1.0 Project (JP15dm0207001 to M.M.); AMED-Brain/Minds 2.0 Project (JP23wm0625001 to M.M.); Grant-in-Aid for Transformative Research Areas (B) from the JSPS (JP20H05775 to M.M.); Grant-in-Aid for Transformative Research Areas (A) from the JSPS (JP24H02313 to M.M.); Grant from Toray Science Foundation (to M.M.); the Uehara Memorial Foundation, Overseas Post-doctoral Fellowship (to Y.O).

Conflicts of interest

The authors declare that they have no conflicts of interest.

Author notes

†

contributed equally to this work.