Model of multiple synfire chains explains cortical spatio-temporal spike patterns

It has been postulated that information processing in the brain is based on precise temporal correlation of neural activity across populations of neurons. In a recent study we found spatio-temporal spike patterns in experimental recordings from monkey motor cortex, and here we study if those could be explained by a synfire chain (SFC) like model. The model is composed of groups of neurons connected in feed-forward manner from one group to the next with high convergence and divergence. When activated, e.g., by a current pulse to the first group, spiking activity in the SFC is synchronous within neurons of the same group and propagates from group to group. When a few neurons from different groups are recorded from such an SFC, and the SFC is repeatedly activated, we would find a spatio-temporal spike pattern repeating across trials. Here, we take the statistics of the STPs found in the experimental data from 20 sessions as a reference to compare to a simulated network. Distributions of the data we take into account include 1) the pattern sizes, i.e. the number of neurons involved in the patterns, 2) the number of patterns a single neuron is involved in, 3) the durations of the patterns, and 4) the spatial distances of the patterns across the electrode array used to record the data. For the simulations, we embed SFC(s) in an anatomical model of the respective layer of the motor cortex, defined by its height and the density of the neurons. Model parameters are the length of the SFC, the number of neurons per group, the spatial extent of each neuronal group, and the distance between subsequent groups. Given the size and reach of the Utah array electrodes, we derive the probability of recording neurons from the SFC network. An SFC is considered detected if at least two neurons from two different groups are recorded. We find that depending on the model parameters, an embedded SFC can be detected with high probability, despite the massive subsampling of the cortex by the Utah array. Furthermore, to achieve multiple membership of a neuron in different patterns, we embed multiple SFCs that overlap. The fitting of the model to the pattern data constrains the spatial SFC parameters: the chains have to be broadly distributed in space and contain many neurons per group to match the experimental results.

