Mechanisms of generation of membrane potential resonance in a neuron with multiple resonant ionic currents

Neuronal membrane potential resonance (MPR) is associated with subthreshold and network oscillations. A number of voltage-gated ionic currents can contribute to the generation or amplification of MPR, but how the interaction of these currents with linear currents contributes to MPR is not well understood. We explored this in the pacemaker PD neurons of the crab pyloric network. The PD neuron MPR is sensitive to blockers of H- (IH) and calcium-currents (ICa). We used the impedance profile of the biological PD neuron, measured in voltage clamp, to constrain parameter values of a conductance-based model using a genetic algorithm and obtained many optimal parameter combinations. Unlike most cases of MPR, in these optimal models, the values of resonant- (fres) and phasonant- (fϕ = 0) frequencies were almost identical. Taking advantage of this fact, we linked the peak phase of ionic currents to their amplitude, in order to provide a mechanistic explanation the dependence of MPR on the ICa gating variable time constants. Additionally, we found that distinct pairwise correlations between ICa parameters contributed to the maintenance of fres and resonance power (QZ). Measurements of the PD neuron MPR at more hyperpolarized voltages resulted in a reduction of fres but no change in QZ. Constraining the optimal models using these data unmasked a positive correlation between the maximal conductances of IH and ICa. Thus, although IH is not necessary for MPR in this neuron type, it contributes indirectly by constraining the parameters of ICa.

Neuronal membrane potential resonance (MPR) is associated with subthreshold and network oscillations. A number of voltage-gated ionic currents can contribute to the generation or amplification of MPR, but how the interaction of these currents with linear currents contributes to MPR is not well understood. We explored this in the pacemaker PD neurons of the crab pyloric network. The PD neuron MPR is sensitive to blockers of H-(I H ) and calcium-currents (I Ca ). We used the impedance profile of the biological PD neuron, measured in voltage clamp, to constrain parameter values of a conductance-based model using a genetic algorithm and obtained many optimal parameter combinations. Unlike most cases of MPR, in these optimal models, the values of resonant-(f res ) and phasonant-(f ϕ = 0 ) frequencies were almost identical. Taking advantage of this fact, we linked the peak phase of ionic currents to their amplitude, in order to provide a mechanistic explanation the dependence of MPR on the I Ca gating variable time constants. Additionally, we found that distinct pairwise correlations between I Ca parameters contributed to the maintenance of f res and resonance power (Q Z ). Measurements of the PD neuron MPR at more hyperpolarized voltages resulted in a reduction of f res but no change in Q Z . Constraining the optimal models using these data unmasked a positive correlation between the maximal conductances of I H and I Ca . Thus, although I H is not necessary for MPR in this neuron type, it contributes indirectly by constraining the parameters of I Ca .

Author summary
Many neuron types exhibit membrane potential resonance (MPR) in which the neuron produces the largest response to oscillatory input at some preferred (resonant) frequency and, in many systems, the network frequency is correlated with neuronal MPR. MPR is captured by a peak in the impedance vs. frequency curve (Z-profile), which is shaped by a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Introduction
Neuronal network oscillations at characteristic frequency bands emerge from the coordinated activity of the participating neurons. Membrane potential resonance (MPR) is defined as the ability of neurons to exhibit a peak in their voltage response to oscillatory current inputs at a preferred or resonant frequency (f res ) [1]. MPR has been observed in many neuron types such as those in the hippocampus [2][3][4] and entorhinal cortex [2][3][4][5][6], inferior olive [7,8], thalamus [1,9], striatum [10,11], as well as in invertebrate oscillatory networks such as the pyloric network of the crustacean stomatogastric ganglion (STG) [12][13][14]. Neurons may also exhibit phasonance or a zero-phase response, which describes their ability to synchronize with oscillatory inputs at a preferred phasonant frequency (f ϕ = 0 ) [4,[15][16][17][18]. Resonance, phasonance and intrinsic oscillations are related, but are different phenomena as one or more of them may be present in the absence of the others [15,16,18].
Resonant and phasonant frequencies result from a combination of low-and high-pass filter mechanisms produced by the interplay of the neuron's passive properties and one or more ionic currents and their interaction with the oscillatory inputs [1,15,18,19]. The slow resonant currents (or currents having resonant gating variables) oppose voltage changes and act as high-pass filters. They include the hyperpolarization-activated inward current (I H ) and the slow outward potassium current (I M ). On the other hand, the fast amplifying currents (or currents having amplifying gating variables) favor voltage changes and can make MPR more pronounced. They include the persistent sodium current (I NaP ) and the inward rectifying potassium (I Kir ) current. Most previous systematic mechanistic studies have primarily examined models with one resonant and one amplifying current, such as I H and I NaP , respectively [15,[18][19][20]. Currents having both activating and inactivating gating variables (in a multiplicative way) such as the low-threshold calcium current (I Ca ) are not included in this classification, but they are able to produce resonance by mechanisms that are less understood [16,21].
Although a causal relationship between the properties of MPR and network activity has not been established [but see 22], resonant neurons have been implicated in the generation of network oscillations in a given frequency band because the resonant and network frequencies often match up or are correlated. One example is in the hippocampal theta oscillations [23] in which CA1 pyramidal cells exhibit MPR in vitro at theta frequencies of 4-10 Hz [2][3][4]24] (but see [25]). Interestingly, MPR is not constant across the somatodendritic arbor in these neurons profiles. The maximal conductances of ionic currents of neurons in the stomatogastric nervous system vary widely [35][36][37]. We therefore assume that the parameters that determine the Zprofile in the PD neuron vary across animals. Thus, instead of searching for a single model that fit the PD neuron Z-profile, we used a genetic algorithm to capture a collection of parameter sets that fit this Z-profile. To achieve such a fit, we defined a set of ten attributes that characterize the PD neuron Z-profile (e.g., resonant frequency and amplitude) and used a multiobjective evolutionary algorithm [MOEA,38] to obtain a family of models that fit these attributes. We then used this family of optimal models to identify the important biophysical parameters and relationships among these parameters to explain how the PD neuron Z-profile is shaped. We show how the fact that the inactivating calcium current peaks at the same phase as the passive properties, in response to sinusoidal inputs, can explain why resonant and phasonant frequencies are equal. We identify significant pairwise parameter-correlations, which selectively set certain attributes of MPR. We show that, in this neuron, I H does not produce MPR but can extend the dynamic range of I Ca parameters mediating MPR. Furthermore, we identify a subset of models that capture the experimental shift in the resonant frequency with changes in lower bound of voltage oscillation. Finally, we exploit the fact that the resonant and phasonant frequencies are equal for the PD neuron to provide a mechanistic understanding of the effects of the I Ca time constants on the resonant frequency by using phase information. Our results provide a mechanistic understanding for a generic class of neurons that exhibit both resonance and phasonance as the result of the interaction between multiplicative gating variables and complement the studies in [16].

