Spatial representability of neuronal activity

A common approach to interpreting spiking activity is based on identifying the firing fields—regions in physical or configuration spaces that elicit responses of neurons. Common examples include hippocampal place cells that fire at preferred locations in the navigated environment, head direction cells that fire at preferred orientations of the animal’s head, view cells that respond to preferred spots in the visual field, etc. In all these cases, firing fields were discovered empirically, by trial and error. We argue that the existence and a number of properties of the firing fields can be established theoretically, through topological analyses of the neuronal spiking activity. In particular, we use Leray criterion powered by persistent homology theory, Eckhoff conditions and Region Connection Calculus to verify consistency of neuronal responses with a single coherent representation of space.

A common approach to interpreting spiking activity is based on identifying the firing fields-regions in physical or configuration spaces that elicit responses of neurons. Common examples include hippocampal place cells that fire at preferred locations in the navigated environment, head direction cells that fire at preferred orientations of the animal's head, view cells that respond to preferred spots in the visual field, etc. In all these cases, firing fields were discovered empirically, by trial and error. We argue that the existence and a number of properties of the firing fields can be established theoretically, through topological analyses of the neuronal spiking activity. In particular, we use Leray criterion powered by persistent homology theory, Eckhoff conditions and Region Connection Calculus to verify consistency of neuronal responses with a single coherent representation of space.
Physiological mechanisms underlying the brain's ability to process spatial information are discovered by relating parameters of neuronal spiking with characteristics of the external world. In many cases, it is possible to link neuronal activity to geometric or topological aspects of a certain space-either physical or auxiliary. For example, a key insight into neuronal computations implemented by the mammalian hippocampus is due to O'Keefe and Dostrovsky's discovery of a correlation between the firing rate of principal neurons in rodents' hippocampi and the animals' spatial location [1][2][3] . This discovery allowed interpreting these neurons' spiking activity, henceforth called place cells, as representations of spatial domains-their respective place fields (Fig. 1A 4 ) (Throughout the text, terminological definitions are given in italics.). It then became possible to use place field layout in the navigated environment E-the place field map M E -to decode the animal's ongoing location [5][6][7][8] , and even to interpret the place cells' off-line activity during quiescent stages of behavior or in sleep [9][10][11][12][13][14] , which define our current understanding of the hippocampus' contribution to spatial awareness 15-18 . In the 90s, a similar line of arguments was applied to cells discovered in rat's postsubiculum and in other parts of the brain [19][20][21] , which fire at a particular orientation of the animal's head. The angular domains where such head direction cells become active can be viewed as one-dimensional (1D) head direction fields in the circular space of planar directions, S 1 -in direct analogy with the hippocampal place fields in the navigated space (Fig. 1B). The corresponding head direction map, M S 1 , defines the order in which the head direction cells spike during the rat's movements and the role of these cells in spatial orientation [20][21][22] . Recently, place cells and head directions cells were discovered in bats' hippocampi; in contrast with rodents who navigate two-dimensional (2D) surfaces (see however [23][24][25][26] ), bat's voluminous place fields cover three-dimensional (3D) environments and their head direction fields cover 2D tori 27,28 .
The spatial view cells, discovered in the late 90s, activate when a primate is looking at their preferred spots in the environment (Fig. 1C), regardless of the head direction or location [29][30][31] . Correlating these cells' spike timing with the positions of the view fields helped understanding mechanisms of storing and retrieving episodic memories, remembering object locations, etc. 32,33 . The principles of information processing in sensory and somatosensory cortices were also deciphered in terms of receptive fields-domains in sensory spaces, whose stimulation elicits in spiking responses of the corresponding neurons [34][35][36][37][38][39] .
In all these cases, referencing an individual neuron's activity to a particular domain in a suitable representing space X 40 is key for understanding its contribution and for reasoning about functions of neuronal ensembles in terms of the corresponding "maps" [16][17][18] . This raises a natural question: when is a "spatial" interpretation of Implementation. As it turns out, representable simplicial complexes exhibit several characteristic properties that distinguish them among generic simplicial complexes 49,50 . Verifying these properties over biologically relevant 1D, 2D and 3D representing spaces is a tractable problem 49,51 , although exact algorithms for performing such a verification are not known-only in 1D are some methods available [53][54][55][56][57] . Nevertheless, there exist explicit criteria that allow limiting the dimensionality of the representing space X and eliminating manifestly non-representable complexes based on their homologies, combinatorics of simplexes and other intrinsic topological properties, which will be used below.
Specifically, according to the Leray criterion, a complex representable in D dimensions should not contain non-contractible gaps, cavities or other topological defects in dimensionalities higher than (D − 1) 58 . Formally, it is required that the homological groups of and hence its Betti numbers should vanish in these dimensions, b i≥D (�) = 0 . Moreover, the Betti numbers of all the subcomplexes x of , induced by a fraction x of Figure 1. Spatial maps. (A) A simulated place field map of a small ( 1m × 1m ) environment E , similar to the arenas used in typical electrophysiological experiments 66,67 . Dots represent spikes produced by the individual cells (color-coded); their locations mark the rat's position at the time of spiking. The pool of place cell coactivities is schematically represented by a coactivity complex T PC (top right). The navigated trajectory r(t) induces a sequence of activated simplexes-a simplicial path Ŵ ∈ T PC . (B) The head direction cell combinations ignited during navigation induce a coactivity complex T HD (top). The corresponding head direction fields cover a unit circle-the space of directions (bottom). (C) Spatial view cells activate when the primate gazes at their respective preferred domains in the visual field (left). The curves r 1 (t) and r 2 (t) traced by the monkey's gaze induce simplicial paths Ŵ 1 and Ŵ 2 running through the corresponding coactivity complex T VC (right). www.nature.com/scientificreports/ its vertexes should also vanish, b i≥D (� x ) = 0 . In the case of coactivity complexes, such subcomplexes T x ⊆ T have a particularly transparent interpretation: they are the ones generated by x% of the active cells. According to the second criterion, the number of simplexes in all dimensions of must obey Eckhoff 's inequalities-a set of combinatorial relationships discussed in [59][60][61][62] and listed in the Methods" section, where we also briefly detail the Leray criterion 49,58,[62][63][64] . Previous topological studies of the coactivity data were motivated by the Alexandrov-Čech theorem [42][43][44][45] , which posits that the homologies of the nerve complexes produced by the "good" covers (i.e., the ones with contractible overlaps (4), see 65 ), should match the homologies of the underlying space X, H * (N ) = H * (X) , i.e., have the same number of pieces, holes, tunnels, etc. Specifically, this construction was applied to the place cell coactivity complexes, whose representability was presumed [46][47][48] . Persistent homology theory [68][69][70][71][72][73] was used to trace the dynamics of the Betti numbers b i≤D (T ) in physical dimensionalities D ≤ 2 74-79 and D ≤ 3 79,80 , to detect whether and when these numbers match the Betti numbers of the environment, b i≤D (E) , and how this dynamics depends on spiking parameters. It was demonstrated, e.g., that for a wide range of the firing rates and place field sizes referred to as the Learning Region, L(E) , the low-dimensional Betti numbers of T (t) converge to their physical values after a certain period T min , neurobiologically interpreted as the minimal time required to "learn" the topology of the environment ( b i≤D (T (t)) = b i≤D (E) , t ≥ T min 48 ).
Together, these arguments suggest that experimentally discovered representing spaces and firing fields serve as explicit models of the cognitive maps emerging from neuronal activity-a perspective that is currently widely accepted in neuroscience. However, this view requires verification, since the empirically identified firing fields may be contextual offshoots or projections from some higher-dimensional constructs-in the words of H. Eichenbaum, "hippocampal representations are maps of cognition, not maps of physical space" 81 . The way of addressing this question is straightforward: if the spiking activity is intrinsically spatial, i.e., if neurons represent spatial domains, then the coactivity complexes generated by the corresponding neuronal ensembles should be representable-an explicit property that can be confirmed or refuted using Leray, Eckhoff and other criteria. In the following, we apply these criteria to several types of neuronal activity, both simulated and experimentally recorded, and discuss the results.

