Research ArticleResearch Article: New Research, Neuronal Excitability

Neuromodulation Reduces Interindividual Variability of Neuronal Output

Anna C. Schneider, Omar Itani, Dirk Bucher and Farzan Nadim
eNeuro 19 July 2022, 9 (4) ENEURO.0166-22.2022; https://doi.org/10.1523/ENEURO.0166-22.2022

Abstract

In similar states, neural circuits produce similar outputs across individuals despite substantial interindividual variability in neuronal ionic conductances and synapses. Circuit states are largely shaped by neuromodulators that tune ionic conductances. It is therefore possible that, in addition to producing flexible circuit output, neuromodulators also contribute to output similarity despite varying ion channel expression. We studied whether neuromodulation at saturating concentrations can increase the output similarity of a single identified neuron across individual animals. Using the lateral pyloric (LP) neuron of the crab stomatogastric ganglion, we compared the variability of f–I (frequency–current) curves and rebound properties in the presence of neuropeptides. The two neuropeptides we used converge to activate the same target current, which increases neuronal excitability. Output variability was lower in the presence of the neuropeptides, regardless of whether the neuropeptides significantly changed the mean of the corresponding parameter or not. However, the addition of the second neuropeptide did not add further to the reduction of variability. With a family of computational LP-like models, we explored how increased excitability and target variability contribute to output similarity and found two mechanisms: saturation of the responses and a differential increase in baseline activity. Saturation alone can reduce the interindividual variability only if the population shares a similar ceiling for the responses. In contrast, the reduction of variability due to the increase in baseline activity is independent of ceiling effects.

  • bursting neuron
  • central pattern generator
  • stomatogastric
  • variability

Significance Statement

The activity of single neurons and neural circuits can be very similar across individuals although the ionic currents underlying activity are variable. The mechanisms that compensate for the underlying variability and promote output similarity are poorly understood but may involve neuromodulation. Using an identified neuron, we show that neuropeptide modulation of excitability can reduce interindividual variability of response properties at a single-neuron level in two ways. First, the neuropeptide increases baseline excitability in a differential manner, resulting in similar response thresholds. Second, the neuropeptide increases excitability toward a shared saturation level, promoting similar maximal firing rates across individuals. Such tuning of neuronal excitability could be an important mechanism compensating for interindividual variability of ion channel expression.

Introduction

Under similar behavioral conditions across individual animals, neural circuits often produce very similar outputs. Output similarity is observed in all parts of the central nervous system, from spinal cord circuits (Pearson and Rossignol, 1991; Masino and Fetcho, 2005) to large networks involved in learning and memory (Wang, 2010; Howe et al., 2011; Buzsáki and Wang, 2012). Interestingly, interindividual similarity of circuit output is also present in the absence of sensory feedback, as has been observed in isolated invertebrate neural circuits (Hughes and Wiersma, 1960; Selverston and Moulins, 1985; Harris-Warrick and Marder, 1991; Büschges et al., 1995; Marder and Bucher, 2001; Bucher et al., 2005; Mulloney et al., 2006; Wenning et al., 2018). Such similarity is remarkable because, in any neuron type, ionic conductances underlying activity are quite variable. Even in a neuron that exists in only a single copy in each animal, as is common in small invertebrate circuits, the expression of ion channel conductances and their corresponding mRNA levels can vary several-fold, yet the activity of that neuron in the circuit remains remarkably similar (Prinz et al., 2004; Schulz et al., 2006; Goaillard et al., 2009; Tobin et al., 2009; Ransdell et al., 2013a; Marder et al., 2015; Tran et al., 2019). This observation raises the question of how neurons and circuits produce similar outputs despite variable components. A possible explanation is that different voltage-gated conductances are coregulated in a compensatory manner to give rise to similar neuronal excitability across individuals (MacLean et al., 2003, 2005; Khorkova and Golowasch, 2007; Harris-Warrick and Johnson, 2010; Temporal et al., 2012). Alternatively, output similarity may be an emergent property of the full circuit, as compensation for intrinsic variability may include the variability of synaptic strengths or constraints on synaptic current trajectories (Prinz et al., 2004; Anwar et al., 2022).

Output similarity is observed across individuals in comparable behavioral or circuit states, but circuit output is flexible and can change substantially depending on the behavioral needs of the animal. One mechanism providing such flexibility is neuromodulation by a variety of transmitters and hormones that modify properties of ion channels, synaptic release mechanisms, and other circuit components (Marder and Weimann, 1992; Katz et al., 1994; Parker, 2000; Johnson et al., 2011; Bargmann, 2012; Marder, 2012; Nadim and Bucher, 2014). This poses an interesting problem for output similarity. Neuromodulator receptor expression itself can show substantial interindividual variability (Garcia et al., 2015), suggesting that, in addition to the variability of expression, ionic conductances are also subject to variability of neuromodulator effects.

In the pyloric circuit of the stomatogastric ganglion (STG), conductance variability is substantial in both unmodulated and modulated states (Schulz et al., 2006; Khorkova and Golowasch, 2007; Golowasch, 2014; Anwar et al., 2022), indicating that receptor variability does not simply compensate for target ion channel variability. However, neuromodulatory effects interact with the variability of intrinsic neuronal properties, which may lead to reduction in output variability. In addition, output similarity may increase if a modulator increases neuronal excitability toward a saturation level that is shared across individuals. Voltage-gated ionic currents are nonlinear and, in any neuron type, interact in a complex but specific manner to produce neuronal excitability (Goaillard et al., 2009). It is therefore possible that, despite their varying levels, ionic currents interact in a manner that produces consistent output across individuals. The flexibility provided by neuromodulation may fine-tune these interactions to enhance neuronal output similarity. Finally, circuits are under the influence of several neuromodulators at all times, and while some neuromodulators have divergent cellular effects, others converge to modify the properties of some of the same ion channels or synapses (Bucher and Marder, 2013; Nadim and Bucher, 2014). Convergent comodulation may contribute to interindividual output similarity.

We focused on whether neuromodulation can increase output similarity at the level of an isolated neuron and, if so, whether comodulation by convergent neuromodulators enhances this similarity. We addressed these questions in the lateral pyloric (LP) neuron of the STG, which exists as a single copy in each animal. We first examined whether modulator-activated current levels in the LP neuron show similar levels of variability as other ionic currents. We then measured excitability, postinhibitory rebound properties, and response to periodic input in the synaptically isolated LP neuron. Subsequently, we compared the interindividual variability of these measures of activity in control conditions and in the presence of one or two (convergent) peptide neuromodulators. We used the experimental data to construct families of computational LP model neurons, with the same variability range as the data, to understand the mechanisms through which neuromodulation can influence interindividual variability at the single-cell level.

Materials and Methods

Experimental preparation

All experiments were performed on adult male Jonah crabs, Cancer borealis. Animals were obtained from local seafood stores in Newark, NJ, and kept unfed in tanks at 10–13°C. Crabs were placed in ice for at least 30 min for anesthetization before dissections. The stomatogastric nervous system (STNS) was dissected and pinned dorsal side up in a Sylgard (Ellsworth Adhesives) lined Petri dish. The dorsal sheath of the STG was removed with fine tungsten pins. During experiments, the STG was constantly perfused in a petroleum jelly well with 10–13°C saline at a flow rate of ∼4 ml/min. C. borealis saline contained the following (in mm): 440 NaCl, 26 MgCl2, 13 CaCl2, 11 KCl, 10 Tris base, and 5 maleic acid, buffered to pH 7.4.

The majority of chemical synapses in the pyloric circuit are graded, inhibitory, and glutamatergic, and can be blocked with picrotoxin (PTX). PTX (Sigma-Aldrich) was dissolved in DMSO (Thermo Fisher Scientific) and stored as a 10−2 m stock solution at 4°C. The PTX stock solution was diluted in saline to a final concentration of 10−5 m immediately before use. During voltage-clamp (VC) experiments, neurons were prevented from spiking by blocking sodium channels with tetrodotoxin (TTX). TTX (Alomone Labs) was dissolved in distilled water as a 10−4 m stock solution and stored at 4°C. The TTX stock solution was diluted in saline to a final concentration of 10−7 m immediately before use. Proctolin (Proc) and crustacean cardioactive peptide (CCAP; both custom synthesized by RS Synthesis) were dissolved individually in distilled water and stored as 10−3 m aliquots at −20°C. The neuropeptide stock solutions were diluted in TTX or PTX saline immediately before use to a final concentration of 10−6 M Proc and 5 × 10−7 m Proc + 5 × 10−7 m CCAP (henceforth referred to as Proc + CCAP) so that the total neuromodulator concentration in the single-modulator and double-modulator saline was always the same. All chemicals were always bath-applied to the STG.

Electrophysiology

The pyloric rhythm was recorded extracellularly with stainless steel pin electrodes inserted into petroleum jelly wells around sections of the lateral ventricular nerve (lvn), pyloric dilator nerve (pdn), and pyloric nerve. Extracellular electrodes were connected to a differential AC amplifier (model 1700, A-M Systems). The LP and pyloric dilator (PD) neurons were identified by matching their intracellular recorded activity to the extracellularly recorded pyloric rhythm.

We used two-electrode VC to measure synaptic and voltage-gated currents, and current clamp (CC) to measure excitability. Intracellular electrodes were made from thin-walled borosilicate capillaries with filaments, pulled to a sharp tip using a Flaming-Brown P-97 puller (Sutter Instrument), and filled with 0.6 m K2SO4 + 20 mm KCl (resistance: 20–25 MΩ). Intracellular signals were amplified using Axoclamp 900A amplifiers (Molecular Devices). All recordings were digitized at 5 kHz (Digidata 1440A, Molecular Devices) and recorded with Clampex 10.6 (Molecular Devices).

After cell identification, we removed all neuromodulatory inputs by transecting the stomatogastric nerve (stn; decentralization; Fig. 1). Unless otherwise indicated, we ran each of the protocols in decentralized control, Proc saline, Proc + CCAP saline, and after washing out. We did not wash between single and dual neuromodulator applications. Durations for wash in were 10 min and for wash out were 10–30 min before running a set of VC or CC protocols. In five experiments, we used sham wash in/wash out without any modulators in the saline to control for potential rundown of the LP neuron over time. Before and after each set of protocols, we measured the input resistance of the LP neuron by injecting currents pulses of −1 nA for 500 ms (10 sweeps). Experiments were discarded if input resistance was <5 MΩ.

Figure 1.
Figure 1.

