Dimensionality reduction and recurrence analysis reveal hidden structures of striatal pathological states

1 Introduction

The basic mechanisms of brain functions like perception, memory, attention, motor programs, emotions, and decision-making are now being studied using a variety of experimental techniques and theoretical frameworks (Rolls, 2016). Different pieces of knowledge are determined from each experimental/theoretical configuration. Realizing what objective data each one produces and putting them all together into a coherent “big picture” are both difficult tasks. Recently developed technologies for numerous simultaneous recordings and the computing power to evaluate them have led to a controversy over approaches that attempt to comprehend multicellular recordings and neuronal populations without sacrificing or omitting single cell resolution. The discovery that brain neurons do not act alone but rather collaborate to form groupings known as neuronal ensembles, which exhibit spatiotemporal coactivation, is significant (Yuste, 2015; Lara-González et al., 2022). When neurons in an ensemble are engaged in spontaneous, stimulated, diseased, or task-related activity, they fire in a coordinated manner (Ikegaya et al., 2004; Harris et al., 2011; Pérez-Ortega et al., 2016; Hamm et al., 2017; Sheng et al., 2019; Siniscalchi et al., 2019). The neural networks with emergent populational features have connections made between neuron groups rather than between individual neurons (Figure 1; Hopfield, 1982; Buzsáki, 2010; Kampa et al., 2011; Semedo et al., 2019; Rossi-Pool et al., 2021).

www.frontiersin.org

Figure 1. An illustration from Hebb that shows how recurrent transitions between neuronal ensembles might be (represented as “systems” in Hebb’s postulate). It proposes the possibility that ensemble C can serve as a link between two ensemble chains, thus participating in both, being the neuronal substrate for associative learning (Hebb, 1949). In addition, a modular theory can be inferred: each ensemble has inputs and outputs and performs a task step. Alternating activity between them performs the complete task or routine. Alternative pathways are possible to link different procedures. Tools described in the present article can show that this idea can be demonstrated with identified neuronal ensembles (cf., Figure 5).

www.frontiersin.org

Figure 2. Analysis of calcium imaging experiments. (A) Calcium imaging movies can be acquired from neuron cultures (top left), brain slices (bottom left) or in vivo (right). With the aid of a molecule that fluoresces by binding calcium flowing into the activated neurons, the objective is to examine neurons that are present in the field of view. (B) An illustration of a video frame taken in a region of the brain where the fluorescent protein GCaMP6f is expressed in the cells. Regions of interest (ROI) are created from image sequences like this one to collect calcium signals. (C) Top: An illustration of calcium signals extracted from ROIs. Graphed as a raster plot at the bottom, neuronal activity inferred from calcium signals where each row represents the activity of a single neuron. Several techniques exist for creating raster plots from neural activity (Theis et al., 2016). (D) Once a raster plot has been created, an analysis pipeline is proposed and further described in the present work.

www.frontiersin.org

Figure 3. Identification of neuronal ensembles from calcium imaging experiments. (A) Raster plot. Each row corresponds to the activity of a neuron during an experiment. Black periods indicate moments of inferred electrical neuronal activity (events) where neurons have a high probability of firing action potentials. Dozens of neurons can be monitored simultaneously for several minutes, allowing the analysis of phenomena that can only be recorded in neuronal populations. (B) Adjacency matrix of neuronal activity (VA) determined with UMAP, neurons that are close together in the high-dimensional space display a functional connection (Fröhlich, 2016). (C) Adjacency matrix sorted after identifying neuronal communities: neurons show a significant number of connections within the same ensemble and significantly fewer connections with neurons from other ensembles. (D) Sorted raster plot after reordering the neurons according to their communities. These communities exhibit the properties attributed to neuronal ensembles (Buzsáki, 2010). Notice clusters of neurons (colored) alternating their activity following temporal sequences. (E) Low-dimension UMAP projection of active neuronal matrices identified in the experiment. Colors denote the same ensembles present in panel (D). Neuronal ensembles occupy different positions in space after dimensional reduction. (F) Circular visualization of neuronal ensembles and their functional connections is obtained from the adjacency matrix B and colored as in panel (D). Neurons present most of their connections with neurons of their own ensemble, and links among ensembles may represent the connectivity that enables the alternating activity during temporal sequences: how an ensemble is turned off when another is turned on. (G) Bidimensional projections in low dimension UMAP space to better appreciate communities’ separation from different angles.

www.frontiersin.org

