Ensemble synchronization in the reassembly of Hydra's nervous system
Creators: name: Lovas, Jonathan
name: Yuste, Rafael
Publisher: name: Zenodo
Year: 2021

description: Although much is known about how the structure of the nervous system develops, it is still unclear how its functional modularity arises. A dream experiment would be to observe the entire development of a nervous system, correlating the emergence of functional units with their associated behaviors. This is possible in the cnidarian <i>Hydra vulgaris</i>, which, after its complete dissociation into individual cells, can reassemble itself back together into a normal animal. We used calcium imaging to monitor the complete neuronal activity of dissociated <i>Hydra </i>as they re-aggregated over several days. Initially uncoordinated neuronal activity became synchronized into coactive neuronal ensembles. These local modules then synchronized with others, building larger functional ensembles that eventually extended throughout the entire reaggregate, generating neuronal rhythms similar to those of intact animals. Global synchronization was not due to neurite outgrowth but to strengthening of functional connections between ensembles. We conclude that <i>Hydra's</i> nervous system achieves its functional reassembly through the hierarchical modularity of neuronal ensembles.
description: Lovas, J., Yuste, R. (2021) Ensemble synchronization in the reassembly of <i>Hydra's</i> nervous system. <i>Current Biology</i> <b>Data</b> <b>Aggregates – </b> Aggregate Raw Data – Raw .tiff files of aggregates at 8, 24, 48, and 72 hours used for analysis. Analysis – Results of analysis with randomization used to generate figures for the paper Variables of interest: -'spikes' or 'spikesclean': Raster plots of detected spikes (time x neurons). -'normcospikeweights': Weighted adjacency matrix of neuron coactivity (neurons x neurons). -'fsclean': Data structure with (x,y) positions (rows 1 and 2), raw track intensity (row 4), detected spikes (row 5), and count of frames of inactivity (row 6), for every neuron in the experiment (columns) by frame (track field number). -'A': Preceding a number, binary adjacency matrix at particular significance threshold. -'M': Community membership at plateau of maximum modularity. -'Msweep': Commumity membership at gamma values from 0.6:0.01:3. -'Msuball': Community membership of Louvain subdivision of M. -'avgxy': Average (x,y) position of neuron across experiment. Aggregate AVIs - .avi files of each individual experiment Tracking – ICY tracking of experiments used in analysis (Excel Output 2 Files). Includes manual linking of disjoined tracks of the same neuron. <b>Dissociated Cells – </b> DissociatedCellsRawData – Raw .tiff files of experiments with dissociated cells DissociatedCellsAVIs - .avi files of each individual experiment DissociatedCellsTracking – ICY tracking of experiments used in analysis (Excel Output 2). Includes manual linking of disjoined tracks of the same neuron. <b>Intact Animals – </b> Analysis – Annotation of CBs and RPs from intact animals and aggregates. IntactRawData – Raw .tiff files of imaging of intact animals. IntactTracking – Tracking of the experiments used to generate fluorescence traces (Fig. 2). <b>Neurite Outgrowth – </b> 8hr – Images of analyzed neurites at 8hrs in neuronal GFP dilution experiments. 24hr – Images of analyzed neurites at 24hrs in neuronal GFP dilution experiments. 48hr – Images of analyzed neurites at 48hrs in neuronal GFP dilution experiments. 72hr – Images of analyzed neurites at 72hrs in neuronal GFP dilution experiments. Structure Analysis – Annotation of neurite structure. <b>2021_Lovas_Code – </b> Steps for reproduction: Preprocess all raw data as described in the manuscript. Track with ICY using parameters described in the manuscript (or use included files). Set paths and run data through 'preprocessX' script. Manually correct errors in tracking (or use included files). Set paths and run data through 'linktracksaftercuration'. Iterate 4 and 5 until an acceptable level of error is reached. Set paths and run 'cleanlinking'. Set paths and run 'generatesignificantnetworks', and 'brainconnectivitytoolboxanalysis'. Set paths and run 'fig3_synchronizationwithdistance', 'fig4_networkanalysis', 'fig4_Bk_Pcloud', 'fig4_Ck_Pcloud', 'figS1_trackingstats', 'figS2_firingstats', 'figS3_nullvalidation', 'figS4_networkanalysis_cut20pct'. Set paths and run 'louvaingammasweep'. Set paths and run 'fig5_louvaingraphs', 'fig5_maxQmodularityandsubclustering', 'fig5_plotcomsandsubcoms', 'figS4_gammasweepplots'. Set paths and run 'neuritetracing' and 'figS1_neuritetracing' (or use included files). Set paths and conservative neurite distance parameters and run 'fig6_avalanches', 'fig6_avalanchesonlouvain', 'fig6_communityrasters_subcoms', 'figS6_spatialoverlap'. Set liberal distance parameters and run 'fig6_avalanches' and 'fig6_avalanchesonlouvain' to generate Figure S6. Annotate contraction bursts and rhythmic potentials in 72hr aggregates and intact animals (or use included files). Set paths and run 'fig2_72vsintact'. Funding provided by: Defense Advanced Research Projects Agency<br>Crossref Funder Registry ID: http://dx.doi.org/10.13039/100000185<br>Award Number: HR0011-17C-0026Funding provided by: National Science Foundation<br>Crossref Funder Registry ID: http://dx.doi.org/10.13039/100000001<br>Award Number: CRCNS 1822550Funding provided by: Vannevar Bush Faculty Award*<br>Crossref Funder Registry ID: <br>Award Number: ONR N000142012828Funding provided by: National Institutes of Health<br>Crossref Funder Registry ID: http://dx.doi.org/10.13039/100000002<br>Award Number: 5 T32 GM 8798-18Funding provided by: Vannevar Bush Faculty Award<br>Crossref Funder Registry ID: <br>Award Number: ONR N000142012828
description: <b>STAR Methods</b> <b>RESOURCE AVAILABILITY</b> <b>Lead contact</b> Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Rafael Yuste (rafaelyuste@columbia.edu). <b>Materials Availability</b> This study did not generate new unique reagents. <b>Data and code availability</b> The data and codes generated or analyzed during this study are available on Dryad. <b>EXPERIMENTAL MODELS AND SUBJECT DETAILS</b> <b><i>Hydra </i></b> <i>Hydra vulgaris </i>(Strain AEP) expressing GCaMP6s in the interstitial stem cell lineage under control of an actin promoter were cultured in the dark at 18<sup>o</sup>C in standard media as previously described <sup>4,59,60</sup>. Freshly hatched <i>Artemia</i> nauplii were fed to the animals every other day. <b>METHOD DETAILS</b> <b>Hydra Dissociation/Reaggregation</b> For experiments, the dissociation and reaggregation of <i>Hydra</i> were performed as previously described<sup>6</sup> with minor modifications. In brief, after 48hrs starvation, 40-50 medium-sized <i>Hydra</i> were collected in a 3.5cm dish, washed thoroughly with room-temperature <i>Hydra</i> medium (RTHM) five times, washed quickly with ice-cold deionized water five times, and resuspended in ice-cold Gierer dissociation media (DM) - prepared from frozen concentrate - in a 1.5mL tube with one DM exchange. <i>Hydra</i> were incubated at 4<sup>o</sup>C for 2 hours, with media exchanges every ½ hour after vigorous resuspension of the animals with an unpolished Pasteur pipette. After 2 hours, DM was exchanged again and <i>Hydra</i> were triturated until the media became cloudy, indicating the start of tissue dissociation. After this, DM was exchanged one last time, and <i>Hydra</i> were mechanically dissociated with a new unpolished Pasteur pipette firmly at a constant rate avoiding bubbles ~75-125 times until approximately one-half to one third of the volume of intact tissue remained. Tissue fragments were allowed to sediment 90 seconds and the top ¾ of the cell suspension was gently collected. Cells were counted on a hemocytometer and seeded in 0.5mL microcentrifuge tubes in numbers according to the requirements of the experiment, in this case 9k cells. We note that while most cells were single in suspension, a small percentage of cells (&lt;10%) remained in tiny aggregates of ~2-10 cells. Due to their small size, filtration through standard filters was unable to break these aggregates, assuming the cells were not aggregating after complete dissociation in the time before cell counting. Cells were allowed to sediment on ice for 5 minutes, then centrifuged at 4<sup>o</sup>C for 4 minutes at 300g in an Eppendorf 5430 R centrifuge. From here, pellets in 0.5mL microcentrifuge tubes were placed at 18<sup>o</sup>C and incubated for 2 hours, 3cm above a bed of melting ice to slow the temperature transition. Pellets were transferred to 800µL dissociation media after two hours. At 8 hours after dissociation, media was diluted 50% with 800µL <i>Hydra</i> media at 18<sup>o</sup>C, and again with a 50% HM replacement at 22 hours, 46 hours, and 70 hours. <i>Hydra</i> in all experiments survived the dissociation procedure and developed into normal <i>Hydra</i>. <b>Imaging</b> <i>Hydra</i> were imaged as previously described,<sup>4</sup> using 100μm spacers between coverslips for the imaging sessions. Analyzed images were acquired on an Olympus IX70 inverted fluorescence microscope equipped with a GFP filter set, 10x Plan Apo air objective, and a Hamamatsu C9913-EMCCD camera using Micromanager.<sup>61</sup> This setup was illuminated by a mercury arc lamp. Images were acquired at 100% illumination intensity, at 250ms exposure, 1x binning, for 7200 frames (~42.1 minutes with CCD charge transfer time). Intact animals were imaged on a Leica M165 FC fluorescence dissection microscope as previously described.<sup>4</sup> <b>QUANTIFICATION AND STATISTICAL ANALYSIS</b> <b>Image Pre-Processing</b> Images were pre-processed in FiJi. For tracking, either 2 (Intact Animals; 48, 72 hour aggregates) or 10 (Dissociated Cells; 8, 24 hour aggregates) frames were summed, background was subtracted with a 10 pixel rolling ball radius, and a 3-pixel diameter median filter was applied. After tracking, all 7200 frames of each experiment were similarly background subtracted and median filtered for extraction of continuous track fluorescence, interpolating (x,y) positions across 10 or 2 frames, depending on the time point. <b>Cell Tracking and Neuron Identification</b> Cells were tracked using the open-source cell tracking software ICY,<sup>11,62</sup> using parameters tuned to the unique requirements of experiments across time points (Table 1). Spatial coordinates of each track were imported into MATLAB and used to extract fluorescence intensity from pre-processed images with the maximum values along each axis of a 7x7 pixel square around the detected track. The deformable hydrostatic skeleton of <i>Hydra</i> makes tracking cells over large movements difficult. While this problem is remedied to some extent in reaggregating <i>Hydra</i>, especially at early time points when there is not much movement, inaccuracies in cell tracking still exist (Figure S1). Even under static conditions of fluorescence, cell tracking is an unresolved problem and manual curation is necessary to approach flawless results.<sup>63</sup> To ensure accuracy in tracking across experiments, misaligned or broken tracks were manually corrected, and inaccuracies in neuron identification (e.g. autofluorescent spots passing under tissue) were excluded (Video 1). The final analysis was conducted on frames 571-5700 of the original 7200 frames acquired at 2.85<i>hz</i> for several reasons. The dropped time at the start of the experiments allowed the animals to equilibrate to the intense illumination. In addition, this period of time offset some of the effects of photobleaching, as autofluorescent spots on the aggregate bleached rapidly with intense illumination, and GCaMP began to dim toward the end of the imaging session. Finally, the tracking software ICY did not handle the ends of some tracks well, and cutting time removed some bouts of track activity that weren't picked up at either end of the experiment. <b>Manual Annotation of Animal-Wide Activity</b> The two main forms of neural activity in <i>Hydra</i> – contraction bursts and rhythmic potentials – were manually annotated in 72 hour aggregates and intact animals for analysis. We began analysis on the first contraction burst that followed a period of quiescence and ended analysis on the last complete contraction burst followed by a period of quiescence within the 30 minute period of analysis. Given that a subset of contraction burst neurons sometimes fires before the remaining CB neurons are recruited to the burst, we began our annotation of pulses of contraction burst activity with the first firing of the dominant, early subset of neurons in the burst. We note that the analysis of rhythmic potentials was conducted on the Rhythmic Potential 1 circuitry, as the activity of Rhythmic Potential 2 was very rare, if present at all, in the 30 minute imaging sessions.<sup>4</sup> <b>Spike Extraction</b> To account for the effects of mild photobleaching on fluorescence intensity traces, aggregates were segmented in MATLAB and changes in the overall fluorescence intensity across the entire segmented aggregate were subtracted from fluorescence intensity traces using an X-frame sliding window (Figure S2), after smoothing raw fluorescence traces with a 2<sup>nd</sup>-degree polynomial Savitzky-Golay filter. From here, MATLAB's 'findpeaks' algorithm was used to detect both the maximum of the filtered fluorescence intensity as well as the maximum and minima of the local first derivative of fluorescence intensity of the peaks, using a threshold of standard deviation and wider sliding averages exclude peaks in the signal generated by noise or broad sub-spike-threshold calcium fluctuations. 'Spikes' were defined around the maximum and minimum of the first derivative of fluorescence intensity to account for differences in spiking kinetics of different classes of neurons (Figure S2). From here, two additional thresholds were applied to compensate for the effects of photobleaching on the intensity of tracks throughout the experiment, and to exclude noisy autofluorescent spots from analysis. The final analysis was conducted on tracks that were active at least 1/3 of the 30 minute experiment, firing at a frequency of at least 1/180/p(active)<i>hz</i> (10 spikes). <b>Data Analysis </b> Spike rasters were used to generate adjacency matrices for analysis. To account for the differences in the duration of spiking events of individual neurons at the imaging framerate, correlations between neurons were based on the normalized coactivity (a final value between zero and one, if they never or always fired together, respectively). Studies of functional connectivity did not take the distance between neurons into account. Studies of neuron distribution within communities and waves of neuronal activity that permeate the system took into account the fact that Hydra's neuron's limited number of neurites don't project long distances in most cases.<sup>4,52</sup> To provide a rough proxy for structural connectivity in analysis of propagation of activity through the network, we used conservative estimates: the average neurite length per neuron of 48 hour aggregates (58μm; Figure S5, S1), multiplied by 2 to allow connections between distal neurites (115μm; Figures 6, S1), as well as a more liberal estimate: the average length of the longest neurite per neuron, multiplied by 3, to account for any potential incomplete labeling in neurons of aggregates or dim distal extremities of neurites (185μm; Figure S6, S1). To determine 'significant functional connections' for edges in metrics calling for unweighted, undirected networks, for each pair of nodes in the network we circularly shifted the spike trains at uniformly distributed random intervals 1000 times to reveal a probability distribution for the observed coactivity between the neurons. To avoid error propagation within a fairly liberal null model thresholded at a standard a standard α=0.05, we used the most stringent threshold that would allow analysis of our chosen network properties to limit the chance addition of spuriously significant edges to the network (e.g. significant network α = 0.005 vs. α=0.05). We note that significant coactivity is presumably mediated by chains of physical connections/synapses if this distance was longer than twice our observed neurite lengths. As controls to understand effects of network size and connection density, we compared our network metrics to binary networks randomized by circularly shifting adjacency matrices along both axes 20 times to generate p=0.05 density-preserving null distributions for our metrics (Figure 4a,b,d,e, dashed lines, Table 4a,c,e,g). As a control for the inevitable removal of nodes from the developing network, we also analyzed our models with random subsets of 20% of the nodes removed, approximately the number lost at each successive time point (18.21+/-5.8%; Figure S3, Table 9). To uncover community structure in our network models, we use the Louvain method of modularity maximization. This greedy method is fast and precise, and performs better than other methods when presented with network modules comprised of different numbers of nodes.<sup>64</sup> In addition to other graph theoretic analysis in the present work, Louvain consensus community clustering makes use of the Brain Connectivity Toolbox.<sup>19</sup> In brief, the method heuristically optimizes modularity, a measure of the density of edges inside communities to edges outside communities, by first optimizing modularity locally on all nodes, then grouping these small communities into a node and reiterating.<sup>30</sup> The Louvain method is particularly appropriate for studies of hierarchical modularity, as the iterative method inherently reflects any containment hierarchy in the network structure. To reconcile variability in detected communities due to the random initiation of nodes at the start of each iteration of the algorithm, we use the final consensus communities after 1000 runs of the algorithm, clustering the final consensus matrix of agreement values reflecting similar partitioning between nodes with τ=0 to ensure as many nodes as possible were clustered (Figure S4a).<sup>31</sup> With consensus communities in place, to ensure optimal partitioning with the method we then varied the resolution parameter γ, a modification to the modularity objective function which increasingly penalizes out-of-community connections to tune the density of detected communities. <sup>33-35,65,66</sup> In essence, higher gamma values favor detection of smaller, more densely connected communities. Given the stochastic initiation of the algorithm, to ensure we were looking at a stable partition of the weighted network we used the 'plateau method', sweeping the resolution parameter γ of the modularity objective function Q from 0.6 to 3 (e.g. Figure S4b), using communities from the plateau of identical partitions across at least 3 increasing γ increments of 0.01 that gave the largest value of Q.<sup>32,33</sup> Sub-partitioning of the modules of this initial partition was performed with γ=1. Analysis of trends in the values γ sweep in Figure S4b,c was performed with plateaus of length &gt;=0.02 γ increments to provide more data points for analysis. Waves of activity, or 'neuronal avalanches,' depending on the field, were determined using the metrics of structural connectivity defined above. Power-law plots where created using the MATLAB 'plplot' function.<sup>21</sup> Waves of activity were superimposed on Louvain communities using custom MATLAB code (supplemental material). <b>Statistics </b> We tested the significance of synchronization observed in our experiments against the null hypothesis that the synchronous firing of neurons is governed by a random process.<sup>67</sup> Randomized networks were used to establish confidence intervals to determine significant connectivity in analyses. Due to the variable length of spiking of <i>Hydra</i> neurons, random networks were created by randomly circularly shifting individual neuron tracks around the length of the experiment at uniformly distributed random locations, repeated 1000 times for each experiment. For tests of the significance of network structure, the weights or binary values of the original adjacency matrix were randomized 20 times to generate p=0.05 null distributions. The significance of data plotted in power law form was assessed using the non-parametric Wilcoxon Rank-Sum test of significance. Data presented in bar graphs was assessed using the non-parametric Mann-Whitney-Wilcoxon U-test. Bar graphs were generated using custom MATLAB code, and violin plots were generated using a modification of the 'violinplot' function in MATLAB (Bastian Bechtold, 2016).