Isolating the LP neuron. A, Schematic of the STNS (top) and the simplified pyloric network with electrodes to indicate which neurons were recorded from (bottom). Modulatory projection neurons are located in the commissural ganglia (CoG) and esophageal ganglion (OG) and project via the stn to the neurons in the STG. Neurons of the pyloric network are located in the STG. The network consists of a pacemaker group (one AB and two PD neurons), and several follower neurons, of which only the LP neuron provides direct chemical feedback to the pacemaker group. Chemical inhibitory synapses and their transmitters (glutamate or acetylcholine) are shown as circles, electrical coupling is depicted with resistor symbols. B–D, Extracellular recordings of the lvn and pdn, which carry axons of both LP (largest unit on lvn) and PD neurons (mid-sized units on lvn), or only of PD neurons, respectively, and intracellular recordings of LP and the two copies of PD. B, When the preparation is intact (all intrinsic neuromodulators present), the LP neuron receives strong periodic inhibition from the pacemaker group. C, After decentralization (removal of intrinsic neuromodulators by transecting the stn), the pyloric rhythm deteriorates but synaptic connections are still functional as illustrated by LP inhibition during PD bursts, and IPSPs in PD for each LP spike. D, After the addition of 10−5 m PTX, glutamatergic synapses between the pacemaker group and LP are blocked. The cholinergic synapse from PD to LP is still functional. When PDs are hyperpolarized (traces clipped), this synapse is silenced. Recordings in all three panels are from the same experiment and, with exception of the inset in B, on the same voltage scale.

Measurement of ionic currents

For VC experiments, we applied 10−7 m TTX to block transient sodium currents, and then measured several ionic currents in the LP neuron in control and Proc saline. These currents included delayed rectifier and calcium-activated potassium currents measured together as the high-threshold potassium current (IHTK; Khorkova and Golowasch, 2007), the transient potassium current (IA), the hyperpolarization-activated inward current (Ih), the synaptic current from the PD to LP synapse (Isyn), and the modulator-activated inward current (IMI; Golowasch and Marder, 1992b; Fig. 2A). For these measurements, the PD neuron was always held at −60 mV, except during the Isyn measurements. For the total potassium currents (IK), the LP neuron was held at −80 mV and stepped for 500 ms from −50 to 20 mV in 10 mV increments. IHTK was measured in the same way, except that LP was held at −40 mV to inactivate IA. IA was calculated as the difference current between IK and IHTK. To measure Ih, we stepped LP from −40 to −120 mV for 10 s. For Isyn, the LP neuron was held at −50 mV and the presynaptic PD neuron was stepped from −70 to 0 mV in 10 mV increments and a 500 ms step duration. To measure the noninactivating IMI, the voltage of the LP neuron was ramped from −80 to 20 mV (75 mV/s) and back to −80 mV. The I–V curve was then calculated as the difference between the current responses to the negative ramp in Proc and control saline (Schneider et al., 2021).

Figure 2.
Figure 2.

Variability of modulated components is not different from nonmodulated components across individuals. A, Example voltage-clamp recordings for nonmodulated and modulated currents. Current traces for IHTK [transient (t) portion indicated by black arrow, persistent (p) by gray arrow], IA (with total potassium currents IK as gray overlay; IA is the difference current between all potassium currents IK and IHTK); Ih; the synaptic current Isyn from PD to LP; and the neuromodulator-activated, voltage-gated current IMI are shown with the corresponding voltage-clamp protocol as the inset. B, From these voltage traces, the features indicated with white markers were extracted: the maximum transient current for both potassium currents and additional the persistent current for IHTK, the amplitude of the current at the end of the hyperpolarizing voltage step for Ih; the scaling factor for the sigmoid fit for Isyn; and the maximum inward current of the current fit for IMI. C, The distribution of these parameters across experiments are shown. Each dot represents one experiment, and red lines mark the mean. Black is for control (decentralized), blue is for Proc. Transient IHTK is shown in bold colors, persistent in transparent colors. Since IMI is calculated as a difference current (Proc – ctrl), there are no data in control. D, Coefficients of variation (SD/mean) are in the same range for the nonmodulated currents (IHTK, IA, Ih) and modulated currents (Isyn, IMI).

We used the maximum current level of each current in each experiment to calculate the variability of that component. For potassium currents, we used the value in response to the voltage step to +20 mV (Fig. 2B, open circles) for all further calculations. IA was calculated as the difference in current between IK (Fig. 2A, gray traces) and IHTK. For IHTK, we measured the transient and persistent currents separately, indicated by the arrows “t” and “p” in Figure 2A. For Ih, we calculated the current amplitude between the beginning of the hyperpolarizing step and the end, indicated by the arrow in Figure 2B. For Isyn, we first fitted the I–V data from the presynaptic voltage and average postsynaptic current during each presynaptic voltage step with a logistic sigmoid function, as follows: f(V)=a1 + expV−V1/2k, (1)where a is the maximum postsynaptic current, V1/2 is the voltage at midpoint, and k is the slope factor at V1/2. Lower bounds were set at −10 nA for a and −80 mV for V1/2. We used the value of a as the maximum synaptic current, indicated by the open circle in Figure 2B.

For IMI , we first calculated the difference current between Proc and control. From this, we separated the I–V curves for the positive and negative voltage ramp because the positive ramp activates an additional inactivating component, IMI-T (Schneider et al., 2021). The I–V curves of the negative ramps were fitted with the following function: f(V)=a(V−b)1 + expV−V1/2k, (2)with the following bounds: 0 ≤ a, 0 ≤ b (mV) ≤ 40, −40 ≤ V1/2 (mV) ≤0, and 0.1 ≤ k (mV) ≤ 20. The maximum current level for IMI was the minimum negative current obtained from the fit function, indicated by the open circle in Figure 2B.

The parameters we used (Fig. 2B, white markers) were as follows: for the potassium currents, the maximum current value at the voltage step to 20 mV; for Ih, the maximum amplitude during the step to −120 mV; for Isyn, the scaling factor of the sigmoid fit of the mean postsynaptic current at each voltage step; and for IMI, the maximum inward current of the fitted current. Since IMI is activated only in the presence of neuropeptides and is calculated as a difference current, we are unable to provide control data for this current.

Measurements of neuronal excitability

To measure the variability of excitability, we chemically isolated the LP neuron with 10−5 m PTX before running the stimulation protocols, which blocks inhibitory glutamatergic synapses from the pacemaker anterior burster (AB) and follower PY (pyloric) neurons to the LP neuron, as well as all synapses from the LP neuron to its targets (Bidaut, 1980). After 20 min of wash in, when IPSPs from LP to PD were effectively blocked, one or both PD neurons were hyperpolarized by −5 nA current injection or voltage clamped at a holding potential of −90 mV to silence the cholinergic synapse from PD to LP (Martinez et al., 2019; Fig. 1A,D).

To measure the excitability of the LP neuron in the form of frequency–current (f–I) curves, we injected increasing current levels from 0 to 5 nA in 0.5 nA increments, with 5 s duration, followed by the reversed protocol with decreasing current levels from 5 to 0 nA to check for hysteresis.

We measured the rebound properties of the LP neuron in two ways. To examine rebound properties following a long hyperpolarization period, we hyperpolarized the LP neuron with a −5 nA DC current for 10 s, followed by a 10 s interval of no current injection. This protocol was repeated for five sweeps and allowed us to measure the complete histogram of the burst structure of the LP neuron on rebound from hyperpolarization. To measure the steady-state rebound properties of the LP neuron as experienced during normal pyloric activity, we periodically hyperpolarized the LP neuron with 20 current pulses, −5 nA, 1 s on and 1 s off. We empirically determined that this 0.5 Hz periodic stimulation reliably resulted in rebound spiking in all modulatory conditions. Higher-frequency stimulations did not reliably produce spiking at steady state in the LP neuron in control or wash and were therefore not analyzed.

To compare the parameters of the f–I curves, we calculated the average instantaneous spike frequencies at each level of current injection and fit the experimental measurements with the following power function: f(I)=a(I−I0)b, (3)where a is a scaling factor, I0 is the current level that first elicited spikes, and the power b was set to 0 ≤ b ≤ 1 to limit the f–I curve to a sublinear function (Fig. 3Ai). To measure hysteresis of the f–I curve between depolarizing and repolarizing current injections (Fig. 3Aii), we calculated the average instantaneous spike frequency in the range of 2–4 nA and divided that value for the depolarizing current injections by that of the repolarizing current injections. This converted the hysteresis data to a ratio scale to use the coefficient of variation (CV) as a measure for variability. Ratios >1 mean that spike frequencies were greater during the current injections with depolarizing increments.

Figure 3.
Figure 3.

Neuromodulation changes excitability and reduces interindividual variability. A, Example of f – I relationships in four modulatory conditions (black: control; blue: Proc; magenta: Proc + CCAP; gray: wash). Left panel: Intracellular recording at four levels of increasing current injection in control condition. Right panel: Instantaneous spike frequencies at each current level in different neuromodulatory states, shown as the average with SD, were fitted with a power function. The left inset shows how the fit parameters influence the appearance of the curve. B, Distribution of fit parameters. Each dot represents the value from one experiment, and red lines indicate means. Application of one (blue, Proc) or two (magenta, Proc + CCAP) neuromodulators significantly changed parameters (asterisks). C, Application of one or two neuromodulators reduced the variability of the fit parameters compared with control (black). D, f – I curves showed hysteresis depending on in creasing (filled circles) or decreasing (open circles) levels of current injection. Only frequencies between 2 and 4 nA current injection (bold symbols) were used to calculate hysteresis as the ratio between increasing and decreasing current levels. E, Distribution of hysteresis across experiments. Application of one or two neuromodulators significantly changed hysteresis (asterisks). F, Application of one or two neuromodulators reduced the variability of hysteresis compared with control (black). Removing the outlier in control (i.e., the maximum value), for both hysteresis and I0, did not change the statistical significances. Control data for this figure are shown in Extended Data Figure 3-1. Raw data for B and E are provided in Extended Data Figure 3-2.

Figure 3-1

f–I relationships with sham applications of neuromodulators. One example experiment. A, Increasing current application (indicated by arrows). B, Decreasing current application (indicated by arrows). Download Figure 3-1, TIF file.

Figure 3-2

Raw data shown in Figure 3, B and E. Download Figure 3-2, XLS file.

To analyze rebound properties following a long hyperpolarizing DC current (Fig. 4), we considered all five sweeps because we did not observe transient effects across subsequent sweeps. For the rebound in response to periodic hyperpolarizing pulses (Fig. 5), we only considered the last 10 of the 20 pulses where the response of the LP neuron had reached steady state. In both cases, we calculated the latency to the first spike and the rebound spike structure for the selected sweeps (or pulses) in relation to the end of the hyperpolarizing current step. For the spike structure, we fitted the cumulative spike count histogram with the sigmoid function in Equation 1, in which a is now the maximum number of spikes per sweep or pulse and t1/2, instead of V1/2 , is the time of the midpoint relative to the end of the current injection. Additionally, for the steady-state rebound, we approximated the time constant by identifying the pulse number at which the latency dropped to 63% of its total drop value.

