Worm maintenance

C. elegans were stored in the dark, and only minimal light was used when transferring worms or mounting worms for experiments. Strains generated in this study (Extended Data Fig. 1a) have been deposited in the Caenorhabditis Genetics Center (CGC), University of Minnesota, for public distribution. Hermaphrodites were used in this study.

Transgenics

We generated a transgenic worm for interrogating signal propagation, TWISP (AML462), which has been described in more detail previously²³. This strain expresses the calcium indicator GCaMP6s in the nucleus of each neuron; a purple-light-sensitive optogenetic protein system (GUR-3 and PRDX-2) in each neuron; and multiple fluorophores of various colours from the NeuroPAL²⁷ system, also in the nucleus of neurons. We also used a QF-hGR drug-inducible gene-expression strategy to turn on the gene expression of optogenetic actuators only later in development. To create this strain, we first generated an intermediate strain, AML456, by injecting a plasmid mix (75 ng μl⁻¹ pAS3-5xQUAS::Δ pes-10P::AI::gur-3G::unc-54 + 75 ng μl⁻¹ pAS3-5xQUAS::Δ pes-10P::AI::prdx-2G::unc-54 + 35 ng μl⁻¹ pAS-3-rab-3P::AI::QF+GR::unc-54 + 100 ng μl⁻¹ unc-122::GFP) into CZ20310 worms, followed by UV integration and six outcrosses^56,57. The intermediate strain, AML456, was then crossed into the pan-neuronal GCaMP6s calcium-imaging strain, with NeuroPAL, AML320 (refs. ^23,27,58).

Animals exhibited decreased average locomotion compared to the WT (mean speeds of 0.03 mm s⁻¹ off drug and 0.02 mm s⁻¹ on drug compared to the mean of 0.15 mm s⁻¹ in WT animals²³), as expected for NeuroPAL GCaMP6s strains, which are also reported to be overall less active (around 0.09 mm s⁻¹ during only forward locomotion)²⁷.

An unc-31-mutant background with defects in the dense-core-vesicle-release pathway was used to diminish wireless signalling⁵³. We created an unc-31-knockout version of our functional connectivity strain by performing CRISPR–Cas9-mediated genome editing on AML462 using a single-strand oligodeoxynucleotide (ssODN)-based homology-dependent repair strategy⁵⁹. This approach resulted in strain AML508 (unc-31 (wtf502) IV; otIs669 (NeuroPAL) V 14x; wtfIs145 (30 ng μl⁻¹ pBX + 30 ng μl⁻¹ rab-3::his-24::GCaMP6s::unc-54); wtfIs348 (75 ng μl⁻¹ pAS3-5xQUAS::Δ pes-10P::AI::gur-3G::unc-54 + 75 ng μl⁻¹ pAS3-5xQUAS::Δ pes-10P::AI::prdx-2G::unc-54 + 35 ng μl⁻¹ pAS-3-rab-3P::QF+GR::unc-54 + 100 ng μl⁻¹ unc-122::GFP)).

CRISPR–Cas-9 editing was carried out as follows. Protospacer adjacent motif (PAM) sites (denoted in upper case) were selected in the first intron (gagcuucgcaauguugacucCGG) and the last intron (augguacauuggguccguggCGG) of the unc-31 gene (ZK897.1a.1) to delete 12,476 out of 13,169 bp (including the 5′ and 3′ untranslated regions) and 18 out of 20 exons from the genomic locus, while adding 6 bp (GGTACC) for the Kpn-I restriction site (Extended Data Fig. 1b). Alt-R S.p. Cas9 Nuclease V3, Alt-R-single guide RNA (sgRNA) and Alt-R homology-directed repair (HDR)-ODN were used (IDT). We introduced the Kpn-I restriction site, denoted in upper case (gacccagcgaagcaaggatattgaaaacataagtacccttgttgttgtgtGGTACCccacggacccaatgtaccatattttacgagaaatttataatgttcagg) into our repair oligonucleotide to screen and confirm the deletion by PCR followed by restriction digestion. sgRNA and HDR ssODNs were also synthesized for the dpy-10 gene as a reporter, as described previously⁵⁹. An injection mix was prepared by sequentially adding Alt-R S.p. Cas9 Nuclease V3 (1 μl of 10 μg μl⁻¹), 0.25 μl of 1 M KCL, 0.375 μl of 200 mM HEPES (pH 7.4), sgRNAs for unc-31 (1 μl each for both sites) and 0.75 μl for dpy-10 from a stock of 100 μM, ssODNs (1 μl for unc-31 and 0.5 μl for dpy-10 from a stock of 25 μM) and nuclease-free water to a final volume of 10 μl in a PCR tube, kept on ice. The injection mix was then incubated at 37 °C for 15 min before it was injected into the germline of AML462 worms. Progenies from plates showing roller or dumpy phenotypes in the F₁ generation after injection were individually propagated and screened by PCR and Kpn-I digestion to confirm deletion. Single-worm PCR was carried out using GXL-PRIME STAR taq-Polymerase (Takara Bio) and the Kpn-1-HF restriction enzyme (NEB). Worms without a roller or dumpy phenotype and homozygous for deletion were confirmed by Sanger sequencing fragment analysis.