It has been postulated that information processing in the brain is based on precise temporal correlation of neural activity across populations of neurons. In a recent study we found spatio-temporal spike patterns in experimental recordings from monkey motor cortex, and here we study if those could be explained by a synfire chain (SFC) like model. The model is composed of groups of neurons connected in feed-forward manner from one group to the next with high convergence and divergence. When activated, e.g., by a current pulse to the first group, spiking activity in the SFC is synchronous within neurons of the same group and propagates from group to group. When a few neurons from different groups are recorded from such an SFC, and the SFC is repeatedly activated, we would find a spatio-temporal spike pattern repeating across trials. Here, we take the statistics of the STPs found in the experimental data from 20 sessions as a reference to compare to a simulated network. Distributions of the data we take into account include 1) the pattern sizes, i.e. the number of neurons involved in the patterns, 2) the number of patterns a single neuron is involved in, 3) the durations of the patterns, and 4) the spatial distances of the patterns across the electrode array used to record the data. For the simulations, we embed SFC(s) in an anatomical model of the respective layer of the motor cortex, defined by its height and the density of the neurons. Model parameters are the length of the SFC, the number of neurons per group, the spatial extent of each neuronal group, and the distance between subsequent groups. Given the size and reach of the Utah array electrodes, we derive the probability of recording neurons from the SFC network. An SFC is considered detected if at least two neurons from two different groups are recorded. We find that depending on the model parameters, an embedded SFC can be detected with high probability, despite the massive subsampling of the cortex by the Utah array. Furthermore, to achieve multiple membership of a neuron in different patterns, we embed multiple SFCs that overlap. The fitting of the model to the pattern data constrains the spatial SFC parameters: the chains have to be broadly distributed in space and contain many neurons per group to match the experimental results. mechanisms that are different from those of the classical SFC model of [1]. Similar to 31 the original SFC model, they often rely on successively activated groups of 32 synchronously firing neurons [21][22][23][24][25][26]. In this work, we therefore refer to these models as 33 SFC-type models, even though the mechanisms triggering synchronous firing in groups 34 of neurons are different from the classical SFC model, the connectivity structure is not 35 always strictly feed-forward, and the group sizes can be substantially smaller (of the 36 order of 10) if nonlinear dendritic integration (dendritic action potentials) is accounted 37 for [21][22][23]26]. Other models predict the emergence of robust precise firing patterns 38 without synchronous firing, but with specific temporal delays, such as the "braids" 39 model [14], or the polychronization networks [27]. 40 Temporal and spatio-temporal patterns are observed across several species and brain 41 areas, and at different temporal resolutions ranging from milliseconds to 42 seconds [16,[28][29][30][31][32][33][34][35][36][37][38][39]. In a recent study we found spatio-temporal spike patterns (STPs) in 43 experimental data recorded from primary motor cortex (M1) of two monkeys [40]. Here, 44 we investigate if these patterns can be explained by SFCs, given the recording setup 45 with a 10x10 Utah array. Our assumption is that an STP results from recordings of at 46 least two different neurons that take part in the SFC (as sketched in Fig 1B). We 47 conceptually embed SFC(s) into a cortical volume corresponding to the region recorded 48 from by the Utah array: an area of 4 × 4 mm 2 , layer 2/3 and its respective neuron 49 density. Additionally, the sensitivity of the individual electrodes of the Utah array are 50 considered for estimating the neurons recorded from in the model. 51 The distributions of the number of neurons involved in an STP, the temporal extents 52 of the STPs, the spatial cortical distance of neurons involved in STPs, and the 53 involvement of a single neuron in different STPs are the references for comparison. As 54 shown in the modeling results (Section 3.1.2), a single SFC is likely to be detected for a 55 large range of model parameters. However, to account for the STP characteristics found 56 in the experimental data, the model has to include multiple SFCs within the assumed 57 piece of cortex. Given a conservative estimate of how many chains exist in the recorded 58 volume, the detection of SFCs in form of STPs is very likely despite the massive 59 subsampling by the Utah array (Section 3.1.3). Additionally, the fact that individual 60 neurons occur in different patterns requires that these SFCs overlap and include the 61 same neurons (Section 2.3). The final results of the simulations (Section 3.2.4) 62 correspond well to the experimental data, which leads to the conclusion that STPs are 63 well explained by multiple SFCs. 64 2 Methods 65 2.1 Experiment and data 66 We are interested in understanding whether STPs detected in experimental data may be 67 explained by the SFC model. Thus, we examine spatio-temporal spike patterns from density in layer 2/3 is 18 % larger than the average density across all layers [42,43,45]. 80 In a separate study [40], we analyzed 20 experimental sessions (10 per monkey), 81 recorded in different days over the time span of months. For monkey N, the recordings 82 were performed in 2014, from June to July; whereas, for monkey L, recordings started in 83 October 2010 and finished in February 2011. Each session has a duration of 84 approximately 15 minutes and consists of around 120 successful trials. The data from 85 each session were spike sorted independently using the Plexon spike sorter (Plexon Inc, 86 Dallas, Texas, USA, version 3.3), retaining only the single unit activities (SUAs) with a 87 signal-to-noise ratio ≥ 2.5, and with an average firing rate of < 70 Hz across the trials. 88 Artifacts consisting of hyper-synchronous spikes at sampling resolution occurring across 89 electrodes were assumed to be cross-talk artifacts and were therefore detected and 90 removed in a preprocessing step [31]. Over all sessions, the average number of SUAs per 91 session is 107, however, it differs across monkeys: higher for monkey N (mean 92 =143 ± 14.82 units), and lower for monkey L (mean =70.5 ± 13.93 units). The average 93 number of SUAs per electrode is 1.1 averaged over the two monkeys. 94 The two monkeys, L and N, were trained to self-initiate trials by pressing a start 95 button. After a waiting period of 400 ms, a visual cue (yellow LED) was shown to the 96 monkey (waiting signal). The monkey was instructed to wait again for 400 ms, until it 97 was presented to another visual cue, lit for 300 ms (from CUE-ON to CUE-OFF). The 98 task consisted in reaching and grasping an object with the indicated grip type and force 99 level. The grip can be either be a precision grip or a side grip. The CUE-ON signal 100 contained the grip information. The monkey then waited for 1, 000 ms, eventually 101 receiving the GO-SIGNAL, which contained the information on the amount of force to 102 exert to pull the object towards itself. The behavioral conditions were selected 103 randomly for each trial. The requested grip and force had to be maintained for 500 ms, 104 and if it was performed correctly, the monkey received a reward. Further details on the 105 experimental setup can be found in [41,47] and analysis results on the same data are 106 presented in [31,48]. 107 2.2 Analysis approach for spatio-temporal spike patterns in 108 the experimental data 109 We used the SPADE [49][50][51] method to detect STPs in the parallel spike trains of these 110 recordings. The method consists of three successive steps. The first is the detection of 111 all putative patterns in the data at a certain temporal resolution, during which the 112 putative patterns are stored along with the information on when and how often they 113 occur. This is achieved by applying a Frequent Itemset Mining algorithm [52,53]. The 114 second step is the statistical evaluation of the significance of the patterns detected in 115 the first step, under the null-hypothesis of mutual independence of spike trains given 116 their firing rate (co-)modulations [54]. The third step is a conditional test performed on 117 all significant patterns, in order to remove patterns arising from the overlap of true 118 pattern spikes and chance spikes. SPADE outputs all patterns labeled as significant, i.e. 119 STPs, by the statistical tests. Each pattern is returned together with the neuron ids 120 involved in the patterns, the lags between spikes, the times of the occurrences, and the 121 p-value. Details of the method and on its implementation can be found in [51]. The 122 temporal resolution of the detected patterns is given as an input to the SPADE method, 123 and is here fixed to b = 5 ms. Moreover, we fix the maximal duration of STPs to 60 ms. 124 Each recording session contains four trial types (corresponding to the combinations 125 of the two force types and two grip types). For the pattern analysis, we segmented the 126 trials of about 4 s into six 500 ms-long behavioral epochs in each single trial; start, cue, 127 early delay, late delay, movement and reward. The respective epochs of the same trials 128 of the same session are concatenated and analyzed as an independent data set, and thus 129 in total yield (4 trial types × 6 epochs) data sets per session (see [31] for details on the 130 data segmentation into trial epochs). The data sets of each session are analyzed 131 separately, since the probability to record the same neuron(s) across sessions is 132 unlikely [55]. 133

134
We analyzed the data recorded in the delayed reaching-and-grasping task for STPs 135 using SPADE. We detected 119 patterns in total (monkey N: 61, monkey L: 58) in 20 136 sessions. The statistics of the results pooled over all sessions are displayed in Fig 2. 137 STPs in different epochs are not identical, but specific to the behavior. However, for 138 this study we pool across trial types, epochs, and sessions. Typically we find around 6 139 patterns per session (mean = 5.95 ± 3.26). Patterns are formed by spikes emitted by 140 different neurons, and contain spikes of 2-6 neurons (mean = 2.9 ± 0.93 neurons), but 141 the majority of patterns consist of two and three spikes (Fig 2A). Another statistic is 142 the number of different STPs a single neuron (within a single session) is involved in 143 ( Fig 2B). Most STP neurons participate in only 1 pattern, but many also in 2-3 144 patterns, and in a few cases, a single neuron can partake in up to 9 STPs [40,56]. Histogram of the number of STPs a single neuron appears in (multiple membership). (C) Histogram of time lags between neighboring pattern spikes, i.e. only the time difference between subsequent pattern spikes are taken into account. The temporal resolution of the histogram coincides with the resolution of the SPADE analysis (5ms) . The maximum time lag corresponds to the maximal pattern duration (here 60ms). (D) Histogram of the euclidean distances between pattern spikes. Only subsequent pattern spikes are taken into account. The entries for each spatial distance are weighted by the occurrences of the respective distance on the Utah array electrode grid. All histograms are normalized such that the sum of all entries is 1.
Regarding the temporal delays between the spikes in STPs, we observe that the 146 whole range of allowed delays (0 to 60 ms) is covered for both monkeys (Fig 2C). We do 147 not see a particular tendency towards preferred delays or oscillatory activity. However, 148 the results show a slight preference for temporal delays around 10 to 30 ms. We also 149 calculate the euclidean spatial distances between STP neurons of the ones forming 150 successive spikes. The position of the respective neuron is given by the electrode the 151 neuron is recorded from. The distribution of these distances is shown in Fig 2D). This 152 distribution is normalized by the number of distances that exist on the array. The 153 resulting distribution is rather flat, with a peak around 3, 600 µm, however, we tend to 154 ignore that peak since the larger the distances, the fewer number of samples are 155 available which makes the statistics less reliable. 156