Results
The PD neuron produces 1 Hz bursting oscillations with a slow-wave span of approximately -60mV to -30mV (Fig 1a). Driving the neuron through this voltage range with a ZAP function in voltage clamp (Fig 1b top panel) produces a minimum (arrow in Fig 1b bottom panel) in the amplitude of the current response (Fig 1b). The input frequency at which this minimum occurs corresponds to the resonant frequency: a peak in the Z-profile (f res , Z max ; Fig 1c1). The value of f res was 0.86 ± 0.05Hz producing Z max values of 10.23 ± 0.51 MΩ (N = 18; Fig 1d). The φ-profile shows a phasonant frequency f φ = 0 = 0.81 ± 0.05Hz, which in most cases matched f res (Fig 1c2). The PD neuron had a Q Z of 2.77 ± 0.71 MΩ and Λ ½ of 0.53 ± 0.04 Hz. Across preparations, Q Z showed considerable variability, whereas f res , Λ ½ , and f φ = 0 were relatively consistent (Fig 1d). The corresponding median values for f res , Q Z, Λ ½ , and f φ = 0 were 0.83 Hz, 2.77 MO, 0.5 Hz, 0.79 Hz, respectively.
To obtain model parameter combinations constrained by the PD neuron Z-and φ -profiles, we generated a population of models using an NSGA-II algorithm (see Methods). The attributes of a single PD neuron Z-and φ -profiles (Fig 2, filled red circles) constrained the optimization of the parameter values. This resulted in a population of~9000 sets of parameters ("optimal" dataset). All models in the optimal dataset captured the attributes of Z and φ to within 5% of the target (light blue lines in Fig 2), with the exception of φ max , which may be due to the anatomical structure of the PD neuron, a property that is omitted in our single-compartment model, or due to additional ionic currents, such as the potassium A current, which are not included in our model [16,39].

The generation of MPR by the interaction of two resonant voltage-gated currents
To understand how Z is generated by the dynamics of individual ionic currents at different voltages and frequencies, we examined the amplitude and kinetics of ionic currents. In voltage clamp, Z is shaped by active voltage-gated currents, interacting with the passive leak and capacitive currents, in response to the voltage inputs. To understand the contribution of different ionic currents, we measured these currents in response to a constant frequency sine wave voltage inputs (Fig 3a inset) at three frequencies: 0.1Hz, 1Hz (f res ) and 4Hz (Fig 3). For Membrane potential resonance MPR of the PD neuron was measured in voltage clamp. a. During ongoing activity, the PD neuron shows a slow-wave voltage waveform ranging approximately between -60 and -30 mV. b. The membrane potential (V zap ) and the injected current (I PD ) were recorded when the PD neuron was voltage-clamped using a ZAP function between -60 and -30mV and sweeping frequencies between 0.1 and 4 Hz. The arrowhead indicates resonance, where the current amplitude is minimal and Z is maximal. c. The impedance amplitude Z(f) (c1) and phase φ(f) (c2) profiles of the PD neuron recorded in 18 preparations. The cross bars show the mean and SEM of f res and Z max (c1) and f φ = 0 (c2). The shaded region indicates the 95% confidence interval. d. The range of three Z(f) attributes f res , Q Z , and Λ 1/2 and one φ(f) attribute f φ = 0 . Each attribute was normalized to the median of its distribution for cross comparison. CoV is the coefficient of variation.  The Z(f) (a) and φ(f) (b) profiles of 500 randomly selected models from the optimal dataset (light blue curves) are compared to the target neuron's impedance profiles (red circles). All attributes (except φ max ) were captured to within 5% accuracy. The values of the biological target impedance amplitude attributes (in Hz, MΩ) were: (f 0 , Z 0 ) = (0. 1, 8.2), (f res , Z max ) = (1, 13.7), (0.4, 11.65), (2.5, 11.65) and (4, 9.6). The target impedance phase attributes (in Hz, rad) were: (0.1, 0), (f φmax , φ max ) = (0.4, 0.5), (f φ = 0 , 0) = (1.05, 0), (2, -4), (f φmin , φ min ) = (4, -0.4). these frequencies, we plotted the steady-state current as a function of voltage (Fig 3b-3d left) and normalized time (or cycle phase = time x frequency; Fig 3b-3d right). At 0.1 Hz, the amplitudes of I H and I L +I Cm sets I total at low (~-60 mV) and high (~-30 mV) voltages, respectively (Fig 3b left). Since I H deactivation is slow, it also contributes to I total at high voltages (Fig 3b right). At 1 Hz (= f res ), I H still sets the minimum of the total current, but, because of its slow kinetics, its steady-state dynamics are mostly linear (Fig 3c left). However, now I Ca peaks in phase (Fig 3c right) with the passive I L + I Cm at high voltages, thus producing a smaller I total (magenta bar in Fig 3c). The values of I H at 4 Hz are not much different from 1 Hz (Fig 3d). However, I Ca peaks at a much later phase (Fig 3d right), which does not allow it to compensate for I L + I Cm at high voltages, thus resulting in a larger I total (magenta bar in Fig 3d). Note that at 1 Hz, the total current peaks at a cycle phase close to 0.5, thus implying that that the f res and f φ = 0 are very close or equal (Fig 3c right). Although Fig 3 shows the results for only one model in the optimal dataset, these results remain nearly identical for all models in the optimal dataset. The standard deviation of the currents measured, including the total current was never above 0.15 nA over all models. The inset in Fig 3c shows one standard deviation around the mean for the data shown in the right panel, calculated for 200 randomly selected models.
An important collective property of the models we found is that the two frequencies, f res and f φ = 0 coincide (Fig 4a and 4b). We analyzed the experimental data, and confirmed that the coincidence of MPR and phasonance frequencies also occurs in the biological system (Fig 4b inset). This is typically not the case for neuronal models (and for dynamical systems in general), not even for linear systems [18][19][20], with the exception of the harmonic oscillator. However, the fact that it occurs in this system, allows us to use the current vs. cycle phase (current-phase) diagrams to understand the dependence of f res and f φ = 0 on the model parameters (Fig 4c).
The current-phase diagrams are depicted as in Fig 3b-3d, as graphs of I total , I L and I Ca as a function of the cycle phase for each given specific input frequency (Fig 4c). We do not show I H and I Cm in this plot, because at frequencies near f res they do not change much with input frequency. Note that I L is independent of the input frequency (five panels in Fig 4c) because it precisely tracks the input voltage.
In voltage clamp, f φ = 0 = 1Hz is where I total is at its minimum amplitude exactly at cycle phase 0.5, coinciding with the peak of the input voltage (Fig 4c, middle). The fact that I L precisely tracks the input voltage imposes a constraint on the shapes of I Ca and I total . Therefore, by necessity, if the I Ca trough occurs for a cycle phase below 0.5, the I total peak must occur for a cycle phase above 0.5 (Fig 4c, top two panels) and vice versa (Fig 4c, bottom two panels). This is shown by the slope of the line joining the peaks of I total and I Ca and, at f res this line is approximately vertical (Fig 4c middle panel).
We use this tool to explain the dependence of the Z-profile on the time constants t Ca m (Fig 5a) and t Ca h (Fig 5b). The corresponding current-phase diagrams are presented in Fig 5c  and 5d, respectively. In each panel we present the current-phase diagrams for f at 1 Hz (= f res when the parameter is at 100%; middle) and f = f res (sides) when f res is different from 1Hz.
To understand the dependence of Z on changes in t Ca m and t Ca h we have to primarily explain the dependence of the two attributes Z max and f res on these parameters. While f res has a similar monotonic dependence on t Ca m and t Ca h (as these parameters increase, f res decreases), Z max has the opposite dependence on t Ca m and t Ca h . The opposite dependence of Z max on t Ca m and t Ca h is a straightforward consequence of the opposite feedback effects (positive for t Ca m and negative for t Ca h ) that these parameters exert on I Ca . An increase in t Ca m (for fixed values of t Ca h ) results in a smaller I Ca in response to a given voltage clamp input. Because I Ca is smaller and negative, this leads to an increase in I total and a decrease in Z at all frequencies. Similarly, an increase in t Ca h (for fixed values of t Ca m ) results in a larger I Ca , leading to a decrease in I total and an increase in Z. For a fixed value of the input frequency f (e.g. f = 1 Hz as in Fig 5), for Z max to decrease as t Ca m increases (Fig 5a), the cycle phase of peak I Ca is delayed thereby subtracting less from I L on  Interaction of resonant currents and membrane potential resonance the depolarizing phase. This leads to I total to phase advance relative to I L (Fig 5c) and causes f res to decrease. Similarly, for Z max to increase as t Ca h increases (Fig 5b), I Ca has to peak later in the cycle thereby subtracting less from I L on the depolarizing phase, which causes I total to peak earlier in the cycle, which in turn causes the bar also to swing from the left to the right (Fig 5d). Therefore, f res decreases.