Figure 4.
Figure 4.

Neuromodulation changes general rebound properties and reduces interindividual variability. A, Spike raster and intracellular recording of the last of the five sweeps of one example experiment in four modulatory conditions (black: control; blue: Proc; magenta: Proc + CCAP; gray: wash). B, Latency was measured as the time from the end of the hyperpolarizing current injection to the first spike, averaged across all five sweeps. Sweeps of the intracellular recording are shaded in gray, from dark to light for sweeps 1–5. C, Spike histogram (Ci, 200 ms bin size) and sigmoid fit to the cumulative spike histogram (Cii). Dots indicate sigmoid midpoint. A–C are from the same experiment. D, Parameter distribution for latency and fit parameters; t1/2 is relative to the end of the current injection (A, t0). Dots represent individual experiments; the red line indicates the mean. Application of neuromodulators significantly (asterisks) changes most parameters (n.s.: ANOVA not significant). E, Variability of all parameters is reduced in the presence of neuromodulators. Control data for this figure are shown in Extended Data Figure 4-1. Raw data for D are provided in Extended Data Figure 4-2.

Figure 4-1

Rebound with sham applications of neuromodulators. One example experiment. A, Spike raster and corresponding intracellular recordings for all five sweeps. B, Spike histogram and sigmoid fit to the cumulative spike histogram. Dots indicate sigmoid midpoint. Download Figure 4-1, TIF file.

Figure 4-2

Raw data shown in Figure 4D. Download Figure 4-2, XLS file.

Figure 5.
Figure 5.

Neuromodulation changes steady-state rebound properties and reduces interindividual variability. A, Example intracellular recording of the LP neuron with 20 cycles of periodic inhibition. B, Spike raster and intracellular recording of all 20 sweeps of the same experiment in all four modulatory conditions (black: control; blue: Proc; magenta: Proc + CCAP; gray: wash). C, During the first sweeps the latency to the first spike successively decreased; therefore, we only included the last 10 sweeps (steady state) in the further analysis. The open circles indicate the sweep at which the latency was reduced to 63% of its total range. D, Spike histogram (Di, 100 ms bin size) and sigmoid fit to the cumulative spike histogram (Dii). Dots indicate sigmoid midpoint. A–D are from the same experiment. E, Parameter distribution for latency and fit parameters; t1/2 is relative to the end of current injection (B, t0). Dots represent individual experiments; the red line indicates the mean. Application of neuromodulators significantly changes parameters. F, Variability of all parameters is reduced in the presence of neuromodulators. Control data for this figure are shown in Extended Data Figure 5-1. Raw data for E are provided in Extended Data Figure 5-2.

Figure 5-1

Periodic rebound with sham applications of neuromodulators. One example experiment. A, Spike raster and corresponding intracellular recordings for all 20 sweeps. B, Spike histogram and sigmoid fit to the cumulative spike histogram. Dots indicate sigmoid midpoint. Download Figure 5-1, TIF file.

Figure 5-2

Raw data of the graph shown in Figure 5E. Download Figure 5-2, XLS file.

Data analysis