Figure 4. Determining significant coactivity peaks. (A) Raster plot and histogram of coactivity. Neuronal activity exhibits a unique space-time structure. (B) Plot of a representative surrogate neuronal raster obtained by maintaining the number of neurons, time, and number of firings of each one, but randomly permuting their activity on the x-axis. The coactivity histogram (bottom) shows a different pattern with respect to that observed in experimental data. However, note that the activity graph (at right) remains intact, revealing that the same level of neuronal activity is present. This class of shuffled raster plots the probability of executing a type 1 error as evaluated with a corresponding Runs test. (C) Another class of surrogate raster plot maintains neuronal activity with the same restrictions as in panel (B), but with the same distribution of intervals between active moments after the random permutation. The pattern of coactivity determined through this process is like the one observed in the experiment. This yields a surrogate raster that is hard to distinguish from the experimental one. This surrogate raster plot is used to evaluate type 2 error with a corresponding Runs test. (D) Once shown that the experimental coactivity signal is not a product of chance, the significant values of coactivity are determined with a sliding window equal to “n” standard deviations (regularly n ≥ 2). A dynamic threshold is used to capture significant coactivity peaks. Asterisks indicate the time periods when the coactivity exceeds the threshold value.

www.frontiersin.org

Figure 5. Dynamics of neuronal ensembles in different experimental conditions. (A) Significant activity of neuronal ensembles was identified in a raster plot obtained in the striatum under the control condition in the presence of N-methyl-D-aspartate (NMDA) to pharmacologically evoke the activation of network components; this brain nucleus is commonly very silent when there is no stimulus (spontaneous firing: Lara-González et al., 2019). The coactivation peaks are shown by darker colors, while the coactivation peaks that are not significant are indicated by paler colors for the same ensembles of neurons. (B) A low-dimensional UMAP projection of the significant peaks of coactivity (VP) that identify the significant coactive matrices from the raster plot. The colors denote identified neuronal ensembles, each in a specific niche in the UMAP space, confirming differences in the activation patterns of these neuronal groups (separation with the same color denotes the same neurons with different activation patterns). (C) Directed graph of neuronal ensemble temporal transitions; each node represents an ensemble. Arrows Edges are functional connections between them: color indicates the origin and arrowheads the destiny; the thickness of the edges represents the number of times a transition was carried out in the observed data. Some transitions can be identified in the raster plot even if the sorting procedure (see the Section 3.4) did not show them together. (D) Projections on the combinations of UMAP planes of the column vectors shown in section B. (E–H) Ensemble coactivation patterns in the decorticated striatum. Background neuronal activity is considerably lower than in control conditions and compared with the pathological states below (cf. histograms of % cell activity at right). Although numerous groups can be separated in section F, scarce temporal sequences (arrows in G) form a ring structure with many unidirectional connections. (I–L) Coactivation patterns of neuronal ensembles identified in striatal tissue depleted of dopamine (6-OHDA Parkinsonian model). There are fewer connections (K) between ensembles that are activated at regular times (I). A more frequent transition appears (K; blue to orange). (M–O) Patterns of significant coactivation of neuronal ensembles identified in the striatum under the model of L-DOPA induced dyskinesia (LID): ensembles are projected in the low-dimensional space (N), as well as transitions (M,P), with several recurrences and high background activity (% activity histogram at right; note scale change).

Hebb (1949) proposed the neuronal ensembles theory as cell assemblies, and many research teams have recently refined it to propose them as the fundamental nervous system processing units (Figure 1; Buzsáki, 2010; Carrillo-Reid, 2021; Grillner, 2020, 2021; Lara-González et al., 2022). Although they may be referred to by different names, they commonly share several characteristics (Carrillo-Reid, 2021) and can be identified depending on context, i.e., responding together to stimulus (e.g., sensory inputs), causing an output (e.g., behavior), or representing experimental conditions (e.g., pathological states). It is considered that emergent properties of neural networks are a consequence of the spatiotemporal dynamics of coactive sets of neurons that alternate their activity following recurrent sequences that resemble computational routines or algorithmic processes that, together with a “modular” brain theory and analytical techniques, may generate a large-scale understanding of brain functions (Buzsáki, 2010). Here we expose a pipeline that may serve to ask whether neuronal samples of histological scale (dozens of neurons), commonly used for medical diagnosis, are representative of underlying ensembles of a larger scale by showing instructive and heuristic changes under different experimental conditions.