Parameter constraints and pairwise correlations
Previous studies have shown that stable network output can be produced by widely variable ion channel and synaptic parameters [37,40]. Our biological data, similarly, showed that many of the Z-and φ-profile attributes, such as f res , Λ ½ and f φ = 0 are relatively stable across different PD neurons whereas Q Z shows the most variability (Fig 1d). To determine whether the Z-and φ-profile attributes constrain ionic current parameters, we examined the variability of the model parameters in the optimal dataset. We found that some parameters were more constrained while others were widely variable, as measured by the coefficient of variation (CoV ;   Fig 6a). Parameters showing large CoVs were " ; those showing small CoVs were " g L and the time constant of activation of I H and I Ca and half-activation voltage of , " g L (in increasing order of CoV value). A small CoV value implies that the parameter is tightly constrained in order to produce the proper Z-and φ-profiles.
A number of studies have indicated that the large variability in ion channel parameters is counter-balanced by paired linear covariation of these parameters [36,37,[41][42][43]. Considering the large variability, we identified parameter pairs that co-varied (Fig 6b). For this, we carried out a permutation test for the Pearson's correlation coefficients, followed by a Student's t-test on the regression slopes, to identify significant correlations between pairs of parameters (see Methods). The strongest correlations were between the following parameters: " Fig 6b).
1=2 was fixed at -70 mV, using data from experimental measurements in crab [44] (see Methods). However, we also repeated the MOEA with V H m 1 1=2 set to -96 mV, as reported in lobster experiments [45], and found that all correlations observed with the former g H in the optimal models without qualitatively changing the other findings.
In particular, we found that the " g Ca À V Ca h 1 1=2 correlation appeared nonlinear, but there were strong and distinct linear correlations in the two regions " g Ca > 0.05 μS (low " g Ca ) and " g Ca < 0.05 μS (high " g Ca ; Fig 6c). To ensure that our partitioning of the population into different levels of " g Ca was valid, we ran the MOEA two additional times, each time using only the mean , and t Ca m for either the low or the high " g Ca values. These optimal models consistently separated into two non-overlapping model parameters, consistent with the low and high " g Ca models in Fig 6c. We examined if the low and high " g Ca models separated or showed distinct correlations in the remaining parameters. The two groups produced non-overlapping subsets of model graph. We calculated the Pearson's correlation coefficient for each pair of parameters in the low and high " g Ca groups and tested for significance as before (see Table 1). We found that only the high " g Ca group showed a significant t Ca m À t Ca h and " g H À t Ca h correlations (Table 1). Additionally, both low and high " g Ca groups showed the follow- Furthermore, when we ran the MOEA on models where " g H was set to 0, the only optimal models obtained fell within a narrow range of the high " g Ca group (S2 Fig), which is consistent with the distribution of high " g Ca models in the " g H À " g Ca panel of Fig 6d.  Fig 6. The optimal models show variability in individual and pairs of parameters. a. The range of parameters for all optimal models (~9000). Each parameter is normalized by its median value for cross comparison. The median values were . Three representative optimal model parameter sets are shown (cyan, orange, purple solid line segments) indicating that widely different parameter combinations can produce the biological Z(f) and φ(f). CoV is coefficient of variation. b. Pairwise relationships among parameters of all optimal models (black dots). The range of parameter space was sampled within the prescribed limits given to the optimization routine, shown by including the sampled non-optimal models (grey). Permutation test showed significant pairwise correlations (green highlighted boxes with linear fits shown as green lines). c. Optimal models could be separated into two highly significant linear fits (green lines) in g Ca À V Ca h 1 1=2 according to whether g Ca < 0.05 (red; Low g Ca ) or g Ca > 0.05 (cyan; High g Ca ). d. All pairwise relationships, separated on the low or high g Ca (colors as in panel c). Green boxes are the same as in b. https://doi.org/10.1371/journal.pcbi.1005565.g006 Decreasing the lower bound of voltage oscillations influences the measured f res and Z max The lower voltage range of the PD bursting oscillation is strongly influenced by the inhibitory synaptic input from the lateral pyloric neuron (LP), and previous work has shown that f res in the PD neuron is influenced by the minimum of the voltage oscillation (V low ) [14]. In order to explore which subset of our optimal models faithfully reproduce the influence of the minimum voltage range, we measured the Z-profile when V low was changed from -60 to -70 mV (Fig 7a).
To explore whether the shift in f res as a function of V low could be captured by either low or high " g Ca models, we measured the shift in f res and Z max , when V low was changed from -60mV to -70mV. We found that f res decreased by 0.24±0.03 Hz and Z max increased by 5.2±0.6 MO for high " g Ca models, whereas f res decreased by 0.07±0.02Hz and Z max decreased by 2.6±0.2MO for low " g Ca models (Fig 7b, right panel). Therefore, neither model group reproduced the experimental changes in the Z-profile, specifically, a decrease in f res and no change in Z max .
We consequently filtered the full optimal dataset (black dots Fig 7c) to find a subset of models that reproduced the change in f res and Z max (to within 5% of the representative experimental Z(f) shown in Fig 7a) when V low was decreased to -70mV. Of the~9000 models in the population, we found~1000 models that produced the desired change. Interestingly, the resulting models showed a trade-off in values for " g Ca and V Ca h 1 1=2 parameters that showed little overlap with the low and high " g Ca model groups (Fig 7c). To understand why this particular group (which we will term intermediate " g Ca ) produced small changes in Z max when V low was decreased, we plotted the current-voltage relationships for I Ca , I H , I Ca +I H and I total for V low = -60 and -70 mV, measured at f = 1Hz (f res at V low = -60mV) and compared these models with the low and high " g Ca models. For V low = -60mV, the ionic currents behaved similarly for all model groups and I total was maximal at -30mV (magenta curve in Fig 7d1-7d3), indicating the similarity of all models in the optimal dataset. However, when V low was at -70mV revealed differences in peak I Ca , without affecting the peak amplitude of I H across the different " g Ca groups (Fig 7e1-7e3). The differences in peak I Ca accounted for most of the changes in I total across the different " g Ca groups. The Z max values for intermediate " g Ca models reproduced the small shift seen in experiments because I Ca was at the correct level at high voltages (-30 mV) when V low was at -70mV (Fig 7e3). The other two groups did not produce appropriate Z max for V low = -70mV because either I Ca was too small  (and hence I total too large), resulting in a smaller Z max (Fig 7e1) or vice versa (Fig 7e2). It was also clear that the more negative voltages allowed for an increase in I H levels and therefore larger contribution to the total current. With V low at -70mV, not only was there a larger peak amplitude of I H at the lower voltages, but the current at positive voltages also increased because of the very slow deactivation rate. Consequently, I H did not fully turn off when I Ca peaked, so that it also contributes to shaping the upper envelope of the total current. I H kinetics were different across the groups (Fig 7e1-7e3). Taken together with the fact that when I H was removed produced only parameter values with very high " g Ca and very low V Ca h 1 1=2 (S1 Fig), these data suggest that I H could extend the range of I Ca parameters over which MPR could be generated through compensation for variable levels of I H .
The I Ca in low " g Ca models was too small when V low was -70 mV, because the low conductance did not allow for a significant contribution from the additional de-inactivation (considering the higher V Ca h 1 1=2 in this group) and therefore the peak current did not increase enough. Consequently, the contribution of I H at low voltages was greater than that of I Ca at higher voltages (Fig 7e2). Conversely, in the high " g Ca group, V Ca h 1 1=2 was more negative and so many more channels were available for de-inactivation and the contribution of I Ca at higher voltages was much larger than that of I H at low voltages (Fig 7e3). These findings suggest that the balance between these two currents, that shape the lower and upper envelope of the total current response to voltage inputs, is necessary to produce the appropriate shift in f res without influencing Z max significantly.  Fig 7g). Limiting the optimal models to the intermediate " g Ca group also revealed a correlation in the " g Ca À " g H parameters (R 2 = 0.79; p < 0.001; Fig 7h). This new correlation may be produced by the balance of the amplitudes of I H and I Ca at the lower and higher voltages, respectively. f res and Q z are maintained by distinct pairwise correlations To determine if any of the MPR attributes were sensitive to the correlations, we ran a 2D sensitivity analysis on a random subset of 50 models. We tested for significant difference in sensitivity across low, intermediate and high levels of " g Ca . In particular, we tested for significant sensitivity of f res and Q Z when parameters were co-varied in directions parallel (L k ) or perpendicular (L ┴ ) to their respective population correlation lines.
We first examined whether f res and Q Z were sensitive to t Ca m À t Ca h for both high (Fig 8a1), low (Fig 8a2), and intermediate " g Ca (Fig 8a3) when parameters were moved along L k and L ┴ (blue and green line ; Fig 8a1- ANOVA; N = 50, p < 0.001), which had a positive sensitivity (Fig 8b). This result indicates that the correlation did a better job at maintaining the value of f res when the value of " g Ca is intermediate or high. For all " g Ca groups, we found that there was a significant interaction between the Z attribute and direction (2-way RM ANOVA; F(1, 49) = 853.52, p < 0.001). When carrying out a pairwise comparison for each direction within an attribute, we found a significant difference in sensitivity between L k and L ┴ for f res (t(93.57) = 28.251, p<0.001). Similarly, for all " g Ca groups, there was a significant difference in sensitivity between L k and L ┴ for Q Z (t(93.57) = -8.294, p<0.001). Because the difference between L k and L ┴ for Q Z was negative, these results suggest that the t Ca m À t Ca h correlation determines f res and not Q Z (Fig 8b). We next examined whether f res and Q Z were sensitive to the " g Ca À V Ca h 1 1=2 correlation for the three model groups (Fig 9a1-9a3). For all " g Ca groups, we found that there was a significant interaction between the Z attribute and direction (2-way RM ANOVA; F(1, 49) = 1262.73.2, p < 0.001). When carrying out a pairwise comparison for each direction within an attribute, we found a significant difference in sensitivity between L k and L ┴ for f res (t(95.18) = 10.10, p<0.001). Similarly, for all " g Ca groups, we found a significant difference in sensitivity between L k and L ┴ for Q Z (t(95.18) = -35.62, p<0.001). Therefore, these results suggest that the " g Ca À V Ca h 1 1=2 correlation determines Q Z and not f res (Fig 9b). Finally, we tested the sensitivity of f res and Q Z to the " g Ca À " g H correlation in the intermediate " g Ca group (Fig 10a). We found that there was a significant interaction between the Z attribute and direction (2-way RM ANOVA; F(1, 11.12) = 2236.2, p < 0.001). When carrying out pairwise comparisons between directions for each attribute, we found there was a significant difference in f res sensitivity between L k and L ┴ (t(93.93) = 2.65, p = 0.0095; Fig 10). Although the sensitivity of Q Z was not 0 for L k , the difference in sensitivity values between L k and L ┴ was also significantly different (t(93.93) = 62.157, p < 0.0001; Fig 10b). These results suggest that, when V low is at -70 mV, for this subset of models to shift f res with only small shifts in Z max , " g H and " g Ca values must be balanced. It may be possible that the Q Z sensitivity is not closer to zero , which is also negatively correlated with " g Ca , should decrease too to compensate for changes in Q Z .

Discussion
Many neuron types exhibit membrane potential resonance (MPR) in response to oscillatory inputs. Several studies have shown that the resonant frequency of individual neurons is correlated with the frequency of the network in which they are embedded [2,6,12,14,22,46]. Moreover, networks of resonant neurons have been proposed to generate more robust network oscillations than neurons with low-pass filter properties [27,28]. In several cases, the underlying nonlinearities and time scales that shape the Z-profile also shape specific properties of the spiking activity patterns, thus leading to a link between the subthreshold and suprathreshold voltage responses [25,47].
Previous work in the crustacean stomatogastric pyloric network has shown that the resonance frequency of the pyloric pacemaker PD neurons is correlated with the pyloric network frequency and is sensitive to blockers of both I H and I Ca [12][13][14]. However, it was not clear how these voltage-gated ionic currents and the passive properties could interact to generate MPR in the PD neurons. Previous modeling work showed that these currents participate in the generation of resonance in CA1 pyramidal neurons [16,17]. However, due to the differences in I Ca time constants, the interaction between its activating and inactivating gating variables did not produce phasonance in CA1 pyramidal neurons, while it does in PD neurons. On a more general level, it is not well understood how the nonlinear properties of ionic currents affect their interplay. Previous studies have shown such interactions may lead to unexpected results, which are not captured by the corresponding linearizations [16][17][18][19]. This complexity is expected to increase when two currents with resonant components are involved [16,48]. We therefore set out to investigate the biophysical mechanism underlying such interactions by using a combined experimental and computational approach and the biological PD neuron as a case study. The two PD neurons are electrically coupled to the pacemaker anterior burster neuron in the pyloric network and their MPR directly influences the network frequency through this electrical coupling [22]. Consequently, our findings have a direct bearing on how the pyloric network frequency is controlled.
Many studies of biophysical models have explored the parameter space using a brute-force technique, by sampling the parameters on a grid [40,49]. Although this technique provides a rather exhaustive sampling of the parameter space, using a fine grid on a large number of free  7) were changed along a line parallel (k, blue) to the correlation line (black) or along a perpendicular line ( j À À , grey). For each model and each line, k or j À À , we fit a line to the relative change in either f res or Q Z as a function of the relative change in g Ca . b. The sensitivity values of f res or Q Z to k or j À À are shown for the three groups. c. Impedance profiles showing how Q Z changes when the parameters vary along a line parallel (blue) or perpendicular (grey) to the g Ca À g H correlation line in one optimal model. Arrows show the direction of the movement of Z max and f res for the change in parameters along k or j À À . parameters could lead to combinatorial explosion and result in a prohibitive number of simulations. On the other hand, a sparse sampling may miss "good" solutions. A multi-objective evolutionary algorithm (MOEA) can generate multiple trade-off solutions in a single run and can handle large parameter spaces very well. In contrast to a brute-force approach, the MOEA can potentially cover a much larger range with possibly hundreds of values [38]. One disadvantage of the MOEA is that, as the number of objectives increases, the search may miss a large portion of the parameter space. This occurs because randomly generated members often tend to be just as good as others, which means that the MOEA would run out of room to introduce new solutions in a given generation. To try to overcome this problem, we carefully chose the parameters of the MOEA such as population size, mutation and crossover distribution indices (100, 20 and 20, respectively) and ensured that the sampled population covered the parameter space evenly. Additionally, we ran the MOEA multiple times, each time collecting all the good parameter sets until one has exhausted all regions of the parameter space where good models exist.
In previous work, we and other authors have examined how the additive interaction of ionic currents with resonant and amplifying gating variables shape the Z and φ profiles at both the linear and nonlinear levels of description [6,15,18,20,32,33,50]. However, the role of inactivating currents in the generation of MPR is not so clear. Authors have established that I Ca can generate MPR in the absence of additional ionic currents [21], that the activation variable diminishes the propensity for MPR and the interaction with I H enhances the dynamic range of parameters producing I Ca -mediated resonance [16]. Even so, to date, only a descriptive explanation of how the ionic current parameters affect certain attributes of MPR has been provided, but no study has provided a mechanistic understanding in terms of the parameters of I Ca that go beyond numerical simulations.
Similar to [16], the model we used in this paper involves the interaction between resonant and amplifying components. Specifically, this model includes a calcium current with both activation (amplifying) and inactivation (resonant) gating variables, and an H-current with a single activation (resonant) gate. Since I H and I Ca shape the lower and upper envelopes of the voltage response to current inputs, respectively [12], given the appropriate voltage-dependence and kinetics of the currents both could play equal roles at different voltage ranges. In fact, either I Ca inactivation or I H is capable of producing MPR [2,21]. In CA1 pyramidal neurons, the differences in Z profiles are due to the passive properties and the kinetics of I H [4]. It is possible that the kinetic parameters of I H and I Ca are tuned so that they contribute nearly equally to shaping the envelopes of the voltage-clamp current.
By tracking the current response to sinusoidal voltage inputs at various frequencies, we found that the f res and f φ = 0 are driven by the peak phase of I Ca , and that f res and f φ = 0 are nearly equal because of the phase matching of I Ca with I L . This is not always the case for neuronal models, and dynamical systems in general, not even for linear models, except for the harmonic oscillator [18][19][20]. In fact, as we mentioned above, this is not the case for the I Ca model used in [16], although our results on the I Ca inactivation time constant are consistent with that study. In these models phase advance for low input frequencies required the presence of I H . The underlying mechanisms are still under investigation and are beyond the scope of this paper. However, the fact that it occurs was crucial to develop a method to investigate the dependence of the resonant properties, particularly the dependence of the f res on the I Ca time constants, using phase information. To date, no other analytical method is available to understand the mechanisms underlying this type of phenomenon in voltage clamp. The tools we developed are applicable to other neuron types for which f res is equal to or has a functional relationship with f ϕ = 0 . However, the conditions under which such a functional relationship exists still needs to be investigated.
Linear correlations between biophysical parameters of the same or different currents have been reported [37] and may be important in preserving the activity of the model neuron and its subthreshold impedance profile attributes [41]. Previous studies examined combinations of parameters in populations of multi-compartment conductance-based models fit to electrophysiological data [16,51] and found only weak pairwise correlations suggesting that the correlations do not arise from electrophysiological constraints. In contrast, constraining the parameters of the ionic currents found to be essential for MPR in PD neuron by MPR attributes, we observed strong correlations underlying parameters when the Z and φ were constrained by the experimental data. We found that constraining the model parameters by f res produced a correlation between the values of time constants of I Ca among the population of 9000 optimal parameter sets. Furthermore, running a 2D sensitivity analysis confirmed that the time constants were constrained so that the effect of making inactivation slower was compensated for by making activation faster to maintain f res constant.
The optimal model parameter sets showed a nonlinear co-variation relationship between the " g Ca and half-inactivation voltage of I Ca . However, the models could be divided into two groups, low and high " g Ca in each of which this co-variation was close to linear. Interestingly, although I Ca alone was the primary current underlying MPR, in the absence of I H (with " g H ¼ 0) the models were restricted to the high " g Ca group. A 2D sensitivity analysis showed that co-varying parameters in each groups along their respective correlation lines preserved Q Z without affecting f res , indicating that each group requires distinct changes in one parameter to compensate for effects of changes in the other. Local sensitivity analysis showed that changes in V Ca h 1 1=2 had opposite effects on f res between high and low " g Ca groups. Increasing V Ca h 1 1=2 decreased f res in high " g Ca models but increased it in low " g Ca models. A previous modeling study has found that changes in V Ca h 1 1=2 greatly influenced the amplitude of MPR with little effect on post-inhibitory rebound in thalamic neurons [21]. It would be interesting to verify whether the mechanisms that generate MPR overlap with those that contribute to post-inhibitory rebound properties.
Previous work in our lab has shown that the voltage range of oscillations significantly affects f res [13]. Here we show that decreasing, V low , the lower bound of the oscillation voltage of the PD neuron, from -60 to -70 mV, significantly shifted f res to smaller values without affecting Z max . Within our optimal model parameter sets, we obtained a set of~1000 models in the intermediate " g Ca range that produced a similar shift in f res but no change in Z max . Because V low greatly affects both I Ca inactivation and I H activation, this indicated a potential interaction between these two currents. In fact, we found that because I H and I Ca are activated preferentially in different voltage ranges, their amplitudes needed to be balanced to keep Z max unchanged when V low was decreased. If the ratio of I H to I Ca amplitudes is incorrect, then Z will amplify (for high " g Ca models) or attenuate (for low " g Ca models). The intermediate " g Ca models also showed a stronger t Ca m À t Ca h correlation, which may be important in matching the phase of I Ca with that of I L . This group also showed a strong " g H À " g Ca correlation, which may provide a mechanism for controlling the changes in I H amplitude at more negative voltage with similar changes in I Ca amplitude at more positive voltages.
In contrast to the findings of Rathour and Narayanan [16], in our optimal models the I H amplitude was not different across the groups with different I Ca properties. However, since I Ca and I H are differentially modulated [45,52], their functional role may overlap when their voltage thresholds and time constants are shifted by neuromodulation. Therefore, we expect that under certain neuromodulatory contexts, I H may play more of an active role in the generation of MPR. A similar effect of two ionic currents on resonance has been observed in the hippocampal pyramidal cells that participate in the theta rhythm, in which two currents, the slow potassium M-current and I H , were found to operate at the depolarized and hyperpolarized membrane potentials respectively to generate theta-resonance [2].
In general, variability of ionic current expression in any specific neuron type should lead to great variability in network output. Yet, network output in general, and specifically the output of the crustacean pyloric network is remarkably stable across animals [30,53,54]. Our results suggest that in oscillatory networks the interaction among ionic currents in an individual neuron may be tuned in a way that the variability of the output is reduced in response to oscillatory inputs. Although our computational study may provide some insight into how such stability is achieved, it also indicates a need for additional mathematical analysis to elucidate the underlying mechanisms.

Methods Electrophysiology
The stomatogastric nervous system of adult male crabs (Cancer borealis) was dissected using standard protocols as in previous studies [14]. After dissection, the entire nervous system including the commissural ganglia, the esophageal ganglion, the stomatogastric ganglion (STG) and the nerves connecting these ganglia, and motor nerves were pinned down in a 100mm Petri dish coated with clear silicone gel, Sylgard 186 (Dow Corning). The STG was desheathed to expose the PD neurons for impalement. During the experiment, the dish was perfused with fresh crab saline maintained at 10-13˚C. After impalement with sharp electrodes, the PD neuron was identified by matching intracellular voltage activity with extracellular action potentials on the motor nerves. After identifying the PD neuron with the first electrode, a second electrode was used to impale the same neuron in preparation for two-electrode voltage clamp. Voltage clamp experiments were done in the presence of 10 −7 M tetrodotoxin (TTX; Biotium) superfusion to remove the neuromodulatory inputs from central projection neurons (decentralization) and to stop spiking activity [13,14]. Intracellular electrodes were prepared by using the Flaming-Brown micropipette puller (P97; Sutter Instruments) and filled with 0.6M K 2 SO 4 and 0.02M KCl. For the microelectrode used for current injection and voltage recording, the resistance was, respectively, 10-15MΩ and 25-35MΩ. Extracellular recording from the motor nerves was carried out using a differential AC amplifier model 1700 (A-M Systems) and intracellular recordings were done with an Axoclamp 2B amplifier (Molecular Devices).

Measuring the Z-profile
During their ongoing activity, the PD neurons produce bursting oscillations with a frequency of~1 Hz and slow-wave activity in the range of -60 to -30 mV. Activity in the PD neuron is abolished by decentralization. The decentralized PD neuron shows MPR in response to ZAP current injection when the current drives the PD membrane voltage to oscillate between -60mV and -30mV, which is similar to the slow-wave oscillation amplitude during ongoing activity [12]. The MPR profiles are not significantly different when measured in current clamp and voltage clamp [14]. Since the MPR depends on the dynamics of voltage-gated ionic currents, it will also depend on the range and shape of the voltage oscillation. Therefore, to examine how Z(f) in a given voltage range constrains the properties of voltage-gated currents and how factors that affect the voltage range change MPR, we measured Z(f) in voltage clamp [10].
To measure the Z-profile, the PD neuron was voltage clamped with a sweeping-frequency sinusoidal impedance amplitude profile (ZAP) function [55] and the injected current was measured [14]. To increase the sampling duration of lower frequencies as compared to the larger ones, a logarithmic ZAP function was used: The amplitude of the ZAP function was adjusted to range between -60 and -30 mV (v 0 = -45 mV, v 1 = 15 mV) and the waveform ranged through frequencies of f lo = 0.1 to f hi = 4 Hz over a total duration T = 100 s. Each ZAP waveform was preceded by three cycles of sinusoidal input at f lo which smoothly transitioned into the ZAP waveform. The total waveform duration was therefore 130 s.
Impedance is a complex number consisting of amplitude and phase. To measure impedance amplitude, we calculated the ratio of the voltage and current amplitudes as a function of frequency and henceforth impedance amplitude will be referred to as Z(f). To measure φ Z (f), we measured the time difference between the peaks of the voltage clamp ZAP and the measured clamp current. One can also measure Z(f) by taking the ratio of the Fourier transforms of voltage and current. However, spectral leakage, caused by taking the FFT of the ZAP function and the nonlinear response, often resulted in a low signal-to-noise ratio and therefore in inaccurate estimates of impedance. Such cases would lead to less accurate polynomial fits compared to the cycle-to-cycle method described above and we therefore limited our analysis to the cycle-to-cycle method.
Because the average Z-profile may not be a realistic representation of a biological neuron, we used the attributes of Z and φ measurements from a single PD neuron as our target. We characterized attributes of Z into five objective functions used for fitting by specifying five points of the profile (Fig 11a). These five points were: • (f 0 , Z 0 ), where Z 0 = Z(f 0 ) and f 0 = 0.1 Hz, • (f res , Z max ), thereby capturing Qz = Z max -Z 0 , • (f 1 , Z(f 1 )) where f 1 = 4 Hz, • The two frequencies at which Z = Z 0 + Q Z / 2. Pinning the profile to these points captures the frequency bandwidth Λ ½ which is the frequency range for which f > Z 0 + Q Z / 2 (Fig 1a).
We also constructed five objective functions to capture the attributes of φ(f) at five points (Fig 11b): • (f φmin , φ min ) where φ min is the maximum phase delay, • (2 Hz, φ f = 2 ) capturing the phase at 2Hz.

Single-compartment model
We used a single-compartment biophysical conductance-based model containing only those currents implicated in shaping Z and φ [12]. We performed simulations in voltage clamp and measured the current as: where I Cm is the capacitive current (C dV dt in nA), C m is set to 1 nF and I L is the voltage-independent leak current in nA. The voltage-dependent currents I curr (I Ca or I H ) in nA are given by where V is the ZAP voltage input (see below), m curr is the activation gating variable, h curr is the inactivation gating variable, " g curr is the maximal conductance in μS, E curr is the reversal potential in mV, and p and q are non-negative integers. For I Ca , p = 3, q = 1 and, for I H , p = 1 and q = 0. The generic equation that governs the dynamics of the gating variables is: where x = m curr or h curr , and The sign of the slope factor (k x ) determines whether the sigmoid is an increasing (negative) or decreasing (positive) function of V, and V x is the midpoint of the sigmoid.
A total of 8 free model parameters were defined (Table 2), which were optimized in light of the objective functions introduced above, to yield a good fit to the Z-profile attributes as described below. The slope factors k x of the sigmoid functions m Ca 1 ðVÞ, h Ca 1 ðVÞ, and m h 1 ðVÞ were fixed at -8 mV, 6 mV, and -7 mV, respectively. V H m 1 1=2 was fixed at -70 mV, using data from experimental measurements in crab [44]. The voltage-dependent time constant for I H was also taken from [44] to be where the range of t H m is given in Table 2.

Fitting models to experimental data
Computational neuroscience optimization problems have used a number of methods, such as the "brute-force" exploration of the parameter space [51] and genetic algorithms [56]. However, the brute-force method is computationally prohibitive for an 8-dimensional model parameter space, which would require potentially very fine sampling to find optimal models. [57]. We used an MOEA (evolutionary optimization) to identify optimal sets of model parameters constrained by experimental Z and φ attributes. MOEAs are computationally efficient at handling high-dimensional parameter spaces and other studies have used them to search for parameters constrained by other types of electrophysiological activity [57] Evolutionary optimization finds solutions by minimizing a set of functions called objective functions, or simply objectives, subject to certain constraints. In our problem, each objective represents the Euclidean distance between the target and the model attributes of Z and φ. When optimizing multiple (potentially conflicting) objectives, MOEA will find a set of solutions that constitute trade-offs in objective scores. For instance, an optimal parameter set may include solutions that are optimal in f res but not in Q z or vice versa and a range of solutions in between that result from the trade-offs in both objectives. In this paper, we used the non-dominated sorting genetic algorithm II (NSGA-II) [38,58] to find optimal solutions, which utilizes concepts of non-dominance and elitism, shown to be critical in solving multi-objective optimization problems [58]. Solution x 1 is said to dominate solution x 2 if it is closer to the target Z(f) and φ(f) profiles in at least one attribute (e.g., f res ) and is no worse in any other attributes (e.g., Q Z , Z 0 , etc.).
NSGA-II begins with a population of 100 parameter combinations created at random within pre-determined lower and upper limits ( Table 2). The objective values for each parameter combination are calculated and ordered according to dominance. First, the highest rank is assigned to all of the non-dominated, trade-off solutions. From the remaining set of parameters, NSGA-II selects the second set of trade-off solutions. This process continues until there are no more parameter combinations to rank. Genetic operators such as binary tournament selection, crossover, and mutation form a child population. A combination of the parent and child parameter sets form the population used in the next generation of NSGA-II [38,58]. NSGA-II favors those parameter combinations-among solutions non-dominating with respect to one another-that come from less crowded parts of the parameter search space (i.e., with fewer similar, in the sense of fitness function values, solutions), thus increasing the diversity of the population. The crowding distance metric is used to promote large spread in the solution space [38]. We ran NSGA-II multiple times (3-5 times, until the mean values of the distributions of optimal parameters was stable) each time for 200 generations with a population size of 100, and pooled the solutions at the end of each run to form a combined population of~9000 parameter combinations. The algorithm stopped when no additional distinct parameter combinations were found. The Z and φ values associated with the optimal parameter sets match the target features (objectives) defining Z and φ to within 5% accuracy.
To test whether two parameters were significantly correlated in the population of 9000 PD models, we calculated the Pearson's correlation coefficients for each pair of parameters and used a permutation test to determine the number of times the calculated correlation coefficient (using a random subset of 20 models). The p-value was given as the fraction of R-values for the permuted vectors greater than the R-value for the original data [51]. We also used a t-test to determine whether the calculated slope of the linear fit differed significantly from zero, which gave us identical results. We repeated both procedures 2000 times, each time with a random subset of 20 models and calculated the percentage of times we obtained a p-value < 0.01.

Sensitivity analysis
We assessed how the values of f res and Q Z depend on changes in parameter values by performing a sensitivity analysis as in [59]. We split the model parameters into two categories: additive, for the voltage-midpoints of activation and inactivation functions, and multiplicative, for the maximal conductances and time constants. We changed the parameters one at a time and fit the relative change in the resonance attributes as a linear function of the relative parameter change. We changed the multiplicative parameters on a logarithmic scale to characterize parameters with both low and high sensitivity.
Multiplicative parameters were varied as p n+1 = exp(±Δp n ) p 0 with Δp n = 0.001 Ã 1.15 n and the sign indicating whether the parameter was increased or decreased. To ensure approximate linearity, we added points to the fit until the R 2 value fell below 0.98. The sensitivity was defined as the slope of this linear fit (Fig 12). For example, if a resonance attribute has a sensitivity of 1 to a parameter, then a 2-fold change in the parameter results in a 2-fold change in the attribute. We changed additive parameters by ±0.5 mV.
We assessed the sensitivity of f res and Q z to parameter pairs (p 1 and p 2 ) that were correlated. We first fit a line through the correlated values in the p 1 -p 2 space. We then shifted this line to pass through a subset of 50 random points in p 1 -p 2 space, resulting in a family of parallel lines, L k . For each point, we also produced a line perpendicular to a line L ┴ . For each model, we performed a sensitivity analysis as before but used the linear fit equation L k or L ┴ to calculate value of p 2 . We fit the relative change in the Z(f) attribute as a linear function of the correlated change in p 1 and p 2 . We used the slope of the linear fit to represent the sensitivity. We used a 2-and 3-way repeated measures ANOVA and the lsmeans function in R to perform pairwise comparisons of means in testing for significant differences between each group of g Ca , each direction, L k and L ┴ , and between each Z attribute, f res and Q Z .
For each model, we solved a system of three differential equations for m H , m Ca and h Ca (voltage was clamped). All simulations were performed using the modified Euler method [60] with a time step of 0.2 ms. The simulation code, impedance calculations, and MOEA were written in C++. MATLAB (The MathWorks) and R were used to perform statistical analyses. 1=2 space shown for all models (grey dots) and those without I H (blue dots). We removed I H by setting " g H = 0, and ran the MOEA multiple times using the same Z-and φ-profiles to constrain the I Ca parameters. A linear fit (green) shows that, when " g H = 0, the relationship between " g Ca À V Each model parameter was changed from the optimal value (origin) in both directions on a logarithmic scale to characterize parameter sensitivity. The slope of a linear fit of the relative change in the Z(f) attribute and the parameter was measured as sensitivity. The parameter was changed until the fit was no longer linear (R 2 <0.98).