Introduction

Early postnatal development in rodents is a period of intense circuit wiring and remodelling at several scales through various major processes that include neuronal growth, synaptogenesis, apoptosis, migration, rise of intracortical connectivity, functional maturation of inhibitory synapses and the disappearance of transient connectivity schemes ^1–4. Critical periods for various sensory systems open and close during this time, further highlighting the profound reshaping of cortical networks ⁴. This developmental period of remodelling is crucial for establishing the functional architecture of the mature cortex. Most importantly, all of these developmental processes are activity-dependent, with collective dynamics playing a critical role in the proper integration of neurons into functional networks ⁵. These network dynamics sequentially unfold while relying on different mechanisms and circuits for their generation ^1,5. The emergence of functional brain circuits during development is therefore a precisely timed choreography, the timing of which is inherently tied to the age of the organism under study.

However, developmental age lacks precision due to significant variations in physical characteristics and growth patterns, even among offspring from the same genetic lineage (Fig. 1). Hence, longitudinal imaging from the same animal across days is the optimal solution to better capture the evolution of circuit dynamics during development. Additionally, developmental variability extends beyond the organism level, becoming even more pronounced at the single neuron level (Fig. 1). While population-level descriptions of cortical circuit development are crucial, individual neurons exhibit unique developmental trajectories rooted in their specific origin, birth timing and cellular identity ¹. In addition, the singular dynamics of sparse individual neurons can matter, as demonstrated for the rare hub cells ^6,7. Thus, in order to fully understand the circuit basis of cortical development in health and disease, it is crucial to track neuronal activity at both population and individual cell levels in the growing brain of the same animal. This hurdle spans both experimental methodology and data analysis.

Figure 1:

Significant progress has been made over the recent years to develop innovative solutions to meet this challenge on the experimental side, using two-photon calcium imaging. Technical improvements included modified head plates and specialized surgical and care protocols ^6,8–12. Still, in many cases, different neurons from the same mouse were recorded at different ages, or the tracking of individual cells could only be achieved through visual inspection and manual annotation of relatively sparse, often genetically-subscribed neuronal populations expressing a calcium fluorescent reporter. The gold standard for developmental circuit analysis would be automated tracking of densely labeled neuronal populations, enabling efficient longitudinal monitoring while eliminating the burden of manual cell tracking. An algorithm capable of accurate automatic tracking of the same neurons across multiple days of brain development is thus essential.

Methods for tracking the activity of large populations of cells across days have been successfully developed and deployed in the adult brain ^13–16. In adult circuits, cell tracking is indeed relatively straightforward due to their structural stability, with minimal tissue growth, negligible morphological changes, and constant numbers of neurons. In contrast, the first weeks of mouse neocortical development are characterized by rapid and critical developmental changes spanning several scales from synapses, to single neurons and networks. These include extensive brain growth ^17,18, substantial morphological changes ¹⁹ as well as changes to cell numbers due to programmed cell death ¹⁰. This makes tracking the same cells across sessions substantially more difficult during development compared to the adult brain.

To overcome the technical limitations detailed above, we developed an experimental protocol using chronic calcium imaging in mice and a novel cell tracking algorithm (Track2p), specifically tailored to development. This allowed us to track the activity of large, densely labelled, populations of neurons during early postnatal development, which has not been possible before. The algorithm overcomes the challenges of cell tracking during brain growth by applying sequential registration and cell matching steps, using the spatial overlap of cells on adjacent recordings as a matching criterion. Track2p is freely available as an open-source package with an interactive GUI, enabling researchers to analyze longitudinal calcium imaging data from both developing and mature circuits.

Applying the algorithm to a dataset recorded during the second postnatal week in mouse barrel cortex, a critical period for the formation of topographic maps in that region, yielded hundreds of identified neurons tracked across all days. Assessing the quality of the algorithm by benchmarking it on a newly generated ground truth dataset showed high tracking performance ²⁰. Leveraging ground truth benchmarking, we demonstrate that explicitly accounting for developmental processes, such as brain growth, are critical for accurately tracking cells during postnatal development. Our work thus shows that chronic calcium imaging and cell tracking using Track2p can be used to monitor the changing physiological properties of large populations of matched neurons during early cortical development. We demonstrate that the statistics of activity patterns in the tracked population display two periods of stability, with a critical transition point around postnatal day 11 (P11), marking the emergence of a stable behavioral state representation.

Results

Cell tracking using image registration and ‘overlap-based’ matching

Tracking neuronal activity across multiple days presents unique challenges due to the dynamic nature of brain development. We developed a novel tracking algorithm, called Track2p. As in other tracking algorithms for calcium imaging data ^13,16,21, the final goal of Track2p is to follow individual cells across sessions allowing the user to compare their functional properties in downstream analyses. To achieve this goal, the algorithm takes as input a set of preprocessed recordings (using Suite2p ²³), each consisting of a set of regions of interest (ROIs i.e. putative neurons detected based on activity, see Methods) and their respective calcium fluorescence traces, as well as a mean image of the field of view (FOV). Briefly, the algorithm aims to match ROIs in any given pair of consecutive sessions based on their spatial overlap, assuming that the more the two overlap in anatomical space, the more likely they correspond to the same neuron. Due to developmental processes such as brain growth and other experimental factors, it is necessary to account for day-to-day changes that occur between the two recordings before computing spatial overlaps. This is achieved by performing affine image registration on the mean FOV images between consecutive days.

We apply the registration and spatial matching iteratively, starting with the first pair of sessions (s₀ and s₁) as follows (Fig. 2A):

Figure 2:

Firstly, we estimate the spatial transformation between s₀ and s₁ using affine image registration (Fig. 2B, the transformation is denoted as T). We employ affine transformation, since it can account both for rigid transformations (rotations and translations arising from minor mismatches in FOV alignment across experiments) as well as scaling and shearing (mostly due to brain growth). The changes across the two consecutive recordings are approximated as the transformation registering the mean FOV from session s₁ (green in Fig. 2B) to the mean FOV of s₀ serving as reference (red in Fig. 2B).

