Quantifying network behavior in the rat prefrontal cortex

The question of how consciousness and behavior arise from neural activity is fundamental to understanding the brain, and to improving the diagnosis and treatment of neurological and psychiatric disorders. There is significant murine and primate literature on how behavior is related to the electrophysiological activity of the medial prefrontal cortex and its role in working memory processes such as planning and decision-making. Existing experimental designs, specifically the rodent spike train and local field potential recordings during the T-maze alternation task, have insufficient statistical power to unravel the complex processes of the prefrontal cortex. We therefore examined the theoretical limitations of such experiments, providing concrete guidelines for robust and reproducible science. To approach these theoretical limits, we applied dynamic time warping and associated statistical tests to data from neuron spike trains and local field potentials. The goal was to quantify neural network synchronicity and the correlation of neuroelectrophysiology with rat behavior. The results show the statistical limitations of existing data, and the fact that making meaningful comparison between dynamic time warping with traditional Fourier and wavelet analysis is impossible until larger and cleaner datasets are available.


Introduction 37
Hodgkin and Huxley's model of neuron action potentials derived from experiments with 38 the squid giant axon 1 was the seminal event in neurophysiology 2-4 . The central dogma 39 of modern neuroscience is that neuron electrochemical activity and connectivity at the 40 microscopic level can provide a clear understanding of complex behaviors. Thus, 41 measuring electrical signals can be sufficient to integrate neural networks' chemical 42 signaling and activity. 43 The visual cortex, and its relationship to visual processing, is one of the better-44 understood neural networks because there are significant spatial correlations in signals 45 that result from local connections between neurons. These spatial autocorrelations 46 result in a high signal-to-noise ratio due to the strong intensity of the electric field 47 produced by local synchronized neuron firing. Additionally, the neuroelectrophysiology 48 of individual cells in the visual pathway has allowed mathematical modeling of the 49 performance of specific neurons 5,6 . The spatial autocorrelations of visual neural activity 50 are sufficiently high that low-resolution, indirect measurements of neuron activity, such 51 as local blood flow 7-9 and extracranial electrodes 10 , are sufficient to reconstruct 52 perceived images. The organization of the visual cortex may be explained by the fact 53 that images projected on the retina are spatially correlated, and therefore biological 54 neural networks have adapted to take advantage of these correlations. 55 From an evolutionary perspective, the visual cortex is an ancient structure, and 56 therefore the neural networks involved in visual processing have optimized over time. In 57 contrast, the prefrontal cortex as part of the neocortex is one of the least developed 58 brain regions, especially through evolutionary analysis of its size in mammals 11,12 . 59 Unlike the visual cortex, the prefrontal cortex does not seem to possess highly local and 60 regular spatial organization 13 . Because of the lack of a regular organization of the 61 prefrontal cortex and the consequent low intensity of electrical signals, 62 electrophysiology is one of the only methods with adequate sensitivity to probe the fine 63 structure of the prefrontal neural network. Unfortunately, such measurements are highly 64 invasive, and experiments are only feasible on laboratory animals such as rats [14][15][16] . To 65 find the relationship between neuron firing and behavior, methods are required to 66 analyze these complex electrical measurements. 67 The characteristics of this neuron network constrain the firing patterns of neurons. 68 Therefore, the methods are required to capture key network characteristics from 69 recorded spike trains and local field potentials (LFPs). Our hypothesis is that these 70 network characteristics are embodied by the synchronicity of neuron firing: neurons that 71 work in concert will be highly (anti-)correlated, while neurons that process independent 72 data will not be correlated. We chose to investigate neuron synchronicity using the 73 dynamic time warping (DTW) method. DTW is an efficient, non-parametric approach to 74 determine the best alignment of two time series, such that the overall shapes are 75 matched. Compared to classical Fourier and wavelet analysis, DTW gives a global 76 alignment of the data that is robust to local variations in timing. We predicted that DTW 77 could discern if two neurons fire synchronously or asynchronously. 78 Thus, we analyzed electrophysiology data collected from rats performing the T-maze 79 task, which evaluates working memory. We use previously published spike train and 80 LFP recordings taken from the rat medial prefrontal cortex (mPFC) 14-16 , along with new 81 data from three additional rats (rats A, B, and C) recorded using published methods 16 . 82 We investigated if DTW would be useful in elucidating the connection between neural 83 activity and behavior. 84