One multicellular recording technique that makes it possible to observe neuronal ensembles in culture, brain slices, or in vivo brain preparations from different animal models (Carrillo-Reid et al., 2008, 2009, 2016; Ahrens et al., 2013; Lock et al., 2015; Pérez-Ortega et al., 2016, 2021; Serrano-Reyes et al., 2020) is calcium imaging by using fluorescent indicators of neuronal activity, including the expression of genetically engineered calcium sensors in clusters of neurons via viral transfections or from birth using transgenic animals (Lin and Schnitzer, 2016; Wang et al., 2019). The purpose is to obtain videos of neuronal activity with single cell resolution in different experimental contexts (Yang and Yuste, 2017). The proposed pipeline uses the information from a video acquired through a calcium imaging experiment. However, it could be adapted to other multicellular recording techniques. Figure 2 depicts several calcium imaging setups as well as the method for locating regions of interest (neurons) and computing the fluorescence signals indicative of intracellular calcium levels. The final objective is to obtain the spatiotemporal information of the neuronal activity present in the field of view that can go from the histological to the mesoscale level, including thousands of neurons (Stringer et al., 2021). Instead of using “raw” calcium signals, which may show misleading relationships when evaluated as simultaneous recordings, spike timings estimated from calcium signals should be employed (Carrillo-Reid et al., 2008; Pérez-Ortega et al., 2016; Theis et al., 2016; Serrano-Reyes et al., 2020). Once this is done, analytical methods for identifying and studying neuronal ensembles are diverse, and a consensus has not yet been reached (Wenzel and Hamm, 2021). Different methods have been used in our laboratory with consistent and qualitatively similar outcomes, including hierarchical clustering, dimensionality reduction, similarity indices, and correlated activity (e.g., Carrillo-Reid et al., 2008, 2009, 2011; Jáidar et al., 2010; Pérez-Ortega et al., 2016; Serrano-Reyes et al., 2020; Duhne et al., 2021), indicating that, despite the constraints of calcium imaging, ensemble identification and activity sequences may provide robust data. But dimensionality reduction techniques are getting better all the time. Recent applications of the uniform manifold approximation and projection (UMAP) have been used in several contexts (e.g., Becht et al., 2019). Combined with modularity algorithms (Newman, 2006), it may become an optimal technique to identify communities in a population, called neuronal ensembles in the present context. Numerous techniques have identified the Parkinsonian circuit as characterized by exhibiting neuronal hyperactivity in the striatum as compared to the control, accompanied by a highly recurrent ensemble that greatly monopolizes the microcircuit at histological scale: the circuit gets stuck by this ensemble, as a “metaphor” of what happens with patients that cannot move or have difficulties doing it (Jáidar et al., 2010, 2019; Plata et al., 2013; Pérez-Ortega et al., 2016; Lara-González et al., 2019). Perhaps other disease states may also show a “fingerprint” when observed at the ensemble level (Fornito et al., 2016). The proposed pipeline allows us to describe disease by analyzing neuronal populations in ex-vivo brain slices.

2 Materials and equipment

The current work utilizes the database created by earlier laboratory work, which contains all information regarding experimental protocols and experimental subjects (Pérez-Ortega et al., 2016). A total of 37 ex-vivo brain slices from different mice were utilized in n = 12 control, 11 decorticated, 7 parkinsonian, and 7 dyskinetic experiments. The movies were captured at a rate of 4 frames per second. The derivative criteria was applied to identify the periods of neuronal activity. The code in the python language and instructions to follow the analysis pipeline can be found at. Throughout this work, we use the UMAP v.0.5 Python implementation (McInnes et al., 2018), Brain Connectivity Toolbox for Python v.0.5.2 (Rubinov and Sporns, 2010), and the PyRQA tool to perform recurrence analysis in a massively parallel manner using the OpenCL framework (Rawald et al., 2017).

3 Methods 3.1 General considerations

Neuronal activity can be stored in a variety of ways; practical data structures rely on the temporal resolution and recording period (Paninski and Cunningham, 2018). The acquisition rates in calcium imaging experiments that last several minutes enable the construction of a brain activity monitoring matrix. There is previous work comparing actual with inferred neuronal activity from calcium imaging (Pérez-Ortega et al., 2016; Theis et al., 2016; García-Vilchis et al., 2018; Serrano-Reyes et al., 2020; Calderón et al., 2022). This matrix, usually called a raster plot (R), is a binary matrix of [N × F], where the y-axis denotes the number of neurons, N, whose individual activities comprise the rows (row vectors along the x-axis), and the x-axis represents the number of movie frames, F, which shows neurons that activate simultaneously at any given time (column vectors along the y-axis). Thus, R contains all the information necessary to study the population behavior through the simultaneous recordings of several neurons. It serves as the basis for figuring out the sequences of neuronal ensembles that alternate their activity during the experiment.

3.2 Determining the representative graph of neuronal activity using the uniform manifold approximation and projection

Uniform manifold approximation and projection is a dimensionality reduction algorithm that assumes that data samples are uniformly distributed in a topological high dimensional space. It learns the data manifold and then projects it into a lower dimensional space. UMAP accomplishes this goal through two main processes. First, it builds a graph connecting the nearest neighbors of each data point; this is achieved by choosing the distances between the points across the manifold, assuming they are uniformly distributed and connected to at least one other point. The next step for UMAP is to project or map the graph to a lower dimensional space. In this space, it is sought that the distances in the manifold do not vary with respect to the global coordinate system. Once this is achieved, the algorithm can start looking for a good low-dimensional representation by minimizing a cost function (Cross-Entropy) whose goal is to find the optimal weights of the connections. When this is finished, an array of the coordinates of each point in the specified data sample is depicted in a space of lower dimension, keeping the original structure as similar as possible (McInnes et al., 2018). In the raster plot R there are two vectors that may be prone to this dimensionality reduction: those in rows [1 × F] that represent the activity of individual neurons over time (VA) and those in columns [N × 1] that represent the coactive population at each instant of the experiment (VP). To identify and visualize neuronal ensembles and their temporal sequences of activation under different experimental conditions, we combine the UMAP methodology and graph theory algorithms on both the vectors of neuronal activity (VA, Figure 3) and the vectors of coactivation (VP, subsequently Figures 4, 5). Vectors of neuronal activity (VA) help to identify groups of neurons that coactivate with similar spatiotemporal patterns (Pérez-Ortega et al., 2016), that is, neurons that belong to the different neuronal ensembles that can be identified in each experiment (Figure 3A; Yuste, 2015).

