Intricate response dynamics enhances stimulus discrimination in the resource-limited C. elegans chemosensory system

Abstract

Background

Sensory systems evolved intricate designs to accurately encode perplexing environments. However, this encoding task may become particularly challenging for animals harboring a small number of sensory neurons. Here, we studied how the compact resource-limited chemosensory system of Caenorhabditis elegans uniquely encodes a range of chemical stimuli.

Results

We find that each stimulus is encoded using a small and unique subset of neurons, where only a portion of the encoding neurons sense the stimulus directly, and the rest are recruited via inter-neuronal communication. Furthermore, while most neurons show stereotypical response dynamics, some neurons exhibit versatile dynamics that are either stimulus specific or network-activity dependent. Notably, it is the collective dynamics of all responding neurons which provides valuable information that ultimately enhances stimulus identification, particularly when required to discriminate between closely related stimuli.

Conclusions

Together, these findings demonstrate how a compact and resource-limited chemosensory system can efficiently encode and discriminate a diverse range of chemical stimuli.

Background

Living organisms critically rely on chemical signals. These signals direct fundamental behaviors such as locating food sources and mating partners, or avoiding toxins and predators. Sensory systems therefore evolved to differentiate between the multitude of chemical stimuli to allow animals to form an accurate representation of the environment [1, 2].

In vertebrates, chemosensation is segregated into the olfactory and the gustatory modalities. Each of these modalities relays the information through several neural layers, and this distributed information is then integrated in deeper cortical areas [3]. Within the olfactory system itself, accurate identification of odorants is performed already in the glomeruli and relies on both the population coding and the timing in which each olfactory neuron was activated [4, 5].

Invertebrates, particularly those with small nervous systems, have shallower networks with limited neural layers. As such, chemosensory information is likely to be integrated, even partially, already at the level of the sensory layer [6,7,8,9,10,11]. But how can such sensory systems with limited sensory resources uniquely encode many different stimuli?

For example, the nervous system of Caenorhabditis elegans (C. elegans) nematodes consists of 302 neurons with well-established synaptic connections [12,13,14]. The main chemosensory organ, the amphid, is situated anteriorly and includes 12 bilaterally symmetrical pairs of neurons. Many of these chemosensory neurons respond to a variety of chemical stimuli, including olfactory and gustatory cues [15,16,17,18,19]. Indeed, single cell RNA-seq data indicate that each neuron expresses an array of chemosensory receptors [20, 21].

The C. elegans connectome shows a high degree of inter-connectivity within the sensory layer, suggesting that sensory coding may be shaped by neural communication among the sensory neurons or via feedback from the downstream interneurons [22]. While some of the neurons respond to the stimulus directly (primary neurons), others respond through recruitment via chemical synapses, electrical gap junction, or humorally via neuropeptide signaling. This unique representation poses even greater coding limitations since only a smaller subset, consisting of the primary neurons only, can uniquely encode the target stimulus.

Ample studies quantified neural responses to a myriad of stimuli, focusing mainly on the sensory population coding and the response magnitudes [15,16,17, 19, 23]. These studies revealed that the sensory system employs a hierarchical sparse coding scheme, whereby some neurons respond to a wide range of stimuli, while others are more selective. In addition to population coding, response dynamics was also shown to be carrying important functional information. For example, pulsatile activity in the sensory neuron AWA underlies an efficient chemotactic navigation that allows animals to reach attractive cues faster [24, 25].

In this study, we comprehensively analyzed with cellular resolution how the chemosensory system of C. elegans responds to and codes various chemical stimuli. We reveal that stimuli are encoded using a small, usually bilaterally symmetric, subset of neurons, where primary neurons consist of only 2–4 neurons on average. Interestingly, some neurons possess rich response dynamics that is stimulus specific or network-communication dependent. This rich response dynamics repertoire significantly improves discernment between similar neural population coding, thus endowing limited sensory systems with the capacity to encode a larger amount of information.

