Neuromodulation Reduces Interindividual Variability of Neuronal Output

Visual Abstract


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;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(MacLean et al., , 2005Khorkova 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 . 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.

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 1 5 Â 10 À7 M CCAP (henceforth referred to as Proc 1 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 1 20 mM KCl (resistance: 20-25 MV). 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 1 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 MV.

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. 2A). 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 120 mV (Fig. 2B, open circles) for all further calculations. I A was calculated as the difference in current between I K ( Fig. 2A, 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 2A. For I h , we calculated the current amplitude between the beginning of the hyperpolarizing step and the end, indicated by the arrow in Figure 2B. 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: where 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 2B.
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: with the following bounds: 0 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 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 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 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. 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. Higherfrequency 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: where 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. 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.
S(x) denotes the logistic sigmoid function 1/(1 1 exp(-x)). Time constants are in milliseconds. 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.
Statistical tests were performed with SigmaPlot (version 12.0; SyStat Software) or custom-written MATLAB scripts. Significance was assumed at a = 0.05. We performed oneway 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: 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 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 1 is the steady-state value of m or h and t 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: where 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: where I Ca denotes the calcium current, [Ca] 1 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 mm 3 /s, p = 0.0467 (cm · ms) À1 , vol = 6.49 mm 3 , [Ca] 1 = 0.02 mM, and t 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 , g A , g h , g KCa , g MI ) and calcium permeability (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 , g Na , and g K ) and the inactivation time constant (t h ) for I Na . 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 integrateand-fire (LIF) model. The change in voltage with respect to time is given as follows: where a spike occurs when voltage exceeds 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 (t 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 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 integrateand-fire model is available at https://github.com/fnadim/ Schneider_et_al_variabilty_single_neuron and as Extended Data 1.

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 , and synaptic currents 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.  2A,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. 2C), 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. 2D). 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 1 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, 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 I HTK [transient (t) portion indicated by black arrow, persistent (p) by gray arrow], I A (with total potassium currents I K as gray overlay; I A is the difference current between all potassium currents I K and I HTK ); I h ; the synaptic current I syn from PD to LP; and the neuromodulator-activated, voltage-gated current I MI 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 I HTK , the amplitude of the current at the end of the hyperpolarizing voltage step for I h ; the scaling factor for the sigmoid fit for I syn ; and the maximum inward current of the current fit for I MI . 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 I HTK is shown in bold colors, persistent in transparent colors. Since I MI is calculated as a difference current (Procctrl), there are no data in control. D, Coefficients of variation (SD/mean) are in the same range for the nonmodulated currents (I HTK , I A , I h ) and modulated currents (I syn , I MI ).
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. 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 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 1 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 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 I 0 , did not change the statistical significances. Control data for this figure are shown in Extended Data 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 I h  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 1 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 1 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 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. 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 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. 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 1 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 t 1/2 or k. As with our previous measurements, there was no difference for any of the parameters between Proc and Proc 1 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 t and t 1/2 , SD decreased more than the mean, which resulted in smaller CVs in Proc and Proc 1 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 "N" is the number of animals, results from tests for normality and equal variance is given as pass/fail, the test statistics are F for ANOVA and H for ANOVA on ranks, and "res" is the residuals. Significant p-values (p 0.05) are printed in bold. Values for ctrl, Proc, Proc 1 CCAP, and wash are means for ANOVA, and medians for ANOVA on ranks. Latency was bounded by experimental design between 0 and 10; therefore, we logarithmically transformed the data to a normal distribution and did the statistics on the log-transformed data. "N" is the number of animals, results from tests for normality and equal variance is given as pass/fail, the test statistics are F for ANOVA and H for ANOVA on ranks, and "res" is the residuals. Significant p-values (p 0.05) are printed in bold. Values for ctrl, Proc, Proc 1 CCAP, and wash are means for ANOVA, and medians for ANOVA on ranks.
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 g MI , the maximal conductance of I MI .
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 I MI 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 (D g MI ) to the population of models (Fig. 6Aiii, middle). The resulting fits of the cumulative spike histograms are shown in Figure 6B.
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 D g MI added to the population. This suggests that increasing excitability by adding D 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 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 D g MI , we allowed D g MI to be a normal distribution with a moderate mean of D g MI (0.025) and variances ranging from 0 to 0.1 (s g MI ; Fig. 6Aiii, bottom), resulting in the cumulative histogram fits shown in Figure 6C. We found that the effect of increasing s g MI increased the variability of rebound firing compared with adding fixed D 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 (D 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 D g MI values (n = 85).
The most obvious change in the f-I curves was a left shift with increasing D g MI (Fig. 7Aii). 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. 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 D g MI values. However, we were unable to add larger amounts of D g MI because models stopped generating spikes due to depolarization block. With increased D g MI , variability decreased only for a and b, but not for I 0 (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 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 D g MIÀL . The "N" is the number of animals, results from tests for normality and equal variance are given as pass/fail, the test statistics are F for ANOVA and H for ANOVA on ranks, and "res" is the residuals. Significant p-values (p 0.05) are printed in bold. Values for ctrl, Proc, Proc 1 CCAP, and wash are means for ANOVA, and medians for ANOVA on ranks. Latency was bounded between 0 and 1 by experimental design. Therefore, we transformed the data to a normal distribution by multiplying by p and calculating the arctangent. 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 D 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 (D g MI , middle) or a variable amount of g MI with a fixed mean (s 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 (D 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 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.

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 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.
continued 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 (s g MI ). Individual dots represent values from an individual LP model, blue bars indicate the mean. 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 21 dependence, as currents both shape and depend on voltage trajectories and Ca 21 fluctuations. In fact, I MI seems to be partially carried by Ca 21 ions and influence internal Ca 21 levels (Zhao et al., 2011;Gray et al., 2017;Schneider et al., 2021) and may therefore impact both voltage-and Ca 21 -dependent currents. Notably, there are substantial differences in how sensitive different neuronal activity attributes are to the variability of different currents , suggesting that not all variability has the same functional impact. In addition to their maximal conductances, voltage and Ca 21 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 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 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-proteincoupled 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.