A comprehensive functional analysis of the chemosensory system

To systematically study how the compact chemosensory system of C. elegans worms encodes various chemical stimuli, we imaged activity from virtually all of the chemosensory neurons. For this, we used a transgenic strain expressing the genetically encoded calcium indicator GCaMP in all amphid sensory neurons (Fig. 1A, B). Individual animals were inserted into a custom-made microfluidic device [26], and neuronal activity was measured in response to diverse olfactory and gustatory stimuli, representing both attractive and repulsive agents: isoamyl-alcohol (IAA), diacetyl (DA), sodium chloride (NaCl), hyperosmotic (1 M) glycerol (Gly), quinine (Quin), and sodium dodecyl sulfate (SDS), where IAA, DA, and NaCl are attractive cues, while glycerol, quinine, and SDS are repellents. For all conditions, we assayed neural activity for both the presentation (ON step) and the removal (OFF step) of the stimulus (Fig. 1C). To verify that the delivery, or the removal, of the stimulus was temporally accurate, we added a fluorescent dye (rhodamine) to the stimulus. We therefore also assayed neural responses to the buffer supplemented with rhodamine only (control group, Fig. 1C). In general, rhodamine alone elicited minimal responses and all our statistical analyses account for these background-level responses (see the “Methods” section).

Fig. 1

Functional dynamics of the C. elegans chemosensory system in response to a variety of chemical stimuli. A A confocal image of the right side of the amphid organ. Imaging was done using the strain azrIs280 [osm-6::GCaMP3, osm-6::mCherry-NLS] [29, 30]. Red, nuclear mCherry; green, cytoplasmic GCaMP. Neuron identification relies on known anatomical position. B Visualization of the amphid nuclei segmented from A. Fluorescence intensity was measured from the segmented spheres. Right side, purple; left side, blue. C Mean neural dynamics of individual neurons following stimulus presentation and removal. White dashed lines indicate ON/OFF steps. Note that the ASK, ASH, and ADL neurons respond to blue light, hence the activity at the start of the imaging period. Conditions tested: control (n = 7); DA, diacetyl 10−4(n = 23); IAA, isoamyl alcohol 10−4 (n = 18); NaCl, sodium chloride 50 mM (n = 26); glycerol 1 M (n = 7); Quinine 5 mM (n = 11); SDS 0.1% (n = 12). A fluorescent red dye (500 nM rhodamine) was added to the stimuli to verify accurate stimulus switch. The control condition consisted of switching between buffer and buffer + dye. Responses observed in the control condition served as the baseline responses for neurons that may have responded to the dye only. The AWC pair is sorted by activation strength in each worm and is marked AWCs (strong) and AWC.w (weak)

To simultaneously image all of the neurons, we used a confocal system equipped with a fast-resonating scanner that allowed imaging the entire brain volume at 2 Hz (30–40 slices, at a 0.6–0.7 μm Z-resolution), providing the necessary spatiotemporal resolution to reliably extract activity from individual sensory neurons [16, 17, 23, 27,28,29,30]. These acquisition settings, coupled with our analysis pipeline (see the “Methods” section), allowed tracking and measuring activity of all chemosensory neurons from both the right and left lateral sides (22 in total), excluding only AFDL/R which were often below detection levels (Fig. 1A, B).

The populations of responding neurons per each stimulus were generally in line with previous reports (Fig. 1C). For example, the AWC-type neurons responded to the removal of most stimuli in an OFF-step response manner [31]. Similarly, ASH, known as polymodal aversive neurons [32], responded upon encountering noxious stimuli, such as the hyperosmotic solution of 1 M glycerol, and SDS. The ASEL and ASER neurons responded to the addition and removal of NaCl, respectively [18]. Consistent with previous reports [15, 16], we observed a functional hierarchy in the chemosensory network. This hierarchy is reflected in that some neurons are general responders (AWC, AWB), responding to most or all of the tested stimuli, whereas other neurons are more selective to particular stimuli (e.g., ASG, AWA). All of the examined neuron classes responded to at least one of the tested conditions. Together, these findings indicate that individual neurons in this transgenic reporter strain are functionally intact and that our automated analysis system reliably segments and identifies individual target neurons to extract accurate dynamic responses to various stimuli (Fig. 1C).

Lateral symmetric neurons generally show highly correlated activity

