## 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.

## 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 MgCl_{2}, 13 CaCl_{2}, 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 K_{2}SO_{4} + 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Ω.

#### 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 (*I*_{HTK}; Khorkova and Golowasch, 2007), the transient potassium current (*I*_{A}), the hyperpolarization-activated inward current (*I*_{h}), the synaptic current from the PD to LP synapse (*I*_{syn}), and the modulator-activated inward current (*I*_{MI}; Golowasch and Marder, 1992b; Fig. 2*A*). For these measurements, the PD neuron was always held at −60 mV, except during the *I*_{syn} measurements. For the total potassium currents (*I*_{K}), the LP neuron was held at −80 mV and stepped for 500 ms from −50 to 20 mV in 10 mV increments. *I*_{HTK} was measured in the same way, except that LP was held at −40 mV to inactivate *I*_{A}. *I*_{A} was calculated as the difference current between *I*_{K} and *I*_{HTK}. To measure *I*_{h}, we stepped LP from −40 to −120 mV for 10 s. For *I*_{syn}, 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 *I*_{MI}, 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).

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. 2*B*, open circles) for all further calculations. *I*_{A} was calculated as the difference in current between *I*_{K} (Fig. 2*A*, gray traces) and *I*_{HTK}. For *I*_{HTK}, we measured the transient and persistent currents separately, indicated by the arrows “t” and “p” in Figure 2*A*. For *I*_{h}, we calculated the current amplitude between the beginning of the hyperpolarizing step and the end, indicated by the arrow in Figure 2*B*. For *I*_{syn}, 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:
*a* is the maximum postsynaptic current, *V*_{1/2} is the voltage at midpoint, and *k* is the slope factor at *V*_{1/2}. Lower bounds were set at −10 nA for *a* and −80 mV for *V*_{1/2}. We used the value of *a* as the maximum synaptic current, indicated by the open circle in Figure 2*B*.

For *I*_{MI} , 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, *I*_{MI-T} (Schneider et al., 2021). The *I–V* curves of the negative ramps were fitted with the following function:
*a*, 0 ≤ *b* (mV) ≤ 40, −40 ≤ *V*_{1/2} (mV) ≤0, and 0.1 ≤ *k* (mV) ≤ 20. The maximum current level for *I*_{MI} was the minimum negative current obtained from the fit function, indicated by the open circle in Figure 2*B*.

The parameters we used (Fig. 2*B*, white markers) were as follows: for the potassium currents, the maximum current value at the voltage step to 20 mV; for *I*_{h}, the maximum amplitude during the step to −120 mV; for *I*_{syn}, the scaling factor of the sigmoid fit of the mean postsynaptic current at each voltage step; and for *I*_{MI}, the maximum inward current of the fitted current. Since *I*_{MI} 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. 1*A*,*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:
*a* is a scaling factor, *I*_{0} 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. 3*Ai*). To measure hysteresis of the *f–I* curve between depolarizing and repolarizing current injections (Fig. 3*Aii*), 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-1

*f–I* relationships with sham applications of neuromodulators. One example experiment. ** A**, Increasing current application (indicated by arrows).

**, Decreasing current application (indicated by arrows). Download Figure 3-1, TIF file.**

*B*### 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 *t*_{1/2}, instead of *V*_{1/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-1

Rebound with sham applications of neuromodulators. One example experiment. ** A**, Spike raster and corresponding intracellular recordings for all five sweeps.

**, Spike histogram and sigmoid fit to the cumulative spike histogram. Dots indicate sigmoid midpoint. Download Figure 4-1, TIF file.**

*B*### Figure 4-2

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

### Figure 5-1

Periodic rebound with sham applications of neuromodulators. One example experiment. ** A**, Spike raster and corresponding intracellular recordings for all 20 sweeps.

**, Spike histogram and sigmoid fit to the cumulative spike histogram. Dots indicate sigmoid midpoint. Download Figure 5-1, TIF file.**

*B*### Figure 5-2