All data were analyzed using custom scripts written in MATLAB (releases R2018a and R2020b; MathWorks). Clampex files were imported with abfload (version 1.4.0.0; https://www.mathworks.com/matlabcentral/fileexchange/6190-abfload).

Statistical tests were performed with SigmaPlot (version 12.0; SyStat Software) or custom-written MATLAB scripts. Significance was assumed at α = 0.05. We performed one-way ANOVA when data passed the Shapiro–Wilk normality test and Levene’s equal variance test, or one-way ANOVA on ranks if at least one of the tests failed. We used Tukey’s test for post hoc multiple comparisons if the data had equal variance, and Dunn’s test if the variance was not equal. ANOVA results are listed in Tables 3, 4, and 5. We did not use paired tests because four of our single-cell variability experiments did not spike in control condition and were replaced with the control condition of four of the five sham experiments in which we never applied any modulators. Variability was calculated as CV (SD normalized to the absolute mean) and adjusted for the sample size (Haldane, 1955) if data were on a ratio scale, or as SD if data were on an interval scale. One exception is latency, where we also use SD as a variability metric. To calculate CV for latency, we must use the onset of the stimulus as a reference time point. Since we have stimulus durations of 10 s for most rebound experiments, the CVs in this case would be a scaled version of the SD.

The data presented are from 11 experiments for component (ionic current) variability, 16 experiments for single-cell (neuronal) variability, and 5 experiments with no modulator application.

LP model structure and activity

The LP neuron model was built as described in the study by Schneider et al. (2021), and implemented and run using NEURON (Carnevale and Hines, 2006) in Python. In brief, the model is composed of two coupled compartments, one representing the soma/neurite, and the other representing the axon. Spiking is generated in the axon, which has Ileak, INa, and IK, while the soma produced bursting activity with Ileak, ICa, IKCa, IA, Ih, and IMI. Calcium accumulation was tracked, where ICa contributed to intracellular calcium concentration ([Ca]in), and [Ca]in influenced both ICa and IKCa as described in Table 1. In addition to ICa, all currents were modeled as standard Hodgkin–Huxley type currents with the following general form: Ix=g¯xmxphxq(V−Ex), (4) dzdt=z∞(V)−zτz(V), (5)where x is the current type, g¯x is the maximal conductance, p and q are, respectively, the activation (m) and inactivation (h) variable (non-negative integer) exponents, and Ex is the reversal potential. Reversal potentials of axonal currents are Eleak = −55 mV, ENa = 70 mV, and EK = −80 mV. The reversal potentials of somatic currents are Eleak = −50 mV, EKCa = −80 mV, EA = −80 mV, Eh = −20 mV, and EMI = −10 mV. The change of activation and inactivation terms is given by Equation 5, where z∞ is the steady-state value of m or h and τz is the corresponding time constant. The parameters for the model are provided in Table 1. The calcium current was modeled using the Goldman–Hodgkin–Katz formalism, as follows: ICa=PCamCa3hCaFζ([Ca]oute−ζ−[Ca]ine−ζ−1), (6) ζ=zCavol⋅FRTV, (7)where PCa is the total permeability of the current, mCa and hCa are activation and inactivation variables given by Equation 5, vol is the volume of the microdomain influencing the current, F is Faraday’s constant, R is the universal gas constant, and T is temperature. [Ca] is the calcium concentration outside (out) and inside (in) the cell. The internal calcium concentration is given by the following: d[Ca]indt=[Ca]∞−[Ca]inτCa−P1zCaF⋅vol⋅PICa, (8)where ICa denotes the calcium current, [Ca]∞ denotes the steady-state calcium concentration inside the cell, P1 is the maximal per cluster permeability of ICa, and P is the total permeability over all clusters of interest. The parameters of ICa are P1 = 0.014 (cm · ms)−1, T = 283.15 K, [Ca]out = 13 mm, P1 = 1.1675 μm3/s, p = 0.0467 (cm · ms)−1, vol = 6.49 μm3, [Ca]∞ = 0.02 mm, and τCa = 25 ms.

View this table:
Table 1

Model parameters for the LP-like model

The difference between the LP model used here compared with the model in Schneider et al. (2021) is that, here, the kinetics of ICa, IA, and Ih were tuned to capture rebound firing, as seen in the biological data of this study. The protocol for rebound firing in the model is the same as in the biological experiments, where −5 nA is injected into the LP model soma for 10 s on and 10 s off. The f–I curves were generated by injecting 5 s current steps into the LP soma. The current step amplitudes ranged from −10 to 5.5 nA in increments of 0.5 nA. This larger range was chosen to better fit the f–I power function curve parameter I0. Cumulative spike histograms of model rebound experiments, and f–I power function curves of model f–I experiments were fit using scipy.optimize (Virtanen et al., 2020).

Generating a family of LP models

To address the questions of the effects of modulation on population variability, families of models were generated such that their activity summary statistics matched those of the experimental observations in the control condition (decentralized). This was done by generating a pool of candidate models using the simulation-based inference toolbox SBI (Gonçalves et al., 2020; Tejero-Cantero et al., 2020), and then selecting the final family of models according to criteria constrained by experimental results. For the rebound firing case, the summary statistics matched were the latency and power function coefficients (a, t1/2, k). The model parameters estimated were the somatic maximal conductances (gleak, g¯A , g¯h , g¯KCa , g¯MI ) and calcium permeability (PCa). The selection criteria imposed on the rebound firing population were that rebound latency must be >0.1 and <1.5 s. This resulted in a family of models (n = 198) whose summary statistics matched the experimentally measured firing rebound statistics in the decentralized control condition.

A similar approach was used to generate a family of LP models whose f–I curve fit parameters matched those of experimental data in the control (decentralized) condition (n = 85). In this case, we used the f–I power function fit parameters a, b, and I0 as summary statistics and estimated the same model parameters as in the rebound case, as well as of axonal maximal conductances (gleak, g¯Na , and g¯K ) and the inactivation time constant (τh) for INa. The selection criteria used were the exclusions of models that went into depolarization block in the range of additional g¯MI added, and models where the f–I power function fit parameter a was <6.

Population of linear integrate-and-fire models

To check whether our observations could be explained from first principles, we used a standard linear integrate-and-fire (LIF) model. The change in voltage with respect to time is given as follows: Cmdvdt=Iapp−gleak(v−Eleak)−gMI(v−EMI), (9)where a spike occurs when voltage exceeds vth, after which the voltage is immediately set to vreset.

Cm is the membrane capacitance, and Iapp is the amount of DC current injected. The ionic currents included were Ileak and IMI-L, a linearized version of IMI. Here, Ileak is a standard leak current with reversal potential Eleak = −60 mV and conductance gleak, and IMI-L is a positive leak current with reversal potential EMI = 10 mV and conductance gMI-L. All values used are given in Table 2.

View this table:
Table 2

Parameter values used for LIF simulations

At steady state, the frequency of the LIF neuron is given as follows: f=1τmlnvss−vresetvss−vth. (10)

The membrane time constant (τm) is as follows: τm=Cmgleak + gMI-L. (11)

Steady-state voltage (vss) is as follows: vss=Iapp + gleakEleak + gMI-LEMIgleak + gMI-L. (12)

A family of LIF model was constructed where gleak was sampled from a uniform distribution while all other parameters were kept fixed. To increase excitability, we added fixed amounts of g¯MI-L to the family of LIF models. The frequencies of the models were calculated using Equation 10 and fit to Equation 3, and the measures of variance were calculated in the same way as for the biological data.

Data availability

Source code for the LP model and linear integrate-and-fire model is available at https://github.com/fnadim/Schneider_et_al_variabilty_single_neuron and as Extended Data 1.

Extended Data 1

Supplementary Model code. Download Extended Data 1, ZIP file.

Results

Interindividual variability is not different between nonmodulated and modulated ionic currents

Most neurons express a multitude of ionic currents and receive synaptic input from many neurons. A modulatory input often activates or modifies the levels of a subset of these ionic currents. Substantial interindividual variability has been shown in identified pyloric neurons both for unmodulated voltage-gated currents (Schulz et al., 2006; Hamood and Marder, 2014), modulator-activated currents (Goaillard et al., 2009), and synaptic currents (Goaillard et al., 2009; Anwar et al., 2022). As a first step, we wanted to confirm the substantial interindividual variability of ionic currents found in the LP neuron and that neuropeptide neuromodulation does not systematically change component variability. We compared the variability of ionic currents that are targeted by proctolin and currents that are not, by measuring these currents in the same preparations. The nonmodulated currents that we measured included IHTK (transient and persistent portions), IA and Ih, and the modulated currents included Isyn and IMI (Fig. 2A,B; Golowasch and Marder, 1992b; Li et al., 2018). All currents except IMI were measured twice, once in control and once in Proc (Fig. 2C), whereas IMI was measured as a difference current between control and Proc (see Materials and Methods). For each current, the CV (SD divided by the absolute mean) was calculated across individuals (Fig. 2D). As expected, CV values were larger for the currents that had a small overall magnitude (Ih and IMI), likely because of increased measurement error, and therefore larger SD, for small currents. Because CV is the ratio of SD to the mean, changes in CV can be because of the following three mechanisms: (1) if the mean is constant but the SD is increasing, CV would increase; (2) if SD is constant but the mean increases, CV would decrease; and (3) if both mean and SD are changing, changes in CV would depend on the ratio changes of both. Generally, these measurements showed that CV values were similar between modulated and nonmodulated currents, and between control and Proc. Hence, the addition of Proc did not systematically reduce current variability, and the variability of modulated currents is in the same range as that of nonmodulated currents. Consequently, any change in variability of neuronal response properties caused by neuropeptide modulation is not simply because of altered component variability, but must arise from the interactions of ionic currents.

Neuromodulation reduces interindividual variability of the f–I relationship

The f–I relationship of a neuron describes the gain of its input–output function and is an established measure of neuronal excitability (Ermentrout, 1996; Skinner, 2013). We therefore measured the variability of the f–I relationship of the LP neuron across preparations (Fig. 3). We generated f–I curves by applying different levels of DC current to the LP neuron and measuring its firing frequency at each current level. To compare the f–I relationships across preparations, we used the f–I data from the increasing current steps and fit this curve with a power function for each modulatory condition (control, Proc, Proc + CCAP, wash; Fig. 3A). We fitted the f–I curves to obtain three parameters that capture the relationship and compared the fit parameters across all modulatory conditions (Fig. 3B).

Application of neuromodulators consistently altered the f–I curve of the LP neuron (Table 3, one-way ANOVA results). In contrast, sham neuromodulator application (see Materials and Methods) did not change the f–I relationship (Extended Data Fig. 3-1, example). Neuromodulators significantly increased the scaling factor a, indicating that the spike frequency of the LP neuron increased at any applied current level in the presence of neuromodulators. This is a well known phenomenon and is expected from modulators that activate inward currents (Heckmann et al., 2005). The exponent b decreased with modulators, indicating that the firing frequency started saturating at lower current levels, as higher values of b (closer to 1) indicate a more linear f–I relationship. Finally, modulator application lowered the firing threshold, as I0, the current level at which LP first spiked, decreased significantly, again an expected and well known phenomenon (Binder et al., 1993).

View this table:
Table 3

One-way ANOVA results for f–I experiments

To compare interindividual variability of the parameters under different modulatory conditions, we calculated the CV for parameter a, and the SD for parameters b and I0 (Fig. 3C; see Materials and Methods). For all parameters, variability was highest in control and decreased with modulator application. After the wash, variability of all parameters except for I0 increased. For all parameters, variability was reduced in the presence of modulators. However, the coapplication of Proc and CCAP did not consistently decrease variability more than Proc on its own. Since the SD for a was similar among the four conditions, the decrease in CV was because of the increased mean in Proc and Proc + CCAP compared with control and wash.

The f–I relationship with decreasing current steps was not identical to that with increasing steps (Fig. 3D), indicating hysteresis (Lee and Heckman, 1998). To quantify this hysteresis, we used the average instantaneous frequencies for current step values between 2 and 4 nA and measured the ratio of this measurement for increasing over decreasing steps (Fig. 3E, Table 3). This ratio was >1 in all conditions meaning that, at the same current level, the LP neuron firing frequency was always larger when increasing the applied current than decreasing it. With modulators, the ratio was closer to 1, indicating a reduction in hysteresis. As for the f–I curve fit parameters, modulator application reduced the interindividual variability of hysteresis, as seen in the smaller CV values (Fig. 3F). The decrease in CV together with the decrease in the respective means indicated a strong contribution of the reduction of SD. Raw data for these calculations is available in Extended Data Figure 3-2.

Neuromodulation reduces interindividual variability of rebound properties

In the intact pyloric network, the LP neuron generates bursts of spikes when it rebounds from inhibition by the pacemaker neurons. To assess rebound properties, we recorded the responses after release from a 10 s hyperpolarizing current step of −5 nA under different neuromodulatory conditions (Fig. 4A). We determined the latencies to the first rebound spike (Fig. 4B) and generated histograms of the rebound spike trains for the 10 s window after hyperpolarization (Fig. 4Ci). From the histograms, we fitted the cumulative spike counts with a sigmoid to extract parameters that can be used for comparisons of variability (Fig. 4Cii). Sham neuromodulator application (see Materials and Methods) did not change the rebound properties (Extended Data Fig. 4-1, example).

The LP neuron was usually silent in the absence of modulators. Upon hyperpolarization, we observed a slow voltage sag, which is known to be due to the activation of Ih (Fig. 4A; Golowasch and Marder, 1992a). Upon rebound, the LP neuron fired a train of spikes that terminated within the 10 s recording interval. In the presence of Proc or Proc + CCAP, however, LP was more depolarized and tonically firing without external current input. Upon hyperpolarization, it stopped spiking and produced a similar sag as in the control saline. In all conditions, the LP neuron produced a rebound after hyperpolarization and then its spike rate tapered off to its baseline level (Fig. 4A). The LP neuron produced its rebound spiking with a brief latency (Fig. 4B). In each condition, this latency to the first spike was similar for all sweeps. However, the application of Proc + CCAP significantly reduced the latency compared with control (Fig. 4D, Table 4). The spike histogram reflected that the spike rate first increased after the release from hyperpolarization and then tapered off (Fig. 4Ci). The total spike number, measured as the sigmoid fit parameter a, significantly increased in the presence of neuromodulators (Fig. 4D, Table 4), consistent with the increased excitability we observed in the f–I curves. The t1/2, relative to the release from hyperpolarization, occurred later with Proc than in control, but k, the steepness of the cumulative spike histogram, was not significantly changed with any of the modulators (Fig. 4D, Table 4). However, the variability of all parameters was always reduced in the presence of neuromodulators, indicated by lower values of SD and CV, regardless of any changes in the parameter (Fig. 4E). Variability of all parameters returned toward control levels in wash. For rebound latency, the decreased CV was because of a larger decrease in SD than the mean. For a, the SD increased, which was compensated for by the larger increase in the mean. For t1/2, the mean was increased and the SD was decreased, resulting in a smaller CV. Raw data for these calculations are available in Extended Data Figure 4-2.

View this table:
Table 4

One-way ANOVA results for general rebound experiments

View this table:
Table 5

One-way ANOVA results for steady-state rebound experiments

Neuromodulation reduces interindividual variability of periodic rebound bursts at steady state

In the intact circuit, the LP neuron generates bursts of action potentials on rebound from regular periodic inhibition by the pyloric pacemaker neurons. To mimic the response of the LP neuron to this periodic input, we hyperpolarized the neuron with 20 −5 nA current pulses, applied periodically with a 1 s on/1 s off protocol. The LP neuron responded to the first few pulses with fewer spikes and longer latencies before the responses stabilized to a steady state (Fig. 5A–C).

We analyzed rebound latencies and spike histograms only from the last 10 of the 20 cycles, in the same manner as shown in Figure 4. The steady-state rebound latency was smaller in Proc + CCAP than in control (Fig. 5E, Table 5). The only significant difference for spike count parameters was an increase in the total spike count, a, in the presence of modulators (Table 5), which is consistent with the results of Figure 4, but there was no modulator-induced change in t1/2 or k. As with our previous measurements, there was no difference for any of the parameters between Proc and Proc + CCAP application. Sham neuromodulator application (see Materials and Methods) did not change the rebound properties (Extended Data Fig. 5-1, example).

Although the application of neuropeptides did not always result in significant differences of the mean values, once again the variability for any of the measured parameters was always lower in the presence of modulators, and this variability increased after wash (Fig. 5F). For τ and t1/2, SD decreased more than the mean, which resulted in smaller CVs in Proc and Proc + CCAP. For a, the increase in mean compensated for the increase in SD so that CV decreased. Furthermore, the reduction of the CV of k was because of both an increase in mean and a decrease of SD. Raw data for these calculations is available as Extended Data Figure 5-2.

The influence of modulatory currents on rebound properties in a family of LP model neurons

As shown in the experimental results, peptide modulation is sufficient to reduce the interindividual variability of f–I curves and rebound properties at the single-cell level. In the STG, peptide modulators are known to activate an inward current, IMI, and to modify synaptic properties (Thirumalai et al., 2006; Zhao et al., 2011; Li et al., 2018). Generally, IMI increases the excitability of a neuron (Zhao et al., 2010). Given that the experiments were conducted on a synaptically isolated LP neuron, we used modeling to examine whether the activation of IMI was sufficient to decrease interindividual variability. To examine the contribution of IMI to the reduction of interindividual variability, we built a family of 198 LP model neurons (see Materials and Methods). In these models, the addition of peptides was modeled by increasing the level of g¯MI , the maximal conductance of IMI.

We tuned the kinetics of the LP model neurons to allow for rebound spiking as in the biological LP neuron (Fig. 6Ai). The models captured the increase and tapering of spike frequency on rebound from hyperpolarization and produced tonic spiking when high levels of IMI were activated. The population of LP model neurons was selected for similar rebound properties as measured in the biological data (Fig. 6Aii) and had a baseline normal distribution of low levels of g¯MI (Fig. 6Aiii, top). To test whether the interindividual variability of rebound decreases in response to modulation, we added fixed amounts of g¯MI (Δg¯MI ) to the population of models (Fig. 6Aiii, middle). The resulting fits of the cumulative spike histograms are shown in Figure 6B.

Figure 6.
Figure 6.

Increasing g¯MI in a family of LP models reduces the variability of rebound parameters. A, A family of 198 LP models was tuned to the rebound statistics from the biological experiments. Ai, Example voltage traces and spike histogram of one of the rebound LP models at different levels of added Δg¯MI . Aii, Paired plots of the rebound parameters from the biological experiments (black) and the family of LP models (orange). Aiii, The baseline g¯MI distribution (top) was shifted by either adding a fixed amount of g¯MI (Δg¯MI , middle) or a variable amount of g¯MI with a fixed mean (σg¯MI , bottom). B, Fits to the cumulative spike histograms (top row), the fit parameters and latency (middle row), and the corresponding measure of variability (bottom row) when adding increasing levels of g¯MI with a fixed distribution (Δg¯MI ). Individual dots represent values from an individual LP model, blue bars indicate the mean. An asterisk above a CV bar indicates that the CV for this group is significantly different from the CVs of all other groups. C, Fits to the cumulative spike histograms (top row), the fit parameters and latency (middle row), and the corresponding measure of variability (bottom row) when adding variable levels of g¯MI with a fixed mean (σg¯MI ). Individual dots represent values from an individual LP model, blue bars indicate the mean.

The metrics of variability were measured as in the biological experiments (see Materials and Methods). We observed a monotonic decrease of variability with increasing amounts of Δg¯MI added to the population. This suggests that increasing excitability by adding Δg¯MI is sufficient to decrease the interindividual variability of rebound firing.

The actions of modulators can be inherently variable in magnitude, as supported by reports of variable measurements of IMI across individuals and variable receptor expression (Garcia et al., 2015). To address whether variability in the activation of IMI would modify the effect seen in the case of adding fixed amounts of Δg¯MI , we allowed Δg¯MI to be a normal distribution with a moderate mean of Δg¯MI (0.025) and variances ranging from 0 to 0.1 (σg¯MI ; Fig. 6Aiii, bottom), resulting in the cumulative histogram fits shown in Figure 6C. We found that the effect of increasing σg¯MI increased the variability of rebound firing compared with adding fixed Δg¯MI . This suggests that the statistics of the g¯MI distribution activated will influence the effect observed on the interindividual variability of rebound firing, where, for a given mean value of g¯MI , a narrower distribution will allow for a greater reduction of variability than one with a larger variance.

Is the reduction of variability only due to an increase in excitability?

Generally, the increase of g¯MI decreased the variability of rebound parameters in our family of LP model neurons. However, f–I relationships capture a different aspect of neuron output than rebound properties. Therefore, we tuned a family of LP model neurons to have f–I statistics similar to those observed in the biological experiments and added fixed amounts of g¯MI (Δg¯MI ) to the baseline distribution, as we did for the model rebound experiments (Fig. 7Ai). We only included models that were reliably producing spikes at all Δg¯MI values (n = 85).

Figure 7.
Figure 7.

Increasing gMI in a family of LP models reduces the variability of f–I parameters. A, Schematic of the right shift of g¯MI to increase excitability (Ai) and f–I curves of a family of 85 LP models for selected Δg¯MI values (Aii). B, Fit parameters (top row) and corresponding variability measures (bottom row). Dots represent individual experiments; blue bars indicate the mean.

The most obvious change in the f–I curves was a left shift with increasing Δg¯MI (Fig. 7Aii). This shift simply reflects a lower threshold for excitability (i.e., a more negative I0). In addition, the means of scaling factor a increased and those of b decreased; Fig. 7B), indicating, respectively, an increase in the maximum spike rate and an overall reduction in gain. In contrast to the rebound parameters, the f–I parameters did not seem to approach a saturating value across the range of Δg¯MI values. However, we were unable to add larger amounts of Δg¯MI because models stopped generating spikes due to depolarization block. With increased Δg¯MI , variability decreased only for a and b, but not for I0 (Fig. 7B). Even so, the reduction in variability was far less than what we observed for the model rebound parameters.