85
In the data we used from the three published studies (Table 1), rat neuron spike trains  86 and/or LFPs were recorded as the rats performed in T-maze (Figure 1-A). The exact 87 experimental details for the three studies differed slightly 14-16 , but their combined 88 scientific goal was to analyze the link between neuron recordings and working memory 89 performance in the T-maze task. We focused on the four-second window of time 90 centered on the moment at which the rat leaves the T-intersection (decision box in 91 Figure 1-A), presumably making its choice. This window contained two seconds of pre-92 decision neural activity (e.g., cognition and decision-making process) and two seconds 93 of post-decision neural activity (e.g., evaluation of reward or lack of reward). Our 94 hypothesis was that neural activity was significantly affected by the independent 95 variables of (1) the timing of the recording relative to the choice, (2) the correctness of 96 the choice, and (3) the sampled neurons (Table 2). 97 We used the DTW method on the neuron spike trains from all three studies and from 98 the local field potentials for Ito et al. 14 and Yang et al. 16 (Figure 1-  We defined an event as a single traversal of the T-maze by a rat. We defined a run as 111 all events recorded on a single day for a single rat. In the Ito 14 and Stout 15 studies, 112 electrodes were moved between rats, and, therefore, a different set of neurons were 113 recorded between runs. We set the level of statistical significance to α =0.05, so that we 114 expect a false positive rate of α . 115 Neural network firing is consistent across a single run 116 Our first question was: when analyzing d crit , can we pool events from a single run? The 117 null hypothesis was that during a single run, we may assume that the electrode 118 placement is unchanged and that the rat basal activity is also unchanged between 119 recordings in the same run, as measured by d crit , The alternative hypothesis is that 120 specific factors affect d crit between events in the run, such as the categorical variables in 121 Table 2. 122 First, we analyzed the normality of d crit with the Shapiro-Wilk test (Supplemental Tables  123  1 and 2). For neuron spike trains, we found that in 96.7% (60 out of 62) of runs, the 124 samples were not normally-distributed with p < 0.05, and for LFPs, we found that 80.0% 125 (20 out of 25) were not normally-distributed. Comparing spike trains to LFP with 126 Boschloo's exact test, these fractions were not significantly different (p=0.615). We 127 performed subgroup analysis on the categorical variables (1) correctness and (2) timing 128 and we found that this proportion was unaffected by further stratification by these 129 variables. Therefore, we concluded that the Kruskal-Wallis non-parametric test was 130 more appropriate than classical one-way analysis of variance (ANOVA), due to violation 131 of the assumption of normally-distributed data and due to the limited number of events 132 per run. 133 we also performed ANOVA and found that 11.2% (7 out of 62) of the spike train runs 137 displayed significant event-to-event variance, regardless of stratification; the 138 corresponding proportion for LFPs was 20% (5 out of 25). We tested if Kruskal-Wallis 139 gave different results as compared to ANOVA using Boschloo's exact test, and found 140 that neither the spike train data (p=0.238) nor the LFPs (p=0.378) were significantly 141 different. 142 Since few of the runs (within our 5% margin of error for false positives from had a significant event-to-event variance in the spike trains, we concluded that we may 144 perform pooled analysis of d crit for the spike trains of events across a single run. 145 Boschloo's exact test on the LFPs also demonstrated a non-significant deviation from 146 the 5% margin of error (p=0.413 for Kruskal-Wallis, p=0.140 for ANOVA). 147 Neural network structure is revealed by the DTW matrix 148 Instead of characterizing network connectivity with a single number (d crit ), we used all 149 the information available in the DTW matrix to determine if firing patterns were 150 consistent, using non-parametric Mantel tests. We performed all pairwise comparisons 151 for events in the same rat to ask the question: given a single run, is there a difference in 152 Pearson correlation coefficient of before vs after depending on if we compare (1) before 153 vs after of the same event to (2) before vs after of different events? We performed 154 pairwise Mantel test comparisons among events in the same run. We performed 155 Kruskal-Wallis and Kolmogorov-Smirnov tests to determine if these distributions of the 156 correlation coefficient and p-value were different (Supplemental Figures 2 and 3,  157 Supplemental Table 3). Overall, we found that there was no significant difference event-158 to-event for the before vs after DTW matrix correlation. 159 When we stratified by correctness (Supplemental Table 3 Figures 6-8), the ROC curve 207 and AUC vary for specific rats/runs (Supplemental Figure 9), and in some cases 208 demonstrate high predictive value. 209 For the Yang study 16 , we were able to obtain the exact geometry of the electrode. 210 Therefore, we were able to stratify DTW distances by the physical distance of the 211 neurons measured (Figure 3). For this stratification, we were unable to determine a 212 consistent trend. Rats B and C demonstrate decreased DTW distance (increased 213 synchronicity) in both neural spike trains and LFPs, whereas the remaining rats 214 demonstrated slightly increasing DTW distance (decreased synchronicity) in both, 215 except for rat S. For rat S, however, the large decrease in DTW distance seen in the 216 spike trains may not be captured in the LFPs, because there was no zero-distance 217 comparison for LFPs, whereas multiple spike trains recorded from the same electrode 218 could be thought of as having zero distance. 219 Discussion 220 There are two major issues with attempting to measure a specific population of neurons. 221 First, precise surgical implantation of electrodes is required; second, enough neurons 222 must be sampled to capture overall neural network dynamics. These two issues 223 synergize at our current level of intracranial electrode technology, complicating fine 224 measurements of neurons. While the rat brain is far smaller than a human brain (~1 cm), 225 it still contains an estimated 21 million neurons 18 . Reproducible surgical implantation in 226 a specific area of the rat brain such as the medial prefrontal cortex is therefore highly 227 dependent on the fine-motor skill of the researcher. Additionally, due to the lack of 228 spatial organization of the prefrontal cortex, we postulate that the surgical sampling of 229 neurons is essentially random. Therefore, the chances of measuring the same set of 230 neurons among different rats or after adjusting electrode positioning is extremely low, 231 an issue which can only be mitigated through the measurement of larger populations of 232 neurons. In the Methods, we discuss these issues further with a simple mathematical 233 illustration (Reproducibility of electrode recordings, Figure 1d). 234 In our analysis, we demonstrated that the DTW distance and the computed parameter 235 ݀ c r i t captured some of the mPFC neural network firing dynamics for both the spike trains 236 and LFPs. We found weak evidence that the correctness of rat choice influences the 237 firing dynamics (Supplemental Table 3). However, we expect that many more 238 experiments are needed to confirm the classification of the correctness of a decision 239 based on mPFC activity because our results indicated that firing events could be pooled 240 across the same run, but not across different runs or studies (Supplemental Tables 1  241 and 2). In terms of spatial dependence of firing synchronicity, we found that for specific 242 rats and runs, the classification of DTW distances for recordings taken from the same vs 243 different electrodes was both statistically significant (Supplemental Figures 4, 6-8) and 244 useful as a regressor (Supplemental Figure 9). For the Yang study 16 , the spatial 245 dependence of DTW distances was highly dependent on the specific rat sampled 246 (Figure 3), providing evidence of the lack of spatial organization of the prefrontal cortex. 247 The difficulties we encountered in analyzing spike trains and LFPs can be addressed 248 through changes to the experimental design. The number of neurons sampled is an 249 important characteristic to consider, especially if only a few neurons were recorded and 250 some recordings represented only a single neuron (e.g., Stout 15 ). If the mPFC is 251 composed of separate neural networks which perform specific, modular tasks, it is 252 unlikely that recording a low number of neurons will provide meaningful, interpretable 253 results. Since firing dynamics were consistent across a single run, we posit that multiple 254 recordings of a single neuron population are the best approach to characterizing 255 network firing behavior. The DTW spatial dependence remained consistent across 256 multiple runs for rats A and C (Figure 3), therefore we are reasonably confident that we 257 recorded the same population of neurons across the runs. In contrast, because the 258 electrodes were adjusted after each run, Ito 14 and Stout 15 only recorded a single run per 259 neuron population, which complicates the modelling we must perform. We hypothesize 260 that the additional variation across runs in Yang 16 may be due to changes in the resting 261 cognitive state of the rats. For the data analyzed in this work, it was not possible for us 262 to consistently distinguish between variance due to noisy measurements and variance 263 due to rat behavioral changes. It may be necessary to perform many control runs to 264 characterize in adequate detail all of the resting cognitive states of a specific rat. 265 Crucially, recording a single control run as Yang 16 did for the seven other rats (B, K, Q, 266 R, S, T, X) probably fails to detect within-rat variations in the resting state. To 267 complement the recording of many control runs, recording additional rat behaviors may 268 yield insight as to how the rat resting cognitive state influences the network firing 269 dynamics. 270 To provide a guide for researchers investigating the connection between 271 neuroelectrophysiology and behavior, we asked what did each of the three studies do 272 right and what could be improved? There are three parameters we considered: n neurons 273 (the number of neurons recorded); n network repeat (the number of times the same neural 274 network is recorded); and n rats (the number of biological replicates). We contend that all 275 three of these parameters must be relatively high to make meaningful conclusions, 276 especially if one wishes to use more powerful parametric statistical tests over the non-277 parametric tests we used in this study. 278 Stout 15 recorded too few neurons (n neurons ~1), and changed the neurons recorded 279 across runs (n network repeat = 1) for n rats = 5. Ito 14 recorded a moderate number of neurons 280 (n neurons ~20), but still changed the neurons recorded across runs (n network repeat = 1); n rats 281 = 3. Yang 16 recorded a moderate number of neurons (n neurons ~20), and kept the same 282 neurons recorded across runs (n network repeat = 5), but only performed these replicates for 283 n rats = 2; the other n rats = 7 had (n network repeat =~1). 284 Recent advances provide hope that we may one day understand the relationship 285 between neuron firing and consciousness. Whole brain connectomes have been 286 mapped for simple organisms (C. elegans 19 , Ciona intestinalis larvae 20 , and Platynereis 287 dumerilii 21 ); and recently, all 3,000 neurons and 548,000 synapses of a fruit fly larva 288 (Drosophila melanogaster) have been mapped 22 . Although mammalian brains are 289 orders of magnitude larger than insect brains, whole-brain connectomes for rats may 290 not be out of reach in the coming decades. These connectomes are essential for 291 understanding cognition, but they provide only a static reference. This study highlights 292 the data requirements for the next steps in investigating the function and dynamics of 293 biological neural networks. 294

Methods and Materials 295
Dynamic time warping and statistical tests 296 Mathematical details on dynamic time warping, the Mantel test, and the Kolmogorov-297 Smirnov test are included in the Supplemental Methods. 298

Estimated number of neurons in the rat medial prefrontal cortex 299
Our goal here is to provide a reasonable estimate for the number of neurons that one 300 must record to reconstruct the dynamics of the entire mPFC. We used existing 301 estimates of neuron and synaptic density in the rat prefrontal cortex as starting points 302 for our modeling. We set the total number of neurons in the rat brain 18 to N brain = 2.1 x 303 10 7 with volume 23 V brain = 2,500 mm 3 . We set the volume of the mPFC 23 to V mPFC = 20 304 mm 3 . Assuming a uniform distribution of neuron count throughout the brain, the 305 estimated number of neurons in the mPFC was

Reproducibility of electrode recordings 307
Here, we discuss the problem of reproducibility in neuroelectrophysiology. Because of 308 the apparent lack of spatial organization in the mPFC, when electrodes are implanted, 309 they uniformly sample the estimated N mPFC = 50,000 neurons in the mPFC. Depending 310 on how cooperative and synchronized neuron firing is, we assume that there is a certain 311 number of neural clusters for which the measurement of a single neuron 312 sufficiently captures the behavior of the entire cluster, and that measurement of these 313 clusters correlates with rat behavior. We also assumed that sampling is done with 314 replacement, i.e., n cluster << N mPFC. For reproducibility, one must (1) record the same 315 neural clusters across biological replicates and (2) identify or classify the neural clusters 316 to show that findings in one animal generalize to other animals. 317 (1) The first issue of capturing a certain percentage of the neural clusters in the rat brain 318 corresponds to a classical problem in combinatorial probability known as the coupon 319 collection problem. The mathematical question is: given c = n cluster categories, what is 320 the expected number of samples ܵ which needs to be drawn from those categories so 321 that all categories are represented at least once? In the case of equal probabilities for 322 drawing each category, 323 where H c is the c th harmonic number. Alternatively, we may require that a certain 325 number 1 < k < c out of all the categories is represented 24 , so that Eq. ( 1 ) is a special 326 case of 327 ( 2 ) 328 In both cases, S ~ c log c asymptotically. We have plotted the expected number of 329 clusters we may identify depending on the maximum number of neurons recorded from 330 each of the studies (Figure 1-D). In practice, it is necessary to first determine the 331 number of neural clusters that exist and to choose a number of neurons sampled that 332 can adequately sample all the clusters with high probability. 333 (2) The second issue of identifying neural clusters across different animals relies on the 334 researcher's ability to determine which of the n cluster ! possible matchings is the correct 335 matching between any two rats.    The study to which the rat belongs rat 17 total animals The rat identity run 62 total runs (with more than 1 neuron recorded) The rat identity and the date of recording