Model for spatial embedding of SFCs in a cortical network 157
We aim to assess the detectability of SFCs in the experimental setting presented in 158 Section 2.1. Since the recording was performed with a 10 × 10 electrode Utah array 159 recording with 96 active electrodes arranged at 400 µm distance, we only look at the 160 cortical volume below the area of the electrode array (4 mm×4 mm) for this study. 161 Moreover, we only consider layer 2/3 by setting the height of the simulation volume to 162 its thickness: h SFC Volume = 1.5 mm (see [42][43][44][45]). The volume in which neurons can be 163 detected by the electrodes of the array is confined to spheres with a radius 164 corresponding to the sensitive range (see Section 2.5). The simulation volume and its 165 height h SFC Volume = 1.5 mm corresponds to the volume in which SFCs are embedded in 166 our model. Since the sensitive range of the electrodes is a separate parameter, 167 increasing the volume height effectively increases the subsampling of the neurons in the 168 simulation volume by the Utah array. We choose a default value of ρ = 35, 000 mm −3 169 for the neuron density in our simulation, but we will investigate the impact of varying it 170 in the following sections. In our simulation, we first fill the volume under the Utah 171 array with neurons according to the neuron density ( Fig 3A). Then, we simulate a single 172 SFC realization by positioning a sequence of SFC groups in the space of the cortical 173 area. In particular, the SFC embedding procedure starts by drawing the position of the 174 first group center within the borders of the electrode array from a two-dimensional 175 uniform distribution covering the area of the array. Around the first group center, we 176 select w neurons inside a cylinder with radius r group and height h SFC Volume (Fig 3B).

177
The parameter w controls the number of neurons per group, and r group represents the 178 radius of the cylinder in which all w neurons of a single group reside. For the 179 subsequent group, the center position is drawn from a two-dimensional, symmetrical 180 Gaussian distribution which is centered on the previous group center (Fig 3C). The 181 standard deviation σ group distance of this Gaussian distribution determines the average 182 spatial distance between groups ( d ): The next group of neurons is selected inside the cylinder around the center position 184 analogously to the first group. This procedure is repeated to distribute l groups 185 iteratively, where l represents the length of the SFC, i.e. the total number of groups in 186 the chain (Fig 3D). We allow groups of a chain to be partially or fully outside of the 187 simulation volume (e.g. the final group of the SFC shown in Fig 3D), since there is no 188 plausible reason for SFCs to be constrained to the volume below an electrode array.

189
However, these neurons are not available to be 'recorded' from. The simulated volume is homogeneously filled with neurons according to the local neuron density (here the density is reduced to 10 mm −3 to unclutter the plot). Its size is given by the area covered by the Utah array along the x and y axes, and by the thickness of the cortical layer 2/3 (1, 500 µm) in which the electrode tips are placed in the z axis. To illustrate the relation of the shown spatial scales to the one of an electrode and its sensitive range one electrode (at the corner of the Utah array) is shown in gray on the right side of the plot and its sensitive range is illustrated as a green sphere at its tip. (B) The center of the first synfire group is randomly placed inside the volume. A group is confined by a circle with radius r group in the horizontal plane, and by the thickness of the cortical layer along the third dimension, resulting in a cylinder (not to be confused with a micro-column). w neurons inside the cylinder are randomly assigned to the group. For illustration purposes we show only groups of w = 10 neurons here.
(C) Top down view on the cube showing everything projected onto an x-y plane. The Utah electrodes and their sensitive ranges are illustrated across the whole area by black dots surrounded by white circles, respectively. The center of the second SFC is drawn from a two-dimensional Gaussian centered on the previous group center, with a standard deviation σ group distance . The neurons of the second group are assigned analogously to the first group. (D) This procedure is repeated for l groups, resulting in one SFC being embedded in the simulated volume.
To reduce computational cost and account for participation of one neuron in 191 multiple SFCs we always embed multiple SFCs per simulation. In order to achieve that, 192 we repeat all of the above steps n chains times to embed n chains SFCs in the volume. A 193 neuron can be part of multiple groups and multiple chains. This setup allows us to 194 study multiple trials if we consider single chains one after another. Since an electrode only records neurons whose soma is close to the electrode tip, we 197 need to consider their distance to the electrode [57]. We make the simplifying 198 assumption that for isolating the activity of a single neuron, its cell body needs to lie 199 within a certain radius around the electrode tip, which we refer to as the electrode's 200 sensitive range (r sens ): Here, x, y and z are the spatial coordinates of the neuron and x el , y el and z el are the 202 coordinates of the electrode tip (the tip is where the actual recording takes place). We 203 set the depth of the electrode tips to r sens . This is arbitrary, but z el does not have any 204 impact on the results since the z coordinates do not affect the assignment of neurons to 205 SFC groups (see Section 2.4). Thus, a neuron can be recorded if it lies within a sphere 206 with a radius equal to the sensitive radius of the electrodes (r sens ). The experimental 207 data show that on average 1.1 single units are recorded from a single Utah array 208 electrode, estimated from all 20 sessions (see Section 2.1). Thus, we set the isolation 209 probability of a neuron located in the sensitive sphere around an electrode tip with 210 radius r sens accordingly: with the neuron density ρ in the vicinity of the electrode and the number of single 212 units per electrode in the numerator. Note that this means that independently of the 213 SFC parameters, the constant average number of detected neurons per electrode is 214 always 1.1. This matches the number of isolated single units, which might seem low at 215 first glance. However, to be selected, they also have to be isolated by the spike sorting 216 software, which requires that neurons have to fire often enough, have a high enough 217 signal to noise ratio and have a spike waveform that is clearly distinguishable. This 218 number is thus to be expected to be quite low, up to 10-to 100-fold lower than expected 219 from the neuron density [58]. Assuming a sensitive radius of r sens = 50 µm [57] and a neuron density of ρ = 35, 000 mm −3 [46], there are 18 neurons in the sensitive volume of 221 each electrode of which we can isolate 1.1 per electrode, matching the 10-to 100-fold 222 subsampling observed by Shoham et al. (2006) [58]. See [47] for more detailed 223 information on the applied spike sorting procedure.