The first step is to build a weighted graph, G, where each edge represents the probability that two nodes are functionally connected in our high-dimensional manifold (Fröhlich, 2016; McInnes et al., 2018). To determine connectivity, each VA is considered a sample from a continuous high-dimensional subspace (topology of neural activity). UMAP extends radiuses from each point to connect them by choosing the nth nearest neighbor. The connection of each node with its neighbors allows the local structure of the nodes to be maintained when performing subsequent manipulations (McInnes et al., 2018).

How “near” a particular point is to another is shown by the strength of each connection in the weighted graph. Since each point in this diagram represents a neuron’s vector of activity (VA; Figure 3A), the fact that two points are “near” indicates that their activity patterns are comparable. The result is the adjacency matrix G of size [N × N] that represents the weighted graph of the experiment. The next step is to determine a division of the graph into communities that reveals the existence of neuronal groups with similar activity patterns, that is, neuronal ensembles (Fröhlich, 2016; McInnes et al., 2018). The UMAP parameters are briefly discussed below in order to do this.

3.3 Description of uniform manifold approximation and projection parameters

The UMAP algorithm parameters are used to control the balance between local and global structure in the final projection of the data. The first parameter described is the approximate number of nearest neighbors (N_NEIGHBORS), which is used to construct the initial high-dimensional graph. It restricts the size of the local neighborhood that UMAP will look at when trying to learn the manifold structure of the data. Low values of this parameter will force UMAP to focus more on local structure by restricting the number of neighboring points considered when analyzing high-dimensional data. Small values should be used to capture fine details in the structure of neuronal activity (Becht et al., 2019) because it controls how UMAP balances the global vs. local structure of the data, ensuring fine granularity in building the ensembles. One concern is if this method can “dissolve” the high recurring ensemble discovered during parkinsonism by utilizing other clustering methods. If this occurs, what are the dynamics of the individual components? Do they follow specific types of temporal sequences? Another parameter is the METRIC, which controls how distances in multi-dimensional space are calculated from the input data. UMAP supports a wide variety of metrics, including Euclidean, normalized, angular, those used for binary data, and those based on paired correlations. In agreement with the original ideas of neuronal ensembles (Hebbian correlated activity), we use the metric based on correlations. In this way, the neurons that have functional connections represent those with strongly correlated activity patterns. The next parameter is the minimum distance between points in the low-dimensional space (MIN_DIST), which controls how tight the points are in an identified cluster. Larger values of MIN_DIST pack points more loosely, while smaller values lead to tighter packs. For finer clustering, the algorithm is favored by small values of this parameter (Becht et al., 2019). Essentially, the minimum spread of points can be controlled, thus avoiding scenarios with many points located on top of each other in the lower dimensional embedding. The default value MIN_DIST = 0.1, recommended in the UMAP documentation (see text footnote 2), is used in this paper. The last parameter to consider is N_COMPONENTS, which allows the user to determine the dimensionality of the reduced space where the data will be embedded to allow visualization. The most frequent numbers are 2 or 3, which correspond to the traditional dimensions that are simple to perceive. To generate the weighted graph, N_NEIGHBORS and METRIC are the most critical parameters, whereas MIN_DIST and N_COMPONENTS are crucial for producing the projection onto the low-dimensional UMAP space (McInnes et al., 2018). UMAP is one of the best tools for displaying highly dimensional data, and its success may be attributed to two very significant qualities: greater global structure preservation and more intelligible parameters. A crucial point to bear in mind is that anybody using the pipeline presented here should adjust the parameters (if necessary) in accordance with the characteristics of their data and their goals.

3.4 Identification of neural ensembles by using modularity and graphs to cluster neuronal activity

The relationships between the neurons (nodes) in the raster plot are depicted in the weighted graph G (see Section 3.2). The objective is to identify neuronal ensembles: groups of neurons that coactivate with comparable patterns of activity. Graph theory offers community extraction procedures to locate areas of the graph where groups of neurons (nodes) are highly coupled. The “modularity” algorithm (Newman, 2006) is one mechanism that maximizes the number of connections between elements that are in the same group and minimizes the number of connections between other groups. Modularity positive values indicate the possible existence of a community structure in the graph. Conversely, negative values indicate that a graph cannot be efficiently divided into communities.