Apart from the ASE [18] and the AWC [33] neurons, it is generally assumed that the left and right bilaterally symmetric sensory neurons exhibit similar neural responses [16, 15]. We therefore utilized our ability to simultaneously measure functional responses from both lateral organs and analyzed the correlation between them. For this, we performed a pairwise correlation analysis between all the responding neurons across all six conditions.

Indeed, the lateral right- and left-symmetric neurons showed highly correlated activity dynamics across all conditions (Fig. 2A), and these correlations tended to increase with the response amplitude (Fig. 2B, C). The only exceptions to this were the AWC neurons, which responded asymmetrically in some conditions, and the ASER/L neurons that showed negative or no correlation at all (Fig. 2). The correlation matrices also show that neuron types cluster in a stimulus-specific manner, as clustering varied across the different conditions (Fig. 2A). For example, in response to IAA, the activities of the AWA, AWB, and ASER neurons are correlated with each other and negatively correlated with the activity of the AWC neurons. However, in response to glycerol, activity of AWB and AWC is highly correlated and negatively correlated with the activity of ASER. These results indicate a unique correlation pattern for each condition, providing a “finger-print” of the neuronal representation of a given stimulus. Due to the symmetry in responses, all neuron pairs, aside from AWC and ASE, were grouped for subsequent analyses.

Fig. 2

Activity of the right and left laterally symmetric neurons is highly correlated. Pairwise time-series correlation matrices of the amphid neurons response dynamics. Correlations were first calculated across all neurons of each worm and then averaged over all worms in a condition. Each matrix was sorted using agglomerative hierarchical clustering. Pairs of right and left symmetric neurons are indicated by connecting lines. A Histograms of the correlations of responding left–right neuron pairs. Only pairs with a mean response amplitude above 0.1 were used. B Scatter plots of the mean pair response amplitude vs the pair’s correlation. The ASE and some of the AWC neurons show neutral or negative correlations. C Clustered activity across all conditions reveals the overall lateral functional symmetry. Neuron pairs are sorted by the mean correlation over all conditions. Only the ASER/L pair does not show correlated activity across all conditions

Neural dynamics varies in a stimulus-dependent manner

The strong correlation between the two lateral amphid neurons effectively reduces the number of “coding units” in the system by roughly a half. We therefore asked whether, in addition to the ensemble of responding neurons, stimulus identity could be further signaled by the activation dynamics of specific neurons.

Upon exposure to (or removal of) a stimulus, responding neurons typically show a sharp calcium increase that slowly, over several seconds, decays to baseline levels (Fig. 3, blue). But how stereotypic are these response dynamics? For example, do individual neurons show stereotypic responses regardless of the specific stimulus? Do certain stimuli elicit the same response dynamics in different neurons? To address these questions, we performed a PC analysis on the response traces of all responding neurons across all of the conditions (Fig. 3). The first two principal components combined explain ~ 80% of the variance, and appear to reflect the absolute activity levels before and after presentation of the stimulus (Additional file 1: Fig. S1A-B). However, clustering by the PCs 3–4 (accounting for ~ 10% of the variance) provides a clear separation into three clusters, based on the shape of the response dynamics (Additional file 1: Fig. S2A-B). Most responses (~ 75%) form a single cluster representing the stereotypical response dynamics of a sharp rise in Calcium levels to a narrow peak followed by an exponential-like decrease until resuming baseline levels (Fig. 3, blue). This cluster includes each of the responding neurons in at least one condition, both ON and OFF step responses, and all tested stimuli.

Fig. 3

Activity dynamics varies in a stimulus-dependent manner. PC analysis of neuronal response dynamics. The PCA was performed on individual neuron traces, and each point is the average trace of a single neuron across all worms in a condition projected onto the PC space. K-means clustering revealed the dynamics differences in PCs three and four. K = 3 was chosen based on silhouette scores (Additional file 1: Fig. S2A). The blue cluster represents stereotypical dynamics and includes examples from all neuron classes. The red and green clusters show non-stereotypical dynamics and consist mostly of the ASH, ASI, and AWC neurons in response to specific stimuli. Notably, these three neurons are also represented in the stereotypical-dynamics blue cluster. Inset traces show representative examples of response dynamics for each cluster. Each trace is normalized to its maximal level