Is the reduction of variability only due to an increase in excitability? To address this question from a first principles perspective, we used a family of LIF models (n = 500). Since the frequency for these models does not saturate, we truncated all f–I curves at 20 Hz to fit the power function to the same range of frequency values in each condition. The underlying variability within this family of models comes from the variability of the leak conductance. We expected that, like in the other models, increasing excitability in the LIF models (by increasing IMI-L levels) would result in a left shift of the f–I curves, but the other parameters would not be changed. As expected, I0 became more negative with increased excitability, and the SD remained constant (Fig. 8). However, the scaling factor (a) increased, and the curves became less linear (b decreased). The decrease in b corresponds to the widening of the base of the f–I curves with increasing Δg¯MI-L . The variability of these parameters was reduced in the same way as for the LP model family (Fig. 8). Thus, it appears that the reduction in variability due to increased excitability is a more generic property and is not limited to the LP model neuron.

Figure 8.
Figure 8.

Increasing excitability in a family of LIF models partially reduces the variability of f–I parameters. A, f–I curves for a family of 500 LIF models at three different Δg¯MI-L values. As in the LP models, increasing excitability shifts the curves to the left. B, Fit parameters (top row) and corresponding variability measures (bottom row). Dots represent individual experiments; blue bars indicate the mean.

Discussion

Neural circuit output can show variability across individual animals (Marder and Taylor, 2011; Wenning et al., 2018; Anwar et al., 2022; Gorur-Shandilya et al., 2022), but some attributes of this output must be constrained to provide biologically meaningful function characteristic for a given circuit state. For example, in the pyloric circuit, oscillation frequency can vary substantially across individuals under control conditions, while the relative timing and duty cycles of different neuron types is maintained (Bucher et al., 2005; Goaillard et al., 2009; Anwar et al., 2022). Such interindividual similarity of circuit output attributes has been carefully documented in central pattern-generating circuits, in which bursting neurons maintain a relatively constant phase in each oscillation cycle despite variations in cycle frequency (Grillner, 2006; Mullins et al., 2011; Zhang et al., 2014; LeGal et al., 2017; Martinez et al., 2019). The mechanisms that constrain aspects of circuit output to give rise to interindividual similarity are not well understood, particularly because ionic currents in identified neurons vary substantially across individuals (Schulz et al., 2006; Khorkova and Golowasch, 2007; Golowasch, 2014; Anwar et al., 2022). It is also not clear to what degree output similarity arises at the level of individual neurons or at the level of the whole circuit. Compounding this puzzle is the fact that circuit output is shaped by the actions of neuromodulators, which can also vary across individuals. Here we show that excitatory neuromodulation can increase the interindividual similarity of response properties in an isolated identified neuron. Further work is required to show whether such reduction of variability translates to the circuit output level, where differentially modulated neurons and synapses increase the complexity of neuromodulator actions (Harris-Warrick and Johnson, 2010; Johnson et al., 2011; Oleisky et al., 2020).

Variability of modulator-activated currents

Both peptides that we used in this study converge to activate the same ionic current, IMI, in the LP neuron (Swensen and Marder, 2000, 2001). In the same neuron type, ionic currents, including IMI, can greatly vary across animals (Schulz et al., 2006; Goaillard et al., 2009; Ransdell et al., 2013a,b), but it is possible that IMI has more consistent levels or that it is coordinated with other ionic currents, thereby promoting a similar neural activity. We therefore recorded both unmodulated currents and proctolin-activated currents and found similar levels of variability, indicating that interindividual output similarity is not because of consistent levels of the modulated current or simply because of a reduction in component variability.

Reduction in the variability of response properties must therefore be an emergent property (i.e., it must arise from the way that the currents activated by neuropeptides interact with other currents). Principally, there are two types of current interactions. First, as previous theoretical work has demonstrated, disparate current combinations across individuals can result in similar circuit output and voltage trajectories of individual neurons (Golowasch et al., 2002; Prinz et al., 2004). Indeed, despite substantial interindividual variability, pairs of ionic conductances are often correlated (McAnelly and Zakon, 2000; MacLean et al., 2003; Khorkova and Golowasch, 2007; Schulz et al., 2007; Cao and Oertel, 2011; Amendola et al., 2012; Temporal et al., 2012; Ransdell et al., 2013a; Tran et al., 2019). Such correlations can be independent of activity (MacLean et al., 2005; Schulz et al., 2006) and under neuromodulatory control (Khorkova and Golowasch, 2007; Temporal et al., 2012). Second, ionic currents within a cell have complex and nonlinear interactions through their voltage or Ca2+ dependence, as currents both shape and depend on voltage trajectories and Ca2+ fluctuations. In fact, IMI seems to be partially carried by Ca2+ ions and influence internal Ca2+ levels (Zhao et al., 2011; Gray et al., 2017; Schneider et al., 2021) and may therefore impact both voltage- and Ca2+-dependent currents. Notably, there are substantial differences in how sensitive different neuronal activity attributes are to the variability of different currents (Taylor et al., 2009), suggesting that not all variability has the same functional impact. In addition to their maximal conductances, voltage and Ca2+ dependence of ion channels can also covary, and their dependence on neuromodulators can affect this covariation, with potentially substantial consequences for neuronal activity (Amendola et al., 2012).

While we focused on currents that can be easily measured in the soma, the interindividual variability of fast axonal currents in the LP neuron is unknown. However, since currents, such as IMI, that arise in the neurites can influence neuronal excitability, they can also impact the spiking output that a neuron generates and thus alter activity phases or the number of spikes per burst.

Potential mechanisms for the reduction of variability by proctolin

Our computational models showed that the activation of IMI is sufficient to reduce the interindividual variability of LP activity. Notably, in a family of model LP neurons, adding moderate, but variable, levels of g¯MI increased output similarity compared with low levels of g¯MI with a narrower distribution. Thus, increasing levels of g¯MI can reduce the variability of neural activity even if the g¯MI levels are variable.

Increasing IMI increases excitability because IMI is a regenerative inward current (Zhao et al., 2010). Our models indicated two effects that contribute to reduced variability. (1) Firing frequency increased more in neurons with an initially low spiking activity, which promoted interindividual similarity. This is reminiscent of the finding that the activation of peptidergic modulatory neurons increases pyloric cycle frequency primarily when the baseline frequency is low (Nusbaum and Marder, 1989; Bartos and Nusbaum, 1997). The effect we describe here is independent of saturation, but rather reflects the larger slope and wider base of the f–I curve at low frequencies, as we demonstrated with LIF models that do not saturate (Fig. 8). (2) The other effect was saturation. Neurons have a maximal firing frequency, which is approached when they receive increasing excitation. If maximum firing rates are similar across individuals (i.e., a similar ceiling for firing rates), increasing excitability can reduce variability by reaching this ceiling. Maximum firing rates are constrained by the kinetics of ionic currents, in particular, those that influence the refractory period of the neuron. Similar kinetics would produce similar maximum firing rates. In contrast, highly variable ion channel kinetics across a population may lead to an increase, rather than a decrease, in output variability. Our family of model LP neurons in fact had near-identical ionic current kinetics and were only distinct in the maximal conductances of each current.