Secondly, the computed transformation is applied to the ROIs from session s₁ (green in Fig. 2C) to align them to ROIs from session s₀ (red in Fig. 2C). The amount of spatial overlap after registration (yellow in Fig. 2C, bottom), can indicate the accuracy of the estimated transformation between the two days. Assuming that the transformation is accurately estimated, ROIs corresponding to the same cells display substantial spatial overlap, with some ROIs from one recording also potentially overlapping poorly, if not at all, with any ROI from another.

Thirdly, once the ROIs are aligned, the algorithm proceeds with the matching (Fig. 2D). This is done by computing a spatial similarity metric (intersection over union - IoU) between each ROI from session s₀ and each transformed ROI from session s₁. Matches are then assigned in a globally optimal way by maximising the sum of IoU values across all matches using a linear sum assignment algorithm.

Finally, since two consecutive sessions contain different sets of detected cells (see * in Fig. 2D, bottom) and since ROIs can overlap even if the signal does not come from the same cells (see + in Fig. 2D, bottom) we perform an additional filtering step on the assigned matches. Assuming that the IoU values for putative true and putative false matches come from different distributions, we would expect a bimodal distribution of IoU values across all assigned matches (see histogram in Fig. 2E). To reject the putative false positives, we compute a threshold based on automatic thresholding methods (Otsu’s method ²⁴). This ensures that assigned matches with low spatial similarity are rejected (see + in Fig. 2D and E) while the matches with high similarity are accepted (see x in Fig. 2D and E). This yields the final matching for the first pair of consecutive imaging sessions (s₀ and s₁). In the case of more than 2 recordings, this tracking procedure is then iteratively applied for all consecutive pairs of sessions (s₀ to s₁, s₁ to s₂ … s_N-1 to s_N , and so on), with tracks being extended sequentially from s₀ to s_N based on the identified matches, and terminated if a match was not identified at any particular session (Fig. 5F).

We provide an open-source implementation of the algorithm, combined with a user- friendly graphical user interface (GUI), allowing non-specialist users to run the algorithm and interact with its outputs (see Fig. 2A top, Extended Data Fig. 1 and Supp. Info.: Software). Both the algorithm and the GUI come with a simple installation procedure and extensive documentation to facilitate the ease of use and accessibility of Track2p.

Tracking neurons across days in the early postnatal neocortex

In order to validate our method on experimental data, we next applied Track2p to a longitudinal dataset consisting of daily recordings of the same 720x720 μm FOV throughout the second postnatal week of mouse barrel cortex development (postnatal days 8 (P8) to P14, n=7 imaging sessions from one mouse; for more details see Methods.) In addition to the calcium indicator (GCaMP8m), we virally expressed a sparse anatomical marker (tdTomato) targeting GABAergic neurons using conditional expression in GAD67-Cre mice ²². This dual- labelling strategy facilitated reliable FOV identification across sessions and provided an anatomical reference for Track2p’s initial cross-session registration and tracking (Fig. 3, for comparison with tracking based on GCaMP8m registration see the ‘Benchmarking on manually tracked cells’ section). Before running Track2p, each recording was preprocessed using Suite2p ²³ for active cell detection (segmentation) and calcium fluorescence trace extraction (see Methods). The outputs of Suite2p (ROIs, traces and FOV images) were then used as inputs to the Track2p algorithm.

Figure 3:

The Track2p algorithm appeared to successfully register the mean FOV for each successive pair of imaging sessions (based on visual inspection of Track2p outputs; for the first and last pairs see Fig. 3A and B bottom left; for the equivalent of Fig. 3 across all days see Extended Data Fig. 2A). Registering the ROIs for each pair of consecutive recordings showed a great degree of overlap for the majority of ROIs (yellow area in Fig. 3A and B bottom). We did however observe several ROIs that were only present in one session from the pair, likely due to differences in cell detection and developmental factors such as growth, silencing or cell death. When matching ROIs across sessions using spatial overlap (IoU), we would expect these ROIs to show significantly lower IoU values compared to the cells that were present in both sessions. Indeed, the IoU distribution of assigned matches revealed a bimodal distribution for each pair of imaging sessions, allowing us to use classical histogram thresholding methods ²⁴, clearly separating the distributions of putative true and false matches (Fig. 3A and B bottom right, Extended Data Fig. 2B). Propagating the putative matches yielded a total of 728 ROIs that were tracked across all days in this example mouse (out of 1988 ROIs detected in the first recording, Fig. 3E). We used the subset of cells successfully tracked across all days for all our future analyses.

Plotting ROIs overlaid on top of the mean image of the GCaMP8m imaging channel (Fig. 3D for contours on whole FOV, C for magnified example matches) and visually inspecting the matches using the Track2p GUI indicated excellent tracking with Track2p. These visualisations also revealed the substantial growth of the neocortex during the course of the experiment (FOV area covered by ROIs at P8 compared to P14 in Fig. 3D), with some cells growing out of the FOV (see top and left in Fig. 3D) and the pairwise distances between matched ROIs increasing by approximately 15% (Fig. 3F).

Applying Track2p on an example longitudinal imaging dataset thus demonstrated that it can successfully track activity from a large number of putatively matched neurons throughout early postnatal development in mice, as confirmed by visual inspection. To assess the tracking performance in a more quantitative and objective way, we next benchmarked Track2p’s performance on manually generated ground truth for our specific experimental setting.

Benchmarking on manually tracked cells