224
Sensitive Range of Utah Electrodes In order to estimate the sensitive range of 225 the Utah array electrodes, we have to extrapolate from measurements in a different 226 setting due to the unavailability of directly comparable data. [57] estimated the sensitive 227 range of extracellular electrodes to be 50 µm by comparing simultaneous extracellular 228 and intracellular recordings. The electrodes they used have a significantly higher 229 impedance than Utah array electrodes (400 kΩ to 600 kΩ vs 50 kΩ, cf. [59]), which 230 results in a smaller sensitive range and a higher signal-to-noise ratio, but they were also 231 recording in rat hippocampus which has a ten-fold higher neuron density than macaque 232 M1 [60,61] which we are considering in our experimental setting (cf.Section 2.1). Since 233 these differences have both positive and negative contributions to the sensitive range, 234 we assume a range of 50 µm for most analyses but investigate the effect of smaller (30 235 µm) and larger ranges (70 µm) on the detectability of SFCs to quantify the associated 236 uncertainty. For illustration of the size relationships we refer to Fig 3A, which shows a 237 sketch of an electrode of the Utah array including its sensitive range (green sphere).  Next, we would like to quantify the ability to detect an SFC in such a setting. For 246 that, we assume the dynamics in such an embedded SFC as determined by [3] for the 247 case of stable propagation. Diesmann et al. [3] assumed an all-to-all connectivity from 248 one group to the next and found that at least 80 neurons are required per group to get 249 a stable propagation. If the chain is activated by a strong current pulse in parallel to all 250 neurons of the first group, this first group then exhibits synchronous firing such that all 251 neurons of the next group are also activated synchronously at some delay 252 (corresponding to the time delay from one group to the next). That way, a packet of 253 synchronous activity propagates through the chain. For the detection of an active SFC, 254 we require that at least two neurons from two different groups are recorded, since only 255 spike patterns with delays show propagation of activity along the chain. Thus, we define 256 the SFC detectability as the probability to detect neurons from two or more different 257 groups, cf. Fig 4. The detected neurons correspond to an STP where the time difference 258 between spikes of the different neurons are assumed to be given by the number of SFC 259 groups between the two neurons multiplied by the time to propagate from one group to 260 the next, here assumed as 5 ms. With these assumptions, we get one STP per SFC.

262
In this section, we first present the results of the STP analysis of experimental data, 263 and show the corresponding statistics. Since we assume that the subsampling of an SFC 264 may lead to such spike patterns, we aim to verify that SFCs can be detected given the 265 experimental Utah array recording setup described in Section 2.1 and our spike pattern 266 analysis (Section 2.3). Thus, we come back to our model for spatial embedding of SFCs 267 (see Section 2.4), and verify the influence of model parameters and recording constraints 268 on the probability of observing a single SFC (SFC detectability). Within the parameter 269 regions with good SFC detectability, we compare the predicted spike pattern statistics 270 obtained by multiple simulations of our model to the statistics obtained experimentally. 271 This yields the best-fit parameter set in which our model can explain the spike patterns 272 in experimental data assuming that the underlying network utilizes SFCs for 273 information propagation. 274

SFC detectability 275
Our next question is if one embedded SFC is detectable in 10x10 Utah array recordings. 276 In that context we first aim to answer how the different SFC parameters and recording 277 constraints influence the SFC detectability. We start by fixing all but one parameter and investigate the effect of the varied 280 parameter on the SFC detectability. Therefore, we have to decide on reasonable default 281 values for all parameters. We assume that the temporal delay for the assumed 282 propagation of the SFC activity from one group to the next is b = 5 ms. In the 283 experimental data we found STPs with a maximum temporal delay of 60 ms, which was 284 the maximal allowed duration of the patterns in the analysis. This corresponds to a 285 maximal SFC length of l = 12, since l · b = 12 · 5 ms = 60 ms (see Section 2.3). For the 286 SFC width we choose for a start an intermediate value of w = 1, 300 which is larger 287 than the minimal value w = 100 for stable, robust propagation of synchronous spiking 288 along the SFC [3] but far away from the maximal value of the number of neurons that 289 fit into the smallest assumed group disk r group = 300 µm we consider. The maximal 290 value would be w = π · r 2 group · h SFC Volume · ρ = 14, 844 neurons per group. Getting close 291 to this maximum would also mean that most neurons within r group would be inside the 292 same group of the same SFC, since no other neurons would be available anymore in the 293 volume and thus would not allow for any additional embedding flexibility.

294
The spatial parameters r group and σ group distance are limited by two aspects: very 295 small values would result in SFC groups that fit in between the electrodes' sensitive 296 ranges and in an SFC that does not move spatially, respectively. On the other hand, 297 very large values would result in parts of a group or parts of the chain residing outside 298 the volume below the area of the Utah array. Therefore, we fix both parameters to an 299 intermediate value of r group = σ group distance = 900 µm. For the sensitive range of the 300 electrodes we choose r sens = 50 µm, but we expect this not to have an effect as discussed 301 in the appendix (S1 Appendix: The influence of the sensitive range of the electrodes). 302 The  The SFC width w determines the number of neurons in the SFC groups and 315 increasing w also leads to an increase of the detectability, however with a different slope. 316 That curve does not flatten out, since for example, a higher w for a given l just adds 317 more neurons to the existing groups and thus further increases the detectability 318 ( Fig 5B).

