Unsupervised functional neurocartography of the Aplysia buccal ganglion

A major limitation of large-scale neuronal recordings is the difficulty in locating homologous neurons in different subjects, referred to as the “correspondence” issue. This issue stems, at least in part, from the lack of a unique feature that unequivocally identifies each neuron. One promising approach to this problem is the functional neurocartography framework developed by Frady et al. (2016), in which neurons are identified by a semisupervised machine learning algorithm using a combination of multiple selected features. Here, the framework was adapted to the buccal ganglia of Aplysia. Multiple features were derived from neuronal activity during motor pattern generation, responses to peripheral nerve stimulation, and the spatial properties of each cell. The feature set was optimized based on its potential usefulness in discriminating neurons from each other, and then used to match putatively homologous neurons across subjects with the functional neurocartography software. An alternative matching method based on a cyclic matching algorithm was also developed, which allows for unsupervised extraction of groups of neurons and automated selection of high-quality matches. This improvement enabled unsupervised implementation of the machine learning algorithm, thereby enhancing scalability of the analysis. This study paves the way for investigating the roles of both well-characterized and previously uncharacterized neurons in Aplysia, and helps to adapt the neurocartography framework to other systems.

well-characterized and previously uncharacterized neurons in Aplysia, and helps to adapt 60 the neurocartography framework to other systems.

INTRODUCTION
scalability of the framework to accommodate a larger number of subjects. These 122 adaptations and enhancements will facilitate the future use of neurocartography for 123 analyses of the Aplysia buccal ganglion, as well as other nervous systems.

125
Spatial and functional features 126 In order to build an information-rich feature-space, we first collected data, defined 127 many spatial and functional features of each neuron, and assessed the usefulness of were used to record activity in buccal nerves 1-3 (cBn1, cBn2, cBn3, iBn2, iBn3-where 139 c and i denote whether the nerve was contralateral or ipsilateral to the imaged 140 hemiganglion), the esophageal nerve 2 (iEn2) and the radular nerve (iRn). Activity in 141 these nerves is correlated with specific phases of the animal's feeding behavior (i.e., 142 protraction, retraction, post-retraction; Cropper et al., 2004;Morton and Chiel, 1993b). 143 The relationship between the activity of individual neurons and these phases, in turn, 144 provides potentially useful information about the behavioral role of each cell (e.g., Morton 145 and Chiel, 1993a).   Neurons were classified as bursting or non-bursting using previously described bursting-positive neuron, including the firing rate during bursts (frequency), and the start 174 time and duration of burst-like activity during BMPs (Fig. 2B). Start time and duration were extracted after normalizing the duration of each phase of a motor pattern, whereas 176 frequency was calculated prior to normalization.

177
Neurons within the buccal ganglia may send processes through one or more 178 peripheral nerves (Evans and Cropper, 1998;McManus et al., 2014;Morton and Chiel, 179 1993b; Scott et al., 1991), as well as receive synaptic input from cells whose processes 180 pass through these nerves (Nargeot et al., 1999). Nerve stimulation can therefore activate 181 neurons both directly and indirectly, thereby providing a potential rich source of 182 information on many cells. VSD imaging was used to capture the responses of all neurons 183 to stimulation of each of the recorded nerves (Fig. 3). The amplitude of neuronal activity 184 following nerve stimulation yielded seven nerve response features, one for each of the 185 seven stimulated nerves (cBn1-3, iBn2-3, iEn2, iRn).  window than the 25 ms window for the amplitude parameter was used for improved image quality).

199
Finally, the spatial features for all imaged neurons were parametrized. The  Feature assessment 207 We next evaluated what subset of features would form the optimal feature-space 208 for the purpose of identifying putatively homologous neurons. The quality of a feature 209 depends primarily on three factors: i) how much information it provides; ii) how much of 210 that information is unique relative to other features (non-redundancy); and iii) the extent 211 of its variability within and across subjects (consistency). Thus, a high-quality feature must 212 offer non-redundant information that is sufficiently consistent to allow neurons to be 213 distinguished from one another.