In order to benchmark Track2p, expert annotators manually tracked cells across sessions (similarly as in the Cell Tracking Challenge ^20,25). Benchmarking was performed based on a dataset from 3 mouse pups imaged under the same experimental conditions as described in the previous section (including the dataset shown in Fig. 3). For each experiment, we chose 64 homogeneously distributed ROIs detected on the first day and tracked them based on visual inspection across consecutive sessions (for more details see Methods and Extended Data Fig. 3). This left us with on average 20 neurons per experiment that we were able to manually track across all days (‘GT’ in Fig 4.). We then proceeded to compare these to cell tracks identified by Track2p. For the purposes of evaluation we used a fully automatic tracking procedure (without manual curation and with default Track2p parameters).

Figure 4:

Different metrics exist for evaluating cell tracking. Since we are aiming to track cells across all days, a robust cell tracking metric should reward perfect track matches with ground truth while imposing penalties for missed or incorrectly identified tracks. For this reason we used a biologically inspired ‘complete tracks’ quality metric, which corresponds to the F1 score for completely reconstructed full tracks (we refer to this value as ‘CT’ according to ^20,25). A CT score of 1 would correspond to perfect tracking, while a CT score of 0 would correspond to no matches or a large proportion of wrong matches (see Fig. 4A and Methods for more details). Additionally, to assess tracking performance over time, we quantified the proportion of reconstructed ground truth tracks across progressively shorter time intervals (‘Prop. correct in Fig. 4C-F, see Methods).

Track2p was evaluated under three different scenarios, and compared to the performance of a widely used algorithm developed for longitudinal tracking in the adult brain (CellReg ²¹, see Methods). The initial evaluation was done using day-to-day registration based on a sparse anatomical marker (as in Fig. 3). Calculating the CT score in this condition showed remarkably high tracking performance for all datasets (see ‘Anatomical’ in Fig. 4B, mean CT = 0.93). To test whether a sparse marker was strictly necessary for successful tracking, we next ran the algorithm using the mean image of the GCaMP8m channel as a comparison. Interestingly, this evaluation showed a similar performance (see ‘Functional’ in Fig. 4B, mean CT = 0.91), indicating that dense calcium indicator labelling alone can provide comparable tracking performance, eliminating the requirement for sparse anatomical labeling.

To show the importance of a method tailored to the growing brain, we next compared the baseline tracking performance with two alternative conditions: firstly, by using Track2p without explicitly accounting for day-to-day growth, performing a rigid instead of affine registration of consecutive recordings; and secondly, by comparing Track2p to CellReg tracking ²¹. In both scenarios, tracking performance significantly deteriorated compared to our baseline method (see ‘Rigid’ and ‘CellReg’ in Fig. 4B, mean CT = 0.22 and CT = 0 respectively). Interestingly however, these methods still managed to correctly reconstruct a portion of tracks across shorter age spans in certain cases (Fig. 4E and F), with the performance dropping significantly for longer tracks in comparison to the baseline condition (Fig. 4C and D).

These comparisons demonstrate Track2p’s robust cell tracking performance in the developing brain, maintaining high accuracy over extended periods of development. However, tracking performance certainly depends on the type of experimental data (age, model system, brain area, imaging parameters, FOV alignment etc.), hence we also provide additional resources and documentation, allowing users to benchmark Track2p for their particular use case.

Development of firing statistics from hundreds of tracked neurons across postnatal development

Having validated the performance of cell tracking, we next analyzed the development of functional properties of the tracked population of neurons across the second postnatal week of mouse barrel cortex development. For this, we used a full dataset of 6 mice imaged daily for a minimum of 6 consecutive days within the second postnatal week (P7 to P14, see Fig. 5B for summary, we used this dataset for all subsequent analyses). On average 526 (± 190 std) neurons per mouse were successfully tracked across all days using Track2p. During the course of the experiment, the weight of a mouse increased on average by 55 % (± 18% std), with the pairwise distance between neurons increasing on average by 15 % (± 7% std) (Fig. 5C and Extended Data Fig. 5F). Weight was found to vary significantly across mice, as early as the first experimental day, and to correlate with brain growth (r=0.92, Extended Data Fig. 5E). Our robust tracking capability provides a comprehensive view of the evolution of neuronal dynamics within single mice, while also allowing us to take into account the heterogeneity of the developmental timelines across mice.

Figure 5:

The early postnatal period of cortical development studied here is characterized by the transition from synchronous calcium events recruiting many neurons to progressively more decorrelated population dynamics ^26,27. In order to visualise these global changes in our dataset of longitudinally imaged mice, we first plotted rasterplots of calcium fluorescence traces as a function of time for all tracked neurons on each imaging day (Fig. 5A for P8 and P14, all days in Extended data Fig. 4A). Visual inspection of these rasterplots clearly indicates the disappearance of recurring periods of highly synchronous activity from around P11 onwards as well as a global increase in single-neuron activity rates. Accordingly, quantifying calcium event rates as a function of age indicated a gradual increase in the mean and a widening of the distribution with age (Fig. 5D and Extended Data Fig. 5A and G, all statistical comparisons were done between an ‘early’ (≤ P11) and ‘late’ (> P11) set of recordings). This evolution also signalled the transition towards long-tailed firing rate distributions, which are ubiquitous in adult brain circuits ²⁸.

Next, we analyzed the distributions of pairwise correlations between activity traces from all longitudinally tracked neurons as a function of mouse age (Fig. 5D and Extended Data Fig 5B and H). Consistent with the disappearance of highly synchronous network events observed in the rasterplots, we found a significant decrease in mean pairwise correlations, indicating a progressive decorrelation of neuronal activity (Fig. 5E). The spatial distribution of pairwise correlations also evolved across days from highly correlated locally to more broadly distributed (Fig. 5F), suggesting the gradual breakdown of spatially clustered assemblies ⁵.