Results

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⁻⁴(n = 23); IAA, isoamyl alcohol 10⁻⁴ (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 AWC^s (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⁻⁴, 10⁻⁵ and 10⁻⁶ concentrations. *p < 10⁻³ (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.⁻⁴ (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.

Discussion

We studied how a compact chemosensory system, consisting of limited neural resources, encodes a variety of chemical cues. Using C. elegans as a model system, we found that animals use a small set of neurons to encode each of the stimuli, where ~ 2–4 sensory neurons are the primary sensors of the stimulus, and a few additional sensory neurons are recruited via inter-neuronal signaling. Interestingly, while most neurons show a stereotypic response dynamics of sharp increase in calcium levels followed by a slow decay, some neurons show variable dynamics that depend on either stimulus identity or inter-neuronal communication. These fine response dynamics significantly improve stimulus identification, effectively enhancing the coding capacity of the compact sensory system.

To analyze a wide space of possible sensory responses, we employed a variety of stimuli that represent both olfactory and gustatory cues, some of which are attractive while others are repulsive (Fig. 1). Overall, our results agreed with previous works in regard to response profiles of the sensory system to various chemical stimuli [15,16,17,18, 31, 34, 35]. In addition, most symmetric neuron pairs showed highly correlated activity, with the only exceptions being the AWC and the ASE neuron pairs (Fig. 2). However, it is plausible that testing additional stimuli will reveal more neurons with differential bilateral functionality.

Our results recapitulate a previous report showing that the neurons respond in a hierarchical manner, where some are broadly tuned (e.g., AWC) to respond to all stimuli, whereas others are more finely tuned and respond to specific stimuli only [15]. A similar principle was also observed when analyzing responses to a wide array of olfactory stimuli [16], though some of the broadly and narrowly tuned neurons differed between this and our study. This difference could be due to the nature of the stimuli used in each study. For example, the ASI neurons were broadly activated in our study mostly in response to gustatory stimuli. This may explain why this neuron was not detected as a broad responder when assaying responses to volatile cues [16]. Thus, the segregation to broadly and narrowly tuned neurons may heavily depend on the stimuli used. Nevertheless, in both studies, AWC is consistently identified as a broadly responding neuron.

While we find a small fraction of the chemosensory neurons to respond to each of the stimuli, Lin et al. found that a higher fraction of the sensory neurons is activated when presented with volatiles only [16]. The lower number of responding neurons in our analyses may be because we separated the responses to ON and OFF steps since some of the neurons are classically on-step responders (e.g., AWA) while others are consistently off-step responders (e.g., AWC). Furthermore, we considered the two bilateral AWC and ASE neurons as distinct entities, analysis that further contributed to the sparser coding conclusion. Slight differences in imaging protocols between studies could also result in different individual responses. Conditions such as feeding state, light adaptation, and ambient temperature could all prime the sensory neurons in various ways. However, these differences are unlikely to cause dramatic changes in the fundamental principles of information processing in the network, as evidenced by the replication of such principles across studies [15, 16].

Analyses of mutant strains revealed that the encoding ensemble is determined by the expansion of the responding subset through recruitment of secondary responders, mainly through neurotransmission (Fig. 4, Fig. 5C, D). In contrast, neuropeptides play a minor role in the encoding process (Fig. 5E, F) and primarily modulate the activity of the responding ensemble (Additional file 1: Fig. S5). Therefore, primary/secondary assignments were determined based on the influence of synaptic input on a neuron’s response. Primary neurons were defined as neurons that showed a significant response to the stimulus and were not significantly impacted in the unc-13 mutant. Secondary neurons were those that displayed a significant difference between the WT and mutant responses. Of note, these assignments were based on the step where the main response occurred. For example, AWA’s direct response to DA occurs in the ON step. The change in activity observed in the OFF step reflects the return to baseline following the removal of the stimulus. Therefore, only the ON step was used to define AWA as a primary responder to DA. While many of the effects we describe are robust, here, too, the choice of the precise protocols and statistical analyses, as well as background conditions and experimental variability, might affect the assignments of primary/secondary roles to specific neurons.

An additional intern-neural signaling route involves electrical gap junctions. For example, ASH responses are regulated non-cell-autonomously by gap junctions [38]. C. elegans contains 25 genes coding for innexins [39,40,41], the invertebrate analogs of the vertebrate connexins that make up electrical gap junctions. This substantially complicates analysis of their involvement in recruitment of secondary responders. Nevertheless, if gap junctions also play a role in activating secondary neural responders, then the actual number of primary neurons might be lower than the 2–4 neurons based on our findings (Fig. 4B). This possibility further lowers the already limited coding capacity of the sensory system and underscores the importance of inter-neural communication to fine tune response dynamics that improve stimulus identification.

The sparse representation in the chemosensory system presents a combinatorial problem for the animal, where potentially thousands of different chemical cues are encoded by a small subset of the sensory neurons. This problem is exacerbated by the further division of responding neurons into primary (directly sensing) and secondary (network recruited) neurons, as revealed by the inter-neuronal signaling mutants. This observation bears some resemblance to the “primacy coding” model proposed for mammalian olfactory systems, whereby the identity of a stimulus is encoded by a small set of early-responding glomeruli in a concentration-independent manner [4, 5]. According to this model, coarse stimulus identification occurs first, with fine-tuning of odor identity and concentration being mediated by higher-latency glomeruli. While the relatively slow calcium dynamics does not allow us to resolve response latency, it is an intriguing possibility that a similar “primacy code” occurs in C. elegans, with the broad categorization of the stimulus being determined by primary responding neurons, followed by secondary-responder mediated fine-tuning.

Furthermore, our analyses showed that the relative contribution of the response dynamics to stimulus identification differs between neurons (Fig. 6). It is appealing to consider this feature as another variation of a “primacy code,” where the initial response amplitude provides the coarse stimulus identification, and the longer-scale temporal dynamics of the same neuron assists the secondary responders in its fine-tuning. Indeed, one of the three “slow” and response-rich neurons (Fig. 3), ASI, performed the worst of all responding neurons (excluding ASK and ADL, Fig. 6), for encoding stimuli based on response amplitude, but was markedly better based on dynamics. The two other neurons with variable responses, ASH and AWC, generally display fast initial responses, rendering the dynamics less informative. Notably, these temporal features may be particularly important in C. elegans as they can propagate from the sensory neurons across the network to direct matching behavioral outputs [24, 42].

Activity of secondary responding sensory neurons relies on lateral signaling within the sensory layer or through feedback from downstream interneurons. The importance of such signaling was demonstrated, for example, in the improvement of the signal-to-noise ratio to support a more robust chemotaxis [22, 31]. Recruitment of unique subsets of secondary neurons could also increase coding capacity, potentially alleviating the combinatorial limitations of a small sensory system and the limited number of primary responding neurons. For example, if two different stimuli activate the same primary responding neurons, the secondary recruited neurons may differ. This could be due to the fact that each sensory neuron expresses several chemical receptors [20, 21]. Thus, even if two different stimuli elicit a response in the same sensory neurons, distinct intracellular signaling paths may lead to two different synaptic outputs that will eventually recruit different secondary responders. This in turn generates stimulus-specific codes that allow the animal to discriminate between the two stimuli, despite the fact that they are sensed by the same sensory neurons.

We propose an additional strategy by which the coding capacity could be increased. Our results show that the sensory neurons exhibit a variety of possible calcium dynamics determined by both extrinsic (stimulus) and intrinsic (inter-neuronal signaling) factors. Attempting to predict stimulus identity based on single-neuron responses revealed that these dynamics carry additional information that may help differentiate between stimuli. This may be particularly useful when confronted with a large stimulus space composed of many closely related compounds or different concentrations of the same stimulus. Indeed, it is possible that neurons employ faster, more intricate dynamics than are observable with a relatively slow calcium indicator. C. elegans neurons have long been thought to only employ graded potentials. However, calcium-mediated action potentials have been described in several neurons, including the AWA sensory neurons [43]. In AWA, calcium transients correlate with spike trains suggesting that different activation frequencies may code for different stimuli, thus possibly further increasing the coding capacity.

The encoding ensembles raise several behavioral predictions that would be interesting to test. For example, stimuli with highly similar responding ensembles should be harder to distinguish for WT animals. However, the more differences in the neurons’ roles (i.e., primary responders of one stimulus are secondary responders of the other), the easier the stimuli should be to separate in the absence of network recruitment. For instance, IAA and DA share many of the responding neurons and, depending on concentrations, are similarly attractive to the worms [29, 44]. Based on our results, as the AWA neurons are primary responders for DA, but secondary for IAA, we would expect a reduction of network signaling to AWA to shift the preference towards DA. Similarly, predictions can be made for neurons with known behavioral roles. For example, ASH is the main nociceptive neuron, and its activation causes immediate backward movement [45]. We can therefore expect disruptions in network signaling to ASH to affect aversion behaviors. Our results suggest that the elimination of network input to ASH would alter the behavioral response to Gly, but not to SDS.

Conclusions

Together, this study reveals the principles by which a chemosensory system that is limited in its sensory resources may uniquely encode a large repertoire of chemical stimuli. Similar strategies may have evolved in higher brain systems where a more elaborate sensory system is required to uniquely code and discriminate between a greater space of different stimuli.

Methods

Strains

ZAS280 In[osm-6::GCaMP3, osm-6::ceNLS-mCherry-2xSV40NLS] [29].

ZAS325 is a cross between ZAS280 and unc-31(e928) [30].

ZAS371 is a cross between ZAS280 and unc-13(s69).

Worm cultivation

All worms were grown on NGM plates seeded with OP 50 and kept at 20 °C. Age synchronization was performed by bleaching. All experiments were done 3 days post the bleach, using young adult worms.

Imaging the chemosensory system

Worms were starved for 20 min on un-seeded NGM plates and then inserted into a microfluidic chip [26] where they were partially anesthetized with 10-mM levamisole and left to habituate for 10 more minutes. Recordings lasted for 2.5 min starting with 30 s of light habituation. ON/OFF steps lasted for 1 min each, with stimulus presentations occurring at 30 and stimulus removal at 90 s. Since ASH, ADL, and ASK showed strong responses to light itself, for the ON step responses of these neurons, we used recordings collected from a second round of ON steps at 150 s from the start of imaging.

Stimuli used were diacetyl 10⁻⁴ (DA), isoamyl alcohol 10⁻⁴ (IAA), NaCl 50 mM (NaCl), glycerol 1 M (Gly), SDS 0.1% (SDS), and quinine 5 mM (Quin).

Imaging was performed on a Nikon A1R + confocal laser scanning microscope with a water immersion × 40 (1.15NA) objective at ~ 1.5 volumes/second and with Z-axis intervals of 0.5–0.8 μm. To account for slightly varying acquisition rates, all traces were interpolated to 2 Hz. The system was controlled by the Nikon NIS-elements software.

Neuron identification and signal extraction

Neurons were detected and tracked based on the nuclear mCherry signal using an algorithm developed by [46]. Neuronal identities were determined visually based on their anatomical positions. Neurons that could not be unambiguously identified were removed from analysis. The left–right assignment of the AWC^ON/AWC^OFF neurons is random in each worm [47], and our reporter strain does not allow to differentiate between them. To maintain the functional distinction, the AWC pair in each worm was sorted by activation strength. Where relevant, the neurons are marked AWC^s (strong) and AWC^w (weak).

Calcium traces were extracted from within 0.9 of the radius of the originally segmented neural sphere, in an effort to reduce cross-reads from nearby neurons. In cases where cross-reads persisted, the signal was read from the half of the sphere facing away from the signal donor. Signal intensity was normalized by baseline activity which was defined as the lowest value of a 10-frame running average. All data analyses were performed using in-house MATLAB, Python, and R scripts.

Time-trace correlations

To calculate the correlations in neuron dynamics, pairwise correlations were calculated for all neurons in each worm and then averaged across all worms in a given condition. Overall correlations were determined by averaging all of the worms in all conditions. Each matrix was then sorted using agglomerative hierarchical clustering with Euclidean distance as the distance metric. Due to the strong correlations between bi-laterally symmetrical neurons, in the statistical analysis of all neurons (excluding the AWC and ASE pairs), the left and right sides were pooled together, unless specified otherwise. The ASJ neurons were pooled despite displaying relatively low bilateral correlations because their asymmetry stems from a tendency to respond with stochastic pulses (Supplementary Fig. 5), rather than a bilateral distinction. The AWC neurons within each worm were sorted by overall activity levels and designated AWC^s (strong) and AWC^w (weak).

Calculating response magnitudes

To determine response magnitudes, average baseline activity of 10 s prior to the step was subtracted from the maximum activity in the 7 s following the step. To account for activity caused by changes in the flow direction and the fluorescent dye, the respective control response was subtracted from the stimulus response. Significance of changes in WT amplitudes were tested against zero using one-sample t-test, and mutants were tested against the WT using two-sample t-test. All p-values were adjusted for multiple comparisons using FDR.

For classification, missing data was imputed as described in [16], using a matrix completion algorithm based on minimization of the nuclear norm provided by [48]. (https://github.com/udellgroup/Codes-of-FGSR-for-effecient-low-rank-matrix-recovery).

Neuron activity dynamics

For trace dynamics analysis, only the neurons that were determined to respond (and therefore displayed response dynamics to analyze) were used. The response traces of each step were aligned to the maximum first derivative value around the stimulus presentation and z-normalized. Each trace lasted from 20 s before the stimulus presentation to 50 s after. PCA was performed on the individual traces, and the mean trace was calculated per neuron, per stimulus and projected onto the PC space. Clustering was performed using K-means, with the number of clusters chosen using silhouette scores.

Stimulus identity classification

A random forest classifier (with 100 trees and a depth of 4) was used to predict the identity of a stimulus from neuronal responses. ON and OFF step responses of all neurons from each worm were concatenated into 26-dimensional vectors, which were then randomly divided into training (80%) and test (20%) sets. For single-neuron classification, each point consisted of the ON and OFF step responses of one neuron. For cross-validation, the model was trained 10 times on different train-test sets, and the average F-score across all 10 trials was calculated.

The contribution of response amplitudes and trace features to classification by neuron or stimulus was calculated by taking the mean F-score across stimuli/neurons respectively. For example, the success of the classification of diacetyl was calculated as the mean F-score of diacetyl across all individual neurons. Similarly, success of classification by AWA was calculated as the mean F-score of AWA across all stimuli. F-scores were used to evaluate the performance of the classifier because of the imbalanced sample sizes in the data.

Classification by dynamics

For the purpose of classification, the dimensionality of individual step responses was reduced to 22 trace features using the Catch22 algorithm [37]. These features were then z-normalized and PCA was performed. The first 3 principal components were used as additional features for the random forest classifier (https://github.com/DynamicsAndNeuralSystems/catch22).

Resampling data

The data provided in [16] contain neuron activity traces grouped by the stimulus presented. Our classifier, however, requires the data to be arranged as “complete” worms, where each observation (worm) consists of the activity features of all of the amphid neurons. These “synthetic” worms were constructed by independently and randomly sampling each neuron from all of the responses of that specific neuron to a specific stimulus, e.g., a synthetic diacetyl-sensing worm was composed of a random AWA response to diacetyl, an independently randomly chosen AWB response to diacetyl, etc. Where trace features were also used, all features of a neuron were sampled from the same response trace (i.e., the amplitude and trace feature PCs were taken together, rather than each being independently sampled). The model was cross-validated 10 times with different train-test (80%/20%) splits. Because the resampling of the data equalized the sample sizes between conditions, accuracy was used as the evaluation metric.