To cross-validate GUR-3/PRDX-2-evoked behaviour responses, we generated the transgenic strain AML546 by injecting a plasmid mix (40 ng μl⁻¹ pAS3-rig-3P::AI::gur-3G::SL2::tagRFP::unc-54 + 40 ng μl⁻¹ pAS3-rig-3P::AI::prdx-2G::SL2::tagBFP::unc-54) into N2 worms to generate a transient transgenic line expressing GUR-3/PRDX-2 in AVA neurons.

Cross-validation of GUR-3/PRDX-2-evoked behaviour

Optogenetic activation of AVA neurons using traditional channelrhodopsins (for example, Chrimson) leads to reversals^45,60. We used worms expressing GUR-3/PRDX-2 in AVA neurons (AML564) to show that GUR-3/PRDX-2 elicits a similar behavioural response. We illuminated freely moving worms with blue light from an LED (peaked at 480 nm, 2.3 mW mm⁻²) for 45 s. We compared the number of onsets of reversals in that period of time with a control in which only dim white light was present, as well as with the results of the same assay performed on N2 worms. Animals with GUR-3/PRDX-2 in AVA (n = 11 animals) exhibited more blue-light-evoked reversals per minute than did WT animals (n = 8 animals) (Extended Data Fig. 2h).

Dexamethasone treatment

To increase the expression of optogenetic proteins while avoiding arrested development, longer generation time and lethality, a drug-inducible gene-expression strategy was used. Dexamethasone (dex) activates QF-hGR to temporally control the expression of downstream targets⁶¹, in this case the optogenetic proteins in the functional connectivity imaging strains AML462 and AML508. Dex-NGM plates were prepared by adding 200 μM of dex in dimethyl sulfoxide (DMSO) just before pouring the plate. For dex treatment, L2/L3 worms were transferred to overnight-seeded dex-NGM plates and further grown until worms were ready for imaging. More details of the dex treatment are provided below.

We prepared stock solution of 100 mM dex by dissolving 1 g dexamethasone (D1756, Sigma-Aldrich) in 25.5 ml DMSO (D8418, Sigma-Aldrich). Stocks were then filter-sterilized, aliquoted, wrapped in foil to prevent light and stored at −80 °C until needed. The 200-μM dex-NGM plates were made by adding 2 ml of 100 mM dex stock in 1 l NGM-agar medium, while stirring, 5 min before pouring the plate. Dex plates were stored at 4 °C for up to a month until needed.

Preparation of worms for imaging

Worms were individually mounted on 10% agarose pads prepared with M9 buffer and immobilized using 2 μl of 100-nm polystyrene beads solution and 2 μl of levamisole (500 μM stock). This concentration of levamisole, after dilution in the polystyrene bead solution and the agarose pad water, largely immobilized the worm while still allowing it to slightly move, especially before placing the coverslip. Pharyngeal pumping was observed during imaging.

Overview of the imaging strategy

We combined whole-brain calcium imaging through spinning disk single-photon confocal microscopy^62,63 with two-photon⁶⁴ targeted optogenetic stimulation⁶⁵, each with their own remote focusing system, to measure and manipulate neural activity in an immobilized animal (Fig. 1a). We performed calcium imaging, with excitation light at a wavelength and intensity that does not elicit photoactivation of GUR-3/PRDX-2 (ref. ⁶⁶) (Extended Data Fig. 2b). We also used genetically encoded fluorophores from NeuroPAL expressed in each neuron²⁷ to identify neurons consistently across animals (Fig. 1c).

Multi-channel imaging and neural identification