319
The spatial SFC parameters r group and σ group distance have quite different effects on 320 the detectability. The spatial SFC group size r group hardly affects the detectability 321 ( Fig 5C). The increased area covered by an SFC group is counteracted by the decreasing 322 density of SFC neurons inside that area, such that the detectability only increases very 323 slowly for larger radii. That slow increase is however consistently monotonous. In 324 contrast to that, increasing the spatial SFC inter-group distance σ group distance decreases 325 the detectability (Fig 5D), since a higher inter-group distance results in a higher chance 326 for the SFC to leave the simulation volume. The supplementary figure (S2 Figure: 327 Impact of SFCs leaving the simulation volume) supports this argument; if the boundary 328 conditions for the recorded area/volume were periodically continued there would be no 329 decay of SFC detectability. 330 The sensitive range of the electrodes r sens has no effect on the SFC detectability in 331 our model (Fig 5E). This seems rather counter-intuitive, but the sensitive range is 332 implicitly included in the number of SUAs that can be recorded from an electrode 333 (Eq 3). This parameter has to be fixed in our model since we have to select the number 334 of neurons that can be recorded by an electrode. In order to do this we have to set a 335 maximal distance r sens from the electrode tip. The exact value we choose does not 336 matter since it eventually cancels out as discussed in the appendix (S1 Appendix: The 337 influence of the sensitive range of the electrodes). 338 We vary the neuron density between ρ = 20, 000 − 50, 000 mm −3 to account for the 339 measurement uncertainty and the variance between measurements in different 340 studies [42,43,45,46]. The neuron density has a straight-forward effect on the SFC 341 detectability. Increasing the neuron density results in more neurons in the simulation 342 volume but decreases the ratio of neurons within one SFC to total neurons. This results 343 in a decreasing detectability (Fig 5F). In the previous section, we derived how single parameters affect the SFC detectability. 346 For this, we fixed all but one parameter and investigated the effect of the varied 347 parameter. In order to check for possible co-dependencies, we will now vary all 348 parameters at once and calculate the SFC detectability for parameter combinations 349 across our parameter space.  Increasing the group width w (along the y-axis within each quadrant), the SFC 361 detectability increases. If in addition the SFC length l is increased (across the 362 quadrants along the y-axis) the SFC detectability increases even more. Intuitively, the 363 more neurons belong to an SFC, the more likely it is that two neurons are detected from 364 different groups. Increasing the SFC group displacement σ group distance (along the x-axis 365 within the quadrants) lowers the SFC detectability for a given fixed group disk radius 366 r group , whereas a change of the group disk radius r group (across the quadrants along the 367 x-axis) does not lead to big differences.

368
In summary, the SFC detectability is highest for high l and w and low group 369 distances σ group distance . The SFC detectability reaches a maximum of 0.975 for l = 48, 370 w = 2, 100, σ group distance = 100 µm and r group = [1, 300, 1, 700, 2, 100] µm. This means 371 that for such values the detectability of an SFC is very likely, and we should observe a 372 spatio-temporal pattern of at least two neurons involved. These parameter 373 configurations with large groups are actually feasible and lead to good agreement with 374 experimental data as will be seen in the following sections. In the previous section we estimated the detectability of a single SFC which lies in the 377 cortical volume where the Utah array is recording from. However, assuming that there 378 is more than one SFC in the recorded cortical volume, the probability of detecting an 379 SFC should be larger. As explained in Section 2.6 each SFC results in one STP if at 380 least two participating neurons are recorded and cannot result in multiple STPs since 381 they would simply be merged into a larger STP in the spike pattern detection process. 382 Thus also the probability to detect spatio-temporal spike patterns should be larger if The volume in which we can measure SFCs is given by the area of the Utah array in 393 the horizontal plane (A array ) times the SFC volume height (h SFC volume ) along the Thus for our parameter ranges that would result in 7 to 560 chains. The probability of 403 detecting at least one SFC when multiple SFCs are present is binomially distributed: where n SFCs is the number of SFCs in the cortical volume and p is the probability to 405 detect a single SFC (the SFC detectability). The required number of SFCs in the 406 cortical volume to detect at least one SFC with a probability equal to or greater than α 407 is obtained by plugging P (k SFC ≥ 1) = α into Eq 7 and solving for n SFCs : For an example parameter set with l = 12 groups, w = 100, r group = 900 µm,  In order to compare the spike pattern statistics predicted by our model to the statistics 443 observed in the experiment, we need a way to quantify their similarity. We will compare 444 the distributions of the four measures presented in Fig 11 for   Euclidean distance between two such histograms is always constrained to 0 ≤ d hist ≤ 2. 461 This ensures that we can get a combined measure to which all histograms contribute 462 similarly (independently of their number of bins) by simply summing all distances. 463 We will use the Euclidean distances d pattern size , d multiple membership , d group lags and 464 d spatial distance and the total distance d total in the following.  Fig 8A). For very short chains not enough large patterns are found to 487 match the experimental data. Making the chains longer alleviates this issue. This effect 488 also saturates around l = 12 since due to the 60 ms duration limit we cannot indefinitely 489 make an STP larger by making the corresponding SFC longer. The other histograms do 490 not show strong effects when varying l since they only depend on the total number of 491 chains and the spatial spread of SFCs as we will see in the following paragraphs.

492
We have already seen that increasing the SFC width results in an increase of and d spatial distance with increasing r group . This is probably due to the fact that a small 502 group radius coincides to an increase of synchronous patterns of neurons recorded at the 503 same electrode, since it becomes more likely to detect multiple neurons of a group at the 504 same electrode. However, synchronous patterns and short spatial distances between 505 pattern spikes are not prominently observed in the experimental data ( Fig 2C). These 506 effects are much weaker than the ones observed when varying the other parameters, but 507 since they have the same sign they will add up when we look at the total distance.

508
The s.d. σ group distance of the two-dimensional Gaussian distribution for the distance 509 between successive SFC groups has a twofold effect on the histogram distances ( Fig 8D). 510 Increasing σ group distance increases the chance that later groups of an SFC leave the 511 simulation volume. This reduces the chance of finding large patterns and the chance of 512 finding long lags between pattern spikes, which is reflected in an increase of d pattern size 513 and d group lags , respectively. However, d spatial distance decreases for larger σ group distance . 514 For a small s.d. there is a bias towards short spatial distances between pattern spikes, 515 which is not what we observe in the experimental data. Since this parameter has 516 opposite effects on multiple distances, we expect to find an optimal value which depends 517 on the parameter set, and we expect the total distance to increase in both directions 518 from that optimal value.