Thus, our models show two potential strategies for the system to reduce population variability: allowing for the output attributes to approach a ceiling or shifting these attributes away from low values (i.e., raising the floor).

Combined application of CCAP and proctolin does not additionally reduce variability

Neuropeptides mostly act through specific G-protein-coupled receptors (Brain and Cox, 2006; Huang and Thathiah, 2015), and in the STG, the subcellular pathways activated by the receptors of excitatory neuropeptides converge downstream to activate IMI (Swensen and Marder, 2001; Garcia et al., 2015). In the LP neuron, the combined activation of IMI by coapplication of CCAP and proctolin at concentrations >10−7 m is simply additive up to saturation (Li et al., 2018). We did not find an additional reduction of variability by coapplication of proctolin and CCAP compared with proctolin alone. Most likely, both 10−6 m proctolin and the combination of 5 × 10−7 m for each proctolin and CCAP were close to saturating (i.e., activated all IMI channels in the LP neuron so that the modulator-mediated reduction in variability was similar in both cases). Interestingly, the combined activation of IMI by the coapplication of CCAP and proctolin is sublinear when at least one of them is applied at a lower concentration (Li et al., 2018).

How would neuromodulation at nonsaturating concentrations affect output similarity? If two neuromodulators additively converge to activate a downstream target, but act through independently varying receptors, then their combined actions would result in less variability in the activation levels of that target. This is because adding two independently varying quantities has a smaller variability than either quantity, due to signal averaging. However, as mentioned above, the interaction of peptide modulators in activating IMI can be nonlinear (Li et al., 2018), which could result in a more complex interaction in how comodulation influences output variability.

Potential pitfalls in interpreting variability metrics

To our knowledge, few attempts have been made to quantify interindividual variability or output similarity (Wenning et al., 2018). SD is not dimensionless and therefore impossible to compare, for example, between spike number and current amplitude. Furthermore, SDs of the same unit cannot be compared when the means are very different. In contrast, the CV is dimensionless and scales the SD to the mean, which allows for the direct comparison of variability across different parameters. However, the CV is only valid for data on a ratio scale and therefore could not be used for all parameters. Furthermore, CV can be sensitive to parameters with small means: experimental error introduces some variability in any dataset that does not scale with the parameter of interest. In such cases, dividing SD by a small mean value results in a CV that is most likely an overestimation of the population variability. It is therefore possible that the interindividual variability of small ionic currents in the LP neuron, such as Ih and IMI, is influenced by measurement errors and therefore the biological variability is smaller than what we reported.

In our analysis, both the SD and the CV yield only a single value across animals, which makes statistical comparisons difficult. This is different from experiments in which intraindividual variability is compared between different conditions (Arancillo et al., 2015; Boele et al., 2018). In those experiments, a distribution of CVs exists for each condition (one CV per animal per condition), which can be compared with common statistical tests. In contrast, in our experiments, we quantified the interindividual variability with a single value. Unfortunately, available tests to compare the CV equality of single values require data to be normally distributed and only have sufficient power for CV values <0.5 (Sokal and Braumann, 1980; Feltz and Miller, 1996; Krishnamoorthy and Lee, 2014). The parameters and their CVs in our datasets did not always meet those requirements, which is why we refrained from using these tests. However, given that both measures of variability (CV and SD) consistently decreased in the presence of neuromodulators and usually recovered after washing, we are confident that neuromodulators did reduce the interindividual variability in the LP neuron.

Acknowledgments

Acknowledgments: We thank Elizabeth Cronin for sharing data and for helpful suggestions and comments. We also thank Horacio Rotstein for suggesting the use of reduced models to explore the reduction of variability. In addition, we thank Jorge Golowasch, Smita More-Potdar, Rodrigo Pena, and Astrid Prinz for helpful suggestions and comments.

Footnotes

  • The authors declare no competing financial interests.

  • This work was supported by National Institutes of Health Grant MH-060605 to F.N. and D.B.; and Deutsche Forschungsgemeinschaft Grant SCHN 1594/1-1 to A.C.S.

  • Received April 26, 2022.
  • Revision received May 27, 2022.
  • Accepted June 6, 2022.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.