We finally turned to characterising the dominant population patterns of neural activity using principal component analysis (PCA). Interestingly, we observed an increase in the number of components required to explain a fixed amount of variance in the neural data across days, suggesting a developmental increase in dimensionality (Fig. 5G), consistent with the statistics of spontaneous activity described in the adult brain ²⁹. Notably, the plots of summary statistics for our entire dataset, indicated a clear outlier, consistent across all quantifications (Fig. 5D, E, F and G). This outlier mouse displayed a similar but delayed developmental trajectory, compared to the other mice. Interestingly, it also showed the lowest initial weight and a less pronounced growth in weight and cortical size (Fig. 5C), likely indicating a lower maturation stage at the onset of the experiments, although the contribution of other experimental factors cannot be excluded given the invasiveness of imaging surgery.

To fully leverage our longitudinal approach, we next turned to comparing the functional properties of the same neurons across development, which is only possible when tracking cells across days.

A marked reorganisation of the functional network structure during the second postnatal week

We first examined the stability of the correlation structure, which we will refer to as ‘Functional connectivity’ (FC). Two alternative possibilities could be envisaged (see Fig. 6A,: (i) FC could be conserved across days, meaning that a highly connected pair on a given day would remain highly connected on the next (Fig. 6A, left); (ii) FC could reorganise, losing the structure from the previous day (Fig. 6A, right). In the first case, the FC values for all pairs would remain similar across the two days (Fig. 6A, bottom left) whereas in the second case they would be different (Fig. 6A, bottom right). To discriminate between the two cases, we defined a ‘FC similarity’ score for any pair of imaging days with matched neurons for a given mouse (Pearson correlation (r) across all neuron pairs, see Fig. 6A, bottom). Such FC similarity metric, encapsulates a variety of possible sources of network changes, from local or long-range connectivity, intrinsic excitability to neuromodulation.

Figure 6:

Quantifying FC similarity across two early sessions for an example animal (P8 to P9, Fig. 6B top left) indicated remarkable stability, almost identical to the FC similarity within a given imaging day session which we took as a reference (computed as FC similarity between the first and second half of a same day recording, Fig 6B, bottom left). Conversely when comparing across more distant developmental ages, we noticed that a large part of the correlation structure was lost, resulting in lower FC similarity (P8 to P14, Fig. 6B, top middle). Interestingly, FC similarity between a pair of later developmental sessions, was again comparable to within-session similarity (P13 to P14 and P14 within, Fig. 6B, right), indicating that the functional network structure was stable, but different from that of earlier ages. Of note, within session stability was consistently lower later than earlier. Quantifying calcium event rate similarity across days indicated a similar pattern (see Extended Data Fig. 6F-J).

To investigate this further we proceeded to compute the FC similarity for all possible combinations of sessions in all mice. Interestingly we observed a sigmoid-like decay of FC similarity when taking P8 as a reference (Fig. 6C). Plotting the full FC similarity matrix for all possible combinations of sessions also indicated two stable FC regimes with seemingly two blocks along the diagonal (Fig. 6D left for example, 6D right for average and 6E for statistical comparison). Interestingly, the sharpness of the transition in FC varied across mice (Extended Data Fig. 6C) and could in part be explained by the weight of the mouse at the onset of the experiment (Fig. 6C, linear part of sigmoid fit). Such variation across mice could reflect either inherent individual differences in developmental processes or experimental factors, where mice with a lower initial weight may have experienced delayed development due to higher sensitivity to the imaging procedure.

Altogether, our analyses indicate that the second postnatal week marks a transition between two stable functional network structures in the barrel cortex of developing mice. Additionally, this shows how tracking cells during development opens new analysis possibilities using self-referencing between neurons, providing new insights into how the developmental choreography unfolds in the same population of neurons across days.

Developmental emergence of stable behavioral state modulation

As a first attempt to investigate the mechanisms underlying this network transition, we examined neural activity in relation to behavioral state. The second postnatal week marks the emergence of active sensation, suggesting that changes in arousal state modulation might drive this network reorganization. In the adult cortex, arousal strongly shapes both global firing rates and neural correlations ^30–33, making it a promising candidate for orchestrating the developmental shift we observed.

To this aim, we sorted the neurons based on the similarity of their activity patterns, using Rastermap ³⁴, and examined the relationship between neuronal activity and behavioral states, indirectly measured using the videos capturing mouse movement (as in ³⁰, see Methods). Interestingly, this sorting revealed a subpopulation of neurons that were highly correlated with the animal’s motion at later, but not earlier developmental stages (Fig. 5A, see Extended Data Fig. 4A for all days). We first quantified this relationship by computing the temporal correlation between the first principal component of neuronal activity with the animal’s motion. Interestingly, this correlation showed a steep increase after P11 for most mice (Fig. 7D). This suggested that global fluctuations in neural activity are more strongly modulated by active movement at later than earlier developmental stages.

Figure 7:

Finally, we asked whether behavioral state could be decoded from population activity dynamics using regression analysis. Despite substantial variability, the same-day decoding performance increased with development (Fig. 7A-C). Since we tracked the same neurons across days, we could also probe the stability of this representation across days, by fitting a model on one day and testing it on all other days (cross-day decoding, Fig. 7E-H, Extended Data Fig. 7B and C^15,35). Interestingly, this analysis showed that, once developed, this representation was indeed stable, allowing for accurate cross-day decoding, with the same neurons showing either consistently positive or negative modulation across days (see Extended Data Fig. 7D for the activity traces of an example cell with a high weight for decoding).

Hence, by combining Track2p longitudinal cell tracking with decoding approaches, one can systematically map the emergence and stability of cortical representations in the developing brain.

Discussion

Here we described Track2P, a cell tracking algorithm that can be used to follow the changing dynamics of hundreds of matched neurons from daily calcium imaging recordings in the growing brain of living mouse pups. At the core of the method lies sequential affine registration followed by cell matching for each subsequent pair of recordings, which lead to a better tracking performance during development when using Track2p compared to other methods. Using this benchmarked approach, we observed a sharp developmental transition from highly synchronized activity to multidimensional, behavior-state dependent neural dynamics. Beyond this key finding, our method opens a much needed tool for investigating the developmental functional trajectories of individual neurons during early postnatal brain development, and their deviations caused by genetic mutations or environmental perturbations.

Developing and benchmarking Track2p showed that it is able to track large populations of neurons despite substantial brain growth and other developmental changes through extended periods of time during development. We suggest this is mostly due to affine transform being a good approximation of the growth processes that occur between two consecutive days of imaging, with the cell matching step allowing to account for slight non- linearities in tissue growth day-by-day. Currently available methods typically register all imaging FOVs to a single reference, which we believe result in the accumulation of non- linearities making registration more difficult and resulting in poorer tracking performance.

From a practical perspective, Track2p reduces manual cell tracking time by transforming a tedious process of manual annotation into a fast and automated procedure. Indeed, considering the time needed for manually tracking all cells in our dataset, we estimate that it would require approximately 5 person-days (120 hours) of work per subject under our experimental conditions, highlighting the necessity of automatic tracking in large scale calcium imaging recordings. To aid future research, we provide Track2p fully open-source, with substantial documentation and an accompanying GUI facilitating the use of the algorithm by users without previous coding knowledge.

Besides developing the algorithm, we also showcase analysis techniques that can be used to gain unique insights from longitudinal calcium recordings, drawing inspiration from previous research studying plasticity, learning and representational drift in the adult brain ^15,35,36. We specifically highlight: 1) quantifications of functional statistics of the tracked population across days (Fig 4); 2) correlation analyses to compare these across days (Fig 5); 3) regression analyses to study the emergence and stability of neural representations across recording days (Fig 6).