Volumetric, multi-channel imaging was performed to capture images of the following fluorophores in the NeuroPAL transgene: mtagBFP2, CyOFP1.5, tagRFP-T and mNeptune2.5 (ref. ²⁷). Light downstream of the same spinning disk unit used for calcium imaging travelled on an alternative light path through channel-specific filters mounted on a mechanical filter wheel, while mechanical shutters alternated illumination with the respective lasers, similar to a previously described method⁵⁸. Channels were as follows: mtagBFP2 was imaged using a 405-nm laser and a Semrock FF01-440/40 emission filter; CyOFP1.5 was imaged using a 505-nm laser and a Semrock 609/54 emission filter; tagRFP-T was imaged using a 561-nm laser and a Semrock 609/54-nm emission filter; and mNeptune2.5 was imaged using a 561-nm laser and a Semrock 732/68-nm emission filter.

After the functional connectivity recording was complete, neuron identities were manually assigned by comparing each neuron’s colour, position and size to a known atlas. Some neurons are particularly hard to identify in NeuroPAL and are therefore absent or less frequently identified in our recordings. Some neurons have dim tagRFP-T expression, which makes it difficult for the neuron segmentation algorithm to find them and, therefore, to extract their calcium activity. These neurons include, for example, AVB, ADF and RID. RID’s distinctive position and its expression of CyOFP allowed us nevertheless to manually target it optogenetically. Neurons in the ventral ganglion are hard to identify because it appears as very crowded when viewed in the most common orientation that worms assume when mounted on a microscope slide. Neurons in the ventral ganglion are therefore sometimes difficult to distinguish from one another, especially for dimmer neurons such as the SIA, SIB and RMF neurons. In our strain, the neurons AWCon and AWCoff were difficult to tell apart on the basis of colour information.

Volumetric image acquisition