References

  1. Amendola J, Woodhouse A, Martin-Eauclaire M-F, Goaillard J-M (2012) Ca2+/cAMP-sensitive covariation of IA and IH voltage dependences tunes rebound firing in dopaminergic neurons. J Neurosci 32:2166–2181. doi:10.1523/JNEUROSCI.5297-11.2012 pmid:22323729
  2. Anwar H, Martinez D, Bucher D, Nadim F (2022) Inter-animal variability in activity phase is constrained by synaptic dynamics in an oscillatory network. eNeuro 11:ENEURO.0027-22.2022. doi:10.1523/ENEURO.0027-22.2022
  3. Arancillo M, White JJ, Lin T, Stay TL, Sillitoe RV (2015) In vivo analysis of Purkinje cell firing properties during postnatal mouse development. J Neurophysiol 113:578–591. doi:10.1152/jn.00586.2014 pmid:25355961
  4. Bargmann CI (2012) Beyond the connectome: how neuromodulators shape neural circuits. BioEssays 34:458–465. doi:10.1002/bies.201100185 pmid:22396302
  5. Bartos M, Nusbaum MP (1997) Intercircuit control of motor pattern modulation by presynaptic inhibition. J Neurosci 17:2247–2256. pmid:9065486
  6. Bidaut M (1980) Pharmacological dissection of pyloric network of the lobster stomatogastric ganglion using picrotoxin. J Neurophysiol 44:1089–1101. doi:10.1152/jn.1980.44.6.1089 pmid:6256507
  7. Binder MD, Heckman CJ, Powers RK (1993) How different afferent inputs control motoneuron discharge and the output of the motoneuron pool. Curr Opin Neurobiol 3:1028–1034. doi:10.1016/0959-4388(93)90177-z pmid:8124071
  8. Boele H-J, Peter S, Ten Brinke MM, Verdonschot L, IJpelaar ACH, Rizopoulos D, Gao Z, Koekkoek SKE, De Zeeuw CI (2018) Impact of parallel fiber to Purkinje cell long-term depression is unmasked in absence of inhibitory input. Sci Adv 4:eaas9426. doi:10.1126/sciadv.aas9426 pmid:30306129
  9. Brain SD, Cox HM (2006) Neuropeptides and their receptors: innovative science providing novel therapeutic targets. Br J Pharmacol 147:S202–S211. doi:10.1038/sj.bjp.0706461
  10. Bucher D, Marder E (2013) SnapShot: neuromodulation. Cell 155:482–482.e1. doi:10.1016/j.cell.2013.09.047 pmid:24120143
  11. Bucher D, Prinz AA, Marder E (2005) Animal-to-animal variability in motor pattern production in adults and during growth. J Neurosci 25:1611–1619. doi:10.1523/JNEUROSCI.3679-04.2005 pmid:15716396
  12. Büschges A, Schmitz J, Bässler U (1995) Rhythmic patterns in the thoracic nerve cord of the stick insect induced by pilocarpine. J Exp Biol 198:435–456. doi:10.1242/jeb.198.2.435 pmid:9318078
  13. Buzsáki G, Wang X-J (2012) Mechanisms of gamma oscillations. Annu Rev Neurosci 35:203–225. doi:10.1146/annurev-neuro-062111-150444 pmid:22443509
  14. Cao X-J, Oertel D (2011) The magnitudes of hyperpolarization-activated and low-voltage-activated potassium currents co-vary in neurons of the ventral cochlear nucleus. J Neurophysiol 106:630–640. doi:10.1152/jn.00015.2010 pmid:21562186
  15. Carnevale NT, Hines ML (2006) The NEURON book. Cambridge: Cambridge UP.
  16. Ermentrout B (1996) Type I membranes, phase resetting curves, and synchrony. Neural Comput 8:979–1001. doi:10.1162/neco.1996.8.5.979 pmid:8697231
  17. Feltz CJ, Miller GE (1996) An asymptotic test for the equality of coefficients of variation from k populations. Stat Med 15:647–658. doi:10.1002/(SICI)1097-0258(19960330)15:6<647::AID-SIM184>3.0.CO;2-P
  18. Garcia VJ, Daur N, Temporal S, Schulz DJ, Bucher D (2015) Neuropeptide receptor transcript expression levels and magnitude of ionic current responses show cell type-specific differences in a small motor circuit. J Neurosci 35:6786–6800. doi:10.1523/JNEUROSCI.0171-15.2015 pmid:25926455
  19. Goaillard J-M, Taylor AL, Schulz DJ, Marder E (2009) Functional consequences of animal-to-animal variation in circuit parameters. Nat Neurosci 12:1424–1430. doi:10.1038/nn.2404 pmid:19838180
  20. Golowasch J (2014) Ionic current variability and functional stability in the nervous system. Bioscience 64:570–580. doi:10.1093/biosci/biu070 pmid:26069342
  21. Golowasch J, Marder E (1992a) Ionic currents of the lateral pyloric neuron of the stomatogastric ganglion of the crab. J Neurophysiol 67:318–331. doi:10.1152/jn.1992.67.2.318 pmid:1373762
  22. Golowasch J, Marder E (1992b) Proctolin activates an inward current whose voltage dependence is modified by extracellular Ca2+. J Neurosci 12:810–817. doi:10.1523/JNEUROSCI.12-03-00810.1992 pmid:1347561
  23. Golowasch J, Goldman MS, Abbott LF, Marder E (2002) Failure of averaging in the construction of a conductance-based neuron model. J Neurophysiol 87:1129–1131. doi:10.1152/jn.00412.2001 pmid:11826077
  24. Gonçalves PJ, Lueckmann J-M, Deistler M, Nonnenmacher M, Öcal K, Bassetto G, Chintaluri C, Podlaski WF, Haddad SA, Vogels TP, Greenberg DS, Macke JH (2020) Training deep neural density estimators to identify mechanistic models of neural dynamics. eLife 9:e56261. doi:10.7554/eLife.56261
  25. Gorur-Shandilya S, Cronin EM, Schneider AC, Haddad SA, Rosenbaum P, Bucher DM, Nadim F, Marder E (2022) Mapping circuit dynamics during function and dysfunction. eLife 11:e76579. doi:10.7554/eLife.76579
  26. Gray M, Daudelin DH, Golowasch J (2017) Activation mechanism of a neuromodulator-gated pacemaker ionic current. J Neurophysiol 118:595–609. doi:10.1152/jn.00743.2016 pmid:28446585
  27. Grillner S (2006) Biological pattern generation: the cellular and computational logic of networks in motion. Neuron 52:751–766. doi:10.1016/j.neuron.2006.11.008 pmid:17145498
  28. Haldane JBS (1955) The measurement of variation. Evolution 9:484–484. doi:10.2307/2405484
  29. Hamood AW, Marder E (2014) Animal-to-animal variability in neuromodulation and circuit function. Cold Spring Harb Symp Quant Biol 79:21–28. doi:10.1101/sqb.2014.79.024828 pmid:25876630
  30. Harris-Warrick RM, Johnson BR (2010) Checks and balances in neuromodulation. Front Behav Neurosci 4:47.
  31. Harris-Warrick RM, Marder E (1991) Modulation of neural networks for behavior. Annu Rev Neurosci 14:39–57. doi:10.1146/annurev.ne.14.030191.000351 pmid:2031576
  32. Heckmann CJ, Gorassini MA, Bennett DJ (2005) Persistent inward currents in motoneuron dendrites: implications for motor output. Muscle Nerve 31:135–156. doi:10.1002/mus.20261 pmid:15736297
  33. Howe MW, Atallah HE, McCool A, Gibson DJ, Graybiel AM (2011) Habit learning is associated with major shifts in frequencies of oscillatory activity and synchronized spike firing in striatum. Proc Natl Acad Sci U|S|A 108:16801–16806. doi:10.1073/pnas.1113158108 pmid:21949388
  34. Huang Y, Thathiah A (2015) Regulation of neuronal communication by G protein-coupled receptors. FEBS Lett 589:1607–1619. doi:10.1016/j.febslet.2015.05.007 pmid:25980603
  35. Hughes GM, Wiersma CAG (1960) The co-ordination of swimmeret movements in the crayfish, Procambarus clarkii (Girard). J Exp Biol 37:657–670. doi:10.1242/jeb.37.4.657
  36. Johnson BR, Brown JM, Kvarta MD, Lu JY, Schneider LR, Nadim F, Harris-Warrick RM (2011) Differential modulation of synaptic strength and timing regulate synaptic efficacy in a motor network. J Neurophysiol 105:293–304. doi:10.1152/jn.00809.2010 pmid:21047938
  37. Katz PS, Getting PA, Frost WN (1994) Dynamic neuromodulation of synaptic strength intrinsic to a central pattern generator circuit. Nature 367:729–731. doi:10.1038/367729a0 pmid:8107867
  38. Khorkova O, Golowasch J (2007) Neuromodulators, not activity, control coordinated expression of ionic currents. J Neurosci 27:8709–8718. doi:10.1523/JNEUROSCI.1274-07.2007 pmid:17687048
  39. Krishnamoorthy K, Lee M (2014) Improved tests for the equality of normal coefficients of variation. Comput Stat 29:215–232. doi:10.1007/s00180-013-0445-2
  40. Lee RH, Heckman CJ (1998) Bistability in spinal motoneurons in vivo: systematic variations in persistent inward currents. J Neurophysiol 80:583–593. doi:10.1152/jn.1998.80.2.583 pmid:9705452
    1. Hooper S,
    2. Büschges A
    LeGal J-P, Dubuc R, Smarandache-Wellmann CR (2017) Coordination of rhythmic movements. In: Neurobiology of motor control: fundamental concepts and new directions (Hooper S, Büschges A, eds), pp 305–340. New York: Wiley.
  42. Li X, Bucher D, Nadim F (2018) Distinct co-modulation rules of synapses and voltage-gated currents coordinate interactions of multiple neuromodulators. J Neurosci 38:8549–8562. doi:10.1523/JNEUROSCI.1117-18.2018 pmid:30126969
  43. MacLean JN, Zhang Y, Johnson BR, Harris-Warrick RM (2003) Activity-independent homeostasis in rhythmically active neurons. Neuron 37:109–120. doi:10.1016/s0896-6273(02)01104-2 pmid:12526777
  44. MacLean JN, Zhang Y, Goeritz ML, Casey R, Oliva R, Guckenheimer J, Harris-Warrick RM (2005) Activity-independent coregulation of IA and Ih in rhythmically active neurons. J Neurophysiol 94:3601–3617. doi:10.1152/jn.00281.2005 pmid:16049145
  45. Marder E (2012) Neuromodulation of neuronal circuits: back to the future. Neuron 76:1–11. doi:10.1016/j.neuron.2012.09.010 pmid:23040802
  46. Marder E, Bucher D (2001) Central pattern generators and the control of rhythmic movements. Curr Biol 11:R986–R996. doi:10.1016/S0960-9822(01)00581-4 pmid:11728329
  47. Marder E, Taylor AL (2011) Multiple models to capture the variability in biological neurons and networks. Nat Neurosci 14:133–138. doi:10.1038/nn.2735 pmid:21270780
    1. Kien J,
    2. McCrohan C,
    3. Winlow W
    Marder E, Weimann JM (1992) Modulatory control of multiple task processing in the stomatogastric nervous system. In: Neurobiology of motor programme selection: new approaches to mechanisms of behavioural choice (Kien J, McCrohan C, Winlow W, eds), pp 3–19. Oxford, UK: Pergamon.
  49. Marder E, Goeritz ML, Otopalik AG (2015) Robust circuit rhythms in small circuits arise from variable circuit components and mechanisms. Curr Opin Neurobiol 31:156–163. doi:10.1016/j.conb.2014.10.012 pmid:25460072
  50. Martinez D, Santin JM, Schulz D, Nadim F (2019) The differential contribution of pacemaker neurons to synaptic transmission in the pyloric network of the Jonah crab, Cancer borealis. J Neurophysiol 122:1623–1633. doi:10.1152/jn.00038.2019 pmid:31411938
  51. Masino MA, Fetcho JR (2005) Fictive swimming motor patterns in wild type and mutant larval zebrafish. J Neurophysiol 93:3177–3188. doi:10.1152/jn.01248.2004 pmid:15673549
  52. McAnelly ML, Zakon HH (2000) Coregulation of voltage-dependent kinetics of Na+ and K+ currents in electric organ. J Neurosci 20:3408–3414. doi:10.1523/JNEUROSCI.20-09-03408.2000
  53. Mullins OJ, Hackett JT, Buchanan JT, Friesen WO (2011) Neuronal control of swimming behavior: comparison of vertebrate and invertebrate model systems. Prog Neurobiol 93:244–269. doi:10.1016/j.pneurobio.2010.11.001 pmid:21093529
  54. Mulloney B, Harness PI, Hall WM (2006) Bursts of information: coordinating interneurons encode multiple parameters of a periodic motor pattern. J Neurophysiol 95:850–861. doi:10.1152/jn.00939.2005 pmid:16236775
  55. Nadim F, Bucher D (2014) Neuromodulation of neurons and synapses. Curr Opin Neurobiol 29:48–56. doi:10.1016/j.conb.2014.05.003 pmid:24907657
  56. Nusbaum MP, Marder E (1989) A modulatory proctolin-containing neuron (MPN). II. State-dependent modulation of rhythmic motor activity. J Neurosci 9:1600–1607. pmid:2566659
  57. Oleisky ER, Stanhope ME, Hull JJ, Christie AE, Dickinson PS (2020) Differential neuropeptide modulation of premotor and motor neurons in the lobster cardiac ganglion. J Neurophysiol 124:1241–1256. doi:10.1152/jn.00089.2020 pmid:32755328
  58. Parker D (2000) Presynaptic and interactive peptidergic modulation of reticulospinal synaptic inputs in the lamprey. J Neurophysiol 83:2497–2507. doi:10.1152/jn.2000.83.5.2497 pmid:10805651
  59. Pearson KG, Rossignol S (1991) Fictive motor patterns in chronic spinal cats. J Neurophysiol 66:1874–1887. doi:10.1152/jn.1991.66.6.1874 pmid:1812222
  60. Prinz AA, Bucher D, Marder E (2004) Similar network activity from disparate circuit parameters. Nat Neurosci 7:1345–1352. doi:10.1038/nn1352
  61. Ransdell JL, Nair SS, Schulz DJ (2013a) Neurons within the same network independently achieve conserved output by differentially balancing variable conductance magnitudes. J Neurosci 33:9950–9956. doi:10.1523/JNEUROSCI.1095-13.2013 pmid:23761890
  62. Ransdell JL, Temporal S, West NL, Leyrer ML, Schulz DJ (2013b) Characterization of inward currents and channels underlying burst activity in motoneurons of crab cardiac ganglion. J Neurophysiol 110:42–54. doi:10.1152/jn.00009.2013 pmid:23576706
  63. Schneider AC, Fox D, Itani O, Golowasch J, Bucher D, Nadim F (2021) Frequency-dependent action of neuromodulation. eNeuro 8:ENEURO.0338-21.2021. doi:10.1523/ENEURO.0338-21.2021
  64. Schulz DJ, Goaillard JM, Marder E (2006) Variable channel expression in identified single and electrically coupled neurons in different animals. Nat Neurosci 9:356–362. doi:10.1038/nn1639 pmid:16444270
  65. Schulz DJ, Goaillard JM, Marder EE (2007) Quantitative expression profiling of identified neurons reveals cell-specific constraints on highly variable levels of gene expression. Proc Natl Acad Sci U|S|A 104:13187–13191. doi:10.1073/pnas.0705827104 pmid:17652510
  66. Selverston AI, Moulins M (1985) Oscillatory neural networks. Annu Rev Physiol 47:29–48. doi:10.1146/annurev.ph.47.030185.000333 pmid:2986532
  67. Skinner FK (2013) Moving beyond type I and type II neuron types. F1000Res 2:19. doi:10.12688/f1000research.2-19.v1 pmid:24358846
  68. Sokal RR, Braumann CA (1980) Significance tests for coefficients of variation and variability profiles. Syst Biol 29:50–66. doi:10.2307/2412626
  69. Swensen AM, Marder E (2000) Multiple peptides converge to activate the same voltage-dependent current in a central pattern-generating circuit. J Neurosci 20:6752–6759. doi:10.1523/JNEUROSCI.20-18-06752.2000
  70. Swensen AM, Marder E (2001) Modulators with convergent cellular actions elicit distinct circuit outputs. J Neurosci 21:4050–4058. pmid:11356892
  71. Taylor AL, Goaillard JM, Marder E (2009) How multiple conductances determine electrophysiological properties in a multicompartment model. J Neurosci 29:5573–5586. doi:10.1523/JNEUROSCI.4438-08.2009 pmid:19403824
  72. Tejero-Cantero A, Boelts J, Deistler M, Lueckmann J-M, Durkan C, Gonçalves PJ, Greenberg DS, Macke JH (2020) sbi: a toolkit for simulation-based inference. J Open Source Softw 5:2505. doi:10.21105/joss.02505
  73. Temporal S, Desai M, Khorkova O, Varghese G, Dai A, Schulz DJ, Golowasch J (2012) Neuromodulation independently determines correlated channel expression and conductance levels in motor neurons of the stomatogastric ganglion. J Neurophysiol 107:718–727. doi:10.1152/jn.00622.2011 pmid:21994267
  74. Thirumalai V, Prinz AA, Johnson CD, Marder E (2006) Red pigment concentrating hormone strongly enhances the strength of the feedback to the pyloric rhythm oscillator but has little effect on pyloric rhythm period. J Neurophysiol 95:1762–1770. doi:10.1152/jn.00764.2005 pmid:16319213
  75. Tobin A-E, Cruz-Bermúdez ND, Marder E, Schulz DJ (2009) Correlations in ion channel mRNA in rhythmically active neurons. PLoS One 4:e6742. doi:10.1371/journal.pone.0006742 pmid:19707591
  76. Tran T, Unal CT, Severin D, Zaborszky L, Rotstein HG, Kirkwood A, Golowasch J (2019) Ionic current correlations are ubiquitous across phyla. Sci Rep 9:1687. doi:10.1038/s41598-018-38405-6 pmid:30737430
  77. Virtanen P, et al. (2020) SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods 17:261–272. doi:10.1038/s41592-019-0686-2 pmid:32015543
  78. Wang X-J (2010) Neurophysiological and computational principles of cortical rhythms in cognition. Physiol Rev 90:1195–1268. doi:10.1152/physrev.00035.2008 pmid:20664082
  79. Wenning A, Norris BJ, Günay C, Kueh D, Calabrese RL (2018) Output variability across animals and levels in a motor system. eLife 7:e31123. doi:10.7554/eLife.31123
  80. Zhang C, Guy RD, Mulloney B, Zhang Q, Lewis TJ (2014) Neural mechanism of optimal limb coordination in crustacean swimming. Proc Natl Acad Sci U|S|A 111:13840–13845. doi:10.1073/pnas.1323208111 pmid:25201976
  81. Zhao S, Golowasch J, Nadim F (2010) Pacemaker neuron and network oscillations depend on a neuromodulator-regulated linear current. Front Behav Neurosci 4:21. doi:10.3389/fnbeh.2010.00021 pmid:20514340
  82. Zhao S, Sheibanie AF, Oh M, Rabbah P, Nadim F (2011) Peptide neuromodulation of synaptic dynamics in an oscillatory network. J Neurosci 31:13991–14004. doi:10.1523/JNEUROSCI.3624-11.2011 pmid:21957260