The algorithm starts by calculating the modularity matrix B of a graph. Which is defined as:

Bij=Aij-kikj2m

Where the values Aij are the elements of the adjacency matrix G determined previously with the help of UMAP (Figure 3B), ki and kj are the degrees of the nodes (i.e., their number of connections) and m=12∑iki is the total number of connections in the graph. Identifying the biggest eigenvalue and determining its corresponding eigenvector comes after this matrix B has been calculated. In accordance with the sign of the eigenvector’s component elements, the graph is then split in half. The process is repeated for each of the parts using the formulation of the generalized modularity matrix B(g):

Bij(g)=Bij-δij∑k∈gBik

Where δij is the Kronecker δ-symbol and Bij(g) is an array [ngx ng] with elements indexed by the labels i, j of vertices within a group g. At each stage, the contribution to the total modularity △Q is calculated through the following equation:

△Q=12(12∑i,j∈gBij(sisj+1)-∑i,j∈gBij)

Where, for a particular division of the graph into two clusters let si = 1 if vertex i belongs to group 1 and si = −1 if it belongs to group 2. If at any stage of the algorithm, a proposed division is found to make a null or negative contribution to total modularity, the corresponding subgraph remains undivided. Numerous clustering proposals are gathered through this technique. To select a particular version, Bruno et al. (2015) iterate over the various clusterings until the same outcomes are found in most of them, producing a “consensus algorithm.” The algorithm finishes when it reaches this stage. There are a variety of iterations of the original method depending on how the beginning circumstances of the network division are put forth and how the modularity function is maximized (Newman, 2006).

Going back to the pipeline, Figure 3B displays the adjacency matrix G of the original raster (Figure 3A). The identical adjacency matrix G is shown in Figure 3C, but the nodes have been rearranged in accordance with the neuronal ensembles that have been identified and the interconnection that exists between them. Thus, Figure 3D displays a raster plot of the same neural activity as Figure 3A, but neuronal ensembles are shown from bottom to top according to the order in which they first appeared. Each ensemble is then given a color, and the sorted raster plot shows the temporal order. The color scheme utilized in the raster plot is retained in Figure 3E, which shows the projection of the activity vectors VA determined using UMAP in a three-dimensional UMAP space. The ensembles are separated from one another in the UMAP space, and the ensembles’ neurons are close to one another. Another method of visualizing the connections between neurons is to create a functional connectivity graph in a circular representation by arranging the colored-labeled neurons in the Figure 3D in a circle. The neurons of the same ensemble are grouped together in Figure 3F, and the size of the nodes corresponds to their number of connections. To see more clearly Figure 3E, projections of the vectors VA in the various planes are displayed in Figure 3G. The results show that UMAP identifies more ensembles than previous clustering methods (see below).

3.5 Determination of neuronal ensembles’ significant coactivity peaks (VP)

Once neuronal ensembles have been identified, it is necessary to show the transitions between them. As shown in Figure 3D, ensembles initiate, reach a peak, and wane. It is necessary to identify the precise instant at which they take part in a temporal ensemble sequence. Figure 4A top, shows an experimental raster plot obtained under control conditions in the striatum (Pérez-Ortega et al., 2016). At the bottom, a plot of neuronal coactivity is illustrated. It is a time series which consists of the sum of active neurons over time. To capture the occurrences in which a particular ensemble participates in ensemble sequences, periods in which neurons in a raster plot R exhibit significant peaks in coactivation should be identified. Time instants are the frames of the video, △t, each with a picture of the tissue at a certain time tn, which corresponds to the columns of the raster plot R (VP). It must be demonstrated that the significant coactivation peaks observed experimentally were not produced randomly. Here, we demonstrate how the “Runs test” —a non-parametric hypothesis test based on the binomial distribution— may be used to determine if a set of data can be explained by a random process (Bradley, 1968). The alternative hypothesis asserts that these values are not the result of a random process, while the null hypothesis claims that coactivity values randomly grow and decrease in accordance with a binomial distribution (Bujang and Sapri, 2018). Considering the two categories of statistical error: type 1 error (rejecting the null hypothesis when it is true) and type 2 error (failing to recognize significance when it is present). The significance level of the test, which is commonly set at α < 0.05, denotes the likelihood of making a type 1 error (for a discussion of this threshold, see Cowles and Davis, 1982; Glantz, 2012; Martinez and Martinez, 2015).

Runs test is important in determining whether a trial outcome is random for subsequent analysis (Bujang and Sapri, 2018) using the sorted neuronal ensemble subrasters (set of vectors VA belonging to neuronal ensembles) for further analysis. The first step in the Runs test is to count the number of runs in the data sequence. A “run” is defined as a series of consecutive increasing or decreasing values with respect to the mean, where the duration of the run is given by the number of these values. In a random data set, the probability that the (i + 1)-th value is greater than or less than the i-th value follows a binomial distribution. The statistic is:

Z=T-T¯sT