214
The quality of the feature-space was quantified using a metric (termed 215 discriminability index) that reflects the percentage of neurons that are distinguishable in 216 a certain feature-space (Frady et al., 2016). Briefly, this index asks whether, upon splitting 217 the data in two halves, a given neuron is closer to itself than to other similar neurons in 218 the feature-space. If that is the case, the neuron is deemed discriminable, that is, 219 distinguishable from other neurons. The discriminability index was computed for    Conversely, the feature-space for rostral preparations consisted of start time and duration 235 burst properties, in addition to nerve responses for iRn, iBn2, iBn3, cBn2, and cBn3.

236
Responses to iEn2 stimulation were small compared to other nerves, and we could not 237 rule out the possibility that these responses were due to noise. Therefore, we   neurons a3, b3, c3, and d3 form best-match pairs (red circles) a3-b3, a3-c3, a3-d3, b3-c3, b3-d3, and c3-d3. The  One drawback of the pairwise approach is that information across all possible pairs 291 needs to be "stitched" (i.e., integrated in some way to obtain unified insights; Frady et al.,  Thus, supervision effort (i.e., the number of matches the analyst has to judge) scales 301 rapidly with the number of animals.

302
Consequently, a more automatic approach for the stitching process was sought.

303
Here, matches were assessed based solely on the consistency of pairwise matching 304 across animals. A set of four neurons was considered successfully matched across 305 subjects when every possible pairing among them was the best possible match for each 306 neuron in pairwise matching (termed "perfect group"). Four was chosen as the minimum 307 majority of 7 caudal experiments or 6 rostral experiments so that the corresponding 308 neurons cannot appear in two separate sub-groups. Furthermore, when two or more such 309 groups of four neurons shared at least one member neuron, these groups were merged  in ganglia e and f, suggesting a failure to find a reliable counterpart. B2, the reference was another neuron 364 from ganglion a. In this example, there were many candidate neurons in all ganglia, implying that matching 365 is unstable across runs. The same analysis was conducted for every neuron from every experiment, and a 366 pair of neurons is considered as a reliable match when they reciprocally demonstrate highly exclusive matches should be selected as correct. Although such data was provided manually by 379 the analyst in the original framework, the copresence results provide that instruction data 380 in the present study. As the high copresence count implies highly consistent matching 381 results across different runs of cyclic matching, it was considered as a proxy of the 382 reliability and used to label the match as correct or not for this analysis.

383
Selected matches were determined using a copresence proportion criterion, 384 computed over 100 repetitions of cyclic matching. For each match, the copresence 385 proportion was the copresence count of the match (e.g., stems in Fig. 6B) divided by the 386 appearances of the reference neuron (e.g., dashed lines in Fig. 6B). Matches that 387 exceeded a given copresence proportion threshold were selected as "correct" for feature-388 space skewing. We compared various skewed and unskewed feature-spaces by using 389 the sparseness of the full copresence profile of each reference neuron as an overall 390 quality metric. Neurons whose matches are highly random across runs are expected to 391 have low sparseness, whereas neurons that have consistent matches will have sparse 392 copresence profiles. In other words, the sparseness is high when the copresence counts 393 converge to one or just a few neurons, and is low when copresence counts distribute

447
The above process was applied to the cyclic matching results obtained after 448 feature skewing. Nineteen groups from the caudal dataset (Fig. 8A-D), and 22 groups 449 from the rostral dataset (Fig. 8E-F) were identified. For comparison, the cross-animal 450 stitching of pairwise matches (i.e., Fig. 5C-D) was repeated using the feature-skewed 451 datasets, which yielded 35 and 18 groups in caudal and rostral datasets, respectively.

452
To compare cyclic groups and pairwise groups, we counted groups that had 453 identical member neurons, as well as those that were not identical, but still shared four or 454 more members. In the caudal dataset, there were five identical groups and six groups 455 with 4+ shared members (Fig. 8A, C). The same analysis was performed using the 456 dataset without feature skewing and obtained one identical group and seven groups with 457 4+ shared members (Fig. 8D). Of note, no seven-member groups were obtained from 458 cyclic matching with the unskewed data. The same comparison of group members was 459 conducted in the rostral dataset (with feature-skewed data), which resulted in four 460 identical groups and three groups with 4+ shared members (Fig. 8E). In addition, the 461 uniqueness index of the group was calculated, which is defined as the mean uniqueness   neurons ranged from small to large, but most appeared to be on the larger end (Fig. 9A). neurons (e.g., ganglia d-g in Fig. 9A), and share many features they can often be 500 identified due to their unique combinations of features (Fig. 9B). This also appears to be 501 the case for many neurons from the rostral surface (Fig. 10A, B), where most putatively 502 identified cells were small neurons that were responsive to stimulation of many nerves.   The original framework of neurocartography was developed using VSD imaging 533 data from the leech ganglion. Even though the framework itself is highly versatile, many 534 functional features are specific to each system. Therefore, the first step of our analysis 535 was establishing the feature-space for the buccal ganglion that contained the optimal set 536 of features. The automated algorithm selected parameters related to BMPs. This is not  Fig. 6B).