Last, we would like to emphasize that this analytical advancement is grounded on optimized experimental procedures to image daily cortical dynamics through a glass window mounted on the developing pup head. Previous studies had described adaptations of the surgical procedure or head fixation to developing pups ⁸, however, by measuring the increase in pairwise distance between tracked neurons as a function of age per animal, we now provide a quantitative metric to estimate the invasiveness of the procedure. Interestingly, we observed a growth rate that closely matched ex vivo quantifications ¹⁸, indicating a minimal impact of our experimental preparation on anatomical brain growth, also confirmed by mouse weight measurements. While we observed a slightly slower weight gain in imaged pups compared to their littermates, further investigation is needed to definitively assess whether repetitive daily imaging might influence cognitive developmental trajectories, for example through maternal separation.

Applying these techniques in our dataset during early postnatal development, we observed a significant and fast change in the nature of neuronal dynamics and state modulation in the barrel cortex, centered on P11, that manifests in several ways: (1) an increase in activity rates; (2) a decrease in pairwise correlations; (3) a shift from locally clustered to widely distributed correlations; (4) an increase in the dimensionality of spontaneous activity; (5) a remapping of the functional pairwise correlation structure; (6) the emergence of a stable representation of spontaneous active movement.

If the first three changes had already been described previously in different cortical areas ^{5,26,27,29,37–40}, the last two could not be observed without the dense longitudinal tracking of the same neurons across several days in the growing brain permitted by Track2P. This multifaceted transition likely reflects several concurrent developmental processes, including changes in sleep architecture and neuromodulatory tone ^38,41–43, maturation and rewiring of local inhibitory circuits ^{6,12,37,40,44–47} or disengagement from peripheral sensory inputs ^38,48. This postnatal spurt is correlated with sparsification, disappearance of the spindle-burst oscillations and increase in dimensionality of the representations ^{26,27,29,37,38,49–52}. Such drastic and fast changes are accompanied by significant behavioral changes, maybe the most salient ones being the change in the nature and duration of sleep as well as the start of active exploration⁵³. The start of behavioural state-modulation and dimensionality increase could reflect the same phenomenon by which spontaneous motor behavior, including facial movements, drives multidimensional brainwide activity in the adult visual cortex ³⁰. This representation is present and stable from P11 onwards, as revealed using our decoding approach. The early absence of representation of spontaneous motor behavior is in agreement with the previously reported lack of reafferent brain activity in response to self-generated wake movement until P11 ^54–56. Thus, although we did not distinguish here between sleep and wake-generated movements, it is likely that our decoding mostly took into account wake behavior, given the long periods of motor activity observed. Hence, we cannot exclude that twitching activity occurring during active sleep could contribute to an earlier representation of spontaneous motion, as reported previously ⁵⁷. More refined behavioral analysis, particularly of sleep-wake transitions and twitches, could provide additional insights.

Therefore, the mid-second postnatal week marks a transition between two stable functional connectivity structures, indicating that globally, early and late functional connectivities differ. However, singular developmental trajectories cannot be excluded. Future studies could examine whether different cell types transition at different times or different paces, for example depending on their time of birth ¹, or whether unique neurons, such as hub cells, could maintain exceptional and stable functional connectivity ⁶. One population of particular interest for tracking are the neurons that contribute the most to behavioral prediction. Indeed, these could form a distinct population of movement-correlated neurons embedded in specific wiring schemes, as recently demonstrated ⁵⁸.

This developmental transition centered around P11, which corresponds roughly to birth in human brain development, spans a period between one and four days, depending on the initial weight of the imaged pup. This transition may not yet be the last step before a mature adult-like network but instead represent a transient state ⁵⁹ preceding the emergence of mature activity patterns. Its timing, correlated with animal weight, suggests it may represent a conserved developmental milestone. Hence, comparative studies across brain regions and species, including potential parallels in human development, would be valuable.

We believe that combining Track2p tracking and the analysis methods described here provides a template for future investigations of more complex developmental phenomena, such as the emergence of sensory representations and cognitive functions. Also, longitudinal imaging uniquely enables investigation of activity-dependent development, including how early activity patterns predict later functional properties and assembly formation of the same cells. Indeed, there is nothing precluding further studies to continue tracking the activity of the same cells until adulthood. Additionally, this approach opens possibilities for targeted manipulation studies to examine how early perturbations affect subsequent circuit development. The ability to track the same cells throughout early postnatal development should thus open doors to entirely new classes of experiments not possible before, providing deeper mechanistic insights into developmental principles and pathologies. This is even more important considering that alterations of developmental trajectories at early postnatal time points are starting to be pointed out as the roots for many developmental disorders ^5,60.