519
Additionally, we inspect how varying the number of SFCs embedded in the 520 simulation volume n chains influences how often single neurons participate in multiple 521 SFCs (Fig 8E). We observe that d multiple membership decreases from low to higher values 522 of n chains , until it reaches a minimum at n chains = 200. However, an even higher n chains 523 increases the chance for neurons to be involved in multiple patterns, since the number of 524 chains is higher but the number of neurons stays the same. Neurons then have to be in 525 more chains on average, which results in them being in more patterns, thus shifting the 526 distribution of number of patterns a neuron appears in to the right (not shown). This 527 measure thus constrains the total number of chains in our model, and we can use it to 528 calculate the optimal number of embedded SFCs for each parameter set, which we will 529 do in the following section. This will essentially allow us to eliminate one SFC 530 parameter. We can do this since all other measures do not depend on n chains . For them, 531 n chains essentially determines the number of samples, but does not change the 532 distribution. We see an increase towards the extreme case n chains = 1 and some noise 533 for very few chains, but then all measures except for d multiple membership stabilize. 534 Finally, in (Fig 8F) we investigate the impact of the neuron density ρ on the 535 histogram distances. d pattern size increases slightly with increasing ρ, which matches our 536 findings in Section 3.1.1, where we showed that increasing the total neuron density ρ 537 decreases the SFC detectability. Since the number of chains and the number of neurons 538 in a chain remains constant, but the total number of neurons increases, the detection 539 probability decreases as the ratio of neurons in one chain to the total number of neurons 540 goes down, thus reducing the chance to record two or more neurons from a chain. Here, 541 this is reflected in an increase in d pattern size since not enough large patterns are found 542 to match the experimental data due to the decreased probability of recording multiple 543 neurons from a chain. Besides this small effect on d pattern size , the neuron density does 544 not show a significant impact on the histogram distances, since it does not affect the 545 spatial distribution of SFC neurons in any way.  Fig 8E. However, this optimal 551 value also depends on l and w, since for a neuron to appear in the same number of 552 SFCs, more chains are required if each chain contains fewer neurons (i.e. lower l or w). 553 In order to verify that we can use d multiple membership to find the optimal number of 554 simulated SFCs n chains for each parameter configuration, we plot d multiple membership as a 555 function of n chains for all combinations of l and w in Fig 9. Varying the combination of l 556 and w corresponds to varying the number of total neurons in an SFC which is simply the 557 product l · w. We can see in Fig 9 that for every parameter combination, we find a clear 558 minimum of d multiple membership which we can use to fix n chains to the corresponding 559 number of simulated SFCs. The optimal values range from 40 for very long and wide 560 chains containing many neurons each to 10, 000 for very narrow and short chains with 561 few neurons per chain. In the following, we will fix n chains for each parameter set using 562 this method and only use the optimal value for all further plots and analyses. The product l · w, which corresponds to the number of neurons per SFC, is denoted by the line color. Along the x-axis, we vary the number of simulated SFCs n chains , and on the y-axis, the histogram distance for the multiple occurrences of neurons in SFCs is shown. Only data points for cases in which patterns have been detected are shown. No patterns were detected in the cases of chains with few neurons (low l · w) for very low values of n chains ∼ 10, which is very far from the optimal values of n chains ∼ 10 4 . (B) Optimal number of chains corresponding to the minimal d multiple membership as a function of l · w, corresponding to the minima of the graphs in A.
We can estimate in how many chains a neuron participates in by dividing the 564 number of "neuron slots" in all chains l · w · n chains by the total number of neurons  Having fixed the optimal number of SFCs , and thus having eliminated one parameter, 572 we can move on to a parameter scan to see how the total histogram distance behaves as 573 a function of the remaining parameters. We showed in Section 3.2.2 that the neuron 574 density does not have a significant impact on the histogram distances. Thus, we will 575 focus on the SFC parameters l, w, r group , and σ group distance and show the results in   Fig 4). The black dot represents the parameter configuration having the lowest total histogram distance.
Across the board, the histogram distance is quite consistently decreasing for 578 increasing values of r group (lighter shades of blue when moving box-by-box to the right). 579 This is expected from Section 3.2.2, where we found that r group has to be high to fit the 580 distribution of time lags between pattern spikes.

581
In most sub-panels of Fig 10, for fixed l and r group , we see a diagonal band of lighter 582 blue and even white shading, corresponding to smaller total histogram distances. This 583 means that higher values of σ group distance require higher values of w for a good fit. This 584 effect has contributions from several single histogram distances, but can be seen most 585 prominently for the pattern size histogram (see Fig 8). High values of w are required 586 such that enough neurons per chain are detected to get large spike patterns as observed 587 in the experimental data (see also Fig 8). However, there is a limit to this since at some 588 point the patterns become larger than the experimental ones. This can be 589 counter-balanced by increasing σ group distance , which decreases the pattern size since 590 parts of the SFC are more likely to be outside of the simulation volume and are thus not 591 detectable and cannot contribute to patterns. A higher SFC length l also increases the 592 chance of larger patterns, which is why the diagonal of lighter shades moves down to 593 lower values of w in the sub-panels for higher values of l. 594 are ruled out since it is very unlikely that enough patterns can be detected in those 596 cases. This corresponds to the parameter sets with a detection probability lower than 597 0.8 in Fig 7. These parameter configurations are indicated by red crosses in Fig 10. The 598 corresponding total histogram distances are also rather high for these parameter sets as 599 can be seen by the darker shades of blue in the respective sub-panels in the figure. 600 We have selected the parameter set with the lowest total histogram distance 601 d total = 0.40 (black dot) and we will look at its pattern statistics in the next section.

602
This best fit of the experimental pattern statistics is reached at l = 24, w = 1, 300, 603 r group = 2, 100 and σ group distance = 500. As discussed in Section 3.2.2, a value of l > 12 604 does make sense here since it helps to increase the number of patterns with long 605 temporal delays by increasing the number of possibilities to observe two neurons that 606 are 12 SFC groups apart. In this parameter configuration with rather high values of w 607 and r group , SFCs contain many neurons per group, and a single group covers most of 608 the area of the Utah array. There is however a quite large region in the parameter space 609 with an almost similarly good fit, so this is not the only possible setting in which the 610 model can explain the experimentally found spike pattern statistics.  (Fig 11). Overall, the distribution of the best-fit parameter setting (in light blue) 616 matches well with the experimental one (in orange) for most of the distributions 617 (Fig 11A-C). The distribution of the Euclidean distance between pattern spikes 618 ( Fig 11D) is very variable due to the finite size of the electrode array. However, both the 619 experimental and the model distribution show entries across most distances (see also 620 Section 3.2).