569
The parametrization of individual neuron pairs through the copresence analysis 570 opened a path to the automation of feature-skewing (Fig. 7), allowing for the entire 571 analytic pipeline, including matching of putatively homologous neurons, to be automated 572 and unsupervised. Thus, the combination of cyclic matching with the copresence analysis 573 is a promising scalable approach for functional neurocartography. It should be noted, 574 however, that this approach would not work well when the number of subjects is small 575 because of the limitation in the number of possible random permutations. In principle, the 576 ability to estimate intersubject variability is dependent upon the number of subjects from 577 which the estimation is performed. Thus, it is to be expected that a reduced number of 578 animals should impair performance of any approach that needs to estimate intersubject 579 variability. It is worthwhile to point out that the original manual framework is well suited 580 for cases where the number of subjects is prohibitively small for the automated approach. Finally, the analytical pipeline may be further refined through its application to a 615 dataset for which the ground-truth is known. We had limited means to assess the 616 accuracy of matching results, primarily due to the lack of ground-truth data. Although  The following nerves were recorded: the ipsilateral and contralateral buccal nerves 2 and 657 3 (iBn2, cBn2, iBn3, cBn3); contralateral buccal nerve 1 (cBn1); ipsilateral radula nerve 658 (iRn); and ipsilateral esophageal nerve 2 (iEn2).

659
The responses of the neurons in the imaged hemiganglion to brief electric shocks 660 to the peripheral nerves of the ganglion were characterized. This procedure was designed 661 to activate neurons directly if a particular neuron had an axonal projection in the 662 stimulated nerve, or indirectly via synaptic drive from other neurons, or a combination of 663 the two. This procedure was repeated over 40 trials for each of the nerves (i.e., iBn2, 664 iBn3, iEn2, iRn, cBn1, cBn2, cBn3). Stimulation was automatically triggered by a WPI 665 Pulsemaster A300 and delivered by a WPI 850A stimulus isolator. Stimulus intensity, 666 pulse duration and frequency were respectively 15 V, 0.5 ms, 2 Hz for caudal and 30 V,  The position of a neuron was defined as the x and y coordinates of the centroid of the 745 pixels included in the ROI. The size was defined as the total number of pixels of the ROI.  The feature parameters were also calculated from the half-split data to assess their 789 discriminability (see the next subsection). Data from spontaneous activity recordings were  positive, the neuron was said to be discriminable in this feature-space.

811
The subset of features that was most beneficial for identifying neurons across 812 experiments should maximize the number of discriminable neurons. The total feature-813 space discriminability was given by the following expression: where H(x) represents the Heaviside step function, used to convert positive 816 discriminabilities into ones, which could then be tallied up to represent the total 817 discriminability of feature-space f. Consequently, the feature-space is maximally 818 informative when D is maximized.

819
To determine the optimal subset of features to maximize the feature-space 820 discriminability, at first, a subtractive assessment of discriminability was employed: 821 starting with the full feature set, each feature was removed one-by-one and 822 discriminability was reassessed. If discriminability decreased with the removal, the feature 823 was considered informative and kept in the feature set. Otherwise, the feature was 824 considered not helpful and removed. Next, because how informative a feature is may be 825 dependent on the redundancy of information with other features, the once-removed 826 features were added back one by one, and discriminability was reassessed. If the added 827 feature increased discriminability, it was retained in the feature space.