Data and code availability

All Track2p code (algorithm and GUI) is available at: https://github.com/juremaj/track2p, with more extensive documentation and demos available at: https://track2p.github.io/home.html. The code used for running downstream analyses and the data used in this study will be made available upon publication.

Methods

Data acquisition

Animals

All experimental procedures were approved by the French ethics committee (Ministère de l’Enseignement Supérieur, de la Recherche et de l’Innovation (MESRI); Comité d’éthique CEEA-014; APAFiS # 30716-2021032215171216 v8 ) and conducted in agreement with the European Council Directive 86/609/EEC. GAD67-Cre mice were kindly donated by Dr. Hannah Monyer (Heidelberg University). Mice were bred and stored in an animal facility with room temperature (RT) and relative humidity maintained at 22 ± 1°C and 50 ± 20%, respectively. Mice were provided ad libitum access to water and food. A total of 6 mice were used in the study, all heterozygous GAD67-Cre transgenics.

Virus injections

We performed viral injections at P0 as previously described ⁶¹. Briefly, we prepared a solution of AAV-hSyn-GCaMP8m and AAV-FLEX-tdTomato (2:1 volumetric ratio, 10¹² genome copies per millilitre; Addgene) with a small volume (10:1) of 0.05% trypan blue (T8154 Sigma) to verify the success of the injection. We then briefly anaesthetised the mouse on ice and injected 2μL of the solution in the right lateral ventricle using a glass micropipette. Pups were left to recover on a heating pad at 37°C before being returned to the dam.

Cranial window surgery

All procedures were performed as in ⁶. Surgeries were performed at P7 for all mice. Briefly, betadine and lidocaine cream were applied topically, covering the area of the intended incision. Isoflurane was used for induction of global anaesthesia and maintained via a nose cone throughout the procedure. A heating pad was used to maintain body temperature. After skin removal a head plate (4mm inner diameter, Luigs and Neumann) was fixed to the part of the skull covering the barrel cortex using glue (SuperGlue3) and Super Bond (DSM Dentaire). A craniotomy was performed within the head plate opening before finally applying a thin layer of Kwik-Sil (WPI) to the surface of the dura and covering it with a glass cover-slip (3mm, Warner Instruments). The cover-slip was last fixed to the headplate and the skull again using glue and Super Bond. Mice were left to recover on a heating pad at 37°C for at least 1 hour before returning them to the home cage.

Chronic 2-photon calcium imaging

Longitudinal two-photon calcium imaging was performed for each mouse for at least 6 consecutive days (see Fig. 5 for details). Imaging was performed using a Bruker (Ultima 2P) microscope with a Coherent Mai-Tai excitation laser (950 nm excitation light) and a 16x Nikon objective (NA 0.8). Before detection, emitted light was split into two optical paths each associated with a filter (red and green, 580-620 and 500-550 nm respectively), allowing us to simultaneously record the GCaMP8m and the tdTomato signal. The acquisition was performed using the Prairie view software. All recordings were performed in layer 2/3 (depth between 100-200 µm from the pial surface) with a 720 x 720 µm field of view and 512 x 512 pixel resolution. Imaging rate was 30 Hz (resonant scanner) and each session lasted 20 minutes. All experiments were performed in the dark, under sensory-minimised conditions, with mice being free to spontaneously run on a non-motorised treadmill (Luigs and Neumann). To facilitate alignment and cell tracking, we kept the alignment of the head mount with respect to the microscope fixed across all sessions for the same mouse. To record from the same region across days we manually aligned the imaging plane in x, y and z to best match the reference images of the red channel (tdTomato) from all previous recording days. During the course of each recording, pups were kept warm by a heating element mounted to the imaging setup. After each imaging session, the pups were returned to their home cage with the dam and their littermates.

Videography

All videography was performed using a Basler camera (Basler ACE2 1920), with an infrared LED light source (ThorLabs 850 nm) pointed at the mouse. Videos were recorded at 30 Hz, with the microscope acquisition acting as a trigger for camera frame acquisition, also allowing for simple synchronisation across the two modalities. All acquisition was done using custom python scripts using the PyPylon wrapper for the Basler camera software (Pylon Camera Software Suite, https://github.com/basler/pypylon).

Processing of calcium imaging data

Calcium imaging data was preprocessed using the Suite2p pipeline, sequentially performing motion correction, ROI detection, signal extraction and spike deconvolution ²³ for each recording separately. Suite2p additionally provides a cell classification feature, providing a probability of a classifier categorising an ROI as being a true cell. We considered all ROIs above the default threshold of 0.5 as true cells. We used baseline corrected fluorescence traces as our dF/F (using the default Suite2p parameters) for all subsequent analyses.

Preprocessing videography

We used the global movements of the mouse as a proxy of its arousal state. Similarly as in ⁶² we quantified these by looking at the pixel-wise difference of consecutive frames in the videography recording, by squaring this difference and summing the values across all pixels. This yielded a scalar value quantifying the motion of the mouse at each time point, which was used for all subsequent analyses.

Cell tracking

Image and ROI registration

As explained in the main text, Track2p aligns the ROIs based on mean FOV image registration. The implementation allows the user to choose which channel to use for computing the transformation (‘anatomical’ or ‘functional’ channel), as well as which types of transformation to estimate (rigid or affine). Unless otherwise stated, we used the ‘anatomical’ (tdTomato in GAD67-Cre expressing cells) channel with affine transformation (referred to as the ‘baseline’ condition). For the purposes of evaluation we additionally ran tracking using the functional channel with affine registration (‘functional’ in Fig.4) and using the anatomical channel with rigid registration (‘rigid’ in Fig.4). Once registering the images, we applied the same transform to all ROIs from the subsequent session that were considered as cells (see main text). All Track2p image registration is done using the itk-elastix toolbox ⁶³.

Cell matching

Cell matching was done in a way to maximise the sum of the intersection over union (IoU) across all matches of ROIs between two sessions. For this reason we defined a cost matrix M with entries corresponding to:

Where r_{s_k, i} is the i-th ROI in session s_k and r’_{s_(k+1), j} is the j-th ROI in the subsequent session s_k+1 after registering it to the coordinate system of s_k. Assigning matches between the two sets of ROIs in this way corresponds to a linear sum assignment problem, with the goal to find a matrix X that yields the optimal assignment cost:

Where X_i,j equals 1 if and only if ROI r_s,i is assigned to r_s+1,i and 0 otherwise. Additionally, since the number of ROIs is not necessarily the same across two sessions, M and X are not necessarily square. In this case, if there are more rows than columns, then not every row needs to be assigned to a column, and vice versa. We solve this problem by using the algorithm described in ⁶⁴ and implemented in SciPy.

Additionally, some of the matches can be putative false positives. To avoid this issue, we use an approach similar to the one described in ²¹, using a statistical threshold to filter matches based on their IoU distribution. To do this we use the Skimage implementation of Otsu’s method applied to the IoU histogram ²⁴.

Generating a ground truth dataset

Since manually tracking all cells would require prohibitive amounts of time (see Discussion), we decided to generate sparse manual annotations, only tracking a subset of all cells from the first recording day onwards. To do this, we took the first recording (s₀), and we defined a grid of 8 x 8 (64) equidistant points over the FOV and, for each point, identified the closest ROI in terms of euclidean distance from the median pixel of the ROI (see Extended Data Fig. 3A). We then manually tracked these 64 ROIs across subsequent days. The manual tracking was done using the Suite2p GUI, by visualising the FOV and masks from two recordings side by side and choosing the matching ROI from the subsequent recording or terminating the track if we could not find a match (see Extended Data Fig. 3B). Only neurons that were detected and tracked across all sessions were taken into account and referred to as our ground truth dataset (‘GT’ in Fig. 4). When comparing the GT to Track2p tracks, we only considered Track2p tracks that originated from one of the 64 ROIs chosen for evaluation initially. All manual tracking was performed blind to the Track2p outputs.

Tracking evaluation metrics

For the evaluation of Track2p under different conditions and comparison to CellReg, we used the ‘Complete tracks’ metric (^20,25, Fig. 3B), defined as:

Where T_rc is the number of perfectly reconstructed tracks, T_c is the number of total tracks identified by Track2p for the 64 s₀ ROIs chosen for manual tracking and T_gt is the number of total tracks in the ground truth dataset (from the same 64 s₀ ROIs). In the case of perfect tracking, CT will be equal to 1 (where all computed tracks are equivalent to the ground truth tracks (T_rc=T_c=T_gt)). Conversely, in the case of failed tracking, the value of CT will be close to 0 (for example when there is a poor match between ground truth and reconstructed tracks (T_rc≈0) or when there are many falsely reconstructed tracks not present in the ground truth (T_c>>T_gt)). The CT metric is mathematically equivalent to an F1 score where true positives are defined as perfectly reconstructed tracks (Fig. 4A row 1), false negatives as tracks from GT without a match in Track2p (Fig. 4A row 2) and false positives as tracks from Track2p without a match in GT (Fig. 4A rows 3, 4 and 5).

In the final step of evaluation we looked at when the algorithms from the original evaluation failed in their tracking of the ground truth neurons. For this we used a metric of the proportion of correctly reconstructed tracks for an increasingly longer number of sessions (s) (referred to as ‘Prop. correct’ in Fig 4C-F):

Where T_gt is equivalent as above and T_rc(s) is equivalent as T_rc but for full tracks up to session s. Since in this case we kept the GT consistent with the initial evaluation, including only the neurons identified across all days, this metric is agnostic to possible false positives at shorter time epochs.

To facilitate evaluation by individual users under different experimental conditions, we provide a helper script (Jupyter notebook) aiding the whole evaluation process by both defining a grid of cells to manually track (see previous section), as well as to compute the tracking quality metrics once the ground truth dataset is completed.

CellReg Tracking

To compare Track2p with CellReg, we ran the MATLAB implementation of CellReg tracking (https://github.com/zivlab/CellReg) according to the provided user manual. We first attempted tracking using the default rigid registration, which yielded an error, terminating the algorithm and indicating that subsequent sessions do not resemble the reference session and suggesting to use non-rigid registration. We then re-launched the algorithm with non-rigid registration and used those results to evaluate the tracking the same way as for Track2p (see Results). Notably, even when running CellReg using non-rigid registration, we noticed that the algorithm did not find any tracks spanning all days (T_rc = 0), which explains the CT score of 0 for all day evaluation.

Functional properties of tracked neurons

Calcium event rates

To quantify the calcium event rate statistics we used SciPy’s peak detection algorithm (Python). For this we first denoised the traces slightly by averaging using a bin size of 10 frames and then proceeded to detect any peaks with a height and prominence of at least one standard deviation. We then calculated the rate as the number of peaks per minute within the recording. When quantifying the stability of the rates across days we computed the Pearson correlation coefficient across all neurons for each possible combination of sessions recorded from the same mouse.

Pairwise correlations and PCA

We quantified pairwise correlations by computing the Pearson correlation between the full traces of all pairs of simultaneously recorded neurons within each session. To quantify the spatial properties of pairwise correlations, we additionally calculated the euclidean distance between ROI centroids for each corresponding pair of neurons. When plotting the pairwise correlations as a function of pairwise distance we averaged in bins of 30 µm. The correlation of neighbouring neurons was estimated as the intercept of an exponential fit to the full data.

When comparing the stability of the correlation structure across days, we refer to this as ‘FC stability’, which we calculate as the Pearson correlation between all pairs of original correlation matrices for a given mouse. When fitting a sigmoid using the first day as the reference, we fix the upper limit of the sigmoid to within-day FC stability for that day and leave the other parameters free. We calculate the inflection point (‘transition age’) and the linear portion (‘transition time’) of the sigmoid using the extrema of the first and fourth derivatives respectively.

For principal component analysis (PCA) we used the scikit-learn implementation (Python).

Decoding

All decoding was done using linear regression with ridge regularisation (ridge regression) to avoid overfitting given the large number of neurons. Ridge regression optimises the weights (β) that minimise the following loss function:

Where in our case y is behavioural data, X is neural data and λ is a regularisation parameter. To choose the optimal λ and to reliably estimate the model performance on same-day decoding we used nested cross-validation (see for example ^65,66), ensuring data efficiency and a splitting between training, validation and test data. Briefly, a grid search for the optimal λ is performed by using the train set to fit models with different λ values and choosing the value corresponding to the model with the best performance on the held out validation set (model selection). This model is then evaluated on the test data that was used neither in model fitting nor in the hyperparameter optimization (model evaluation). To ensure data efficiency, nested cross-validation repeats this procedure for all possible splits of the data, with an outer cross-validation loop used for model evaluation and an inner loop for model selection. We used 5 fold splits for both the inner and outer loops, splits were done on consecutive 2 minute blocks of the recording. For cross-day decoding we fit a new model with the optimal lambda for that given day and evaluated it on all other days. For all decoding analysis we slightly denoised the dF/F as well as the behaviour traces by averaging in bins of 10 consecutive timestamps. All models were implemented and fitted using PyTorch (Python).

Extended Data Figures

Extended Data Figure 1.

Extended Data Figure 2.

Extended Data Figure 3.

Extended Data Figure 4.

Extended Data Figure 5.

Extended Data Figure 6.

Extended Data Figure 7.

Supplementary Information: Software

Interactive visualisation and curation using the Track2p GUI

To facilitate the ease of as well as to visually evaluate and curate the quality of Track2p tracking, we’ve designed an intuitive graphical interface (GUI) that makes it easy for users to launch the Track2p algorithm and interact with the outputs. Track2p GUI allows the users to perform cell tracking as well as generate some preliminary visualisations of cell activity without requiring any programming skills on behalf of the user (Extended Data Fig. 1A). Additionally, we’ve made the algorithm and the GUI simply installable via pip, and provide plenty of documentation on the project’s website: https://track2p.github.io/.

Track2p provides several data visualisation functionalities. Firstly, while running the algorithm, a set of figures is generated that can be used to monitor the progress at each consecutive step, as well as to provide a visual overview of tracking quality, for example based on a visual assessment of registration quality, IoU histograms or zoomed-in images of example tracked cells. The GUI itself provides an interactive view of the FOV overlaid with ROIs tracked across all days (Extended Data Fig. 1B, similar to Fig. 3D). The user can then select a cell from this window that is then displayed in a zoomed in view across all days (Extended Fig. 1B, similar to Fig. 3C). The visualised cell can be selected interactively by selecting it from contours displayed over the whole FOV. Based on these visualisations, users can also assess the tracking performance and, if necessary, manually curate the results by manually identifying cases of incorrect tracking (Extended Fig. 1E). Under conditions shown in this paper manual curation should not be strictly necessary (see Fig. 4 for evaluation without manual curation), however this will depend greatly on the particularities of the dataset used.

Additionally, the GUI also provides activity visualisation, plotting the activity time series of the selected ROI for each recording (Extended Fig. 1C, similar to Extended Data Fig. 4B), while also offering the visualisation of the activity of the whole population across all days by generating raster plots of neuronal activity (Extended Fig. 1F, similar to Fig. 5A and Extended Data Fig. 4A). These visualisations can allow the user to visually inspect changes to activity statistics at the single cell level as well as changes to population level phenomena such as synchronisation (assemblies) and sequences at the level of all tracked neurons. Raster sorting capabilities can also give a visual insight into the stability of these phenomena across recordings. Currently the user can choose between PCA (Principal Component Analysis) or t-SNE (t-distributed Stochastic Neighbor Embedding) techniques to re-sort the rows of the raster based on dimensionality reduction. These methods aim to sort cells with similar fluorescence traces in adjacent rows of the raster plot, visualizing co-varying activity patterns at the population level. Additionally, the user can choose between individual sorting and sorting that is preserved across days. Individual sorting better shows these population phenomena for each specific day, however in this case the correspondence between the cells across days is lost due to independent sorting. On the other hand, across-day sorting allows the user to focus on one day to visualize the stability of these population phenomena and see if they are preserved across other days, since the sorting of the neurons remains consistent across sessions.

Acknowledgements

We thank all the members of the Cossart lab for helpful discussions and constructive feedback. We thank Drs. Lorenzo Fontolan and Michel Picardo for critical feedback on the manuscript. We thank INMED’s animal facility and PBMC technological platform for excellent technical support. This work was supported by the European Research Council under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 951330), by the Fondation Bettencourt Schueller and the Fondation Roger de Spoelberch; R.C. is supported by CNRS. J-C. P. is supported by INSERM. The project leading to this publication has received funding from the French government under the “France 2030” investment plan managed by the French National Research Agency (reference : ANR-16-CONV000X / ANR-17-EURE-0029) and from Excellence Initiative of Aix-Marseille University - A*MIDEX (AMX-19-IET-004).

Additional files

Supplementary Video 1. Example of videography and 2-photon calcium imaging data acquired throughout the second postnatal week. Note the substantial growth of the mouse, the development of the whisker pad and opening of the eyes; and the loss of synchrony and spatial organisation in the calcium imaging data. For visualisation purposes, the videos were manually aligned using Fiji.