Ubiquitous proximity to a critical state for collective neural activity in the CA1 region of freely moving mice
Yi-Ling Chen, Chun-Chung Chen, Yu-Ying Mei, Ning Zhou, Dongchuan Wu, Ting-Kuo Lee
UUbiquitous proximity to a critical state for collectiveneural activity in the CA1 region of freely movingmice
Yi-Ling Chen , Chun-Chung Chen , Yu-Ying Mei , Ning Zhou , Dongchuan Wu , andTing-Kuo Lee Institute of Physics, Academia Sinica, Taipei ,11529, Taiwan Brain Institute, National Tsing Hua University, Hsinchu City,300, Taiwan Institute of Neuroscience, National Yang Ming University, Taipei,112, Taiwan Graduate Institute of Biomedical Science, China Medical University,Taichung, 404, Taiwan iHuman Institute, ShanghaiTech University ,Shanghai, 201210, China Department of Physics, National Sun Yat-sen University, Kaohsiung, 804, Taiwan * [email protected] ABSTRACT
Using miniscope recordings of calcium fluorescence signals in the CA1 region of the hippocampus of mice, we monitor theneural activity of hippocampal regions while the animals are freely moving in an open chamber. Using a data-driven statisticalmodeling approach, the statistical properties of the recorded data are mapped to spin-glass models with pairwise interactions.Considering the parameter space of the model, the observed system is generally near a critical state between two distinctphases. The close proximity to the criticality is found to be robust against different ways of sampling and segmentation ofthe measured data. By independently altering the coupling distribution and the network structure of the statistical model, thenetwork structures are found to be vital to maintain the proximity to the critical state. We further find the observed assignmentof the coupling strengths makes the net coupling at each site more balanced with slight variation, which likely helps themaintenance of the critical state. Network analysis on the connectivity obtained by thresholding the coupling strengths find theconnectivity of the networks to be well described by a random network model. These results are consistent across differentexperiments, sampling and segmentation choices in our analysis.
Introduction
The CA1 region of hippocampus is an area of brain that is important for representing space information of the environment .Notably, place cells, the neurons that respond to a certain location or landmark in the environment, are first discovered in thisregion. Recently developed technologies such as miniscopes allow optical recording of large numbers of neurons in this regionwhile the animal is freely behaving in an experimental environment. Through calcium ( Ca + )-dependent fluorescence, thefiring activities of neurons can be inferred simultaneously. This opens up the possibility for studying population coding in theseregions of the brain.To understand implications of the large volume of data, maximum-entropy modeling has been used to match the statisticalproperties of the observed firing images to that of a pairwise-interaction model of binary spin glass. While only matchingthe first and second order correlation statistics of the spins, such models have been shown to well reproduce the higher-ordercorrelations as well as other properties of the observed data for the neural systems . Generalizing the model system to abroader parameter space, it has also been shown that the original parameter values corresponding to the observed brain imagesare poised near a critical point in the parameter space, consistent with the so-called critical brain hypothesis .To find the best values of model parameters that reproduce the statistical properties of observed data, we use the Boltzmannlearning (BL) method , which amounts to performing gradient descent on the Kullback–Leibler (KL) divergence betweenthe observed and modeled distributions of the binarized system states.For many of the reported cases in statistical modeling of neural dynamics, in the proximity of a critical state seems to be acommon observation for these systems. Here we refer the system to be in the proximity of a critical state simply because thetrue critical state should be considered for an infinitely large system while we only consider 60 to 100 neurons here. Even whenthe system size is finite, many properties are still quite similar to the critical state . Since the recorded neurons are only a verysmall fraction of the total neurons, it is desirable to verify whether the same conclusions would have been reached had thesystems been sampled or segmented differently. Based on the miniscope data from four different experiments, the fitted models1 a r X i v : . [ q - b i o . N C ] F e b orroborate the robustness of the proximity to a critical state under different subsampling and segmentation processes. Therecovered best-fitted parameters of the external fields and the couplings in statistical modeling can reproduce the observed meanand pairwise correlations accurately. Similar distributions of these structural parameters are found under different sampling andsegmentation conditions.Having established the stability of the model parameters, we further analyze the connectivity of the inferred spin networkswith different thresholds on the coupling strength and calculate their network properties, such as, cluster coefficient, averagepath length, and degree distribution. The network turned out to be well-described by a random network model for allexperiments and subsamples when various thresholding criteria of connectivity are applied.The main result is presented in section Results. The structure of the coupling strength parameters and the criticality of thesystems are examined in section Parameters of spin-glass model. In section Analysis of network connectivity. the networkconnectivity is analyzed with different thresholds. It shows the property of a random network model. The discussion is in thesection Discussion. The statistical modeling is described in section Methods. In order to compare the fitted parameters of local field h and coupling J for the four datasets, their distributions are calculatedand shown in Fig. 1A and Fig. 1B, respectively. The detail is described in the Supplementary S1. The mean values and theirstandard deviations are listed in Table 1. Since neurons or the identified regions of interests (ROIs) in these systems in generalfavor silent states, the mean values of the local fields are all negative. The distributions of h are quite broad but mostly havevalues within the same range of -3 to 1 for the four cases. The mean values of J are all very small and it is close to the Gaussiandistribution except the long tails at both ends. Actually, the long tail at the positive end is a bit more extended than the Gaussian.More detailed validation of the parameters is in Supplementary S2 and model predictions in the Supplementary S3. h i D i s t r i b u t i o n A Mouse 1-1SDMouse 2-1SDMouse 3-1SD Mouse 4-1SD J ij D i s t r i b u t i o n B Mouse 1-1SDMouse 2-1SDMouse 3-1SDMouse 4-1SD
Figure 1.
Distributions of best-fitted local field and coupling. A.
The distributions of the best-fitted local field h in the fourmice B. The distributions of the best-fitted coupling J . Table 1.
The means and standard deviations (in parenthesis) of the distributions shown in Fig. 1.dataset local field h coupling J Mouse 1-1SD − . ( ± . ) . ( ± . ) Mouse 2-1SD − . ( ± . ) . ( ± . ) Mouse 3-1SD − . ( ± . ) . ( ± . ) Mouse 4-1SD − . ( ± . ) . ( ± . ) To characterize the thermodynamic properties of the model, we extend the parameter space to include arbitrary T and calculatethe specific heat of the model C v = NT σ ( E ) here σ ( E ) is the variance of energy, as a function of temperature. As commonly seen in earlier studies , the specific heatcurve exhibits a single peak around the original system temperature T = Temperature T Sp e c i f i c h e a t C v Mouse 1-1SDMouse 1-2SDMouse 2-1SDMouse 2-2SDMouse 3-1SDMouse 3-2SDMouse 4-1SDMouse 4-2SD
Figure 2.
Specific-heat curves.
Specific-heat curves for statistical models for the four experiments. The experimental data isat T = .
0. The solid (dashed) curves are obtained with segmentation threshold of 1 and 2 standard deviation (SD) σ (cid:16) L Ca + i (cid:17) for the calcium-dependent fluorescence trace of each ROI i .In the thermodynamic limit where the system has an infinite number of spins , the specific-heat peak can diverge andsignify a critical point of the system. However, as shown in Ref , when there is only a finite number of spins, the specific heatwill have a peak near the critical temperature. The peak becomes broader and peak temperature ( T p ) moves away from the truecritical temperature for systems with fewer and fewer spins. Our result shows T p close to 1, where T = T p of specific heat as the critical temperature.To make sure that the observed critical state is not incidental to our choice of segmentation threshold, it is important torepeat the calculations using 2 times the standard deviation of each calcium trace as the binarization threshold. The resultingspecific-heat curves using both σ (cid:0) [ Ca + ] i (cid:1) (referred to as 1SD from now on) and 2 σ (cid:0) [ Ca + ] i (cid:1) (2SD) as thresholds are shownin Fig. 2. The peak positions of specific heat as listed in Table 2 are all close to T =
1. In addition to the peak temperatureof the specific heat, the critical point of the spin system can also be identified with the maximum of the slope dm / dT ,where the magnetization m is thermodynamic average of the sum of all the spins. The detail about m vs T is discussed in theSupplementary S3 temperature dependence of the magnetization . To compare the critical temperatures obtained from bothcriteria, the temperatures of peak magnetization slope dm / dT are plotted against the temperatures of peak specific heat C v for all calculated cases of various thresholds, subsamples, and mice in Fig. 3. The plot of magnetization slope dm / dT vs.temperature for Mouse 1-1SD is shown in the Supplementary Figure S7. This is the most important result of this work that wefind these critical temperatures within 20% of T = dm / dT and C v are eithergreater or smaller than 1 together. It is well known for this kind of model with a finite number of spins, the peaks of dm / dT and C v will move away from the ideal critical temperature in the same way. While the volume of the acquired data in the experiments may be large, it typically represents only a very small portionof the neural systems. Or, the measurements may not have single-cell resolution. Additionally, the way of segmenting therecordings into discrete states may further discard a significant amount of the information from the measurements. With allthese uncertainties, it requires us to further examine the conclusion of the last section that the data represents a state in theproximity of a critical state. One approach to address this issue is by studying subsamples consists of a subset of ROIs randomlyselected from the original experimental set. Then, each is treated as a new set of data and the calculations are repeated to find .7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7
Peak temperature of dm / dT P e a k t e m p e r a t u r e o f C v Mouse 1-1SDMouse 1-2SDMouse 3-1SDMouse 3-2SDMouse 2-1SDMouse 2-2SDMouse 4-1SDMouse 4-2SDTemperature=1Mouse 1-1SD-10subsamples: 63 ROIsMouse 1-2SD-13subsamples: 63 ROIsMouse 2-1SD-8subsamples: 63ROISMouse 4-1SD-9subsamples: 63ROIsMouse 4-1SD-10subsamples: 79ROIsMouse 4-2SD-10subsamples: 79ROIs
Figure 3.
Peak temperature of dm / dT versus Peak temperature of C v . Critical temperature determined by the peak ofspecific heat (vertical axis) versus that determined by the peak of magnetization slope (horizontal axis) for the four mice underdifferent ways of segmentation and subsampling.
Table 2.
Peak temperature of dm / dT versus peak temperature of C v for four different experiments. Using differentthresholds for segmentation, and some random subsamples of given size from the measured data. For multiple subsamples ofthe same size and from the same data, the mean values of the peak temperatures are shown followed by their standarddeviations in parenthesis. dataset peak of dm / dT peak of C v ( T ) Mouse 1-1SD 0.94 0.96Mouse 1-1SD-10subsamples: 63 ROIs 0 . ( ± . ) . ( ± . ) Mouse 1-2SD 1.03 1.04Mouse 1-2SD-13subsamples: 63 ROIs 1 . ( ± . ) . ( ± . ) Mouse 2-1SD 0.83 0.88Mouse 2-1SD-8subsamples: 63 ROIs 0 . ( ± . ) . ( ± . ) Mouse 2-2SD 1.02 1.03Mouse 3-1SD 1.10 1.03Mouse 3-2SD 1.15 1.09Mouse 4-1SD 0.81 0.84Mouse 4-1SD-9subsamples: 63 ROIs 0 . ( ± . ) . ( ± . ) Mouse 4-1SD-10subsamples: 79 ROIs 0 . ( ± . ) . ( ± . ) Mouse 4-2SD 0.96 0.98Mouse 4-2SD-10subsamples: 79 ROIs 1 ( ± . ) ( ± . ) the best fitted model. Here we consider random subsets of sizes N s = 64, 32, 16, and 8. Each size of the subsamples is repeated16 times to estimate the standard deviation for the specific heat and the peak position as shown in Fig. 4.The resulting peak positions of the specific heat show a weak dependence on the system size within 10% down to N s = for a system of nearest neighbor Ising model, the peak becomes broader and peaktemperature further shifts away as the number of spins is reduced. Here we only show the subsamples of one dataset, Mouse 1with 1SD. Similar result from the Mouse 4 with 1 and 2 SD is shown in the Supplementary Figure S8. The coupling parameters of the spin-glass model obtained from statistical modeling constitute an effective or functionalnetwork structure that can reproduce the observed state properties. One can expect that such network structure may bear somesignificance pertinent to the functional dynamics of the brain. However, before such possibility is pursued, we need also tomake sure the structure is quite generic with little dependence on segmentation and sampling. For all of results reported above,the coupling strength J has much larger influence than the local field h on the statistical properties of the model as long as h is .6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Temperature T Sp e c i f i c h e a t C v original, T peak = 0.96subsamples N s , peak position T peak
64, 0.96 ± 0.0132, 0.95 ± 0.0216, 0.93 ± 0.048, 0.90 ± 0.05
Figure 4.
Specific heat curve of original Mouse 1-1SD and subsamples.
Specific heat as a function of temperature forrandom subsamples of measured ROIs in Mouse 1-1SD experiment. The shaded regions represent the spread of the results for16 samples for a given size of subsamples.within the range shown in Fig. 1A. Hence, only coupling J will be considered in the analysis below. To find out whether the distributions of coupling strength J obtained in Fig. 1B is sensitive to the number of neurons detected,subsets of neurons were randomly selected to form a number of subsamples. Then the model parameters of each subsamplewere calculated by the BL method. In Fig. 5, the coupling-strength distribution of the subsystems is compared to that from thefull dataset and quantified using the Kolmogorov–Smirnov (KS) test (see the Supplementary S4 Kolmogorov–Smirnov test forthe definition of KS statistic) . Similar with the temperature dependence of specific heat, the KS statistic of subsamples showsweak size dependence and remains quite small until the subsample size reaches below 30 as shown in the inset of Fig. 5.The robustness of the coupling strength distribution under the subsampling gives us an encouraging hope that, in ourapproach, it may not be necessary to measure all the neurons but only a relatively small subset of neurons to obtain a similarresult.Besides overall distribution, the value of the coupling strength between any two spins in a subsystem is compared tothe coupling strength between the same spins in the original full system. Since, in principle, two very different datasets areseparately fitted, the Pearson correlation coefficient between them should be very small. However, Fig. 6 shows the coefficientsfor the different subsamples all above 0.6. It is important to note that the subsystems have a very strong correlation with theoriginal much larger system with Pearson coefficients between 0.75 to 0.9 for a system with only 40% neurons left. Detailsabout calculations of Pearson correlation coefficients can be found in the Supplementary S4 Pearson correlation coefficients.Then, the next question is whether having a very similar pair-coupling structure is what allows the subsample systems to remainin the proximity of the critical state as shown in Figs. 3 and 4. To answer the above question, it is important to first understand the unique characteristics of the pair-coupling structure. Thustwo alterations to the coupling structure of the model are considered. Firstly, we reshape the coupling strength distributioninto a Gaussian with the same mean and variance and randomly redraw the values of coupling strength while maintaining therank of strength order for all spin pairs. The resulting specific heat-curves shown in Fig. 7 show a broad variation of peaktemperatures. Unlike the original data shown in dashed line, most subsamples have the peak temperatures of specific heat faraway from T =
1, which is the state related to the recorded data.In the second case, the pair-coupling strength J i j is randomly shuffled. This manipulation destroys the network structurewhile preserving the distribution of coupling strength. All the resulting specific-heat curves as shown in Fig. 8 have their peaksshifted to the right with the average value ¯ T p ≈ .
05. This implies that at T =
1, systems represented by the shuffled J are in thelow temperature phase of their respective models. The two studies in this subsection reveal a couple of important characteristicsof the pair-coupling structure. Randomly assigned coupling strength between pairs, in general, will not lead to a pair-interactionmodel with the peak of specific heat sitting near T =
1. The specific pair-coupling structure obtained from the experiment mustbe maintained to keep the state in the proximity of a critical state. .5 1.0 0.5 0.0 0.5 1.0 1.5 2.0
Coupling strength D i s t r i b u t i o n full systemsubsystems643216 8
20 40 60 80
Subsample size K S s t a t . o f J i , j w/ orig.w/ norm.orig. w/ norm. Figure 5.
Distributions of coupling strength of original Mouse 1-1SD and subsamples.
Distributions of couplingparameter J of the statistical models for the subsampled systems (drawn lines) compared with that for the original system(dashed line). The shaded area is the standard deviation estimated using 16 randomly drawn subsamples for each subsamplesize. The inset shows the Kolmogorov–Smirnov (KS) statistics of the subsamples with original (orig.) and a Gaussian (norm.)of the same mean and variance. The KS statistic between the distribution for the original 79 ROIs of the Mouse 1-1SD systemand a Gaussian distribution of the same mean and variance is marked with a square.
10 20 30 40 50 60
Subsample size P e a r s o n c o rr e l a t i o n c o e ff i c i e n t h i J i , j J i , j in full system J i , j i n s u b s a m p l e Figure 6.
Pearson correlation coefficients of the best-fitted local field and coupling between the original Mouse 1-1SDand subsamples.
Pearson correlation coefficients of the parameters J i j and h i between the models of subsamples and themodel for the original system of Mouse 1-1SD. The coefficients remain above 0.6 down to subsample size 8. The insets showcomparisons of J i j values between the subsamples (vertical axis) and the original (horizontal axis). .5 1.0 1.5 2.0 2.5 3.0 Temperature T Sp e c i f i c h e a t C v Reshape distribution of J i , j to Gaussian original Mouse 1-1SDGaussianized J i , j distribution Figure 7.
Specific-heat curves of the original Mouse 1-1SD and Gaussians.
Specific-heat curves of models with newcoupling strength values randomly drawn from a Gaussian distribution of the same mean and variance as the couplingdistribution of the original model. The rank of coupling-strength order of the spin pairs are preserved when these new strengthvalues are assigned.
Temperature T Sp e c i f i c h e a t C v original Mouse 1-1SDwith randomly shuffled J i , j Figure 8.
Specific-heat curves of original Mouse 1-1SD and shuffled network structures.
Specific-heat curves ofmodels with randomly shuffled network structures from the original statistical model of Mouse 1-1SD. The manipulation keepsthe set of coupling strength values and reassigns them to different spin-pairs.
Net coupling j J i , j R O I d i s t r i b u t i o n observedshuffled Figure 9.
Distribution of net couplings for the ROIs of Mouse 1-1D.
Distribution of net couplings J net i ≡ ∑ j J i j for theROIs in the model of Mouse 1-1SD dataset (shaded area), compared with the average of net ROI coupling distribution ofshuffled networks from the same coupling J i j distribution. Further investigation of the resulting network reveals that the observed network is generally more balanced compared to arandomly shuffled model. That is, if we define the net coupling of an ROI as the sum of the coupling strengths of all linksconnecting to this ROI, J net i = ∑ j J i j , the distribution of net couplings for the observed network is much narrower compared to atypical randomly shuffled samples as shown in Fig. 9.The narrower distribution of net ROI coupling can be observed for all our four experiments. (See the SupplementaryS4 Distributions of net coupling strength of ROIs.) To see how the balance of the coupling strength distribution affects thecriticality of a network, we repeatedly apply selective random swapping of the link strengths of a shuffled network in orderto make its net-ROI-coupling distribution narrower. We find the peak of its specific-heat curve does move towards T = . ± .
28 to1 . ± .
21 when randomly shuffled model is rebalanced to as good as the observation, and to 1 . ± .
25 at the limit of thebalancing update when the net ROI coupling distribution is twice as narrow as the observation. Again, it is not easy to obtain amodel with the peak temperature at T = In the section discussed above, the network structure obtained from the mice data shows several very interesting characteristicsthat are important to the results presented so far. But what kind of network structure is it? Here we will try to answer thisquestion. First, let us define some quantities that are commonly used in classifying network structure. For a connected network G consists of the set of vertices (nodes or ROIs) V , the average local clustering coefficients of the network is given by CC ( G ) = | V | ∑ v ∈ V CC ( v ) , (1)where the local clustering coefficient CC ( v ) of a vertex v is defined asCC ( v ) ≡ N v k v ( k v − ) , (2)where k v is the degree or the number of edges/bonds of v , and N v is the number of connected vertex pairs whose both verticesare connected to v . The path length d ( v , v ) is the minimal number of links that connect the vertices v and v to each other.The average path length L ( G ) of a graph G is the average path lengths over all vertex pairs in the graph .In the model, the coupling J i j has both positive and negative values, hence it is separated into excitatory and inhibitorynetwork with positive and negative J , respectively. To examine the structure more carefully, the properties of excitatory network log (Degree ( k )) l o g ( V e r t e x c o un t ) Mouse 3-1SD: X cl = 0.35Mouse 3-1SD: X cu = 0.29Mouse 4-1SD: X cl = 0.37Mouse 4-1SD: X cu = 0.36 B (79, 0.055) B (79, 0.06) B (96, 0.045) B (96, 0.04) Figure 10.
The log-log graph of degree distributions.
The log-log plot of degree distributions at the critical threshold X cl in the excitatory networks and the critical threshold X cu in the inhibitory networks with their fitted binomial distributions.will be calculated with a lower threshold X l for J i j , so that only bonds or edges with J i j > X l will be considered. For larger X l ,there are few edges or bonds between nodes and few nodes are connected. But as X l decreases, more connections are formingand all the nodes may become connected. Thus the degree distribution, clustering coefficient and path length will all depend onthe lower threshold value X l . There is a critical value X cl that for X l > X cl , the network is disconnected such that not all thenodes are connected in one maximum cluster and we have isolated clusters or nodes. Similarly, properties of inhibitory networkchanges with the upper threshold X u , when only bonds with coupling J i j less than X u are considered. There is also a criticalthreshold X cu that for X u < X cu , the network is disconnected.In a typical random graph model, the degrees follow a binomial distribution B ( k , p ) , where k is the degree and p is theprobability of having an edge. One of the earliest random graph models is the Erd˝os–Rényi model , where the graphs ischaracterized by a parameter p for the probability of finding an edge between a given pair of vertices and there is no correlationbetween any edges. For networks of a large size N , a sharp transition occurs at the threshold p = p c = log N / N : The networkswill almost surely be connected for p > p c and almost surely be disconnected for p < p c .In Fig. 10, the degree distributions at the critical thresholds for all excitatory X cl and inhibitory X cu cases are plotted andcompared with the corresponding binomial distribution.The fraction of bonds between all possible pairs of vertecis in the network is a function of X l ( X u ) for the excitatory(inhibitory) network and it is considered the same as probability p in the Erd˝os–Rényi model. At X cl and X cu , the fraction ofbonds p c is the probability p at the critical thresholds. As seen in Fig. 11, the values of p c for all the excitatory and inhibitorycases of all our experiments are close to the predicted values log N / N from the random network model. The deviation isprobably due to the small size of our networks. While criticality have been observed in neural signals of many studies of brains, it is generally not very clear how universalthese phenomena of critical states are. In the process of mapping the brain dynamic to simple statistical model, many detailsand much information is discarded. Furthermore, even with the ever improving technology in experimental studies of brains,the observed signals generally can only represent a small part of the entire system. In the current study, we verify the robustnessof the proximity to a critical state in the CA1 recordings from four free-moving mice. We use different thresholds and differentways of subsampling the available ROIs into different sizes. The critical temperatures as determined by the peak of specific heatas well as the peak of magnetization slope are summarized in Fig. 3 for all considered cases of thresholding and subsampling ofthe data from the four mice. Majority of the peak temperature values are within 10% of T = .04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 log ( N )/ N C r i t i c a l r a t i o o f b o n d s p c Mouse 1-1SD: X cl = 0.34Mouse 1-1SD: X cu = 0.3Mouse 2-1SD: X cl = 0.47Mouse 2-1SD: X cu = 0.45Mouse 2-2SD: X cl = 0.54Mouse 2-2SD: X cu = 0.55Mouse 3-1SD: X cl = 0.35Mouse 3-1SD: X cu = 0.29Mouse 4-1SD: X cl = 0.37Mouse 4-1SD: X cu = 0.36Mouse 4-2SD: X cl = 0.45Mouse 4-2SD: X cu = 0.41Mouse 1-1SD-10subs-63ROIs: X cl = 0.29Mouse 1-1SD-10subs-63ROIs: X cu = 0.28 Figure 11.
The prediction log ( N ) / N versus the critical fraction of bonds p c . The critical fraction of bonds p c when thenetworks start to break apart from one connected component versus the prediction log ( N ) / N from the random network modelwith N being the system size.subsample sizes down to 30 (see Fig. 5). For subsample size smaller than 30, the KS statistic of a subsample increases quickly.This suggests that a minimal number of ROIs of 30 in experimental measurement may already have the main properties of themodel represented. To understand the significance of the properties of the structure of the coupling network uncovered by thestatistical modeling, the best-fit model is perturbed by independently changing the shape of the distribution and its networkstructure (see Figs. 7 and 8). While this proximity to a critical state seems to be very robust against variations in measurement,random shuffling of the bond coupling strength or assigning a value according to a Gaussian distribution will move the systemfar away from the critical state of the model.We further find that the observed networks are generally more balanced than random shuffles of the same networks. Whilethat can not fully account for the observed proximity to a critical state, making a shuffled network more balanced does help tobring it closer to a criticality. On the other hand, from the application of different thresholding criteria, the connectivity ofthe networks appears well described by the random network model. This suggests that the connectivity may not be specific atthe level of average path length or degree distribution. Thus, the factors contributing to the ubiquity of the critical state maybe more subtle than one would expect. A possible explanation is that the integration of sensory information through layersof neurons that precede the CA1 can effectively act as a renormalization process and filter out the irrelevant interactions orcorrelations leaving the most critical degrees of freedoms to be represented by the activities of the observed neurons at CA1.Facing with the increasing volume of biological or social data, statistical modeling is currently the only quantitative andnontrivial approach that can be generally applied without any domain knowledge. Our results show that the application ofstatistical modeling to in vivo recording of brains can yield model properties that are robust to the arbitrariness and randomnessof sampling and segmentation. Specifically, the network structure of the model is found to be crucial to the critical state. Itwould be interesting to find the basic principles required to construct such a network. It is also interesting to find biologicalconditions for the neurons to be at different critical states or non-critical at all. The male C57BL/6J mice (aged 8–12 weeks) were purchased from National Laboratory Animal Center Inc. (Taipei, Taiwan) andmaintained under control conditions as follows: environmental temperature between 21–25 ◦ C, a relative humidity 60 ± µ L pENN. AAV. CamKII. GCaMP6f.WPRE. SV40 (Addgene, ≥ . × vg/mL) was injected into unilateral dorsal CA1 of thehippocampus (coordination: caudal from bregma -1.94 mm, lateral from bregma -1.25 mm, ventral form dura 1.35 mm) at0.1 µ L/min for 7min in Zoletil-anesthetized mice. Two weeks after virus transduction, a Gradient-Index (GRIN) lens (1.8 m diameter, 4.39 mm length, Edmund Optics Inc.) was implanted into dorsal CA1 (coordination: caudal from bregma − .
94 mm, lateral from bregma − .
25 mm, ventral form dura 1.6 mm). For the relief of pain and excessive innate immuneresponses to the GRIN lens implantation, anesthetized mice were given additional 10 mg/kg carprofen (Sigma-Aldrich), 0.2mg/kg dexamethasone (Sigma-Aldrich) subcutaneously, and 10 mg/kg enrofloxacin (China Chemical & Pharmaceutical Co.,Ltd.) intraperitoneally. After GRIN lens implantation, 0.4 mg/mL enrofloxacin were administrated in the drinking waterfor 7–10 days to prevent bacterial infection. Related surgery information was obtained from the UCLA miniscope website( http://miniscope.org/index.php/Surgery_Protocol ). After 3-week recovery period, we used the integratedminiature endoscope (miniscope V3) to check GCaMP+ regions of interest and then attached the baseplate to the skull of thelightly anesthetized mice with acrylic cement. After baseplate fixation, mice were ready to subsequent behavioral trainings.The 30 × ×
30 cm white square box surrounding with 2 distal visual cues on the contralateral walls was used to examinelocomotor activity of awake-behaving mice. Behavior videos were recorded using Logitech C270 Webcam about 70 cm abovethe arena. After the 3-day habituation, the novel environment-induced hyperactivities were significantly decreased in mice.On the 4th consecutive day, behavior and calcium imaging videos were recorded simultaneously for 5-20 minutes in mice.All behavioral and dynamic calcium imaging data were analyzed by the UCLA Miniscope software based on the constrainednonnegative matrix factorization for microendoscopic data (CNMF-E) (the open-source MATLAB analysis package wasobtained from https://github.com/daharoni/Miniscope_Analysis ). The main data analyzed in the presentstudy included mice positions, the intensity of calcium imaging data (the temporal traces C), and the deconvolution of calciumimaging data (calcium spikes).In the current study, we consider the results from four separate experiments with four different mice. In experiment 1(Mouse 1), 79 ROIs are identified from a 10 min recording. 72 ROI for 10 min. in experiment 2 (Mouse 2), 63 for 20 min.(Mouse 3) in experiment 3 and 96 for 10 min. (Mouse 4) in experiment 4.For the statistical modeling in Methods, we convert the calcium signal trace [ Ca + ] i of each ROI, indexed by i , into binarystates s i = ± θ i = σ (cid:0) [ Ca + ] i (cid:1) , s i = (cid:40) + , [ Ca + ] i > θ i − , otherwise . (3)The collection of image frames of the CA1 recording becomes an ensemble of states for the system of N =
79, 72, 63 and 96spins in the four experiments respectively. Details about transform experimental recordings to binary data can be found in theSupplementary S1 Transform experimental recordings to binary data .
To model the distribution of states for the experimental ensembles, we use an all-to-all pair-interaction spin-glass model similarto earlier works by, e.g., Schneidmanet et al., 2006 . The energy or Hamiltonian of the spin-glass model is given by E ( s ) = − ∑ i h i s i − ∑ i < j J i j s i s j , (4)where spin at site i , s i = ± h i is the local magnetic field at site i , and J i j is the coupling strength between two spins i and j .The second summation is over all N ( N − ) / s = { s i } is given by the Boltzmann distribution P ( s ) = e − β E ( s ) / Z , (5)where the normalization factor Z is the partition function Z = ∑ s e − β E ( s ) , (6)and β = T − gives the inverse temperature of the system. We set T = β = h i and J i j in the model (4) are obtained using a Boltzmann learning algorithm, whichminimizes the KL divergence between the observed and model state distributions, D KL = ∑ s P data ( s ) ln P data ( s ) P model ( s ) . (7) he gradients of the KL divergence in the h i and J i j parameter space are given by ∂ D KL ∂ h i = m data i − m model i , (8) ∂ D KL ∂ J i j = C data i j − C model i j , (9)where m i ≡ (cid:104) s i (cid:105) is the magnetization of spin i and C i j ≡ (cid:10) s i s j (cid:11) is the correlation between spins i and j . The angle bracketsindicate averages over observed state configurations of the experiment system and over (Markov-Chain Monte Carlo, MCMC)simulation-generated configurations for the model system.Various optimization methods based on gradient descent can be applied to find the minimum where the gradients (8) and(9) vanish, leading to m data i = m model i and C data i j = C model i j . However for sizable systems, the convergence of the traditionalBoltzmann learning can be slow or unstable with a fixed learning rate. Also for each step of the BL iteration, an extensiveMCMC run can also be required for reaching desirable accuracies of the magnetization and correlation estimates of the model.We take a multi-staged Boltzmann learning approach in this paper that uses shorter MC runs for earlier BL iterations to savetime at the cost of reduced accuracy and uses longer MC runs at the later stages of the BL to guarantee the precision of theconvergence. Details are given in the the Supplementary S1 The multi-staged Boltzmann learning. Acknowledgments
We are grateful for the funding support by the Ministry of Science and Technology of Taiwan (MOST) to TKL and CCCunder the grant no. 108-2321-B-010-009-MY2, to YLC under grand no. 109-2123-M-001-001, and to DCW under grand no.107-2320-B-039-061-MY3. Also, DCW is supported by the National Health Research Institutes under grant no. NHRI-EX110-10815NI and in part by the China Medical University Hospital under grant no. DMR-108-102.
Supplementary information4 Supplementary Discussion S1: Statistical modelling and multi-staged Boltzmann learn-ing
We analyze data extracted from a genetically-encoded calcium (Ca + ) indicator, GCaNP6f and observe the dynamics of calciumin local CA1 regions of mice that are awake and freely moving in an open box. Unless otherwise mentioned below, we presentdetails of our analysis on the Mouse 1 mouse data while similar results have also been obtained from the other three data sets.For our analysis, the calcium signal is first converted to binary data frames or states using the standard deviation of each ROI(neuron) as a threshold. Supplementary Fig. 12A shows the results for ROI 1, with the blue line representing the raw signalintensity of calcium imaging data at ROI 1 from a video of 17880 images frames (about 10 minutes with 30 frames per second).The green stars show the binary states [-1, 1] scaled up by 100 intensity units representing the [active, silent] states, whichare obtained by applying the threshold of σ (cid:0) [ Ca + ] (cid:1) ≈ . s i of all ROIs (neurons) for all data frames with black or white marking the active s i = + s i = − i is respectively above or below its threshold θ i = σ (cid:0) [ Ca + ] i (cid:1) . We assign a vector of spinvariables s = ( s , . . . , s N ) with s i = i th active ROI and s i = − i th silent ROI in 79 ROIs. In traditional Boltzmann learning (BL) algorithm, the model statistics m model i and C model i j in the gradients (Eqs. (8) and (9) in themain text) are calculated using Markov-chain Monte Carlo (MCMC) method with a fixed large number of steps to ensure theirprecision. An iteration in the BL includes the updates h i ← h i + η (cid:0) m data i − m model i (cid:1) (10) J i j ← J i j + η (cid:0) C data i j − C model i j (cid:1) (11)of the parameters of the model, which is repeated until the model statistics match the data with desirable precision. To speedup this process, instead of using a fixed large number of steps in MCMC throughout all iterations of BL, we use a variablenumber of MCMC steps, adaptive on the size of the difference between the model and the data statistics. Specifically, we usethe root-mean-square difference d rms ≡ C rms + m rms , with Time (17880 frames) I n t e n s i t y A ROI [1]binary version of 100 times [-1,1]threshold=1SD*ROI [1]=63.2
Time (frames) R O I i n d e x [ i ] B Figure 12. A . The calcium intensity of raw data of ROI[1] at 17880 time frames with the threshold= 1SD*ROI[1] at 63.3 (redline) and its binary states (-1,1) times 100 as green points. B . The binary version of 79 ROIs use a threshold to each ROI1SD*ROI[ i ], i = , , , . . . ,
78 in 17880 time frames (10 minutes).
Table 3.
Variable number of MCMC steps in each BL iteration, adaptive to the RMS error in BL d rms Number of MCMC steps T . ≤ d rms . ≤ d rms < .
05 1000000 . ≤ d rms < .
009 10000000 . ≤ d rms < .
005 20000000 . ≤ d rms < .
004 30000000 . ≤ d rms < .
003 4000000 d rms < .
002 5000000 m rms ≡ (cid:115) N ∑ i (cid:0) m data i − m model i (cid:1) (12) C rms ≡ (cid:115) N ( N − ) ∑ i < j (cid:16) C data i j − C model i j (cid:17) (13)to define the distance between the model and the data.We start with h i = J i j = d rms , the MCMC runsfor a variable number of steps T as listed in Supplementary Table 3 and the parameters h i and J i j are updated with rules ofEqs. (10) and (11) using η = .
01 as the learning rate. After 100000 BL iterations, the set of h i and J i j that produced theminimum error d rms is chosen as the best parameters of our model if this minimum d rms is less or equal to 0 . h i and J i j is used as the initial condition for an extra BL run with learning rate η = .
003 and a fixed number of MCMCsteps T = d rms < . η = .
01 and T = d rms . After the convergence of the multi-staged BL method, the best-fitted parameters of local field h and coupling J are first groupedinto histograms. Then the histograms are normalized to have the total area of value 1 to become the distributions shownin Supplementary Figs. 14A and 14B for local field h and coupling J , respectively. The histograms have used the bin size Log (Time (sec)) L o g ( r m s C M ) Mouse 1-1SD using a learning rate 0.01
Multi-staged Boltzmann learningBL using MC steps=10000BL using MC steps=100000BL using MC steps=2000000BL using MC steps=7000000
Figure 13.
The log-log plot of the root mean-square difference versus the computer running time (sec) using multi-staged BL(blue color) and other traditional BL algorithm using different Monte Carlo (MC) steps. Traditional BL algorithms seem to beeasier to get trapped in local minima. Our multi-staged MC method seems to be more efficient and also get consistent resultsfor different runs. h i D i s t r i b u t i o n A Mouse 2-1SDMouse 4-1SD J ij D i s t r i b u t i o n B Mouse 2-1SDMouse 4-1SD
Figure 14. A. Distributions of the best-fitted local fields and B. couplings.determined bymaximum value of h i (and J i j )- minimum value of h i (and J i j )number of bins , (14)where the number of bins is 50. The bin size of local fields is 0.11 and couplings is 0.08 in these cases. We use the all-to-all spin-glass model and the multi-staged Boltzmann learning to get the best-fitted parameters by matching theobserved mean (magnetizations) and pairwise correlation with calculated mean and pairwise correlation. The covariance isdefined byCov ( X , Y ) = ∑ Ni ( X i − X )( Y i − Y ) N (15),where X and Y are the means of X and Y , respectively. The correlation coefficient is r = Cov ( X , Y ) σ X σ Y , where Cov ( X , Y ) is thecovariance matrix and σ X and σ Y are the standard deviations of X and Y , respectively. The comparison of observed with thecalculated magnetizations M i and covariance matrix Cov i j are shown in Supplementary Figs. 15A and 15B, respectively. .90 0.85 0.80 0.75 0.70 0.65 0.60 Data M i M o d e l M i A Data
Cov i , j M o d e l C o v i , j B Figure 15. A. Matching the observed magnetizations with the calculated magnetizations. B. Matching the observedcovariance matrix with the calculated covariance matrix using the data of Mouse 1-1SD.
Data C i , j , k M o d e l C i , j , k r = 0.84 Figure 16.
The observed versus calculated triple correlations using the data of Mouse 1-1SD.
For a more detailed comparison of the collective behaviors in data with calculated populations, the triple correlation is calculatedfor a triplet of spins as given by C i jk = (cid:10) ( s i − (cid:104) s i (cid:105) ) (cid:0) s j − (cid:10) s j (cid:11)(cid:1) ( s k − (cid:104) s k (cid:105) ) (cid:11) . (16)The result of comparing the observed with calculated triple correlations in the case of Mouse 1-1SD is shown in SupplementaryFig. 16 with correlation coefficient r = .
84. We included all possible triples of the observed data of 79 ROIs in 17880 framesand compared with calculated states of 79 ROIs in 17880 frames. K simultaneously active neurons The probability of having K out of N ROIs active simultaneously is shown in Supplementary Fig. 17 for the data Mouse 1-1SD.The agreement between model predictions and the data is reasonable good until K=17. Following Meshulam et al., 2017 , theuncertainty of the data is estimated by using random half of the data. With 100 samples the orange shade in SupplementaryFig. 17 represents the one SD.The probability of K out of N ROIs active simultaneously is shown in Supplementary Fig. 17 for the data Mouse 1-1SD.Following Meshulam et al., 2017 , the uncertainty of the data is estimated by using random half of the data. With 100 samplesthe orange shade in Supplementary Supplementary Fig. 17 represents the one SD. The probability of having all ROIs in silentor s i = − P ( ) = .
021 while the data has 0 . ± . K simultaneously active neurons P ( K ) data Mouse 1-1SD, P(0)=0.015(0.001)model, P(0)=0.021 Figure 17.
The probability of K out of 79 neurons in the population are active simultaneously of the model prediction (bluecurve) and the mean of data prediction (red curve) and the error bars of 1SD using random halves of the data from 100 randomtrials (orange area). T d m / d T Mouse 1-1SD, T peak = 0.94Temperature=1 Figure 18.
Magnetization slope dm / dT versus temperature in the case Mouse 1-1SD. The critical point of the spin glass system can be written as the maximum slope dm / dT (see, e.g., Huang, 1988 ), where themagnetization m is the sum of all the spins, which is equal to the observed mean activity of 79 ROIs in the case of Mouse1-1SD. Supplementary Fig. 18 shows the temperature dependence of magnetization slope dm / dT . The maximum peak valueis situated at T p = .
94 which is very close to T = T = The peak temperatures of the specific heat for subsamples of Mouse 4 data was referred to in Supplementary SupplementaryFig. 3 of the main text. Here we show the detail temperature dependence of specific heat in Supplementary SupplementaryFig. 19. Nine and ten subsamples with 79 and 63 ROIs are randomly selected from the original data with 96 ROIs. Subsamplesstill show their proximity to the critical state. The other three datasets also have similar result.
In Sec. V of the main text, the coupling distributions of the subsamples and original data were compared with the Kolmogorov–Smirnov (KS) test. Here we provide a little bit more information about the KS test and its statistic. But readers are referred to, .6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
Temperature T Sp e c i f i c h e a t C v Mouse 4-1(2)SD T peak = 0.83 ± 0.041SD-10subsamples: 79ROIs, T peak = 0.8 ± 0.061SD, T peak = 0.812SD-10subsamples: 79ROIs, T peak = 1 ± 0.052SD, T peak = 0.96Temperature=1 Figure 19.
Specific heat curve from the statistical models for random subsamples of measured ROIs in Mouse 4-1(2)SDexperiment. The shaded regions show the error bars of one standard deviation of uncertainty.e.g., Sheskin, 2020 for more details. The Kolmogorov–Smirnov statistic between two samples of variable x is defined as D n , m = sup x | F , n ( x ) − F , m ( x ) | , (17)where F , n and F , m are the cumulative distribution functions of the two samples of size n and m , respectively, and sup is thesupremum function. It quantifies the difference between the two distributions. When the measurements of each sample aredrawn independently from an underlying probability distribution, the KS statistic can be related to the probability that the twosamples are drawn from the same distribution. In the current study, we are using the KS statistic merely as a measure of howclose two distributions are to each other. Pearson correlation coefficients between subsamples and original sample are are presented in Fig. 6 of the main text. Herewe provide a little more detail of the calculations. Consider a subsample { s (cid:48) α } , α = , . . . , N S − { s i } , i = , . . . , N − s (cid:48) α = s i α and the condition i α (cid:54) = i β if α (cid:54) = β . The coupling strengths J (cid:48) αβ of the statisticalmodel for the subsample are obtained through fitting the mean and covariance of the size- N S spin-glass system to the statisticsof the subset of ROIs using the BL method. Since some spins of the original system are not included, the coupling strengthbetween two subsample spins α and β is generally different from the coupling between the spins i α and i β in the model of theoriginal size- N system. To find how they are related, we calculate the Pearson correlation coefficient between J (cid:48) αβ and J i α i β over all pairs α , β of the subsamples , ρ = ∑ (cid:104) α , β (cid:105) (cid:16) J (cid:48) αβ − ¯ J (cid:48) (cid:17) (cid:16) J i α i β − ¯ J (cid:17)(cid:114) ∑ (cid:104) α , β (cid:105) (cid:16) J (cid:48) αβ − ¯ J (cid:48) (cid:17) (cid:114) ∑ (cid:104) α , β (cid:105) (cid:16) J i α i β − ¯ J (cid:17) where ¯ J (cid:48) ≡ ∑ (cid:104) α , β (cid:105) J (cid:48) αβ / N S ( N S − ) , and ¯ J ≡ ∑ (cid:104) α , β (cid:105) J i α i β / N S ( N S − ) are the mean coupling strength of subsample bonds inthe subsample and in the original system. For each subsample size N S , the calculation is repeated for 64 random subsamplesand the mean and variance of the Pearson correlation coefficients of their coupling strengths with the original system are usedto plot the line and the shaded area in Fig. 6 of the main paper. The same calculations are also performed for the local field h i with results shown in Fig. 6 as well. The net coupling strength of ROIs in Mouse 1 data set was shown in Fig. 9 of the main text. Below we include the results forthe other three datasets in Supplementary Fig. 20. j J i , j D e n s i t y d i s t r i b u t i o n Mouse 1-1SD
OriginalShuffled 4 2 0 2 4 6Net ROI coupling j J i , j D e n s i t y d i s t r i b u t i o n Mouse 3-1SD
OriginalShuffled4 2 0 2 4 6 8Net ROI coupling j J i , j D e n s i t y d i s t r i b u t i o n Mouse 2-1SD
OriginalShuffled 7.5 5.0 2.5 0.0 2.5 5.0 7.5 10.0Net ROI coupling j J i , j D e n s i t y d i s t r i b u t i o n Mouse 4-1SD
OriginalShuffled
Figure 20.
Net-ROI-coupling-strength distributions of the statistical model for all the four datasets
Author contributions statement
YYM, NZ and DW carried out the mouse experiment and wrote Methods: Experimental setup and data processing. YLC, CCCand TKL wrote the other part of the manuscript. All authors reviewed the manuscript. TKL conceived the idea and guided theproject. YLC, and CCC performed the computations and all the figures in the main text and supplementary information.
Competing interests
The authors declare no competing interests.
References O’Keefe, J. Place units in the hippocampus of the freely moving rat.
Exp. Neurol. , 78–109, DOI: 10.1016/0014-4886(76)90055-8 (1976). Ghosh, K. K. et al.
Miniaturized integration of a fluorescence microscope.
Nat. Methods , 871–878, DOI: 10.1038/nmeth.1694 (2011). Schneidman, E., Berry, M. J., Segev, R. & Bialek, W. Weak pairwise correlations imply strongly correlated network statesin a neural population.
Nature , 1007–1012, DOI: 10.1038/nature04701 (2006). Tkacik, G., Schneidman, E., Berry II, M. J. & Bialek, W. Ising models for networks of real neurons. arXiv:q-bio/0611072 (2006). ArXiv: q-bio/0611072. Meshulam, L., Gauthier, J. L., Brody, C. D., Tank, D. W. & Bialek, W. Collective Behavior of Place and Non-placeNeurons in the Hippocampal Network.
Neuron , 1178–1191.e4, DOI: 10.1016/j.neuron.2017.10.027 (2017). Tkaˇcik, G. et al.
Thermodynamics and signatures of criticality in a network of neurons.
Proc. Natl. Acad. Sci. ,11508–11513, DOI: 10.1073/pnas.1514188112 (2015). Mora, T. & Bialek, W. Are Biological Systems Poised at Criticality?
J. Stat. Phys. , 268–302, DOI: 10.1007/s10955-011-0229-4 (2011). Usher, M., Stemmler, M. & Olami, Z. Dynamic Pattern Formation Leads to $\frac{1}{f}$ Noise in Neural Populations.
Phys. Rev. Lett. , 326–329, DOI: 10.1103/PhysRevLett.74.326 (1995). . Beggs, J. M. The criticality hypothesis: how local cortical networks might optimize information processing.
Philos.Transactions Royal Soc. A: Math. Phys. Eng. Sci. , 329–343, DOI: 10.1098/rsta.2007.2092 (2008).
Ackley, D. H., Hinton, G. E. & Sejnowski, T. J. A Learning Algorithm for Boltzmann Machines.
Cogn. Sci. , 147–169,DOI: 10.1207/s15516709cog0901_7 (1985). Kullback, S. & Leibler, R. A. On Information and Sufficiency.
Annals Math. Stat. , 79–86, DOI: 10.1214/aoms/1177729694 (1951). Stanley, H. E.
Introduction to Phase Transitions and Critical Phenomena (Oxford University Press, 1987).
Erd˝os, P. & Rényi, A. On random graphs I.
Publ. Math. Debrecen , 290–297 (1959). Landau, D. P. Finite-size behavior of the simple-cubic Ising lattice.
Phys. Rev. B , 255–262, DOI: 10.1103/PhysRevB.14.255 (1976). Huang, K.
Statistical Mechanics (John Wiley and Sons (WIE), 1988), 2nd edn.
Wang, X. & Chen, G. Complex networks: Small-world, scale-free and beyond.
Circuits Syst. Mag. IEEE , 6 – 20, DOI:10.1109/MCAS.2003.1228503 (2003). Chen, Q. & Shi, D. The modeling of scale-free networks.
Phys. A: Stat. Mech. its Appl. , 240 – 248, DOI:https://doi.org/10.1016/j.physa.2003.12.014 (2004).
Feller, W.
An Introduction to Probability Theory and Its Applications, Vol. 1 (John Wiley and Sons, Inc., 1968), 3rd edn.
Sheskin, D. J.
Handbook of Parametric and Nonparametric Statistical Procedures, Fifth Edition (CRC Press, 2020).(CRC Press, 2020).