Neural activity was recorded at whole-brain scale and cellular resolution through continuous acquisition of volumetric images in the red and green channels with a spinning disk confocal unit and using LabView software (https://github.com/leiferlab/pump-probe-acquisition/tree/pp), similarly to a previous study⁶⁷, with a few upgrades. The imaging focal plane was scanned through the brain of the worm remotely using an electrically tunable lens (Optotune EL-16-40-TC) instead of moving the objective. The use of remote focusing allowed us to decouple the z-position of the imaging focal plane and that of the optogenetics two-photon spot (described below).

Images were acquired by an sCMOS camera, and each acquired image frame was associated to the focal length of the tunable lens (z-position in the sample) at which it was acquired. To ensure the correct association between frames and z-position, we recorded the analogue signal describing the focal length of the tunable lens at time points synchronous with a trigger pulse output by the camera. By counting the camera triggers from the start of the recording, the z-positions could be associated to the correct frame, bypassing unknown operating-system-mediated latencies between the image stream from the camera and the acquisition of analogue signals.

In addition, real-time ‘pseudo’-segmentation of the neurons (described below) required the ability to separate frames into corresponding volumetric images in real time. Because the z-position was acquired at a low sample rate, splitting of volumes on the basis of finite differences between successive z-positions could lead to errors in assignment at the edge of the z-scan. An analogue OP-AMP-based differentiator was used to independently detect the direction of the z-scan in hardware.

Calcium imaging

Calcium imaging was performed in a single-photon regime with a 505-nm excitation laser through spinning disk confocal microscopy, at 2 vol s⁻¹. For functional connectivity experiments, an intensity of 1.4 mW mm⁻² at the sample plane was used to image GCaMP6s, well below the threshold needed to excite the GUR-3/PRDX-2 optogenetic system²⁴. We note that at this wavelength and intensity, animals exhibited very little spontaneous calcium activity.

For certain analyses (Fig. 6), recordings with ample spontaneous activity were desired. In those cases, we increased the 505-nm intensity sevenfold, to approximately 10 mW mm⁻², and recorded from AML320 strains that lacked exogenous GUR-3/PRDX-2 to avoid potential widespread neural activation. Under these imaging conditions, we observed population-wide slow stereotyped spontaneous oscillatory calcium dynamics, as previously reported^35,68.

Extraction of calcium activity from the images

Calcium activity was extracted from the raw images by using Python libraries implementing optimized versions of a previously described algorithm⁶⁹, available at https://www.github.com/leiferlab/pumpprobe, https://www.github.com/leiferlab/wormdatamodel, https://www.github.com/leiferlab/wormneuronsegmentation-c and https://www.github.com/leiferlab/wormbrain.

The positions of neurons in each acquired volume were determined by computer vision software implemented in C++. This software was greatly optimized to identify neurons in real time, to also enable closed-loop targeting and stimulus delivery (as described in ‘Stimulus delivery and pulsed laser’). Two design choices made this algorithm considerably faster than previous approaches. First, a local maxima search was used instead of a slower watershed-type segmentation. The nuclei of C. elegans neurons are approximately spheres and so they can be identified and separated by a simple local maxima search. Second, we factorized the three-dimensional (3D) local maxima search into multiple two-dimensional (2D) local maxima searches. In fact, any local maximum in a 3D image is also a local maximum in the 2D image in which it is located. Local maxima were therefore first found in each 2D image separately, and then candidate local maxima were discarded or retained by comparing them to their immediate surroundings in the other planes. This makes the algorithm less computationally intensive and fast enough to also be used in real time. We refer to this type of algorithm as ‘pseudo’-segmentation because it finds the centre of neurons without fully describing the extent and boundaries of each neuron.

After neural locations were found in each of the volumetric images, a nonrigid point-set registration algorithm was used to track their locations across time, matching neurons identified in a given 3D image to the neurons identified in a 3D image chosen as reference. Even worms that are mechanically immobilized still move slightly and contract their pharynx, thereby deforming their brain and requiring the tracking of neurons. We implemented in C++ a fast and optimized version of the Dirichelet–Student’s-t mixture model (DSMM)⁷⁰.

Calcium pre-processing

The GCaMP6s intensity extracted from the images undergoes the following pre-processing steps. (1) Missing values are interpolated on the basis of neighbouring time points. Missing values can occur when a neuron cannot be identified in a given volumetric image. (2) Photobleaching is removed by fitting a double exponential to the baseline signal. (3) Outliers more than 5 standard deviations away from the average are removed from each trace. (4) Traces are smoothed using a causal polynomial filtering with a window size of 6.5 s and polynomial order of 1 (Savitzky–Golay filters with windows completely ‘in the past’; for example, obtained with scipy.signal.savgol_coeffs(window_length=13, polyorder=1, pos=12)). This type of filter with the chosen parameters is able to remove noise without smearing the traces in time. Note that when fits are performed (for example, to calculate kernels), they are always performed on the original, non-smoothed traces. (5) Where ΔF/F₀ of responses is used, F₀ is defined as the value of F in a 30-s interval before the stimulation time and ΔF ≡ F − F₀. In Fig. 2a, for example, \({ < \Delta F/{F}_{0} > }_{t}\) refers to the mean over a 30-s post-stimulus window.

Stimulus delivery and pulsed laser

For two-photon optogenetic targeting, we used an optical parametric amplifier (OPA; Light Conversion ORPHEUS) pumped by a femtosecond amplified laser (Light Conversion PHAROS). The output of the OPA was tuned to a wavelength of 850 nm, at a 500 kHz repetition rate. We used temporal focusing to spatially restrict the size of the two-photon excitation spot along the microscope axis. A motorized iris was used to set its lateral size. For temporal focusing, the first-order diffraction from a reflective grating, oriented orthogonally to the microscope axis, was collected (as described previously⁷¹) and travelled through the motorized iris, placed on a plane conjugate to the grating. To arbitrarily position the two-photon excitation spot in the sample volume, the beam then travelled through an electrically tunable lens (Optotune EL-16-40-TC, on a plane conjugate to the objective), to set its position along the microscope axis, and finally was reflected by two galvo-mirrors to set its lateral position. The pulsed beam was then combined with the imaging light path by a dichroic mirror immediately before entering the back of the objective.

Most of the stimuli were delivered automatically by computer control. Real-time computer vision software found the position of the neurons for each volumetric image acquired, using only the tagRFP-T channel. To find neural positions, we used the same pseudo-segmentation algorithm described above. The algorithm found neurons in each 2D frame in around 500 μs as the frames arrived from the camera. In this way, locations for all neurons in a volume were found within a few milliseconds of acquiring the last frame of that volume.

Every 30 s, a random neuron was selected among the neurons found in the current volumetric image, on the basis of only its tagRFP-T signal. After galvo-mirrors and the tunable lens set the position of the two-photon spot on that neuron, a 500-ms (300-ms for the unc-31-mutant strain) train of light pulses was used to optogenetically stimulate that neuron. The duration of stimulus illumination for the unc-31-mutant strain was selected to elicit calcium transients in stimulated neurons with a distribution of amplitudes such that the maximum amplitude was similar to those in WT-background animals, (Extended Data Fig. 2f). The output of the laser was controlled through the external interface to its built-in pulse picker, and the power of the laser at the sample was 1.2 mW at 500 kHz. Neuron identities were assigned to stimulated neurons after the completion of experiments using NeuroPAL²⁷.

To probe the AFD→AIY neural connection, a small set of stimuli used variable pulse durations from 100 ms to 500 ms in steps of 50 ms selected randomly to vary the amount of optogenetic activation of AFD.

In some cases, neurons of interest were too dim to be detected by the real-time software. For those neurons of interest, additional recordings were performed in which the neuron to be stimulated was manually selected on the basis of its colour, size and position. This was the case for certain stimulations of neurons RID and AFD.

Characterization of the size of the two-photon excitation spot

The lateral (xy) size of the two-photon excitation spot was measured with a fluorescent microscope slide, and the axial (z) size was measured using 0.2-nm fluorescent beads (Suncoast Yellow, Bangs Laboratories), by scanning the z-position of the optogenetic spot while maintaining the imaging focal plane fixed (Extended Data Fig. 2a).

We further tested our targeted stimulation in two ways: selective photobleaching and neuronal activation. First, we targeted individual neurons at various depths in the worm’s brain, and we illuminated them with the pulsed laser to induce selective photobleaching of tagRFP-T. Extended Data Fig. 2c,d shows how our two-photon excitation spot selectively targets individual neurons, because it photobleaches tagRFP-T only in the neuron that we decide to target, and not in nearby neurons. To faithfully characterize the spot size, we set the laser power such that the two-photon interaction probability profile of the excitation spot would not saturate the two-photon absorption probability of tagRFP-T. Second, we showed that our excitation spot is restricted along the z-axis by targeting a neuron and observing its calcium activity. When the excitation was directed at the neuron but shifted by 4 μm along z, the neuron showed no activation. By contrast, the neuron showed activation when the spot was correctly positioned on the neuron (Extended Data Fig. 2e). To further show that our stimulation is spatially restricted to an individual neuron more broadly throughout our measurements, we show that stimulations do not elicit responses in most of the close neighbours of the targeted neurons (Extended Data Fig. 2i and Supplementary Information).

Inclusion criteria

Stimulation events were included for further analysis if they evoked a detectable calcium response in the stimulated neuron (autoresponse). A classifier determined whether the response was detected by inspecting whether the amplitude of both the ΔF/F₀ transient and its second derivative exceeded a pair of thresholds. The same threshold values were applied to every animal, strain, neuron and stimulation event, and were originally set to match the human perception of a response above noise. Stimulation events that did not meet both thresholds for a contiguous 4 s were excluded. The RID responses shown in Fig. 4 and Extended Data Fig. 7c are an exception to this policy. RID is visible on the basis of its CyOFP expression, but its tagRFP-T expression is too dim to consistently extract calcium signals. Therefore, in Fig. 4 and Extended Data Fig. 7c (but not in other figures, such as Fig. 2), downstream neurons’ responses to RID stimulation were included even in cases in which it was not possible to extract a calcium-activity trace in RID.

Neuron traces were excluded from analysis if a human was unable to assign an identity or if the imaging time points were absent in a contiguous segment longer than 5% of the response window owing to imaging artefacts or tracking errors. A different policy applies to dim neurons of interest that are not automatically detected by the pseudo-segmentation algorithm in the 3D image used as reference for the point-set registration algorithm. In those cases, we manually added the position of those neurons to the reference 3D image. If these ‘added’ neurons are automatically detected in most of the other 3D images, then a calcium activity trace can be successfully produced by the DSMM nonrigid registration algorithm, and is treated as any other trace. However, if the ‘added’ neurons are too dim to be detected also in the other 3D images and the calcium activity trace cannot be formed for more than 50% of the total time points, the activity trace for those neurons is extracted from the neuron’s position as determined from the position of neighbouring neurons. In the analysis code, we refer to these as ‘matchless’ traces, because the reference neuron is not matched to any detected neuron in the specific 3D image, but its position is just transformed according to the DSMM nonrigid deformation field. In this way, we are able to recover the calcium activity of some neurons whose tagRFP-T expression is otherwise too dim to be reliably detected by the pseudo-segmentation algorithm. Responses to RID stimulation shown in Fig. 4 and Extended Data Fig. 7c are an exception to this policy. In these cases, the activity of any neuron for which there is not a trace for more than 50% of the time points is substituted with the corresponding ‘matchless’ trace, and not just for the manually added neurons. This is important to be able to show responses of neurons such as ADL, which have dim tagRFP-T expression. In the RID-specific case, to exclude responses that become very large solely because of numerical issues in the division by the baseline activity owing to the dim tagRFP-T, we also introduce a threshold excluding ΔF/F > 2.

Kernels were computed only for stimulation-response events for which the automatic classifier detected responses in both the stimulated and the downstream neurons. If the downstream neuron did not show a response, we considered the downstream response to be below the noise level and the kernel to be zero.

Statistical analysis

We used two statistical tests to identify neuron pairs that under our stimulation and imaging conditions can be deemed ‘functionally connected’, ‘functionally non-connected’ or for which we lack the confidence to make either determination. Both tests compare observed calcium transients in each downstream neuron to a null distribution of transients recorded in experiments lacking stimulation.

To determine whether a pair of neurons can be deemed functionally connected, we calculated the probability of observing the measured calcium response in the downstream neuron given no neural stimulation. We used a two-sided Kolmogorov–Smirnov test to compare the distributions of the downstream neuron’s ΔF/F₀ amplitude and its temporal second derivative from all observations of that neuron pair under stimulation to the empirical null distributions taken from control recordings lacking stimulation. P values were calculated separately for ΔF/F₀ and its temporal second derivative, and then combined using Fischer’s method to report a single fused P value for each neuron pair. Finally, to account for the large number of hypotheses tested, a false discovery rate was estimated. From the list of P values, each neuron was assigned a q value using the Storey–Tibshirani method⁴⁰. q values are interpreted as follows: when considering an ensemble of putative functional connections of q values all less than or equal to q_c, an approximately q_c fraction of those connections would have appeared in a recording that lacked any stimulation.

To explicitly test whether a pair of neurons are functionally not connected, taking into account the amplitude of the response, their reliability, the number of observations and multiple hypotheses, we also computed equivalence P_eq and q_eq values. This assesses the confidence of a pair not being connected. We test whether our response is equivalent to what we would expect from our control distribution using the two one-sided t-test (TOST)⁷². We computed P_eq values for ΔF/F₀ and its temporal second derivative for a given pair being equivalent to the control distributions within an \({\epsilon }=1.2{\sigma }_{\Delta F/{F}_{0},{\partial }_{t}^{2}}\). Here, \({\sigma }_{\Delta F/{F}_{0},{\partial }_{t}^{2}}\) is the standard deviation of the corresponding control distribution. We then combined the two P_eq values into a single one with the Fisher method and computed q_eq values using the Storey–Tibshirani method⁴⁰. Note that, different from the regular P values described above, the equivalence test relies on the arbitrary choice of ϵ, which defines when we call two distributions equivalent. We chose a conservative value of ϵ = 1.2σ.

We note that the statistical framework is stringent and a large fraction of measured neuron pairs fail to pass either statistical test.

Measuring path length through the synaptic network

To find the minimum path length between neurons in the anatomical network topology, we proceeded iteratively. We started from the original binary connectome and computed the map of strictly two-hop connections by looking for pairs of neurons that are not connected in the starting connectome (the actual anatomical connectome at the first step) but that are connected through a single intermediate neuron. To generate the strictly three-hop connectome, we repeated this procedure using the binary connectome including direct and two-hop connections, as the starting connectome. This process continued iteratively to generate the strictly n-hop connectome.

In the anatomical connectome (the starting connectome for the first step in the procedure above), a neuron was considered to be directly anatomically connected if the connectomes of any of the four L4 or adult individuals in refs. ¹ and ⁶ contained at least one synaptic contact between them. Note that this is a permissive description of anatomical connections, as it considers even neurons with only a single synaptic contact in only one individual to be connected.

Fitting kernels

Kernels k_ij(t) were defined as the functions to be convolved with the activity ΔF_j of the stimulated neuron to obtain the activity ΔF_i of a responding neuron i, such that \(\Delta {F}_{i}(t)=({k}_{ij}\ast \Delta {F}_{j})(t)\). To fit kernels, each kernel k(t) was parametrized as a sum of convolutions of decaying exponentials

where the indices i, j are omitted for clarity and θ is the Heaviside function. This parametrization is exact for linear systems, and works as a description of causal signal transmission also in nonlinear systems. Note that increasing the number of terms in the successive convolutions does not lead to overfitting, as would occur by increasing the degree of a polynomial. Overfitting could occur by increasing the number of terms in the sum, which in our fitting is constrained to be a maximum of 2. The presence of two terms in the sum allows the kernels to represent signal transmission with saturation (with c₀ and c₁ of opposite signs) and assume a fractional-derivative-like shape.

The convolutions are performed symbolically. The construction of kernels as in equation (1) starts from a symbolically stored, normalized decaying exponential kernel with a factor \(A,A{\gamma }_{0}\theta (t){e}^{-{\gamma }_{0}t}\). Convolutions with normalized exponentials \({\gamma }_{n}\theta (t){e}^{-{\gamma }_{n}t}\) are performed sequentially and symbolically, taking advantage of the fact that successive convolutions of exponentials always produce a sum of functions in the form ∝ θ(t)tⁿe^−γt. Once rules are found to convolve an additional exponential with a function in that form, any number of successive convolution can be performed. These rules are as follows:

1. 1.

   If the initial term is a simple exponential with a given factor (not necessarily just the normalization γ) \({c}_{i}\theta (t){e}^{-{\gamma }_{i}t}\) and γ_i ≠ γ_n, then the convolution is

   $${c}_{i}\theta (t){e}^{-{\gamma }_{i}t}\ast {\gamma }_{n}\theta (t){e}^{-{\gamma }_{n}t}={c}_{\mu }\theta (t){e}^{-{\gamma }_{\mu }t}+{c}_{\nu }\theta (t){e}^{-{\gamma }_{\nu }t},$$

   (2)

   with \({c}_{\mu }=\frac{{c}_{i}{\gamma }_{n}}{{\gamma }_{n}-{\gamma }_{i}},{c}_{\nu }=-\frac{{c}_{i}{\gamma }_{n}}{{\gamma }_{n}-{\gamma }_{i}}\) and γ_μ = γ_i, γ_ν = γ_n.

2. 2.

   If the initial term is a simple exponential and γ_i = γ_n, then

   $${c}_{i}\theta (t){e}^{-{\gamma }_{i}t}\ast {\gamma }_{n}\theta (t){e}^{-{\gamma }_{n}t}={c}_{\mu }\theta (t)t{e}^{-{\gamma }_{\mu }t},$$

   (3)

   with c_μ = c_iγ_i and γ_μ = γ_i.

3. 3.

   If the initial term is a \({c}_{i}\theta (t){t}^{n}{e}^{-{\gamma }_{i}t}\) term and γ_i = γ_μ, then

   $${c}_{i}\theta (t){t}^{n}{e}^{-{\gamma }_{i}t}\ast {\gamma }_{n}\theta (t){e}^{-{\gamma }_{n}t}={c}_{\mu }\theta (t){t}^{n+1}{e}^{-{\gamma }_{\mu }t},$$

   (4)

   with \({c}_{\mu }=\frac{{c}_{i}{\gamma }_{i}}{n+1}\) and γ_μ = γ_i.

4. 4.

   If the initial term is a \({c}_{i}\theta (t){t}^{n}{e}^{-{\gamma }_{i}t}\) term and γ_i ≠ γ_μ, then

   $${c}_{i}\theta (t){t}^{n}{e}^{-{\gamma }_{i}t}\ast {\gamma }_{n}\theta (t){e}^{-{\gamma }_{n}t}={c}_{\mu }\theta (t){t}^{n}{e}^{-{\gamma }_{\mu }t}+{c}_{\nu }(\theta (t){t}^{n-1}{e}^{-{\gamma }_{i}t}\ast \theta (t){e}^{-{\gamma }_{n}t}),$$

   (5)

   where \({c}_{\mu }=\frac{{c}_{i}{\gamma }_{n}}{{\gamma }_{n}-{\gamma }_{i}},{\gamma }_{\mu }={\gamma }_{i}\), and \({c}_{\nu }=-n\frac{{c}_{i}{\gamma }_{n}}{{\gamma }_{n}-{\gamma }_{i}}\).

Additional terms in the sum in equation (1) can be introduced by keeping track of the index m of the summation for every term and selectively convolving new exponentials only with the corresponding terms.

Kernel-based simulations of activity

Using the kernels fitted from our functional data, we can simulate neural activity without making any further assumptions about the dynamical equations of the network of neurons. To compute the response of a neuron i to the stimulation of a neuron j, we simply convolve the kernel k_i,j(t) with the activity ΔF_j(t) induced by the stimulation in neuron j. The activity of the stimulated neuron can be either the experimentally observed activity or an arbitrarily shaped activity introduced for the purposes of simulation.

To compute kernel-derived neural activity correlations (Fig. 6), we completed the following steps. (1) We computed the responses of all the neurons i to the stimulation of a neuron j chosen to drive activity in the network. To compute the responses, for each pair i, j, we used the kernel \({\langle {k}_{i,j}(t)\rangle }_{{\rm{trials}}}\) averaged over multiple trials. For kernel-based analysis, pairs with connections of q > 0.05 were considered not connected. We set the activity ΔF_j(t) in the driving neuron to mimic an empirically observed representative activity transient. (2) We computed the correlation coefficient of the resulting activities. (3) We repeated steps 1 and 2 for a set of driving neurons (all or top-n neurons, as in Fig. 6). (4) For each pair k, l, we took the average of the correlations obtained by driving the set of neurons j in step 3.

Anatomy-derived simulations of activity

Anatomy-derived simulations were performed as described previously⁴⁷. In brief, this simulation approach uses differential equations to model signal transmission through electrical and chemical synapses and includes a nonlinear equation for synaptic activation variables. We injected current in silico into individual neurons and simulated the responses of all the other neurons. Anatomy-derived responses (Fig. 3) of the connection from neuron j to neuron i were computed as the peak of the response of neuron i to the stimulation of j. Anatomy-based predictions of spontaneous correlations in Fig. 6 were calculated analogously to kernel-based predictions.

In one analysis in Fig. 3d, the synapse weights and polarities were allowed to float and were fitted from the functional measurements. In all other cases, synapse weights were taken as the scaled average of three adult connectomes^1,6 and an L4 connectome⁶, and polarities were assigned on the basis of a gene-expression analysis of ligand-gated ionotropic synaptic connections that considered glutamate, acetylcholine and GABA neurotransmitter and receptor expression, as performed in a previous study³⁷ and taken from CeNGEN³⁸ and other sources. Specifically, we used a previously published dataset (S1 data in ref. ³⁷) and aggregated polarities across all members of a cellular subtype (for example, polarities from source AVAL and AVAR were combined). In cases of ambiguous polarities, connections were assumed to be excitatory, as in the previous study³⁷. For other biophysical parameters we chose values commonly used in C. elegans modelling efforts^9,30,47,73.

Characterizing stereotypy of functional connections

To characterize the stereotypy of a neuron pair’s functional connection, its kernels were inspected. A kernel was calculated for every stimulus-response event in which both the upstream and downstream neuron exhibited activity that exceeded a threshold. At least two stimulus-response events that exceeded this threshold were required to calculate their stereotypy. The general strategy for calculating stereotypy was to convolve different kernels with the same stimulus inputs and compare the resulting outputs. The similarity of two outputs is reported as a Pearson’s correlation coefficient. Kernels corresponding to different stimulus-response events of the same pair of neurons were compared with one another round-robin style, one round-robin each for a given input stimulus. For inputs we chose the set of all stimuli delivered to the upstream neuron. The neuron-pairs stereotypy is reported as the average Pearson’s correlation coefficient across all round-robin kernel pairings and across all stimuli.

Rise time of kernels

The rise time of kernels, shown in Fig. 5c and Extended Data Fig. 6d, was defined as the interval between the earliest time at which the value of the kernel was 1/e its peak value and the time of its peak (whether positive or negative). The rise time was zero if the peak of the kernel was at time t = 0. However, saturation of the signal transmission can make kernels appear slower than the connection actually is. For example, the simplest instantaneous connection would be represented by a single decaying exponential in equation (1), which would have its peak at time t = 0. However, if that connection is saturating, a second, opposite-sign term in the sum is needed to fit the kernel. This second term would make the kernel have a later peak, thereby masking the instantaneous nature of the connection. To account for this effect of saturation, we removed terms representing saturation from the kernels and found the rise time of these ‘non-saturating’ kernels.

Screen for purely extrasynaptic-dependent connections

To find candidate purely extrasynaptic-dependent connections, we considered the pairs of neurons that are connected in WT animals (q^WT < 0.05) and non-connected in unc-31 animals (\({q}_{{\rm{eq}}}^{{\rm{unc-31}}}\) < 0.05, with the additional condition q^unc−31 > 0.05 to exclude very small responses that are nonetheless significantly different from the control distribution). We list these connections and provide additional examples in Extended Data Fig. 9.

Using a recent neuropeptide–GPCR interaction screen in C. elegans⁵² and gene-expression data from CeNGEN³⁸, we find putative combinations of neuropeptides and GPCRs that can mediate those connections (Supplementary Table 1). We produced such a list of neuropeptide and GPCR combinations using the Python package Worm Neuro Atlas (https://github.com/francescorandi/wormneuroatlas). In the list, we only include transcripts from CeNGEN detected with the highest confidence (threshold 4), as described previously⁵¹. For each neuron pair, we first searched the CeNGEN database for neuropeptides expressed in the upstream neuron, then identified potential GPCR targets for each neuropeptide using information from previous reports^52,74, and finally went back to the CeNGEN database to find whether the downstream neuron in the pair was among the neurons expressing the specific GPCRs. The existence of potential combinations of neuropeptide and GPCR putatively mediating signalling supports our observation that communication in the candidate neuron pairs that we identify can indeed be mediated extrasynaptically through neuropeptidergic machinery.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.