Two additional clusters represent variable response dynamics, including sustained elevated activity (AWC in IAA/DA OFF step. Figure 3, red), and inhibitory responses with decreased calcium levels with (or without) an initial peak (ASH in Gly OFF step and AWC in IAA ON step, respectively, Fig. 3, green). These variable responses were observed primarily in three neuron classes (AWC, ASH and ASI), suggesting that some neurons possess a larger repertoire of response dynamics than others, possibly providing more nuance in signaling stimulus identity.

Thus, while chemosensory neurons typically respond with very stereotypic activation dynamics, under some conditions, the same neurons exhibit vastly different dynamics. Such alternative responses suggest that sensory neurons may convey different messages depending on the specific stimulus, effectively increasing the information capacity of the sensory layer.

Inter-neuronal communication shapes the sensory response

Sensory responding neurons can be classified as either primary responders, neurons that sense the stimulus directly and independently (e.g., via a dedicated receptor), or secondary responders, neurons that receive significant input from the network that elicits or shapes their response [34, 35]. Importantly, the same neuron can be primary, secondary, or non-responding depending on the stimulus, its concentration, and background conditions. The recruitment of secondary responders may be facilitated by synaptic neurotransmitters, extra-synaptically via secreted neuromodulators/neuropeptides, or through electrical gap junctions. The chemosensory neurons receive all these input types both laterally from sensory-layer neurons and from other neurons, most of which are interneurons [12,13,14, 36].

As internal communication may influence the response dynamics of individual neurons, we set out to discern the degree to which this inter-neuronal signaling shapes such sensory responses. For this, we measured response dynamics in unc-13 and unc-31 mutant strains that are defective in neurotransmitter (synaptic) and neuropeptide (extrasynaptic) release, respectively (Fig. 4A and Additional file 1: Fig. S3).

Fig. 4

The C. elegans chemosensory system relies on extensive inter-neuronal communication. A Changes in neuronal activities in response to ON and OFF steps of the different stimuli. The neural activities are based on the first 7 s after each step. The respective control response of each neuron was subtracted from the stimulus response (see the “Methods” section). Pink asterisks denote significant WT responses (p < 0.05), and black asterisks denote significant differences between the WT and the mutants (p < 0.05). Both sides of each neuron class were pooled, aside from the AWC and the ASE neurons. p-values were obtained using one or two-tailed t-tests corrected for multiple comparisons using FDR. B Schematic representation of primary (red) and secondary (green) responding neurons for each of the tested stimuli as determined by the dependence of the response on synaptic transmission. The ASJ neuron type, whose response was merely modulated (Additional file 1: Fig. S5), is denoted in cyan

Comparing neural responses in these mutant strains to responses in WT animals reveals extensive inter-neuronal communication (Fig. 4A). Over 40% of the neural responses showed altered activity in at least one of the mutants where they were either completely abolished or significantly diminished (marked in asterisks in Fig. 4A). For example, the AWA and AWB neurons responded to IAA presentation in WT worms (pink asterisk) but failed to respond in the unc-13 mutants (black asterisks), suggesting that these neurons are secondary responders to IAA and that they are recruited to respond via neurotransmitter signaling. Examples of the mean traces and response magnitudes of the various neurons across the different conditions and strains are provided in Additional file 1: Fig S4.

Subtler changes were observed in other neurons whose responses were merely modulated rather than completely abolished (Additional file 1: Fig. S5). Examples include the shift in activity of the ASJ neurons to the ON/OFF step (Additional file 1: Fig. S5A-D) and changes in the maintenance of activity throughout the step in AWC and ASH (Additional file 1: Fig. S5E-H). Moreover, modulated activity of some neurons (e.g., ASJ) is stimulus specific and also depends on neuropeptide signaling (Additional file 1: Fig. S5A-D). Dependence of neural activity on internal network signaling and on the specific stimulus may further increase the coding capacity of the sensory layer neurons.

Overall, out of the subset of responding neurons in each condition, only a few (typically 2–4) were unaffected by inter-neuronal communication and can therefore be classified as primary sensory neurons (Fig. 4B). The secondary responders, forming the rest of the encoding ensemble, are recruited by the primary responders via inter-neuronal signaling, presumably to form the unique nuanced response of the specific stimulus.

Stimulus identity can be predicted by network activity