Where T is the observed number of runs, T¯ is the expected number of runs, and sT is the standard deviation of the number of runs. The values of T¯ and sT are calculated as:

T¯=2n1n2n1+n2+1

sT2=2n1n2(2n1n2-n1-n2)(n1+n2)2(n1+n2-1)

Where n1 and n2 denote the number of positive and negative values in the series. The resultant score is compared to the normally distributed, two-tailed confidence interval. It is determined that the alternative hypothesis is correct when the value of the experimental series is higher than that attained by a random series (Mendenhall and Reinmuth, 1982). To apply this method to the data, a dichotomous time series made up of the subrasters coactivity values (VA), where values above the mean are positive and values below the mean are negative, is build. In the case of the Runs test applied to the coactivity of the experimental raster plot of Figure 4A (this test is later applied to Figures 5, 6), a value of Z = 30.05 was determined, which corresponds to p < 0.0001, allowing the null hypothesis to be rejected: the coactivity time series is not a result of chance. Surrogate matrices (Figures 4B,C) must be created while maintaining the same number of neurons, time, acquisition rate, and active frames for each neuron to demonstrate that neither type 1 nor type 2 errors exist. To evaluate the type 1 error, these surrogate matrices are generated by placing in a uniform distribution the instants of activation of each of the neurons during the period of the experiment. Knowing that the null hypothesis should not be rejected, the outcome of test type 1 error should be noted as:

www.frontiersin.org

Figure 6. Recurrence quantification analysis of the firing rate of representative neuronal ensembles. (A–D) Top: recurrence plots of the firing rate of neuronal ensembles of the striatum in the different conditions. At the left and upper ends of the matrices, the traces being analyzed are observed, and for better appreciation, they are also shown at the bottom. Note that, except for matrix B, all conditions exhibit lattice structures, differing in their spacing. (E) The recurrence is higher for the dyskinetic state. (F) The determinism is higher for the parkinsonian condition, suggesting predictability. (G) The divergence is higher for the decorticated condition, suggesting a more chaotic state. (H) The laminarity score is for the parkinsonian circuits, the most regular. (I) The average white vertical line length is higher for the decorticated state, where recurrent times are the longest–almost isolated ensembles. The last row (J–N) shows the comparisons between the samples determined in the complete dataset and essentially confirms what was observed in the representative example quantifications shown in the upper row (E–I). The actual values can be read in the main text.