Synthesis

Reviewing Editor: Arianna Maffei, Stony Brook University

Decisions are customarily a result of the Reviewing Editor and the peer reviewers coming together and discussing their recommendations until a consensus is reached. When revisions are invited, a fact-based synthesis statement explaining their decision and outlining what is needed to prepare a revision will be listed below. The following reviewer(s) agreed to reveal their identity: Ronald Calabrese.

There is agreement between the reviewers that this is overall a well designed study providing an interesting contribution to the field.

The reviewers raised a few concerns about some of the conclusions and provided suggestions to improve clarity.

Additional analysis as suggested by reviewer #2, would be helpful in further supporting the authors’ main points.

The full comments by the reviewers follow below.

Reviewer #1

This is a systematic study of how neuromodulation affects output variability of individual neuronal firing properties using identified and synaptically isolated LP neurons of the stomatogastric pyloric CPG of crabs. The work is carefully executed and statistical approaches are sound and clearly rationalized and explicated. There is clearly a reduction in output variability (by various rigorous measures comparing CV and SD as is appropriate to the variable) with neuromodulator application observed in experiments and this finding is substantiated in an ensemble modeling approach with realistic and leaky integrate and fire neuronal models. The primary conclusion is that increasing excitability by raising baseline firing accounts for the most to the decrease in variability with modulator action. The writing is very clear and the figure very accessible though complicated. The data is convincing.

Minor Suggestions for Revision

1. I was honestly confused by the Abstract’s ending. "With a family of computational LP-like models, we explored how increased excitability and target variability contribute to output similarity and found two mechanisms. (sic) saturation of the responses and a differential increase in baseline activity. Saturation alone can reduce the interindividual variability only if the population shares a similar ceiling for the responses. In contrast, reduction of variability due to the increase in baseline activity is independent of ceiling effect." I could not follow the arguments in the Discussion about saturation, and indeed I did not see the argument framed as similarity of ceiling effects. I suggest you stick to the main point and not emphasize the saturation angle as I don’t see this explored well.

2. I kept wondering about the separation of the spiking mechanism both in time scale and spatially from the measured currents. Do you think the fast spiking system of the axon is less variable than the slow system of the cell body and primary neurite? As far as I can see axonal fast INa has not been measured in LP neurons (maybe there is a fast IK in the axon area?). You measure currents in the slow system and your output measures are largely the fast system. Can you comment on this in discussion?

3. Lines 474-476: Make clear here in the text that excitability in the LIF models is changed by changing IMI.

Reviewer #2

In this article, the authors investigate the changes in variability of neuronal output associated with neuromodulation in a crustacean central pattern generator motorneuron (LP neuron). After evaluating the variability present at the biophysical level by measuring several ion currents (including the neuropeptide-activated current at the center of the study) in voltage-clamp, the authors use different protocols to characterize the neuronal output, including fI curves and more ecological rhythmic rebound stimuli. The authors analyze the variability of each parameter (by calculating either the standard deviation or the coefficient of variation) under control and neuromodulated conditions in order to determine the changes in variability under neuromodulation. Overall, most electrophysiological parameters display reduced variability under neuromodulation. The authors then use a database modeling approach to demonstrate that this reduction in variability is observed even when the neuromodulator-activated current is highly variable from neuron to neuron. Overall, the study is well performed, both at the electrophysiological and modeling level, the data are analyzed correctly and the manuscript is well written. The authors are also conscious of the limits of the statistical analysis.

MAJOR POINTS

As the authors underline it, comparing variability between conditions is to the least a complex question, because variability is quantified by a single parameter (SD or CV in their case), and thus very little statistical argument can be provided to support the conclusions (is a decrease in CV value really a significant decrease?). Notwithstanding this difficulty, and considering the CV as the main factor that has been used by the authors, it would be important to distinguish between several cases that can lead to a decrease in CV under neuromodulated condition. Since the CV is a ratio between SD and the mean, a decrease in CV can be basically a consequence of three different situations: i) a constant SD associated with an increase in the mean ii) a decrease in SD associated with a constant mean and iii) a decrease in the mean associated with a "stronger" decrease in SD. Since the authors use in some cases the SD instead of the CV (particularly when data are not on a ratio scale), I believe distinguishing between these cases is of the utmost importance. I would suggest that, whenever CV is used to estimate variability, the authors comment on the sources of reduction in variability (is the SD decreased? The mean increased? Do both the SD and the mean change in value?). Perhaps considering the range of values (instead of the CV) could also be of interest for some of the electrophysiological parameters analyzed? I would expect the authors to refine a bit their analysis or discussion of the results to distinguish between the sources of the reduction in variability. This supplementary analysis may also help to better define what may be the biological relevance of these different changes in variability under neuromodulation.

Author Response

We thank the reviewers for their helpful comments to improve the manuscript. Our responses are highlighted in red.

Reviewer #1

1. I was honestly confused by the Abstract’s ending. "With a family of computational LP-like models, we explored how increased excitability and target variability contribute to output similarity and found two mechanisms. (sic) saturation of the responses and a differential increase in baseline activity. Saturation alone can reduce the interindividual variability only if the population shares a similar ceiling for the responses. In contrast, reduction of variability due to the increase in baseline activity is independent of ceiling effect." I could not follow the arguments in the Discussion about saturation, and indeed I did not see the argument framed as similarity of ceiling effects. I suggest you stick to the main point and not emphasize the saturation angle as I don’t see this explored well.

The period in the abstract should have been a colon; we corrected this in line 14. We reformatted and extended the corresponding section 6.2 in the Discussion to more clearly delineate which data and arguments belong to the discussion about increasing baseline activity or to the discussion about ceiling effects.

2. I kept wondering about the separation of the spiking mechanism both in time scale and spatially from the measured currents. Do you think the fast spiking system of the axon is less variable than the slow system of the cell body and primary neurite? As far as I can see axonal fast INa has not been measured in LP neurons (maybe there is a fast IK in the axon area?). You measure currents in the slow system and your output measures are largely the fast system. Can you comment on this in discussion?

We added this point to the discussion at the end of section 6.1.

3. Lines 474-476: Make clear here in the text that excitability in the LIF models is changed by changing IMI.

Done.

Reviewer #2

As the authors underline it, comparing variability between conditions is to the least a complex question, because variability is quantified by a single parameter (SD or CV in their case), and thus very little statistical argument can be provided to support the conclusions (is a decrease in CV value really a significant decrease?). Notwithstanding this difficulty, and considering the CV as the main factor that has been used by the authors, it would be important to distinguish between several cases that can lead to a decrease in CV under neuromodulated condition. Since the CV is a ratio between SD and the mean, a decrease in CV can be basically a consequence of three different situations: i) a constant SD associated with an increase in the mean ii) a decrease in SD associated with a constant mean and iii) a decrease in the mean associated with a "stronger" decrease in SD. Since the authors use in some cases the SD instead of the CV (particularly when data are not on a ratio scale), I believe distinguishing between these cases is of the utmost importance. I would suggest that, whenever CV is used to estimate variability, the authors comment on the sources of reduction in variability (is the SD decreased? The mean increased? Do both the SD and the mean change in value?). Perhaps considering the range of values (instead of the CV) could also be of interest for some of the electrophysiological parameters analyzed? I would expect the authors to refine a bit their analysis or discussion of the results to distinguish between the sources of the reduction in variability. This supplementary analysis may also help to better define what may be the biological relevance of these different changes in variability under neuromodulation.

All our figures have panels that show, for each parameter, the distribution of data from individual experiments and the mean across individuals, e.g., Fig 2C. For clarification, we added brief descriptions of the main factor (change in SD, change in mean, or change in both) that contributes to the decrease in CV for the biological experiments. Additionally, we now first list the potential mechanisms for decreasing CV in lines 314-317, and include the raw data for Figures 3-5 as Extended Data.

We also realized that contrary to our descriptions in the Methods and Results, we were showing CV instead of SD for the rebound latencies in Figures 4 and 5. The figures have been corrected and show the SDs now.