Results
Simplicial topology approach. The conventional theory of representability addresses properties of "static" simplicial complexes [49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65] . In contrast, the coactivity complexes are dynamic structures that can be viewed as time-ordered agglomerates of simplexes, restructuring at the moments t 1 < t 2 < t 3 . . ., The exact organization of each complex in the sequence (7) depends on the specifics of the underlying spiking activity, e.g., the initial state of the network, its subsequent dynamics, spiking mechanisms and so forth (in case of the place fields, think of the starting point of navigation, shape of the trajectory, speed, etc.). Thus, verifying representability of these complexes requires testing whether Eckhoff, Leray and other criteria are valid at each moment t.
We constructed coactivity complexes by simulating the rat's navigation through a planar environment E commonly used in electrophysiological experiments (Fig. 1A, see also 66,67 ). The neuronal spikings in this case are generated as responses to the rat's appearances within preconstructed, convex firing domains, e.g., stepping into randomly scattered place fields or facing towards head direction fields centered around randomly chosen preferred angles (see Methods" section and Methods in 48,74,82 ). While the resulting nerve complexes (6) are 2D-representable by design, we inquired whether the corresponding coactivity complexes are also representable, i.e., whether the activity of individual neurons intrinsically represents regions and whether connectivity between these regions is similar to the connectivity between the underlying auxiliary firing fields.
Simulations show that persistent Leray dimensionality D L (above which the spurious loops in T (t) vanish, 83 ) eventually settles at D L = 1 for most complexes, implying that neuronal activity defines a proper planar map. However, this mapping requires time-a Leray period T Lwhich, for the maps populating the learning region L(E) , is typically similar to the learning time T min ( Fig. 2A,B).
Whether a particular value of T L is shorter or longer than the corresponding T min depends on how exactly the coactivity complex is constructed, e.g., whether the simplexes (1) correspond to simultaneously igniting cells groups or assembled from lower-order combinations over an extended period ̟ 84 . Physiologically, the former corresponds to the case when spiking outputs are processed by "coincidence-detector" neurons in the downstream networks and the latter to the case when lower-order coactivities are collected over a certain "spike integration window" ̟-longer than the simultaneity detection timescale w [85][86][87] . Different readout neurons or networks may have different integration periods; to simplify the model, we started by extending the parameter ̟ to the entire navigation period for all cells and cell groups.
The lowest order of coactivity involves spiking cell pairs 88 , which together define a coactivity graph G 89,90 . The cliques ς of this graph produce a clique coactivity complex T ς that generalizes the simplicial coactivity complex T σ , built from simultaneously detected simplexes [75][76][77][78][79][80] . As it turns out, the "coincidence detecting" and the "spike integrating" complexes have different topological dynamics: the former are more likely to start off with a higher www.nature.com/scientificreports/ Leray dimensionality, D L ≥ 3 , that then reduces to D L ≤ 2 ( Fig. 2A), whereas the latter tend to be more stable, lower-dimensional and have shorter Leray and learning times (Fig. 2B).
To test the induced subcomplexes of each T , we selected random subcollections of cells containing x = 50% , x = 33% , x = 25% and x = 20% of the original neuronal ensemble, and found that if the original complex T ≡ T x=1 is representable, then its subcomplexes, T x<1 ⊆ T , typically require less time to pass the Leray criterion, T L (T x<1 ) ≤ T L (T ) , and that for x > 50% the Leray times saturate, T L (T x>0.5 ) ≈ T L (T ) . Thus, the Leray time of the full complex, T L (T ) , can be used as a general estimate of the timescale required to establish representability.
To control the sizes of the coactivity complexes, we used only those periods of each neuron's activity when it fired at least m spikes per coactivity window (i.e., time intervals defining simultaneity of neuronal activity, w ≈ 1/4 secs; for justification of this value see 74,92 ). Additionally, we used only those groups of coactive cells in which pairwise coactivity exceeded a threshold µ (Methods" section). Biologically, these selections correspond to using only the most robustly firing cells and cell assemblies for constructing the coactivity complexes 77 . The results demonstrate that majority of the coactivity complexes T computed for smallest possible m and µ exhibit low persistent Leray dimensionality, D L = 1 , which points at 2D representability of the underlying neuronal activity, with the Leray times T L similar to the corresponding learning times T min (Fig. 2C). We also found that Eckhoff inequalities are typically satisfied throughout the navigation period, i.e., that the Eckhoff criterion does not significantly limit the scope of representable spiking in this case (Fig. 2D).

Region connection calculus (RCC).
An independent perspective on spatial representability is provided by Qualitative Space Representation approach (QSR 93,94 ), which sheds a new light on the dynamics of neuronal maps. From QSR's perspective, a population of cells C = {c 1 , c 2 , . . . , c N } may represent a set of abstract, or formal spatial regions R = {r 1 , r 2 , . . . , r N } , if the relationships between them, as defined by the cells' coactivity, can be consistently actualized in a topological space X by a set of explicit regions, ϒ = {υ 1 , υ 2 , . . . υ N }.
Specifically, regions r i and r j encoded by the cells c i and c j can be: 1. disconnected, DR(r i , r j ) , if c i and c j never cofire; 2. equal, EQ(r i , r j ) , if c i and c j are always active and inactive together; 3. proper part of one another, if c j is active whenever c i is, PP(r i , r j ) , or vice versa, PPi(r i , r j ); 4. partially overlapping, PO(r i , r j ) , if c i and c j are sometimes (but not always) coactive.
These five relations fully capture mereological configurations of regions in a first-order logical calculus known as RCC5 (Fig. 3A 95 ). Using mereological (i.e., parthood-based see 96 and Methods" section), rather than topological, distinctions reflects softness of the firing fields' boundaries: the probabilistic nature of neuronal spiking does not warrant determining whether two regions actually abut each other or not.

com/scientificreports/
A key property of a RCC5-framework defined by spiking neurons-a R 5 schema-is its internal consistency 40 . It may turn out, e.g., that some pairs of cells encode relationships that are impossible to reconcile, e.g., PP(r i , r j ) , DR(r j , r k ) and PO(r i , r k ) . Indeed, if an actual region υ i is contained in υ j then it cannot possibly overlap with a region υ k that is disconnected from υ j . Correspondingly, the neuronal activity that produces such inconsistencies (for the full list see Table 1 in Methods" section) is not representable-not even interpretable in spatial terms. On the other hand, it can be shown that if all triples of relationships are consistent, then R 5 (t) does possess a spatial model, i.e., there exists a set of regions υ i (with no prespecified properties such as convexity, connectivity or dimensionality) that relate to each other as the r i s relate in R 5 [93][94][95][96][97][98][99] .
To verify whether spiking activity is representable in this QSR sense, we constructed an inflating R 5 (t) -schema (an RCC5-framework growing as spiking data accumulates, similar to (2)) for each neuronal ensemble and counted the inconsistent triples of relationships at each moment t. The results show that all R 5 (t)-schemas start off with numerous inconsistencies, which tend to disappear after a certain period T RCC5 that is typically smaller than the Leray time T L (Fig. 3B,C).
The net dynamics of RCC5 relationships is illustrated on Fig. 3C. Note that some of these changes may be attributed to the regions' continuous reshapings or displacements, e.g., two overlapping regions may become disconnected, PO(r i , r j ) → DR(r i , r j ) , is r i moves away from r j , or r i may move into r j , inducing PO(r i , r j ) → PP(r i , r j ) . In contrast, a jump from a disconnect to a containment, without an intermediate partial overlap, e.g., DR(r i , r j ) → PP(r i , r j ) rather than DR(r i , r j ) → PO(r i , r j ) → PP(r i , r j ) , would be a discontinuous, abrupt change. As shown on Fig. 3C, discontinuous transitions are common at the initial stages of navigation, but shortly before T RCC5 they disappear, indicating that the relationships between regions encoded within a sufficiently well-developed R 5 (t) schema evolve in a continuous manner.
These outcomes not only provide an alternative lower-bound estimate for the time required to accumulate data for producing low-dimensional spatial representations, but also help understanding the nature of processes taking place prior to Leray time. In particular, the exuberant initial dynamics, homologically manifested through  Initially, discontinuous events are frequent but shortly before T RCC5 they disappear entirely, leaving the stage to qualitatively continuous sequences. The same barcodes are added in the background, error margins shown by dashed lines. Table 1. RCC5 compositions. Given three regions, x, y and z, and two relationships R 1 (x, y) and R 2 (y, z) , the relationship R 3 (x, z) is not arbitrary. A map is consistent, if every triple of relationships is RCC5-consistent. www.nature.com/scientificreports/ an incipient outburst of spurious loops in the coactivity complexes ( Fig. 2A,B, t < T min ), cannot be interpreted as a mere "settling" of topological fluctuations in the cognitive map-according to Fig. 3B, the RCC5-schema does not form a coherent topological stratum for t < T RCC5 . Rather, the initial disorderly period should be viewed as the time of transition from a nonspatial to a spatial phase, followed by spatial dynamics (for t > T RCC5 ) that involves, inter alia, dimensionality reduction and other restructurings (Fig. 3C).
Current summary. Taken together, these results show that even in the simplest "reactive" model, in which neuronal firings are simulated as responses to regular domains covering a compact space, the low-dimensional representability is not an inherent, but an emergent property. In particular, RCC5-analyses suggest that spatial interpretation of neuronal spiking becomes possible after a finite period. During the times that exceed both T RCC5 and T L , the spiking data can be interpreted in terms of firing fields in a space X of dimensionality higher than the persistent Leray dimensionality of the corresponding coactivity complex, dim(X) >D L (T ).
Physiologically, this implies that the outputs of place cells, head direction cells, view cells, etc., may not be immediately interpretable by the downstream networks as representations of spatial regions-the information required for such inference appears only after a certain "evidence integration. " Correspondingly, the firing field maps constructed according to the standard experimental procedures 1,19,30 also cannot be considered as automatic "proxies" of cognitive maps: such interpretations are appropriate only after the representability of the corresponding coactivity complexes is established. Another principal conclusion is that representability of the spiking activity depends not only on the spiking outputs, but also on how the information carried by these spikes is detected and processed. In particular, spikes integrated over extended periods are likelier to permit a consistent firing field interpretation than spikes counted via coactivity detection. On the technical side, these results imply that an accurate description of the firing fields' plasticity should include possible dimensionality changes 101-103 . Multiply connected place fields. A key simplification used in the simulations described above is that firing fields were modeled as convex regions. While this assumption is valid in some cases 100 , multiply connected firing fields are also commonly observed (Fig. 4A 104,105 ). From our current perspective, the issue is that multiple connectivity of the cover elements (3) may increase the Leray dimensionality of the corresponding nerve complex 49,106 and thus bring additional ambiguity into the analyses. Identification of the firing fields' con-  Fig. 2A, can reach D(T σ ) = 4 if we allow 30% of multiply connected place fields ( 2 − 3 components each). (C) In a clique coactivity complex, the spurious loops in dimensions D = 2 and lower may persist indefinitely, implying either that the firing fields are 3D-representable or that they may be multiply connected. Note that the number of spurious loops in both T σ and in T ς is higher than in the case with convex firing fields ( Fig. 2A,B). (D) The persistence bars computed for the flickering complex F ̟ with spike integration window ̟ = 1 minute, indicate stable mean Leray dimensionality �D(F ̟ )� = 1 , implying that the local charts χ ̟ are planar and hence that the firing fields are two-dimensional. www.nature.com/scientificreports/ nectivity from the spiking data is an elaborate task that requires tedious analyses of the spike trains produced by individual cells or cell groups over periods comparable to the Leray and the learning times 83,109 . To circumvent these difficulties, we reasoned as follows.
Suppose that the spiking activity used to produce a coactivity complex T (t) is generated as a moving agent (animal's body, its head, its gaze) follows a trajectory γ (t) over a space X, covered with stable firing fields υ i , i = 1, . . . , N . Consider a navigation period ̟ that spans over a smaller segment of this trajectory, γ ̟ = {γ (t) : t ∈ ̟ } . If ̟ is sufficiently short, then one would expect γ ̟ to cross at most one component of a typical firing field υ i ; even if γ ̟ meets more than one component of a multiply connected υ i , this property may not manifest itself in the resulting spike trains, i.e., υ i should be effectively simply connected (Fig. 4A). Correspondingly, the Leray dimensionality of the coactivity complex acquired during that period should reflect the dimensionality of a small underlying fragment of X-a local chart χ ̟ -that contains γ ̟ (topologically, χ ̟ (γ ) ∼ = {∪ j υ j |γ ̟ ∩ υ j � = ∅} ). The dimensionality of χ ̟ can then be ascribed to all the contributing υ j s, dim(υ j ) = dim(χ ̟ ).
Further, if the ̟-period is allowed to shift in time, then the segment γ ̟ will also slide along the trajectory γ (t) ; the spikes fired within each t-centered window, ̟ t = [t − ̟/2, t + ̟/2] , will then produce a ̟ t -specific flickering coactivity complex F ̟ (t) ⊆ T (t) , whose topological properties may change with time 78,110,111 . Since F ̟ (t) contains a finite number of elements, it will reconfigure at discrete moments, t 1 , t 2 , . . . , and remain unchanged in-between, F ̟ (t) = F ̟ (t k ) , t ∈ [t k , t k+1 ) . If a given instantaneous configuration F ̟ (t k ) is representable, then its vertexes correspond to the regions comprising the local chart χ ̟ (t k ) , with dimensionality dim(χ ̟ (t k )) ≥D L (F ̟ (t k )) . If two such complexes overlap, F ̟ (t k ) ∩ F ̟ (t l ) � = ∅ (i.e., their vertex sets overlap), then their respective charts also overlap χ ̟ (t k ) ∩ χ ̟ (t k+1 ) � = ∅ , which allows relating their topological properties, including properties of the representing regions.
Clearly, the outcome may depend on how each γ ̟ is embedded into X, the spiking parameters, etc. Moreover, since the Leray dimensionality of the instantaneous complexes can change, so can the dimensionalities of the corresponding local charts: ) . This may seem as a contradiction since the representing space is naturally assumed to be a topological manifold, i.e., all of its local charts, arbitrarily selected, should have the same dimensionality dim(χ ̟ (t)) = dim(X) = D . On the other hand, the deviations of the local dimensionality estimates from a fixed D can be viewed as mere fluctuations caused by occasional contribution of multiply connected firing fields or by other noise sources, e.g., by stochasticity of neuronal spiking 112 . One can hence attempt to discover the true dimensionality of X by evaluating the mean Leray dimensionality of the instantaneous complexes, which physiologically alludes to learning the physical structure of the underlying space from the recurrent information.
Numerical verification of the viability of the proposed approaches can be achieved by simulating multiply connected firing fields and computing homological dynamics of the resulting coactivity complexes. To that end, we randomly added 2 − 3 additional convex components to ∼ 30% of the place fields (Fig. 4A) and simulated the topological dynamics of the corresponding complexes.
The results show that multiple connectivity of the firing fields does indeed increase Leray dimensionality in both the detector and the integrator complexes, T σ (t) and T ς (t) . Moreover, in contrast with the complexes generated off the maps with convex fields, the maps with multiply connected fields tend to produce persistent higher-dimensional loops, notably in the coactivity detecting complexes T σ (t) (compare Figs. 2A and 4B). In the spike integrating clique complex T ς (t) , the Leray dimensionality remains low and may in some cases retain the physical value D L (T ς ) = 1 , although topological loops in dimensions D = 2 and even higher may also appear (Fig. 4C). Thus, multiple firing field connectivity significantly increases the number of spurious 1D holes (by 200-300% ), precluding both types of complexes from assuming the physically expected topological shapes.
Tighter dimensionality estimates can be produced by using shorter spike integration windows ̟ T min and constructing flickering coactivity complexes F ̟ (t) from pairwise coactivities detected over ̟-periods shifting by discrete steps �̟ and yielding an array of windows ̟ 1 , ̟ 2 , ̟ 3 , . . . centered at t k = ̟/2 + (k − 1)�̟ . The specific ̟-values were chosen comparable to the characteristic time required by the rat to run through a small segment of the environment: ̟ ≈ 25-65 s for the arena shown on Figs. 1A and 4A. The Betti numbers for this case were evaluated using zigzag homology theory-a generalization of the persistent homology theory that applies to complexes that can not only grow, but also shrink, break apart, fuse back again, etc. 113,114 . In particular, this approach allows studying how the topological fluctuations in F ̟ (t) affect its Leray dimensionality D L (F ̟ (t)) from moment to moment.
Typical results illustrated on Fig. 4D show that there appears a large number of spurious 0D loops-disconnected pieces-with lifetimes nearly exponentially distributed about the learning periods T ς min , which suggests that fragments of F ̟ (t) appear and disappear at random over such periods. The transient 1D loops also form and decay at ̟-timescale. However, the most important outcome is that the topological dynamics in dimensions D > 1 trivializes-the higher dimensional loops in F ̟ occur very rarely, if ever. These properties are qualitatively unaffected by varying the discretization step �̟ ( ̟/20 �̟ ̟/10 ) or changing the window width ̟ , i.e., the estimates of the mean Leray dimensionality �D L (F ̟ )� = 1 are stable and reveal physical planarity of the representing space.
Verification of the RCC5-consistency of the spiking data produces the same qualitative results as in the case with simply connected firing domains: the R 5 (t)-schemas become consistent after a learning period T RCC5 < T ς min , upon which neuronal activity becomes spatially interpretable, and, by the Leray and Eckhoff arguments, representable in dimensions D ≥ �D L (F ̟ )�. www.nature.com/scientificreports/ Electrophysiological data. We applied the analyses described above to spiking activity recorded in the hippocampus (CA1 area) of rats navigating a linear environment shown on Fig. 5A (for more data description and experimental specifications see 115 ). A typical running session, during which the animal performed 45 − 70 laps between the tips of the track, provided N c 25 simultaneously recorded neurons, allowing to construct small coactivity complexes that quickly become RCC5-consistent, comply with the Eckhoff conditions, and exhibit persistent Leray dimensionality, D L = 0 , with typical persistent Leray time T L ≈ 10 mins (Fig.5B).
The vanishing D L indicates that a linear track illustrated in Fig. 5A is contractible and implies 1D-representability. The latter can also be tested independently via RCC5 analyses, which in this case allows identifying the track's linear structure 116 .
Since some of the hippocampal place fields are multiply connected, we also applied sliding window analyses, adjusting the spike integration period ̟ i to match the duration of the animal's i th run from one end of the track to the other (typically 2 ≤ ̟ i ≤ 10 secs). Computations reveal that the resulting complexes F ̟ i (t) exhibit the same mean Leray dimensionality �D L (F)� = 0 , which is consistent with the persistent D L estimates above. Combining these results produces convergent evidence that in this case the hippocampus does indeed map out a 1D spatial domain (rather than 2D, see 115 ).
The latter conclusion can, in fact, be verified by yet another representability test, which applies only to 1D cases and presumes firing field convexity. The Golumbic-Fishburn (GF) algorithm 53-57 is based on computing a binary index p: the 1D-representable simplicial complexes yield p(�) = +1 , and the non-representable complexes produce p(�) = −1 . For the inflating or flickering coactivity complexes this index becomes time-dependent, p = p(t) , marking the evolution of 1D representability (Methods" section). Applying the GF-algorithm to the inflating coactivity complexes constructed for cells with convex place fields only, we found that 1D representability, p(t) = +1 , appears in about T + ≈ 10 mins, close to the Leray time (Fig. 5B,C), demonstrating consistency with the previously obtained results.
Lastly, we addressed a particular property of the place cell's spiking activity in linear environments-the place fields' directionality: a given place cell may fire during the outbound, but not inbound directions, or vice versa 117 . We verified that the topological dynamics exhibited by the coactivity complexes built using only the outbound or only the inbound activity are very similar to the dynamics of the full (bidirectional) complex T ς (t) (Fig. 5B,C), implying that place cell directionality does not necessarily compromise 1D representability of spiking activity.

Discussion
Topological analyses of the spiking data allow testing whether a given type neuronal activity may arise from a "spatial map, " i.e., whether each neuron's spiking marks a domain similar to a place field, a head direction field, a view field, etc., in a certain low-dimensional space. Thus far, establishing correspondences between neurons and firing fields was based on matching the spike trains with spatial domains empirically, through trial and error 1,19,30 . Here we attempt to address this question in a principled way, through intrinsic analyses of the spiking data, without presuming or referencing ad hoc constructions. A set of hands-off algorithms discussed above allows objective estimates for the dimensionality of a space needed to model the patterns of neuronal firing-a method that is unaffected by technical limitations, experimental ingenuity or complexity (e.g., nonlinearity) of the required firing field arrangements.
To follow the dynamics of the coactivity complexes we extend the conventional approaches of representability theory into the temporal domain, obtaining several complementary time-dependent markers of representability. In particular, we use persistent homologies to extend Leray's theory to the case of inflating simplicial complexes www.nature.com/scientificreports/ and zigzag homologies in the case of flickering simplicial complexes. The latter approach is especially valuable as it allows extracting stable topological information from spiking data that may be generated from the maps with multiply connected firing fields or encumbered by other inherent irregularities, in spirit with the general ideas of topological persistence [68][69][70][71][72][73] . It should also be mentioned that mathematical discussions of the persistent nerve theorem, alternative to ours and more formal, have began to appear 118,119 ; however at this point our studies are independent. A principal observation suggested by our analyses is that representability is a dynamic, emergent property that characterizes current information supplied by the neuronal activity. Moreover, representability depends not only on the amount and the quality of the spiking data itself, but also on the mechanisms used for processing and interpreting this data. Both aspects affect the time required to establish the existence of a representing space and its dimensionality. An implication of this observation is that experimentally constructed firing field maps (place field maps, head direction maps, etc.) cannot be automatically regarded as direct models of cognitive representations of ambient spaces [15][16][17][18][19][20][21] or more general spatial frameworks 120 ; correctness of such interpretations may require more nuanced considerations.

Physiological parameters and constructions.
• Simulated trajectory r(t) = (x(t), y(t)) , used for generating coactivity complexes was obtained by modeling a rat's non-preferential exploratory behavior-navigation without favoring of one segment of the environment E over another (Fig. 1A). The mean speed of about ∼ 20 cm/sec was selected to match experimentally recorded speeds. The direction of the velocity v(t) = (v x (t), v y (t)) defines the "angular trajectory" ϕ(t) = arctan v y (t)/v x (t) that traverses the space of directions, S 1 , allowing to simulate head direction cell activity as the rat explores E 48,74,82 . The simulated navigation period, T = 25 minutes, was selected to match the duration of a typical "running session" in electrophysiological experiments 100 . A shorter spike integration window ̟ ≪ T was used to limit the pool of spiking data for time-localized computations. • Poisson spiking rate of a place cell p depends on the animal's location r(t), where f p is the cell's maximal firing rate and s p defines the size of its place field 101 . A similar formula defines the firing rate of a head direction cell h, h (ϕ) , as a function of the animal's ongoing orientation ϕ , the cell's preferred orientation angle ϕ h , its maximal rate f h and the size of its preferred angular domain s h . In all simulations the firing fields were stable, i.e., the parameters of c and h remained constant.
• Neuronal ensembles produce lognormal distributions of the maximal firing rate amplitudes, f c , and of the firing field sizes, s c 48,121 . We tested about 17, 000 different ensembles, in which the ensemble mean maximal rate f ranged between 4 and 40 Hz for the place cells and between 5 to 35 Hz for the head direction cells. The ensemble mean firing field sizes varied between 10 to 90 cm for the place fields and between 12 • and 36 • degrees for the angular fields. For all ensembles, the firing field centers were randomly scattered over their respective representing spaces. • Multiple Firing Fields were generated by adding two or three randomly scattered auxiliary spiking centers r c ′ , r c ′′ , etc., The maximal firing rates at the auxiliary locations are smaller than the rate at the main location, f c > f c ′ > . . . , as suggested by the experiments 104,105 . • The activity vector of a cell, m c = [m c,1 , . . . , m c,n ], is constructed by binning its spike trains into w = 1/4 seconds long "coactivity windows" 74,92 . Each m c,k specifies how many spikes were fired by c into the k th time bin, n is defined by the duration of navigation, n = ⌊T/w⌋ . High activity periods can be identified by selecting time bins in which the number of fired spikes exceeds an activity threshold m. Topological propaedeutics. Graphs.
• A graph G is defined by its vertices, V = {v 1 , v 2 , . . . , v n } , and a set of edges E that link certain pairs of vertexes. A formal description of a graph is given by its connectivity matrix C(G), with the elements • A coactivity Graph G is built by establishing functional links between cells that exhibit high activity and coactivity ( m c,k ≥ m , m ij ≥ µ see above) 89,90 . (Fig. 6A). www.nature.com/scientificreports/ • Given a graph G, its complement graph G is produced by flipping 0s and 1s in the connectivity matrix C(G), i.e., joining the disconnected vertexes of G and removing the existing edges.
• Topological analyses of simplicial complexes do not address simplexes' shapes and are based entirely on the combinatorics of the vertexes shared by the simplexes. This motivates using abstract simplexes and abstract simplicial complexes that capture the combinatorial structure of κ (d) s without making references to their geometry. Specifically, an abstract 0-simplex is a vertex σ (0) i ≡ v i , an abstract 1-simplex is a pair of vertexes, σ (1) ij = [v i , v j ] ; an abstract 2-simplex is a triple of vertexes, σ (2) ijk = [v i , v j , v k ] , and so forth (Fig. 6C). Thus, abstract complexes may be viewed as multidimensional generalizations of graphs or as abstractions derived from the geometric simplicial complexes.
• Example 2: The combinations of coactive cells define coactivity simplexes (1), which together form a coactivity complex (Fig. 7B). • Example 3. Vertexes of geometric simplexes that form a geometric simplicial complex K define abstract simplexes that form the corresponding abstract simplicial complex (Fig. 7C). • The set of d-dimensional simplexes of a complex forms its (abstract) d-skeleton, sk d (�).
• A clique complex of an undirected graph G is an abstract simplicial complex formed by the cliques (fully interconnected subgraphs) of G 91 , Fig. 6A. Combinatorial properties of cliques are the same as simplexes': a subset of a clique's vertexes form a clique, overlap of two G-cliques is also a clique, ς 1 ∩ ς 2 = ς 3 ∈ G (Fig. 6). Thus, any graph G defines a unique clique complex � (G) . Note, that the 1-skeleton of a clique complex yields its underlying graph, sk 1 (�(G)) = G , but if is not a clique complex, then the clique complex built over its 1-skeleton does not reproduce . • Coactivity complexes used in this study are of two kinds. The first kind is formed by the abstract complexes T σ built from simultaneously coactive cell groups (1). The second kind is formed as the clique complexes of the coactivity graphs G 77,80 . The graph (co)activity thresholds m and µ are used to control the size of the complex T m,µ = T (G m,µ ) : selecting m ≥ 2 , µ = 1 for small maps (i.e., counting cells that produce at least two spikes per time bin w) and m ≥ 2 , µ ≥ 5 for larger maps allows computing the full simplicial complex with dimensionality dim(T m,µ ) ≤ 10 , for which we can numerically apply the Javaplex software 122 .

Mereology.
Mereology is the theory of parthood-relations between whole and part, as well as relations between parts within a whole 96,123 . Mereological level of describing spatial regions is cruder than topological, which also includes contact relationships. Since the latter is generally not captured by spiking data 40,116 , we use mereological RCC5 theory 95 .
• Homological groups are designed to "count pieces" in a space X with suitable coefficients. The key property of these groups is that they remain unchanged-invariant-as X is continuously deformed (see 41,42 for a gentle introduction to the subject). If the coefficients form an algebraic field F, then the homological groups, commonly referred to as the "homologies" of X are simply vector spaces H 0 (X, F), H 1 (X, F), . . . , associated with X (one per dimension). Homologies can be easily computed for spaces whose "pieces" are explicitly defined, e.g., for the simplicial complexes, thus providing a way of identifying their topological structures. In practice, it is easier to use just the dimensionalities of H * s-the Betti numbers b k = dim(H k (X, F)) , to count numbers of connectivity components, cavities, tunnels and other topological features of X in different dimensions 41,42 . For example, if X is the boundary of a hollow triangle (or another noncontractible 1D loop), then β 1 (X) = 1 . If X is 1-dimensional complex, i.e., a graph, then β 1 (X) equals to the number of cycles in X, counted up to topological equivalence. If the triangle is "filled", then it can be continuously contracted into a 0D point; since the latter has no topological structure in dimensions d > 0 , the corresponding Betti numbers also vanish. By the same argument a "filled" tetrahedron has β k>0 = 0 , but if the tetrahedron is hollow, then its boundary, being a 2D noncontractible loop (topologically-a 2D sphere) produces β 2 = 1 , β k>2 = 0 . Similarly, for any d-simplex β k>0 (σ (d) ) = 0 , whereas for its hollow boundary, ∂σ (d) , the Betti numbers are β d−1 (∂σ (d) ) = 1 , β k =0,d−1 (∂σ (d) ) = 0 (Fig. 6). Same results apply to the "abstract" counterparts of all these complexes. Note also that continuous deformations of a 0D point x (a 0D topological loop) amount to "sliding" x inside of a space X that contains x; thus β 0 (X) simply counts such "sliding domains", i.e., the number of connected components in X. As a result, all simplexes and simplicial complexes that consist of one piece have β 0 (X) = 1 . • Persistent homology theory allows tracing the topological structure in a filtered family of simplicial complexes, e.g., describing the topological dynamics of the inflating family (2) [70][71][72][73] . The Betti numbers plotted as function of the filtration parameter (in our case it is time, t) form the barcode, b(T , t) = (b 0 (T , t), b 1 (T , t), . . .) , which provides the exact mathematical meaning to the term "topological shape" used throughout the text. Each bar in b(T , t) can be viewed as the corresponding topological loop's timeline 48,[74][75][76][77]80 . • Zigzag Homology theory allows tracking the Betti numbers of the "flickering" complexes-the ones whose simplexes can not only appear, but also disappear (see 113,114 and Supplement in 110 ). In particular, Zigzag homology techniques allow capturing the times when individual loops appear in the flickering complex, how long they persist, when they disappear, reappear again, etc.

Representability.
A generic algorithm for checking whether a given complex can be built as a nerve of a D-dimensional cover is known only for D = 1 (see below). However, there exist criteria that allow ruling out certain non-representable cases.
• The Leray criterion posits that if a complex is a nerve of a D-dimensional cover with contractible overlaps (4), then its rational homologies in dimensions higher or equal than D should vanish, H i≥D (�, Q) = 0 58 . Moreover, homologies of all the subcomplexes ′ ⊆ , induced by selecting vertex subsets of should also vanish, H i≥D (� ′ , Q) = 0 . These properties can be verified by computing the Betti numbers and verifying that b i≥D (�, Q) = 0 . In practice, it is more convenient to carry out the computations over a finite field, such as Z 2 .
Although the b k (�, Q) numbers may in general differ from the b k (�, Z 2 ) numbers, the latter also have to obey the Leray condition and produce the same Leray dimensionality. As an example, the Leray condition poses www.nature.com/scientificreports/ that the boundary of the triangle is not 1-representable ( β 1 (∂σ (2) , Z 2 ) > 0 ), but the triangle itself may be ( β 1 (σ (2) , Z 2 ) = β 2 (σ (2) , Z 2 ) = 0 ); the boundary of a tetrahedron is not 2-representable ( β 2 (∂σ (3) , Z 2 ) > 0 ), but the tetrahedron may be. • Amenta's theorem connects the Leray dimensionality of a simplicial complex to its Helly number, defined as follows. Let ϒ = {υ 1 , υ 2 , . . . , υ n } be a finite family of regions (Fig. 8). The Helly number h = h(ϒ) of the family is defined to be the maximal number of non-overlapping regions, such that every h − 1 among them overlap. For the corresponding nerve complex N (ϒ) , h = h(N ) = h(ϒ) is the number of vertices of the largest simplicial hole in N (i.e., the dimension of the hole plus 2 49 ). This observation can be used to attribute a Helly number to any simplicial complex , h(�) . From the perspective of representability analyses, a key property of the Helly numbers is that they do not exceed d + 1 for a d-Leray complex 49 . In particular, if the regions υ i ∈ ϒ consist of up to k compact, convex domains in R d , and any intersection υ i 1 ∩ · · · ∩ υ i t also satisfies this property, then h(ϒ) ≤ k(d + 1) 49,106,107 .
. . , f n ) of a simplicial complex is the list of numbers of its k-dimensional simplexes, f k = #{σ i ∈ �| dim(σ ) = k} ("f" is a traditional notation that should not be confused with the firing rates). The h-vector of is defined as where parentheses denote the binomial coefficients. Given the combinatorial decomposition of l, where l k ≥ l k−1 ≥ . . . ≥ l j ≥ j ≥ 1 108 , define the set of numbers with 0 (k) = 0 . Eckhoff 's conjecture 59 , proven in 60 holds that the h-numbers of a d-representable complex must satisfy the following inequalities: which can be verified not only for "static" complexes, but also for the "inflating" (2) and "flickering" complexes, at each step of their evolution. • Qualitative spatial consistency. It can be shown that if the RCC5 relationships among all triples of regions are consistent, then the entire schema R 5 is consistent [93][94][95][96][97][98][99] . The full set of consistent triples is given in the following table.

Definition 1 G(I)
is an interval graph, if each vertex v i ∈ G(I) corresponds to an interval I i ∈ I and a pair of vertexes (v i , v j ) is connected by an edge iff I i and I j intersect.  www.nature.com/scientificreports/ Definition 2 A comparability graph G ⊳ represents an abstract relationship " ⊳ ", if its vertexes v i represent elements of a set, and each link e ij represents a ⊳-related pair, v i ⊳ v j .
An interval graph is hence 1-dimensional skeleton of the nerve of I (Fig. 9A). It can also be verified that the complement of an interval graph is a directed comparability graph G ⊳ (I) , in which the relationship v i ⊳ v j is defined by the order of the overlapping intervals, Definition 3 A directed graph satisfies -property if there are no three vertices v i , v j , v k such that v i , v k are not adjacent, while v i is adjacent to v j and v j is adjacent to v k with the corresponding orientations being e ij and e jk respectively. Definition 4 A comparability graph satisfies ×-property if no four vertexes v i , v j , v k , v l produce disjoint pairs of intervals.
In other words, a situation when v i ⊳ v j and v k ⊳ v l (i.e., the pair of intervals (I i , I j ) overlaps and the pair (I k , I l ) also overlaps), while the remaining pairs remain incompatible, e.g., v j ⋪ v k , v j ⋪ v l , (i.e., I j does not overlap either I k or I l ), etc., does not appear.
Theorem A graph is an interval graph iff its complement is a comparability graph with an order defined by (9), satisfying the ×-property.
This theorem and the definitions motivate the following algorithm for identifying 1D representability of a complex (Fig. 9B): 1. Test whether is a clique complex, i.e., verify whether all (k + 1)−tuples of vertexes v σ = [v i 0 , v i 1 , . . . , v i k ] form a simplex in if and only if each pair of vertexes [v i p , v i q ] ∈ v σ is an edge in its 1-skeleton G = sk 1 (�) . If at least one v σ fails this test, then is not a clique complex and hence not representable. 2. Build the complement G of sk 1 (�) and verify its comparability as follows: i. Choose an edge between v i and v j and define an orientation on it (e.g., e ij = e ji ). If e ij was selected, then search for all vertexes v j ′ that are connected to v j but not to v i (Fig. 9A). If the edge between j and j ′ is not yet oriented, select e j ′ j . If it was already (j ′ j)-oriented, continue on; the opposite, (jj ′ ) -orientation implies that is not representable. If the orientation for new edges cannot be selected, pause the algorithm and dispose of all the edges that have already been oriented. Then pick another unoriented edge and restart the -rule: keep applying it until the process comes to a halt and the next set of edges needs to be removed. Do this until all the edges are serviced and hence oriented.  ii. Verify that no 3-tuple of vertexes (v i , v j , v k ) forms an oriented 3-cycle. If such a cycle exists, is non-representable in 1D. iii. Verify that no triple of vertexes (v i , v j , v k ) is "disconnected, " i.e., given e ij and e jk , there must exist an edge between i and k. If any triple violates this condition, is not representable. Otherwise G = sk 1 (�) is a comparability graph with the order: v i ⊳ v j for each e ij .