Ii={10TypeIerrorismadeTypeIerrorisnotmade

The procedure for determining Ii is repeated for M (1,000) surrogates and the probability to make a type I error is calculated as:

α^=1M∑i=1MIi

This value is an estimate of the significance level of the test for a given critical value. Applied with an exemplary surrogate matrix (e.g., Figure 4B), a value of α^<0.05is determined, suggesting there is not a type 1 error.

To estimate type 2 error, surrogate matrices must not satisfy the null hypothesis, which requires that the raster matrices not be random (e.g., Figure 4C). Consequently, surrogate matrices now conserve the parameters used in the type 1 error test adding the interactivity intervals of each neuron (equivalent to interspike intervals in electrophysiological recordings). With these restrictions, a pseudo-population of neuronal activity like the one observed is determined. The hypothesis test is performed with this surrogate raster matrix, and the value Ii is recorded. In the case of a type 2 error, this value is produced as:

Ii={10TypeIIerrorismadeTypeIIerrorisnotmade

The probability of making the type 2 error after M (1,000) surrogates is:

β^=1M∑i=1MIi

In the case described a value of β^ < 0.0001 is determined, suggesting there is not a type 2 error. Be aware that while the formulae are similar, the surrogate matrix type varies. Once it is determined that our experimental coactivity time series is not randomly generated, its upper extreme values are extracted: the “peaks of coactivity” that can be interpreted as marks in time when column vectors (VP) had significant coactivation of neurons during the experiment.

A sliding window is constructed in the coactivity time series and moved one by one to the right until it reaches the end of our data set in order to achieve these values. This window is built using a value that is expressed as a percentage of the data from the time series’ beginning to end. Typical values include 5, 10, or 20%. For each window, the mean and standard deviation are calculated to yield a local threshold. If it is greater than the mean + two standard deviations, then point is considered to be a significant peak of coactivity. Figure 4D shows the result of applying sliding windows of size 20% to the coactivity plot of Figure 4A. An asterisk is placed at each time where the real value of coactivity exceeds the variable threshold. With this procedure, significant peaks of coactivation (VP) are determined from the subrasters (set of VA vectors). Different thresholds can be used with little change in the outcome.

3.6 Characterizing the activity of neuronal ensembles via recurrence analysis

Recurrence analysis was used to conduct a more quantitative assessment of the ensemble neuronal transition dynamics. This approach allows us to quantify the frequency and length of the neuronal ensemble recurrences. Instead of the suggestions that were previously employed (Carrillo-Reid et al., 2009; Pérez-Ortega et al., 2016; Lara-González et al., 2019; Duhne et al., 2021; Calderón et al., 2022), it presents a new potential of interpretations. In light of this, it is possible to model the behavior of neuronal ensembles as a dynamic system. A dynamical system evolves over time within a state space according to a given rule. This evolution is represented in “phase space.” The state of a system in time t is specified by a point in phase space, and the progression of the system generates an orbit or trajectory through this space. Using recurrence plots, it is possible to demonstrate how neuronal ensembles identified using the methods previously described (see above) behave over time in each diseased condition, leaving a “signature” (Argyris et al., 2015). A recurrence plot is a non-linear data analysis technique that allows visualizing those times in which a state of a dynamic system is repeated, revealing all the moments in which the trajectory of the phase space of a dynamic system visits the same area in the phase space. To build a recurrence analysis, let xi→ be the i-th point on the orbit describing a dynamical system in d-dimensional space, for i = 1, N. The recurrence plot P is a matrix of dots in a [N × N] square, where a dot is placed at (i, j) whenever xj→ is sufficiently close to xi→ (Packard et al., 1980; Eckmann et al., 1987).

Pi,j={1:0:    xi→xi→   ≈≈   xj→xj→   i,j=1,… ,N

In this case, xi→xj→indicates the similarity of a pair of vectors. The matrix P captures a total of N2 binary similarity values. A distance measure is needed to determine the similarity between pairs of vectors. There are several alternatives (Manhattan, Euclidean, maximum distance; Webber and Zbilut, 1994). Here, the Euclidean distance is utilized for simplicity. A neighborhood condition is applied to transform the pairwise similarities into binary values. In the fixed radius condition, the binary values are determined by a threshold ∈. All vectors that lie within the ∈-neighborhood of a query vector xq→ are considered like xq→ (Poincaré, 1890). A strategy is to choose an ∈ threshold based on the density of the recurrence plot (Zbilut et al., 2002). A fixed radius value of ∈ = 1.5 was applied to the data for the current proposal. All the parameters described below were determined from objects in the RQAComputation class of the PyRQA tool (library mentioned in Section 2). The values of the fixed radius (neighborhood requirement) and the metric to measure the similarity between the vectors must be supplied in order to generate such an object (for the implementation see the script titled Figure 6 at see text footnote 1). Statistical comparisons used Mann–Whitney tests with Holm–Sidak post-hoc adjustment.

The parameters that recurrence quantification analysis extracts from recurrence matrices P are: (1) The recurrence point density, or recurrence rate (RR), that is defined as:

RR=1N2∑i,j=1NPi,j

Where Pi, j are the entries of the recurrence matrix P. This parameter quantifies the proportion of recurrence points that are determined with a specified radius. When N →, ∞ RR is the probability that a state recurs to its ∈-neighborhood in phase space. (2) Another measure is based on diagonal lines; it is called determinism (DET). It refers to the portion of recurrence points that form diagonal lines. Only diagonal lines with a length of l ≥ 2 are considered regarding the quantitative analysis:

DET=∑l=2NlPD(l)∑i,j=1Npi,j

Where lPD(l) is the number of points of the recurrence matrix that form diagonal lines of size l. Chaotic signals (aperiodic and presenting sensitivity to initial conditions) yield short diagonal lines, periodic or deterministic signals yield long diagonal lines, and stochastic signals do not show diagonal lines (Webber and Zbilut, 1994). DET parameter is a measure of the order, predictability, or rigidity of the system. (3) Divergence (DIV) is the inverse of the length of the longest diagonal line found in the recurrence plot:

DIV=1max({li}i=1Nl)

Where Nl is the total number of diagonal lines. This parameter is conjectured to be related to the Lyapunov exponent, which estimates the rate at which signal paths diverge. Thus, for larger DIV values, a time series is more chaotic (Eckmann et al., 1987), (4) Laminarity (LAM) captures the number of recurring points that form vertical lines. The equation to determine this parameter is identical to the equation to determine DET, with the exception that, in this instance, the length of the vertical lines is measured rather than the length of the diagonal lines. LAM will decrease if the recurrence plot consists of more single recurrence points than vertical structures. This parameter allows us to investigate how a signal evolves over time through the concept of intermittency. A signal is called intermittent if its temporal evolution appears to be regular for long periods of time and is interrupted from time to time by brief irregular intervals with amplitudes of greater intensity. Intermittency, in this regard, indicates a seamless change from regular periodic behavior to chaotic behavior (Marwan et al., 2007; Argyris et al., 2015). (5) The length of continuous vertical lines formed by values with Pi, j = 0 (average “white” vertical line length W). The length of these lines is equal to the time needed by the system to recur to a previously visited state in such a way that it serves as an estimator for recurrence times (Ngamga et al., 2007).

4 Results 4.1 Vectors of neural activity (VA) are used to determine neuronal ensembles

Hebbian theory (Figure 1; as used by previous authors: Eichenbaum and Davies, 1998) served as the foundation for the proposed pipeline (Figure 2) for determining neuronal ensembles. The analysis of raster plots showing inferred neuronal activity from calcium imaging experiments performed on brain slices was done using the differentiation method (Pérez-Ortega et al., 2016).

The first stage in the process (Figure 2D) is to use UMAP to recognize the neuronal ensembles from a raster plot (Figure 3A). This enables the creation of a representative experiment adjacency matrix (Figure 3B). The ideal split of the matrix into communities that theoretically correspond to the neural ensembles is then determined using the addition of the modularity method (Figure 3C). The detected neural ensembles are then displayed in the original raster plot after it has been sorted (Figure 3D; colored). Every ensemble or collection of neurons creates a time series with recurrence. The same colors used in the sorted raster display are used in Figure 3E to project these neuron groups into a low-dimensional environment. The neurons are then arranged in a circle-shaped representation, with each neuron acting as a node and each edge representing a connection between neurons (Figure 3F). Two-dimensional projections are also provided as a different way to visualize the division of neuron groups (Figure 3G).

Each neuronal ensemble is put through the Runs tests after being divided into its own time series to prove that its activation patterns are not the result of chance (Figure 4; see the Section 3.5). Then, using this data, two analytical procedures are carried out: (1) to examine the significant activation patterns of the sequences of each neuronal ensemble (see below), and (2) to find out whether the activity rates of each neuronal ensemble recur (see below).

4.2 The activation pattern of the significant coactivity peaks (VP) characterizes the pathological conditions in the striatum

In order to highlight potential research directions that address concerns regarding neuronal population dynamics, Figure 5’s left side presents four raster plots (Figures 5A,E,I,M) that were obtained in the striatum under different experimental conditions (Pérez-Ortega et al., 2016). The neuronal ensembles in these striatal microcircuits were determined using previously described methods and are rendered in different colors. A more pronounced hue is used to highlight the significant VP vectors of each neuronal ensemble that were identified via the Runs tests (α^ < 0.05; β^ < 0.0001; in all analyzed raster plots; see the Section 3.5 and Figure 4). Note that not all activity within an ensemble belongs to significant peaks of coactivity, but only precise moments with a low probability of appearing at random are considered. The transitions between the significant coactivity peaks, which occur when one ensemble stops being active and another one starts, define the temporal sequences between neuronal ensembles. The coactivity histogram at the bottom of the raster plots displays colored vertical bars that indicate the intervals in which each neuronal ensemble has a significant VP vector. Each colored vertical vector appears recurrently throughout the course of time, and sequences between different vectors can also be viewed repeatedly (e.g., blue-orange-green-red in Figure 5A). The middle column of Figure 5 shows the projection of these significant coactivation vectors in the low-dimensional UMAP space (R3; Figures 5B,F,J,N; and for a better visualization in different planes: Figures 5D,H,L,P). These projections represent the states of the system and behave as attractors of neuronal activity because, by definition, transitions between them produce trajectories that recur once and again. The right column illustrates graphs showing these state transitions (arrows) and trajectories (cf. Figure 1 showing alternative trajectories), where each node represents a neuronal ensemble (not a neuron as in previous work: Pérez-Ortega et al., 2016; Duhne et al., 2021; Calderón et al., 2022; Figures 5C,G,K,O). Each directed arrow increases its thickness in agreement with the number of times a transition occurs between a given pair of ensembles involved. For better rendering and tracking, the edge color is the same as the origin node, and the destination node is indicated by the arrowhead. Accordingly, directionality indicates the presence of temporal transitions between neuronal ensembles and forms a way of representing the flexibility of a system to codify computations, given that these state transitions with several alternative trajectories (Figure 1) have historically been associated with underlying mechanisms of brain functions (Hebb, 1949; Buzsáki, 2010; Carrillo-Reid, 2021; Lara-González et al., 2022). Next, each experimental condition is described.

The control conditions’ raster plot (with NMDA in the bath since the striatum has low spontaneous activity; Carrillo-Reid et al., 2008; Lara-González et al., 2019) displays seven distinct ensembles (Figure 5A), which are nearly fully segregated from one another in low-dimensional space (Figure 5B) and exhibit alternating activity (denoted by arrows in the state transition graph) that forms temporal sequences and several state transitions throughout time (Figure 5C). Analyzing different control networks in different slices (n = 12), number of ensembles was (mean ± SEM): 7.33 ± 0.51 and the number of transitions: 28.17 ± 2.25. Note that hierarchical clustering showed only three ensembles using the same data (Pérez-Ortega et al., 2016), suggesting that ensembles determined with previous clustering methods can be further subdivided with UMAP and modularity clustering. Background activity around the significant peaks of coactivity is seen as more opaque color, and the whole activity density may be approximated by the histogram of % cellular activity at right. Figure 5B shows the projection of the significant peaks of coactivity (VP) in

留言 (0)

沒有登入
gif