Long-term stability of neuronal ensembles in mouse visual cortex

Coactive neuronal ensembles are found in spontaneous and evoked cortical activity and are thought to participate in the internal representation of memories, perceptions, and mental states. In mouse visual cortex, ensembles can be optogenetically imprinted and are causally related to visual percepts, but it is still unknown how stable they are over time. Using two-photon volumetric microscopy, we performed calcium imaging over several weeks of the same neuronal populations in layer 2/3 of visual cortex of awake mice, tracking over time the activity of the same neurons in response to visual stimuli and under spontaneous activity. Only a small number of neurons remained active across days. Analyzing them, we found both stable ensembles, lasting up to 46 days, and transient ones, observed during only one imaging session. The majority of ensembles in visually-evoked activity were stable, whereas in spontaneous activity similar numbers of stable and transient ensembles were found. Among stable ensembles, more than 60 % of neurons still belonged to the same ensemble even after several weeks. These core ensemble cells had stronger functional connectivity than neurons that stopped belonging to the ensemble. Our results demonstrate that spontaneous and evoked neuronal ensembles can last weeks, providing a neuronal mechanism for the long-lasting representation of perceptual states or memories.

In spontaneous activity, approximately ~50% of ensembles were stable whereas during evoked 47 activity ~70% were stable. We also analyzed neuronal participation within stable ensembles over 48 time, finding stable neurons that belong to the same ensemble through days; neurons that stopped 49 responding or became part of other ensembles; and new neurons that enrolled into an existing 50 ensemble. Functional connectivity analysis revealed that stable neurons were more connected 51 than neurons which were eventually lost. Our results reveal long-term stability over several 52 weeks of ensembles constituted mainly by neurons that are more functionally connected. 53 54 55

57
Experimental and analysis rationale 58 To examine ensembles stability under visually-evoked and spontaneous activity we 59 performed two-photon calcium imaging in layer 2/3 of the visual cortex from transgenic mice 60 (GCaMP6s, n = 4; and GCaMP6f, n = 2) through a cranial window. We head-fixed a mouse in SEM) active neurons during spontaneous and evoked activity imaging periods, respectively. 93 While the number of active neurons was similar in days 1-10, a significant decrease in the 94 number of active neurons occurred in day 43-46, in both spontaneous and evoked activity 95 ( Figure 1H; 51 ± 10 neurons; p = 0.023 and p = 0.013, respectively). This decrease could be due 96 to several factors, including diminishing quality of transgene expression and repeated 97 experimental procedures on the same cortical territory, for example animal manipulation, 98 surgical attachments, microtraumas and laser exposures. On the other hand, the percentage of 99 time that a neuron was active (i.e. the number of frames with activity of a given neuron) did not 100 change, and was 15.5 ± 0.5 % and 15 ± 0.6 % (mean ± SEM), for spontaneous and evoked 101 activity respectively on day 1, with no significance difference (one exception) on the following 102 days ( Figure 1I). This indicated that, on an individual neuron basis, the level of neuronal activity 103 remained similar across days. Thus, if a neuron is active on any given day, this activity is similar. 104 We then explored if active neurons were repeatedly active across days, and analyzed if 105 the neurons that were active on the first session of the first day were also active in later sessions, 106 up to 46 days later ( Figure 1J). A minority of neurons were active across days, and their 107 proportion became reduced over time ( Figure 1K). Even within day 1, only 42 ± 2 % and 37 ± 4 108 % (mean ± SEM) of neurons were active in subsequent imaging sessions, 5 minutes later, during 109 spontaneous and evoked activity (gray on Figure 1K). This low number could partly be 110 explained by the fact that not all active neurons were captured in a 5 minute interval (Figure 1 -111 Figure supplement 1) but also likely are due to the possibility that a significant number of 112 neurons were only active during part of the time. Consistent with this, across days, the 113 percentage of neurons that remained active continued to decrease, with significant decrements on 114 day 10 and 43-46, in both spontaneous and evoked activity (gray on Figure 1K). In spite of this 115 generalized decrease in the number of neurons that remained active, during spontaneous and 116 evoked activity, we found 17 ± 7 and 22 ± 8 (mean ± SEM) neurons on day 1 that were still 117 active 46 days later. We concluded that the neurons that are activated by visual stimuli or 118 spontaneously changed over time, with only a small minority of cells maintaining their responses 119 across weeks. 120 121