Raw data of the graph shown in Figure 5*E*. 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 *I*_{leak}, *I*_{Na}, and *I*_{K}, while the soma produced bursting activity with *I*_{leak}, *I*_{Ca}, *I*_{KCa}, *I*_{A}, *I*_{h}, and *I*_{MI}. Calcium accumulation was tracked, where *I*_{Ca} contributed to intracellular calcium concentration ([Ca]_{in}), and [Ca]_{in} influenced both *I*_{Ca} and *I*_{KCa} as described in Table 1. In addition to *I*_{Ca}, all currents were modeled as standard Hodgkin–Huxley type currents with the following general form:
*x* is the current type, *p* and *q* are, respectively, the activation (*m*) and inactivation (*h*) variable (non-negative integer) exponents, and *E _{x}* is the reversal potential. Reversal potentials of axonal currents are

*E*

_{leak}= −55 mV,

*E*

_{Na}= 70 mV, and

*E*

_{K}= −80 mV. The reversal potentials of somatic currents are

*E*

_{leak}= −50 mV,

*E*

_{KCa}= −80 mV,

*E*

_{A}= −80 mV,

*E*

_{h}= −20 mV, and

*E*

_{MI}= −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:

*P*

_{Ca}is the total permeability of the current,

*m*

_{Ca}and

*h*

_{Ca}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:

*I*

_{Ca}denotes the calcium current, [Ca]

_{∞}denotes the steady-state calcium concentration inside the cell,

*P*

_{1}is the maximal per cluster permeability of

*I*

_{Ca}, and

*P*is the total permeability over all clusters of interest. The parameters of

*I*

_{Ca}are

*P*

_{1}= 0.014 (cm · ms)

^{−1},

*T*=

*283.15 K, [*

*Ca*]

_{out}= 13 mm,

*P*

_{1}= 1.1675 μm

^{3}/s,

*p*=

*0.0467 (cm · ms)*

^{−1}, vol

*=*

*6.49 μm*

^{3}, [Ca]

_{∞}= 0.02 mm, and τ

_{Ca}= 25 ms.

The difference between the LP model used here compared with the model in Schneider et al. (2021) is that, here, the kinetics of *I*_{Ca}, *I*_{A}, and *I*_{h} 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 *I*_{0}. 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*, *t*_{1/2}, *k*). The model parameters estimated were the somatic maximal conductances (*g*_{leak}, *P*_{Ca}). 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 *I*_{0} as summary statistics and estimated the same model parameters as in the rebound case, as well as of axonal maximal conductances (*g*_{leak}, *τ*_{h}) for *I*_{Na}. The selection criteria used were the exclusions of models that went into depolarization block in the range of additional *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:
*v*_{th}, after which the voltage is immediately set to *v*_{reset}.

*C*_{m} is the membrane capacitance, and *I*_{app} is the amount of DC current injected. The ionic currents included were *I*_{leak} and *I*_{MI-L}, a linearized version of *I*_{MI}. Here, *I*_{leak} is a standard leak current with reversal potential *E*_{leak} = −60 mV and conductance *g*_{leak}, and *I*_{MI-L} is a positive leak current with reversal potential *E*_{MI} = 10 mV and conductance *g*_{MI-L}. All values used are given in Table 2.

At steady state, the frequency of the LIF neuron is given as follows:

The membrane time constant (τ_{m}) is as follows:

Steady-state voltage (*v*_{ss}) is as follows:

A family of LIF model was constructed where *g*_{leak} was sampled from a uniform distribution while all other parameters were kept fixed. To increase excitability, we added fixed amounts of