622
In this paper, we asked 1) whether synfire chains can be detected by Utah array 623 recordings, and 2) whether or not the STPs and their statistics which we found in 624 experimental data can be explained by synfire chains embedded below the Utah array. 625 We were able to answer both questions positively, and found a large parameter range for 626 the embedding of multiple SFCs that elicit the same STP statistics as the experimental 627 data.

628
As data reference we used results on spatio-temporal spike patterns from 96-electrode 629 recordings (Utah array) from monkey motor cortex (two monkeys) reported on in detail 630 in a related study [40,56]. The data contained between 80 and 152 simultaneously 631 recorded neurons during motor behavior. Across 20 sessions (10 per monkey) we found 632 STPs and summarized their statistics in four histograms with averaged data from all 633 sessions of both monkeys (Fig 2): A) the number of neurons involved in STPs, B) the 634 number of STPs a single neuron is involved in, C) the maximal temporal duration of 635 STPs, and D) the spatial distances on the Utah array of the spikes involved in STPs. 636 We designed a model of the spatial embedding of one or more SFCs in the cortical 637 area recorded from by the Utah array (motor cortex). The neuron composition of each 638 SFC is parameterized by the length of the chain and the width of the groups. The introduce additional uncertainties. However, we were able to show in the appendix (S1 647 Appendix: The influence of the sensitive range of the electrodes) that the actual 648 sensitive range of the electrodes does not have to be measured and is implicitly included 649 in the number of isolated single units per electrode, a parameter available from 650 experimental data, here 1.1. This is in agreement with [71][72][73] who report between 0.2 651 and 1.2 single units per electrode across multiple subjects and implantation sites. We 652 define a model parameter r sens to select neurons near the electrode tips which are 653 candidates to be detected, but we were able to show that this parameter does not have 654 any influence on the SFC detectability (Section 3.1.1). This is especially advantageous 655 since the sensitive range has not been measured for Utah electrodes and thus underlies a 656 large uncertainty. The neuron density ρ in the cortical layer in which the simulation 657 volume is located constitutes another free parameter. However, we used the most 658 conservative choice, i.e. the highest density (more neurons in total but same number of 659 neurons within one SFC decreases the detectability, cf. Section 3.1.1), for our 660 detectability calculation (cf. Section 3.1.3). Additionally, we showed that for the spike 661 pattern statistics the density is not as critical since it does not affect the spatial 662 distribution of SFC neurons (cf. Section 3.2.2). 663 We constrained the embedded SFCs to the cortical layer 2/3, which the experimental 664 data was recorded from. Using the height of this layer as the height for the SFC volume, 665 h SFC Volume = 1.5 mm, (see [42][43][44][45]) and the neuronal density of ρ = 35, 000 mm −3 , we 666 expect 840, 000 neurons inside this volume. In the experimental data we analyzed, 1.1 667 single units were detected per electrode, which we used as the basis for our neuron whether they can be detected by its electrodes. To this end, we evaluated the 674 probability of recording single neurons, and we developed a detectability measure which 675 requires at least two neurons from different groups of the same SFC to be recorded for 676 the chain to be considered detected. Here, we assumed that every neuron of a chain 677 always fires when the chain is activated, such that one SFC can only result in one STP. 678 Going one step further, we predicted the statistics of the patterns that would be 679 detected in such a setting and fitted the model embedding parameters to the pattern 680 results obtained from the experimental data. We identified a region in the parameter 681 space in which the spike pattern statistics that our model predicts match the 682 experimental results well. The best match was found for a length of l = 24, a group size 683 of w = 1, 300, a spatial radius of each group r group = 2, 100 µm and a group distance of 684 σ group distance = 500 µm, which corresponds to wide SFCs with many neurons per group 685 and a broad spatial distribution. More combinations of l, w, r group and σ group distance 686 provide comparable results, see Section 3.2.4.

687
Since more than one pattern was observed per session of experimental data and since 688 some neurons took part in multiple patterns, we had to embed multiple SFCs at the 689 same time and we had to allow neurons to be in more than one chain. We found an 690 optimal match to the experimental data for broad spatial distributions since STPs in 691 the experiment were detected across all possible spatial distances on the Utah array.

692
As shown in [74,75], SFCs with large group sizes w and dense inter-group 693 connectivity easily undergo an instability where even small fluctuations in the 694 background activity trigger the formation of synchronous bursts of spiking activity.