Neuronal ensembles identification based on functional connectivity
We then explored if neuronal ensembles were present in the dataset. Given that the 123 majority of the neurons did not remain active across different days, to perform the neuronal 124 ensemble analysis we only analyzed neurons that remained active across sessions. To do so, we 125 built binary raster plots of neurons that remained active across sessions ( Figure 2A) and detected 126 ensembles using their functional connectivity (Pérez-Ortega et al., 2016). Specifically, to identify 127 a significant coactivation between every pair of neurons we first generated 1,000 spike raster 128 surrogates by a random circular shift in time of the active frames ( Figure 2B left). Then we 129 tabulated how many times a given pair of neurons were coactive by chance, and used a 95 % 130 threshold on the cumulative probability from surrogate coactivations ( Figure 2B right). A 131 functional connection was assigned if the coactivations were above chance level (p < 0.05). This 132 process was done for every pair of neurons to build a functional network graph ( Figure 2C). We 133 then filtered the raster plot removing the activity without significant coactivation and kept only 134 frames with 3 or more active neurons ( Figure 2D). Then, treating each frame as a vector (1 frame 135 bin = 81 ms), we computed the Jaccard similarity between every pair of vectors ( Figure 2D  Using this approach, we identified an average of 4.57 ± 0.14 and 4.63 ± 0.14 (mean ± 149 SEM) ensembles on spontaneous and evoked activity, with no significant differences between 150 them ( Figure 2 -Figure supplement 1). We then inquired if ensembles were preserved across 151 days. To do so, we defined as the same neuronal ensemble across days if at least 50 % of its 152 neurons were preserved. This 50% criteria was chosen solely as an operational definition of stability. In comparisons made between 2 imaging sessions from either spontaneous or evoked 154 activity, we found that some ensembles were preserved but others were not ( Figure 3A-B). We 155 termed "stable" the ensembles found on day 1 which were preserved in subsequent days, and 156 "transient" all other ensembles. From here on, we focused our analysis on stable ensembles 157 ( Figure 3 -Videos 1-2). On day 1, stable ensembles constituted 54 ± 3 % (mean ± SEM) of all 158 ensembles during spontaneous activity, but interestingly were 72 ± 4 % (mean ± SEM) of 159 ensembles in evoked activity ( Figure 3C). Similar trends were observed on days 2 and 10. On 160 day 43-46, stable ensembles were around ~50 % of all ensembles during spontaneous or evoked 161 activity ( Figure 3C). This suggested that evoked ensembles were more stable across days 1-10, 162 but are similarly stable as spontaneous ones after 43 days. 163 To test if stable ensembles were a merely artifact of the analysis, we shuffled the 164 neuronal activity of the likened session on days 1, 2, 10, 43 and 46. We found 0.8 ± 0.1 (mean ± 165 SEM) stable ensembles from shuffled activity compared with 2.8 ± 0.1 (mean ± SEM) stable 166 ensembles from original data (p = 1x10 -66 ; Figure  We were curious to find if stable and transient neuronal ensembles had different 172 properties. There were no significant differences between stable and transient neuronal 173 ensembles during spontaneous or evoked activity in the number of neurons ( Figure 3D) nor 174 network ensemble density ( Figure 3E). However, ensemble robustness, defined as the product of 175 ensemble duration and similarity of its activity (see Methods), was significantly higher in stable 176 than in transient neuronal ensembles ( Figure 3F). This was the case for both spontaneous and 177 evoked stable ensembles, indicating that stable ensembles had higher consistency in their 178 activity. 179 Finally, we studied neural identity or functional connectivity of stable ensembles. We 180 described single neuron participation on stable ensembles ( Figure 4A-D). Approximately 50 % 181 of neurons belonged to only one stable ensemble (single neurons), and less than 20 % of the 182 neurons belonged to more than one (shared neurons), while the rest of the neurons were not part 183 of any stable ensemble during spontaneous or evoked activity ( Figure 4E). We then inquired what happened to individual neurons of stable ensembles across days. More than 60 % of 185 neurons of a stable ensemble observed on day 1 remained in the same stable ensemble on a 186 second session of day 1, 2, 10 and 43-46, during spontaneous and evoked activity ("maintained" 187 neurons, Figure 4F). The rest of neurons (less than 40 %) changed their participation to another 188 ensemble or stopped participating in detectable ensembles ("lost" neurons). Interestingly, in 189 subsequent sessions, we found "new" neurons joining stable ensembles in a similar proportion as 190 lost neurons. Neither maintained nor lost neurons specifically belonged to one or more than one 191 ensemble ( Figure 4G). However, functional connection density between maintained neurons 192 from the same ensemble (0.71 ± 0.01, mean ± SEM) was significantly higher than density from 193 lost neurons during spontaneous or evoked activity (0.35 ± 0.02, mean ± SEM, p = 7x10 -55 , 194 Figure 4H). Therefore, poor connectivity between lost neurons could be a reason for them not 195 being durable, and high functional connectivity indicates a possible lasting stability. with the hypothesis that they are the same and that sensory stimulus reactivate existing 212 ensembles. 213 Volumetric two-photon calcium imaging allowed carefully tracking the same neurons and 214 identifying neuronal ensembles across days. Nevertheless, our experiments had limitations, as imaging data of animals became noisier across time goes ( Figure 1G), likely due to accumulated 216 experimental issues, which made us difficult to track every neuron with sufficient signal to noise 217 (PSNR > 18 dB, Figure 1H). Nevertheless, neurons tracked for weeks remained at the same level 218 of activity as on the first day ( Figure 1I). Interestingly, even after 5 minutes of recording, only 219 ~50 % of previously active neurons continued to be active. This, and the fact that lost neurons in 220 ensembles are replaced by new ones (Figure 4H enabled. We used two protocols to display in the monitor. The first was in absence of visual 267 stimulation, the monitor was displaying a static blue screen, we used it to record spontaneous 268 activity during 5 min per session. The second protocol was for visual stimulation consisting of a 269 full sinusoidal gratings (100 % contrast, 0.13 cycles/deg, 5 cycles/s) drifting in a single direction 270 per mouse (0º or 270º) presented for 2 s, followed by a random amount between 1 to 5 s of mean 271 luminescence. The visual stimulus is presented 50 times during 5 min per session. We performed 3 consecutive sessions (5 min apart) per protocol per day of experiment. See Figure 1 - Table 1  273 for a detailed sessions recorded per mouse. 274

Volumetric two-photon calcium imaging 275
Imaging experiments were performed from 20-150 days after head-plate procedure. Each mouse 276 was placed on a treadmill with its head fixed under the two-photon microscope (Ultima IV, 277 Bruker). Animals were acclimated to the head restraint for periods between 5 to 15 min for at 278 least 2 days, and exposed to visual stimulation sessions before the recordings presented here. The 279 imaging setup was completely enclosed with blackout fabric to avoid light contamination leaking 280 into the PMT. An imaging laser (Ti:sapphire, λ = 920 nm, Chameleon Ultra II, Coherent) was 281 used to excite a genetically encoded calcium indicator (GCaMP6s or GCaMP6f). The laser beam 282 on the sample (30-60 mW) was controlled by a high-speed resonant galvanometer scanning an 283 XY plane (256x256 pixels) at 17.7 ms (frame period) covering a field of view of 312x312 µm 284 using a 25X objective (NA 1.05, XLPlan N, Olympus). An electrically tunable lens (ETL) was 285 used to change the focus (z-axis) during the recording. We recorded consecutively 3 planes at 286 different depths (-5, 0 and 5 µm from a reference z-axis) waiting 9.3 ms between planes for ETL 287 to stabilize the focus. Thus, we collected three frames, one per depth, every 81 ms for 5 min 288 (single session, 3,704 frames per plane). Imaging was controlled by Prairie View and ETL was 289 synchronized using a DAQ (USB-6008, NI) controlled by a custom-made app on MATLAB 290 (https://www.mathworks.com/matlabcentral/fileexchange/78245-etl-controller-for-volumetric-291 imaging). 292

Recording same neuronal region through days 293
In the first day of experiment, we recorded the vascularization of the pia at 10X and 25X using 294 bright-light microscopy. We fixed the depth to 140 µm from pia to record a reference image 295 (calcium imaging) and a second reference image of the center of the field of view using an extra 296 4X optical zoom (day 1 in Figure 1G). We carefully preserved unchanged the position of the 297 microscope and the base we placed the mice. For following days of recording, we looked for 298 matching the reference image of vascularization at 10X, then 25X. After that, we looked 140 µm 299 depth from pia trying to match the reference image on x and y axis, then we used a 4X extra 300 optical zoom to finely match the second reference image on z axis (day 2, 10 and 43 or 46 in 301 Figure 1G). We performed a volumetric imaging to record 3 planes: one reference plane, one 5 µm above and one 5 µm below in order to emend tilt, which in our system does not go beyond 5 303 µm. We extracted the maximum intensity projection from the 3 planes resulting in a single video 304 for each session ( Figure 1D). 305

Neuronal activity extraction 306
We used a custom-made graphical user interface (GUI) on MATLAB 307 (https://github.com/PerezOrtegaJ/Catrex_GUI) to extract the binary raster activity from every 308 single session video (5 min, 3,704 frames). First, we performed a non-rigid motion correction 309 taking as a reference the mean of the 185 frames (5 %) with less motion artifacts. Then, we 310 searched the ROIs with a modified version of Suite2P algorithm (Pachitariu et al., 2016, Figure  311 1E). ROIs were preserved if they fulfill the following criteria (fixing radius to 4 µm): 0.5*π 312 *radius^2 < area < 4*π*radius^2; roundness > 0.2; perimeter < 3.5*π*radius; eccentricity < 0.9 313 and overlapping < 60 %. Calcium signal from each ROI was extracted measuring the changes in 314 fluorescence with respect of its local neuropil (F raw -F n ) / F n , where F raw is the signal from the 315 ROI and F n is the signal of its local neuropil 10 times ROI radius. ROI local neuropil is not 316 including signal from ROIs if presented within the area. Then we computed the peak signal-to-317 noise ratio PSNR = 20 · log10( max( F raw -F n ) / σ n ), where max represents a maximum 318 function, and σ n represents the standard deviation of the local neuropil. We evaluated the ROIs 319 again keeping them if PSNR > 18 dB. Then we smoothed the calcium signal with a 1 s window 320 average to perform a spike inference using the foopsi algorithm (Friedrich & Paninski, 2016). 321 We binarized the spike inference signal, placing 1 if there was spikes inferred and 0 if not. We 322 placed all binarized signals from every ROI in a N x F raster matrix, where N is the number of 323 active neurons and F the number of frames. This matrix is visualized as a raster plot, where ones 324 in the matrix are the dots representing the active frames of the neurons ( Figure 1F). 325

Tracking neurons across days 326
We computed a rigid and then a non-rigid motion correction between the binary image of the 327 ROIs shape between a single session of day 1 and a single session from day 1, 2, 10 and 43 or 46 328 ( Figure 1J). Then, we looked for the intersection (in pixels) between ROIs of the neurons from 329 two sessions (intersection > 0.5* π *radius^2) and evaluate the Euclidean distance between 330 centroids of the ROIs intersected keeping it if the distance < radius. We used the raster matrix 331 only with the tracked neurons between sessions. 332

Identification of neuronal ensembles based on functional connectivity 333
To analyze neuronal ensembles from raster activity we used a custom-made GUI on MATLAB 334 (https://github.com/PerezOrtegaJ/Neural_Ensemble_Analysis).
Functional connectivity 335 represents the significant coactivity between every pair of neurons from a raster matrix. The 336 number of coactivations Co ab between neuron a and b was computed counting in how many 337 single frames were both neurons simultaneously active. To identify significance, we generated 338 1,000 surrogates of neurons a and b by random circular shifting their activity in time to disrupt 339 their temporal dependency. We counted the number of surrogate coactivations S ab,i in each 340 iteration i, building a cumulative distribution of S ab selecting a threshold T ab of coactivations at 341 the 95%. If the actual number of coactivations Co ab is above threshold T ab we put a functional 342 connection between neuron a and b ( Figure 2B). Doing this with every pair of neurons we got a 343 functional neuronal network, where every node is a neuron and every link represents a 344 significant coactivity between them ( Figure 2C). We used the functional connectivity to rebuild 345 the raster matrix to keep the significant coactivity of the neurons. To do so, we identified the 346 active neurons of every single frame, we looked to their functional connectivity, if a neuron has 347 no connection its activity was removed from that frame. At the end, we also removed the frames 348 with less than 3 coactive neurons ( Figure 2D). Then we computed the Jaccard similarity between 349 all single frames (column vectors) of the rebuilt raster matrix. A hierarchical clustering tree with 350 single linkage was obtained to identify the more similar vectors by keeping the branch with more 351 than 2/3 of Jaccard similarity ( Figure 2E). With the more similar vectors we performed a 352 hierarchical clustering with Ward linkage and grouped based on a contrast index (Beggs & 353 Plenz, 2004). Each group of column vectors is the activity of the neuronal ensemble ( Figure 2F), 354 i.e. the raster matrix E j of an ensemble j of size N x F j , where N is the number of neurons and F j 355 is number of frames where the ensemble j was active. 356

Demixing neuronal ensemble identity 357
The neuronal ensemble activity E j was used to identify the participation of each neuron in the 358 ensemble j. We computed the functional connectivity similarly as described above, but 359 incorporating the correlation of the neurons between the times where the ensemble was active. 360 To do so, we got a binary vector V j representing the times were the ensemble j was active (1) or not (0). Vector V j was of size 1 x F, where F is the number of frames of the session. A Pearson 362 correlation coefficient P j,a between vector V j and the activity of neuron a was computed. Then 363 we got an ensemble weight W j,ab between neurons a and b in ensemble j, which integrates their 364 correlation with the ensemble j and their number of coactivations as follows W j,ab = P j,a · P j,b · 365 Co ab . To identify significance, we generated 1,000 surrogates of neurons a and b shuffling their 366 activity as described before, and assigning randomly a value from the correlation with the 367 ensemble j (P j ). Then we compute the surrogate weight SW j,ab,i in each iteration i, building a 368 cumulative distribution of SW j,ab selecting a threshold TW j,ab of coactivations at the 95%. If the 369 actual ensemble weight SW j,ab is above threshold TW j,ab we put a functional connection between 370 neuron a and b. A neuron is considered to be part of an ensemble if at least had one single 371 functional connection ( Figure 2G). A neuron could be part of more than one neuronal ensemble 372 (shared), only one ensemble (single) or in any ensemble (not participant). 373

Comparing neuronal ensembles between following sessions 374
Taking the first session on day 1 and a second session from day 1, 2, 10, 43 or 43, we got the 375 raster matrix from each session with only tracked neurons (same active neurons in both 376 sessions). Neuronal ensembles were extracted from each raster matrix independently. If there 377 were 50 % or more neurons in an ensemble j from the first session belonging to an ensemble k 378 from the second session we called it "stable" ensemble, otherwise we called it "transient" 379 ensemble. 380

Ensemble measures 381
Ensemble network density, fraction of present functional connections to possible connections 382 within an ensemble. Ensemble robustness, we introduced here as robustness = similarity · 383 activity, where similarity is the average of the Jaccard similarity between every pair of column 384 vectors of the ensemble matrix raster E j , and activity is the fraction of ensemble active frames to 385 the total frames of the session. The higher the value, the higher the robustness. Stability of 386 neurons: given a stable ensemble, a "maintained" neuron participated during the first session on 387 day 1 and during a second session on the following days. A "lost" neuron participated only in the 388 first session but not in a second session, and a "new" neuron did not participate in the first 389 session but participate in a second session. The fraction is based on total neurons in an ensemble 390 from day 1. Promiscuity of neurons: "shared" neurons is the fraction of neurons participating in 391 more than one ensemble; "single" neurons participate in only one ensemble; and "not participant" neurons do not belong to any ensemble. Tuned neurons: we consider a tuned neuron 393 if its number of active frames during visual stimulation was significantly higher than its number 394 of active frames during periods with no visual stimulation (P < 0.5, t-test). 395