### 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 *I*_{HTK} (transient and persistent portions), *I*_{A} and *I*_{h}, and the modulated currents included *I*_{syn} and *I*_{MI} (Fig. 2*A*,*B*; Golowasch and Marder, 1992b; Li et al., 2018). All currents except *I*_{MI} were measured twice, once in control and once in Proc (Fig. 2*C*), whereas *I*_{MI} 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. 2*D*). As expected, CV values were larger for the currents that had a small overall magnitude (*I*_{h} and *I*_{MI}), 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. 3*A*). We fitted the *f–I* curves to obtain three parameters that capture the relationship and compared the fit parameters across all modulatory conditions (Fig. 3*B*).

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 *I*_{0}, the current level at which LP first spiked, decreased significantly, again an expected and well known phenomenon (Binder et al., 1993).

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 *I*_{0} (Fig. 3*C*; 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 *I*_{0} 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. 3*D*), 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. 3*E*, 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. 3*F*). 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. 4*A*). We determined the latencies to the first rebound spike (Fig. 4*B*) and generated histograms of the rebound spike trains for the 10 s window after hyperpolarization (Fig. 4*Ci*). From the histograms, we fitted the cumulative spike counts with a sigmoid to extract parameters that can be used for comparisons of variability (Fig. 4*Cii*). 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 *I*_{h} (Fig. 4*A*; 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. 4*A*). The LP neuron produced its rebound spiking with a brief latency (Fig. 4*B*). 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. 4*D*, Table 4). The spike histogram reflected that the spike rate first increased after the release from hyperpolarization and then tapered off (Fig. 4*Ci*). The total spike number, measured as the sigmoid fit parameter *a*, significantly increased in the presence of neuromodulators (Fig. 4*D*, Table 4), consistent with the increased excitability we observed in the *f–I* curves. The *t*_{1/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. 4*D*, 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. 4*E*). 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 *t*_{1/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.

### 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. 5*A–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. 5*E*, 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 *t*_{1/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. 5*F*). For τ and *t*_{1/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, *I*_{MI}, and to modify synaptic properties (Thirumalai et al., 2006; Zhao et al., 2011; Li et al., 2018). Generally, *I*_{MI} 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 *I*_{MI} was sufficient to decrease interindividual variability. To examine the contribution of *I*_{MI} 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 *I*_{MI}.

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

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 Δ

The actions of modulators can be inherently variable in magnitude, as supported by reports of variable measurements of *I*_{MI} across individuals and variable receptor expression (Garcia et al., 2015). To address whether variability in the activation of *I*_{MI} would modify the effect seen in the case of adding fixed amounts of Δ*Aiii*, bottom), resulting in the cumulative histogram fits shown in Figure 6*C*. We found that the effect of increasing σ

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

Generally, the increase of *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 *Ai*). We only included models that were reliably producing spikes at all Δ*n* = 85).

The most obvious change in the *f–I* curves was a left shift with increasing Δ*Aii*). This shift simply reflects a lower threshold for excitability (i.e., a more negative *I*_{0}). In addition, the means of scaling factor *a* increased and those of *b* decreased; Fig. 7*B*), 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 Δ*a* and *b*, but not for *I*_{0} (Fig. 7*B*). 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 *I*_{MI-L} levels) would result in a left shift of the *f–I* curves, but the other parameters would not be changed. As expected, *I*_{0} 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 Δ

## 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, *I*_{MI}, in the LP neuron (Swensen and Marder, 2000, 2001). In the same neuron type, ionic currents, including *I*_{MI}, can greatly vary across animals (Schulz et al., 2006; Goaillard et al., 2009; Ransdell et al., 2013a,b), but it is possible that *I*_{MI} 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 Ca^{2+} dependence, as currents both shape and depend on voltage trajectories and Ca^{2+} fluctuations. In fact, *I*_{MI} seems to be partially carried by Ca^{2+} ions and influence internal Ca^{2+} levels (Zhao et al., 2011; Gray et al., 2017; Schneider et al., 2021) and may therefore impact both voltage- and Ca^{2+}-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 Ca^{2+} 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 *I*_{MI}, 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 *I*_{MI} is sufficient to reduce the interindividual variability of LP activity. Notably, in a family of model LP neurons, adding moderate, but variable, levels of

Increasing *I*_{MI} increases excitability because *I*_{MI} 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 *I*_{MI} (Swensen and Marder, 2001; Garcia et al., 2015). In the LP neuron, the combined activation of *I*_{MI} 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 *I*_{MI} channels in the LP neuron so that the modulator-mediated reduction in variability was similar in both cases). Interestingly, the combined activation of *I*_{MI} 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 *I*_{MI} 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 *I*_{h} and *I*_{MI}, 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.

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.