If stimulus identity were only signaled by primary responders, the sensory system could face a combinatorial problem in that the variety of distinct environmental stimuli far outnumber the possible combinations of primary responders. To estimate how well population coding discriminates between the various stimuli, we trained a random forest classifier on the peak activities of the 13 chemosensory neurons (9 left–right pairs and the individual AWC and ASE neurons) in response to stimulus presentation and removal, for a total of 26 parameters per observation. The classifier perfectly predicted the identity of the stimulus presented to WT worms, and this prediction accuracy decreased the more neurons were removed from the training set (Fig. 5A and Additional file 1: Fig. S6). This suggests that when considering the entire ensemble of chemosensory neurons, it is very easy to discriminate between the different stimuli.

Fig. 5

Neuronal activity predicts stimulus identity. Confusion matrices of a random forest classifier (100 trees and a depth of 4) trained on response dynamics of WT animals. The classifier was applied to test predictions on WT (A, B), unc-13 (C, D), and unc-31 (E, F) neuron activities. A, C, E Training data contained the responses of the entire network. B, D, F Training data contained single neuron responses to both ON and OFF steps of the stimuli. G, H Confusion matrices for classifying by valence (G) and volatility (H). Numbers next to each row show the stimulus-specific F-scores. The performance scale bar is the same for all panels

Next, we analyzed the relative contribution of individual neurons to coding each of the stimuli by training the classifier on the ON and OFF responses of a single neuron. While the classification accuracy using individual neurons was generally low, neurons varied in the degree and type of classification they allowed for (Fig. 5B and Additional file 1: Fig. S7). For example, activity of the olfactory AWA neurons is sufficient to differentiate between the volatile IAA and DA and the other stimuli; however, as evidenced by the frequent mutual misidentification of the two, it is insufficient to distinguish between the two volatiles. Similarly, using the ASH neurons only, the classifier divided the stimuli into three distinct groups (IAA and DA, Gly and SDS, and NaCl and Quin) but had difficulties differentiating between the stimuli within each group. These results indicate that individual neurons contribute in varying degrees to the signaling of certain stimuli but, when combined, provide sufficient information to accurately identify all of the stimuli in our sample.

We next asked whether the internal communication (neurotransmission, neuropeptide release) in the network is crucial for stimulus identification. For this, we used the classifier trained on neuronal activity recorded from WT worms to predict stimuli based on the neuronal activity of the unc-13 and unc-31 mutants. Tested on unc-31 data, the classifier performed nearly as well as on WT data, suggesting that neuro-peptidergic signaling plays a relatively minor role in stimulus identification (Fig. 5C, D). In contrast, the classifier poorly predicted the stimuli in unc-13 mutants, suggesting accurate coding is heavily reliant on neurotransmission (Fig. 5E, F).

Prompted by the coarse separation of the classifier when relying on single neurons, we asked whether certain neurons are tuned towards specific properties of the stimulus. We therefore divided our stimulus sample by volatility (IAA and DA—volatile, NaCl, Gly, Quin, and SDS—nonvolatile) and valence (IAA, DA, and NaCl—attractive, Gly, Quin, and SDS—aversive) and trained the classifier again using the entire set of chemosensory neurons as well as with single neurons. The classifier performed well on both categories when trained on all chemosensory neurons (Fig. 5G, H). Individual neurons varied widely in their ability to classify stimuli by category, where classification of volatility was best achieved using AWA, ASH, and ASE, whereas classification of valence was best when using the ADF, ASH, and AWB neurons (Additional file 1: Fig. S8-9).

Together, our results suggest that individual neurons can encode specific features of a stimulus and that precise stimulus identification is achieved when combining a small number of responding neurons. Moreover, neurotransmitter, rather than neuropeptide, signaling plays a pivotal role in modulating neural responses to allow stimulus discrimination.

Temporal dynamics improves stimulus discrimination

This far, we showed that considering peak neural dynamics sufficed to accurately identify all stimuli in our data (Fig. 5A). However, this could be due to the diverse nature and the small sample size of the stimuli used herein as well as the large amount of neuronal data collected per trial. But could stimulus discrimination be improved by taking into account activity dynamics, in addition to peak activity?