695
This problem of spontaneous synchronization in broad SFCs can be alleviated by a 696 number of mechanisms, such as by diluting the inter-group connectivity [76], or by 697 accounting for inhibitory feedback [18,77] or synaptic failure [78]. 698 Despite the required large groups in our results, a dilution of inter-group 699 connectivity to 100 to each receiving neuron in the following group guarantees that the 700 activity runs stably through the chain ( [76], and thus our assumption of stably running 701 chains to find the amount of STPs as in the data would still be met. 702 Trengove et al. [18] were able to simulate networks in which 800,000 excitatory and 703 200,000 inhibitory neurons participate in a total of 51.020 SFC groups. Each neuron 704 participates in ∼ 70 chains. In their model, each SFC group activates a random set of 705 inhibitory neurons, which regulates the network activity when a SFC is activated. Our 706 scenario with 840,000 neurons participating on average in 5-30 SFC groups (see 707 Section 3.2.3), depending on the parameter set, is close to the parameters of [18], and 708 thus our setting appears to be feasible in dynamical simulations given that one 709 implements a regulatory mechanism which is as effective as the inhibitory feedback 710 employed in [18].

711
Since the electrodes of the Utah array have the same length, they likely recorded 712 from neurons at the same cortical depth. So we have to assume that the observed 713 phenomena are largely due to horizontal connectivity. The connection probability of two 714 neurons mainly depends on the horizontal distance between their somata and drops off 715 with increasing horizontal distance [79][80][81][82][83]. 80 − 82 % of the projection partners of a 716 neuron in layer 2/3 commonly lie within a distance of 500 µm [79,84,85]. A previous 717 SFC embedding study used a maximal distance of 300 µm between connected neurons 718 as a constraint [76]. Our simulations that match the experimental data showed that 719 large group radii (r group ≥ 900 µm) and large inter-group distances (σ group distance ≥ 500 720 µm) are required. As discussed above, such SFCs would have to be diluted, i.e. not fully 721 connected, to function dynamically. This would help to stay within the anatomical which would also account for the large group radii required in our model. However, a 726 more recent study [84] showed in a dynamic photo-stimulation experiment lateral 727 connections and high response reliability up to a range of 1, 500 µm. These estimates 728 cover the connection distances required for our optimal parameter set and thus 729 constitute another candidate mechanism by which such chains could be realized.

730
Currently, we always fix all parameters to single values within a single simulation.

731
However, a hybrid scenario of broad and narrow spatial profiles could be interesting and 732 more realistic. SFCs with a broad spatial profile could serve information propagation 733 within and between cortical areas, but in addition to this, there could be much more 734 localized SFCs for local computations, which are connected to the grid of broader 735 chains. Such a hybrid scenario, and interlinked chains in general, would increase the 736 chance to find larger patterns. In [11,12], the authors found that information 737 propagation along an SFC can be effectively gated by a second SFC, and, according 738 to [16], the representation of information in spatially propagating structures is 739 advantageous for information binding (cf. [14,15]) across different cortical areas. These 740 findings indicate that SFCs are a solid candidate for information propagation and processing of complex sequences of data.

748
In the same data set that we analyzed for STPs, [90] find a spatial structure of 749 correlated neurons (measured by covariance) which changes across behavioral contexts. 750 The covariances are partly strong within neurons up to few mm distance. This goes well 751 with the fact that we found STPs involving neurons across all possible distances on the 752 4 × 4 mm 2 Utah array. However, it is currently not clear if neurons involved in an active 753 SFC also increase their rates. Grün (private communication, based on data from [91]) 754 showed that neurons involved in an SFC seem to go hand in hand with coherent rate 755 changes as visible in the population raster, however, when neurons are sorted according 756 to their spike times in the SFC, it turns out that these are actually single well timed 757 spikes.

758
A limitation of the SPADE analysis lays in the fact that the method is able to detect 759 only spike patterns repeating exactly (within few ms imprecision) across the data. Thus, 760 it does not allow for temporal imprecision in the spike sequences over a few milliseconds, 761 and does not allow for selective participation of the pattern neurons. In other words, 762 single pattern instances, where some spikes may be missing due to synaptic failure, due 763 to e.g. imprecision in the spike sorting, are not detected by SPADE. However, other 764 approaches relying on different methodologies may detect such "fuzzy patterns" [92][93][94][95][96]. 765 Thus, the results obtained by the analysis of the experimental data may be compared to 766 analysis of other pattern detection methods, perhaps leading to more patterns and/or to 767 somewhat different statistics.

768
In order to keep the model as simple as possible, we did not include the significance 769 testing performed in SPADE [51,54] either. Modeling it more closely could also improve 770 the match of the pattern size distributions, however, it is almost impossible without 771 fully dynamical simulations, since the rates of the neurons and thus the numbers of 772 occurrences of the patterns are important factors required for the significance test.

773
Our model can also be adapted to different recording settings. The electrode 774 positions and sensitive ranges, depth of the recorded cortical layer and the local neuron 775 density however have to be provided. Using other recording techniques and applying our 776 analyses could help to answer remaining and follow-up questions. Neuropixels probes provide a very good spatial resolution orthogonal to the surface of the cortex. Data 778 recorded with such probes could be analyzed and used to model SFCs across multiple 779 cortical layers. A high sampling resolution along the other dimensions is however much 780 more difficult to achieve with such probes since it would require multiple probes to be 781 implanted close to one another.

782
STPs have been detected in many studies across different animals, cortical areas, 783 temporal scales and task conditions. Riehle et al. [29] found synchronous patterns at 784 330 µm to 660 µm distance that would correspond to neurons recorded from a single 785 SFC group. The implied group radius is in accordance with many parameter sets that 786 provide a good fit to the data we analyzed, too. They also found synchronous 787 activations related to behavior of different subsets within a set of three neurons, hinting 788 at multiple participation of neurons in SFCs which we also assume to fit our data. In 789 the same data we analyzed, Torre et al. [31] also found synchronous firing of neurons at 790 distances of 400 µm to 2, 400 µm, implying similar group radii. Prut et al. [28] found 791 spatio-temporal spike patterns with different maximum temporal extents ranging from 0 792 ms to 500 ms within a cortical area. Modeling the latter with SFCs would require very 793 long chains or much larger delays between groups. They even found different inter-spike 794 delays between spikes of the same neurons in different patterns, adding to the evidence 795 for multiple participation of neurons in SFCs. Hemberger et al. [32] found fuzzy spiking 796 sequences which could be triggered by single spikes. The fact that a single spike from 797 one neuron suffices to trigger a reliable spiking sequence also hints at the possibility of 798 very narrow chains with few neurons per group and with strong synapses, which can 799 still be reliable if a single spike is already enough to propagate the activity to the 800 following groups.

801
As follow-ups to this work, SFC spatial embeddings should be complemented by 802 including spiking activity to study the dynamics. A closer modeling of pattern 803 significance testing may also help to match the experimentally found pattern size

833
This isolation probability applies to neurons which are inside the sensitive range of 834 an electrode. So, to get the total number of isolated neurons, we have to multiply this 835 probability by the number of neurons which are inside the sensitive ranges of electrodes. 836 Assuming a homogeneous neuron density, this number exactly corresponds to the 837 denominator in Eq 10: Thus, the total number of isolated neurons is independent of our parameter r sens : 839 N isolated = P 3D (isolation) · N in range = const.
This does not mean that the actual sensitive range of the electrodes does not matter, 840 it is hidden in the 1.1 single units /electrode in Eq 10. Our model parameter r sens is not 841 equal to the actual sensitive range, and it does not have to be, since via Eq 10 we 842 ensure that the overall isolation probability of neurons per electrode in our model 843 matches the experimental data. In order to perform the simulations, we have to assign a 844 value to r sens in order to be able to select neurons close to electrode tips, and we fix 845 r sens = 50 µm.