To utilize time series in the classifier, we used CAnonical Time-series CHaracteristics (CATCH-22) to reduce the dimensionality of the response dynamics of each individual step response [37]. CATCH-22 is a set of diverse time-series analysis methods optimized for classification performance with minimal redundancy. The computed features include the mode of Z-scored distribution, time intervals between extreme values, linear and non-linear autocorrelations, and measures of periodicity. The complete list of 22 features is given in Lubba et al. This reduced each step response trace to 22 features that describe the time series, and allowed us to completely separate the amplitude of the response from the time-dependent dynamics. We then performed a principal component analysis on these features and retrained the classifier by adding the first three components of each neuron (together explaining 55% of the variance, Additional file 1: Fig. S10) as variables. Thus, each dataset consisted of a single neuron’s response represented by eight variables—ON and OFF step response magnitudes and three dynamics features obtained by CATCH-22 and PCA per step, for a total of six trace features. We next compared the performance (F1 score from cross-validation) of the classifiers when trained on the response amplitude only, on the three principal components of the trace dynamics only, and the combined amplitudes and trace dynamics for each individual neuron (Fig. 6A–D).

Fig. 6

Temporal dynamics provides additional information for stimulus identification. A Global F1 scores (across all conditions) for each neuron when considering only the trace dynamics, the amplitudes, or both. B Scatter plot depicting the contribution of the amplitudes and the response dynamics to the overall performance of classification by each neuron as shown in A. Amplitudes and dynamics are expressed as a fraction (relative contribution) of their combined accuracy as shown in the third column of A. C Averaged F1 scores across all neurons for each condition when considering only the trace dynamics, the amplitudes, or both. D Scatter plot depicting the contribution of the amplitudes and the dynamics to the overall performance of classification of each stimulus as shown in C. Amplitudes and dynamics are expressed as a fraction (relative contribution) out of their combined accuracy as shown in the third column of C. E Classifier accuracy scores predicting stimulus identity based on dynamics, amplitudes, and both. The data used herein was obtained from [16] consisting of 11 sensory neurons across 23 different stimuli at 10−4, 10−5 and 10−6 concentrations. *p < 10−3 (one sided t-test, FDR corrected). F Classifier accuracy scores predicting stimulus identity based on dynamics, amplitudes, and both when combining all the data from Lin et al. (irrespective of the specific concentration). *p < 10.−4 (one sided t-test, FDR corrected)

For most neurons, considering both aspects of the response had an additive effect, where performance of the classifier trained on both trace dynamics and amplitudes was better than each by itself (Fig. 6A). However, a classifier that was trained on only the dynamics features of either AWCW or ASI performed better than when trained on amplitudes alone, and combining traces and amplitudes of these neurons resulted in minimal additional improvement (Fig. 6A, B). In contrast, for ASJ and ADF, the amplitude provided the most information, and taking trace features into account did not improve performance further.

We also estimated the contribution of amplitudes and trace dynamics to the representation of each stimulus by taking the mean F1 score of each stimulus over all neurons (Fig. 6C, D). For most stimuli, the effect of combining amplitudes and trace features was additive. Interestingly, overall performance of the classifier tended to be better for the aversive stimuli (Quin, SDS, and Gly) than for the attractive ones (IAA, DA, NaCl), suggesting that the sensory system is more finely tuned to precisely identify noxious stimuli (Fig. 6C).

These results suggest that given a sufficiently diverse stimulus space, the identities of the stimuli can be efficiently encoded using response amplitudes alone (Fig. 4A and Fig. 5A). However, response dynamics of individual neurons carry considerable additional information that could be used to help distinguish between more closely related stimuli, particularly aversive ones.

To better understand the relative contribution of response dynamics to stimulus coding, we analyzed available data that measured the activity of all amphid neurons in response to a panel of 23 different odorants, spanning six chemical classes, each in several concentrations [16]. We first extracted peak activities and repeated the analysis with our classifier.

Overall, our classifier performed similarly to the one described in the paper, reaching comparable prediction accuracy of ~ 70% when trained on response amplitudes alone (Fig. 6E). Training the data using the dynamics yielded a lower accuracy of ~ 40% for each individual concentration. However, combining amplitudes with response dynamics significantly improved the performance of the classifier (for each of the dilutions), thus mirroring the single-neuron classification results obtained in our data (Fig. 6A–D).

We then applied the classifier to the entire dataset (23 odorants at three concentrations, for a total of 69 individual stimuli). Due to the similarity of the population coding to different concentrations of the same stimulus [16], it should be particularly challenging to tell them apart. However, even with such a large number of stimuli, the dynamics significantly added to the overall accuracy compared to the performance when considering amplitudes alone (Fig. 6F).

Taken together, these findings indicate that combining response amplitudes with response dynamics significantly improves stimulus identification by generating unique codes to each stimulus (and its concentration), effectively enhancing the coding capacity of the chemical space across a wide range of concentrations.