Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

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

  • David M. Fox,

    Affiliation Federated Department of Biological Sciences, New Jersey Institute of Technology and Rutgers University, Newark, New Jersey, United States of America

    ⨯
  • Hua-an Tseng,

    Current address: Department of Biomedical Engineering, Boston University, Boston, Massachusetts, United States of America

    Affiliation Federated Department of Biological Sciences, New Jersey Institute of Technology and Rutgers University, Newark, New Jersey, United States of America

    ⨯
  • Tomasz G. Smolinski,

    Affiliation Department of Computer and Information Sciences, Delaware State University, Dover, Delaware, United States of America

    ⨯
  • Horacio G. Rotstein ,

    Contributed equally to this work with: Horacio G. Rotstein, Farzan Nadim

    Affiliations Federated Department of Biological Sciences, New Jersey Institute of Technology and Rutgers University, Newark, New Jersey, United States of America, Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, New Jersey, United States of America

    ⨯
  • Farzan Nadim

    Contributed equally to this work with: Horacio G. Rotstein, Farzan Nadim

    farzan@njit.edu

    Affiliations Federated Department of Biological Sciences, New Jersey Institute of Technology and Rutgers University, Newark, New Jersey, United States of America, Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, New Jersey, United States of America

    ⨯

Abstract

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.

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 the dynamics of voltage-gated ionic currents. Although neuron types can express variable levels of ionic currents, they may have a stable resonant frequency. We used the PD neuron of the crab pyloric network to understand how MPR emerges from the interplay of the biophysical properties of multiple ionic currents, each capable of generating resonance. We show the contribution of an inactivating current at the resonant frequency in terms of interacting time constants. We measured the Z-profile of the PD neuron and explored possible combinations of model parameters that fit this experimentally measured profile. We found that the Z-profile constrains and defines correlations among parameters associated with ionic currents. Furthermore, the resonant frequency and amplitude are sensitive to different parameter sets and can be preserved by co-varying pairs of parameters along their correlation lines. Furthermore, although a resonant current may be present in a neuron, it may not directly contribute to MPR, but constrain the properties of other currents that generate MPR. Finally, constraining model parameters further to those that modify their MPR properties to changes in voltage range produces maximal conductance correlations.

Citation: Fox DM, Tseng H-a, Smolinski TG, Rotstein HG, Nadim F (2017) Mechanisms of generation of membrane potential resonance in a neuron with multiple resonant ionic currents. PLoS Comput Biol 13(6): e1005565. https://doi.org/10.1371/journal.pcbi.1005565

Editor: Koen Vervaeke, University of Oslo, NORWAY

Received: July 29, 2016; Accepted: May 10, 2017; Published: June 5, 2017

Copyright: © 2017 Fox et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: The model neuron (written in MATLAB) and the population of parameter values are available from figshare: DOI 10.6084/m9.figshare.3498473

Funding: This work was supported by National Institutes of Health R01-MH060605 (FN), R01-NS083319 (FN), National Science Foundation DMS1313861 (HGR), National Institutes of Health INBRE 5P20GM103446-15 (TGS), National Science Foundation CREST HRD-1242067 (TGS), and National Science Foundation EPSCoR OIA-1301765 (TGS). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

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 (fres) [1]. MPR has been observed in many neuron types such as those in the hippocampus [2–4] and entorhinal cortex [2–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–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–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 (IH) and the slow outward potassium current (IM). 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 (INaP) and the inward rectifying potassium (IKir) current. Most previous systematic mechanistic studies have primarily examined models with one resonant and one amplifying current, such as IH and INaP, respectively [15, 18–20]. Currents having both activating and inactivating gating variables (in a multiplicative way) such as the low-threshold calcium current (ICa) 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–4, 24] (but see [25]). Interestingly, MPR is not constant across the somatodendritic arbor in these neurons [26]. Hippocampal interneurons also show MPR in vitro, but at gamma frequencies of 30–50 Hz [3, but see 4], and gamma oscillations have been found to be particularly robust in network models containing resonant interneurons [27, 28].

The crab pyloric network produces stable oscillations at a frequency of ~1 Hz, driven by a pacemaker group composed of two neuron types, the anterior burster (AB) and the pyloric dilator (PD), that produce synchronized bursting oscillations through strong electrical-coupling [29]. The PD neuron shows MPR, with fres ~1 Hz that is positively correlated with the pyloric network frequency [12]. Previous work has demonstrated that MPR in this neuron depends on two voltage-gated currents: ICa and IH [12]. Ionic current levels in pyloric neurons are highly variable across animals, even in the same cell type [30]. It is therefore unclear how these currents may interact to produce a stable MPR in the PD neuron and whether this variability persists or is increased or decreased in the presence of oscillatory inputs.

Traditionally, MPR is measured by applying ZAP current injection and recording the amplitude of the voltage response [1, 31]. In some systems, depolarization can increase [32] or decrease [33]) the preferred frequency. Alternatively, resonance is measured by applying ZAP voltage inputs in voltage clamp and recording the amplitude of the total current. Both approaches yield identical results for linear systems, but not necessarily for nonlinear systems. A previous study from our lab using the voltage clamp technique showed that in the PD neuron hyperpolarization decreases both fres and network frequencies [14]. Since MPR results from the outcome of the dynamics of voltage-gated ionic currents activated in different voltage ranges, changing the input voltage amplitude is expected to change fres in an input amplitude-dependent manner. This cannot be captured by linear models in which impedance is independent of the input amplitude. To our knowledge, no study has attempted to understand the ionic mechanisms that produce shifts in fres in response to changes in the voltage range.

Previous studies have explored the generation of MPR by ICa and through the interaction between ICa and IH in hippocampal CA1 pyramidal neurons [16, 17] and thalamic neurons [21], where the resonant and network frequencies are significantly higher than in the crab pyloric network and the ICa time constants are smaller. Based on numerical simulations, these investigations have produced important results about the role of the activating and inactivating gating variables and their respective time constants play in the generation of MPR and the determination of fres. However, a mechanistic understanding of the effects of the interacting time constants and voltage-dependent inactivation that goes beyond simulations is lacking. An important finding for the CA1 pyramidal neurons is that, for physiological time constants, they exhibit resonance, but no phasonance [16]. However, for larger time constants, outside the physiological range for these neurons, they are able to exhibit phasonance. This suggests that PD neurons, which have slower time scale currents, may exhibit resonance and phasonance at comparable frequencies. If so, such a correlation between resonance and phasonance can be used to explain the influence of ionic current parameters.

Our study has two interconnected goals: (i) to understand how the interplay of multiple resonant gating variables shapes the Z- and φ-profiles (impedance amplitude and phase-shift as a function of input frequency) of a biological PD neuron, and (ii) to understand the many ways in which these interactions can occur to produce the same Z-profile in these neurons. For a neuron behaving linearly, e.g., with small subthreshold inputs, this task is somewhat simplified by the fact that linear components are additive. However, neurons are nonlinear and the nonlinear interaction between ionic currents has been shown to produce unexpected results [16, 18, 19].

To achieve these goals we measured and quantified the Z- and φ-profiles of the PD neuron. We then used a single-compartment conductance-based model of Hodgkin-Huxley type [34] that included a passive leak and the two voltage-gated currents IH and ICa to explore what combinations of model parameters can produce the experimentally observed PD neuron Z- and φ-profiles. The maximal conductances of ionic currents of neurons in the stomatogastric nervous system vary widely [35–37]. We therefore assume that the parameters that determine the Z-profile 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 multi-objective 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, IH does not produce MPR but can extend the dynamic range of ICa 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 ICa 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 (fres, Zmax; Fig 1c1). The value of fres was 0.86 ± 0.05Hz producing Zmax 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 fres (Fig 1c2). The PD neuron had a QZ of 2.77 ± 0.71 MΩ and Λ½ of 0.53 ± 0.04 Hz. Across preparations, QZ showed considerable variability, whereas fres, Λ½, and fφ = 0 were relatively consistent (Fig 1d). The corresponding median values for fres, QZ, Λ½, and fφ = 0 were 0.83 Hz, 2.77 MΩ, 0.5 Hz, 0.79 Hz, respectively.

thumbnail
Download:
Fig 1. 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 (Vzap) and the injected current (IPD) 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 fres and Zmax (c1) and fφ = 0 (c2). The shaded region indicates the 95% confidence interval. d. The range of three Z(f) attributes fres, QZ, 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.

https://doi.org/10.1371/journal.pcbi.1005565.g001

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

thumbnail
Download:
Fig 2. Optimal models were fit to the impedance attributes of a single PD neuron.

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: (f0, Z0) = (0.1, 8.2), (fres, Zmax) = (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).

https://doi.org/10.1371/journal.pcbi.1005565.g002

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 (fres) and 4Hz (Fig 3). For 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 IH and IL +ICm sets Itotal at low (~ -60 mV) and high (~ -30 mV) voltages, respectively (Fig 3b left). Since IH deactivation is slow, it also contributes to Itotal at high voltages (Fig 3b right). At 1 Hz (= fres), IH 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 ICa peaks in phase (Fig 3c right) with the passive IL + ICm at high voltages, thus producing a smaller Itotal (magenta bar in Fig 3c). The values of IH at 4 Hz are not much different from 1 Hz (Fig 3d). However, ICa peaks at a much later phase (Fig 3d right), which does not allow it to compensate for IL + ICm at high voltages, thus resulting in a larger Itotal (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 fres 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.

thumbnail
Download:
Fig 3. Passive and voltage-gated currents contribute to the generation of MPR.

a. Z(f) for a random model from the optimal dataset. We measured the steady-state response to sinusoidal voltage inputs (inset) at 0.1 Hz, fres = 1 Hz, and 4 Hz. Voltage-gated (ICa and IH) and passive currents (IL + ICm) are plotted as a function of voltage (left) and normalized time or cycle phase (right) at 0.1 Hz (b), 1 Hz (c), and 4Hz (d). The inset in 5c shows one standard deviation around the mean for the data shown in the right panel, calculated for 200 randomly selected models.

https://doi.org/10.1371/journal.pcbi.1005565.g003

An important collective property of the models we found is that the two frequencies, fres 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–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 fres and fφ = 0 on the model parameters (Fig 4c).

thumbnail
Download:
Fig 4. fres and fφ = 0 of the optimal models are nearly identical.

a. Z(f) (top) and φ(f) (bottom) for a representative optimal model. Green dots indicate fres (top) and fφ = 0 (bottom). b. Histogram showing the difference between fres and fφ = 0 for 500 randomly selected models. A comparison of fres and fφ = 0 of the experimental data of the PD neuron shows a similar distribution (inset, N = 18). (c) Plots of steady-state responses of ICa, IL, and Itotal to sinusoidal voltage inputs at the frequencies marked in panel a shown as a function of normalized time (cycle phase). Dotted vertical line indicates cycle phase 0.5 where the passive currents peak. Solid lines connect the minimum of ICa to the peak of Itotal. The two lines nearly align at fφ = 0.

https://doi.org/10.1371/journal.pcbi.1005565.g004

The current-phase diagrams are depicted as in Fig 3b–3d, as graphs of Itotal, IL and ICa as a function of the cycle phase for each given specific input frequency (Fig 4c). We do not show IH and ICm in this plot, because at frequencies near fres they do not change much with input frequency. Note that IL 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 Itotal 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 IL precisely tracks the input voltage imposes a constraint on the shapes of ICa and Itotal. Therefore, by necessity, if the ICa trough occurs for a cycle phase below 0.5, the Itotal 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 Itotal and ICa and, at fres 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 (Fig 5a) and (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 (= fres when the parameter is at 100%; middle) and f = fres (sides) when fres is different from 1Hz.

thumbnail
Download:
Fig 5. The time constants of ICa activation and inactivation control fres and Zmax.

The Z(f) profiles are plotted for a randomly selected optimal model (green) at different values of (a) and (b). Note that fres of the control (100%) values are at 1 Hz (dashed vertical line). The currents ICa, IL and Itotal plotted as a function of cycle phase at 50% (c1, d1), 100% (c2, d2), and 150% (c3, d3) of the control values of (c) and (d). In each panel of c and d, the currents are shown at 1 Hz (along the dashed lines in a, b) and at fres (filled circles in a, b).

https://doi.org/10.1371/journal.pcbi.1005565.g005

To understand the dependence of Z on changes in and we have to primarily explain the dependence of the two attributes Zmax and fres on these parameters. While fres has a similar monotonic dependence on and (as these parameters increase, fres decreases), Zmax has the opposite dependence on and . The opposite dependence of Zmax on and is a straightforward consequence of the opposite feedback effects (positive for and negative for ) that these parameters exert on ICa. An increase in (for fixed values of ) results in a smaller ICa in response to a given voltage clamp input. Because ICa is smaller and negative, this leads to an increase in Itotal and a decrease in Z at all frequencies. Similarly, an increase in (for fixed values of ) results in a larger ICa, leading to a decrease in Itotal and an increase in Z.

For a fixed value of the input frequency f (e.g. f = 1 Hz as in Fig 5), for Zmax to decrease as increases (Fig 5a), the cycle phase of peak ICa is delayed thereby subtracting less from IL on the depolarizing phase. This leads to Itotal to phase advance relative to IL (Fig 5c) and causes fres to decrease. Similarly, for Zmax to increase as increases (Fig 5b), ICa has to peak later in the cycle thereby subtracting less from IL on the depolarizing phase, which causes Itotal to peak earlier in the cycle, which in turn causes the bar also to swing from the left to the right (Fig 5d). Therefore, fres 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 fres, Λ½ and fφ = 0 are relatively stable across different PD neurons whereas QZ 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 , , , , and ; those showing small CoVs were and the time constant of activation of IH and ICa and half-activation voltage of ICa: , , (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.

thumbnail
Download:
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 according to whether < 0.05 (red; Low ) or > 0.05 (cyan; High ). d. All pairwise relationships, separated on the low or high (colors as in panel c). Green boxes are the same as in b.

https://doi.org/10.1371/journal.pcbi.1005565.g006

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–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: (R = -0.93), (R = 0.73), (R = 0.88), (R = 0.68), (R = -0.82), (R = 0.76), (R = -0.94), and (R = -0.80) (correlations selected with p < 0.01; Fig 6b).

In our experiments, was fixed at -70 mV, using data from experimental measurements in crab [44] (see Methods). However, we also repeated the MOEA with set to -96 mV, as reported in lobster experiments [45], and found that all correlations observed with the former value of remain intact, but simply with a much larger maximal conductance of IH (S1 Fig). In other words, shifting to the left simply results in larger in the optimal models without qualitatively changing the other findings.

In particular, we found that the correlation appeared nonlinear, but there were strong and distinct linear correlations in the two regions > 0.05 μS (low ) and < 0.05 μS (high ; Fig 6c). To ensure that our partitioning of the population into different levels of was valid, we ran the MOEA two additional times, each time using only the mean values of , , , and for either the low or the high values. These optimal models consistently separated into two non-overlapping model parameters, consistent with the low and high models in Fig 6c.

We examined if the low and high models separated or showed distinct correlations in the remaining parameters. The two groups produced non-overlapping subsets of model parameters in the graph. We calculated the Pearson’s correlation coefficient for each pair of parameters in the low and high groups and tested for significance as before (see Table 1). We found that only the high group showed a significant and correlations (Table 1). Additionally, both low and high groups showed the following correlations: , ,, and , . Furthermore, when we ran the MOEA on models where was set to 0, the only optimal models obtained fell within a narrow range of the high group (S2 Fig), which is consistent with the distribution of high models in the panel of Fig 6d.

thumbnail
Download:
Table 1. Statistical p-values obtained using the permutation test of pairwise comparisons for low (lower triangle) and high (upper triangle).

Underlined values are statistically significant (p<0.05).

https://doi.org/10.1371/journal.pcbi.1005565.t001

Decreasing the lower bound of voltage oscillations influences the measured fres and Zmax

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 fres in the PD neuron is influenced by the minimum of the voltage oscillation (Vlow) [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 Vlow was changed from -60 to -70 mV (Fig 7a). Decreasing Vlow significantly decreased fres (by 0.24±0.8Hz), while there was no significant difference in the mean Zmax (-0.15±0.81MΩ) (two-way RM-ANOVA; N = 8, p < 0.001; Fig 7b, left panel).

thumbnail
Download:
Fig 7. The effect of the lower voltage bound Vlow of oscillations on fres and Zmax constrains the optimal models.

a. An example of the change in Z(f) measured in the biological PD neuron for Vlow = -60mV (black line) and Vlow = -70mV (grey line). Inset shows the bounds of voltage clamp inputs in the two cases. b. Shifting Vlow from -60 mV to -70 mV lowers the value of fres measured in the PD neuron significantly, without influencing Zmax (b. Experimental). fres and Zmax values measured in a random subset of optimal model neurons corresponding to low or high values produced the same fres and Zmax values at Vlow = -60mV (black dots), but distinct fres and Zmax values at Vlow = -70mV (low : red dots; high : cyan dots). A subset of optimal models could reproduce the experimental result in which fres shifted to significantly lower values without affecting Zmax. (grey dots). (c) relationship separating out the different groups of models producing different responses to changes in Vlow (colors correspond to b Model panel). Models depicted by grey dots are referred to as intermediate models. (d1-e3) mean voltage-gated ionic currents ICa, IH and ICa+IH and Itotal, shown as a function of voltage for Vlow = -60 mV (d1-d3) and Vlow = -70 mV (e1-e3). Numbers correspond to the location along the as shown in c. f. The intermediate models (grey dots) show a distinct linear correlation. g. Intermediate models (grey dots) show a distinct and tighter correlation compared to all optimal models (black dots). h. Intermediate models (grey dots) show a strong linear correlation that is not observed for all optimal models (black dots).

https://doi.org/10.1371/journal.pcbi.1005565.g007

To explore whether the shift in fres as a function of Vlow could be captured by either low or high models, we measured the shift in fres and Zmax, when Vlow was changed from -60mV to -70mV. We found that fres decreased by 0.24±0.03 Hz and Zmax increased by 5.2±0.6 MΩ for high models, whereas fres decreased by 0.07±0.02Hz and Zmax decreased by 2.6±0.2MΩ for low models (Fig 7b, right panel). Therefore, neither model group reproduced the experimental changes in the Z-profile, specifically, a decrease in fres and no change in Zmax.

We consequently filtered the full optimal dataset (black dots Fig 7c) to find a subset of models that reproduced the change in fres and Zmax (to within 5% of the representative experimental Z(f) shown in Fig 7a) when Vlow 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 and parameters that showed little overlap with the low and high model groups (Fig 7c).

To understand why this particular group (which we will term intermediate ) produced small changes in Zmax when Vlow was decreased, we plotted the current-voltage relationships for ICa, IH, ICa+IH and Itotal for Vlow = -60 and -70 mV, measured at f = 1Hz (fres at Vlow = -60mV) and compared these models with the low and high models. For Vlow = -60mV, the ionic currents behaved similarly for all model groups and Itotal was maximal at -30mV (magenta curve in Fig 7d1–7d3), indicating the similarity of all models in the optimal dataset. However, when Vlow was at -70mV revealed differences in peak ICa, without affecting the peak amplitude of IH across the different groups (Fig 7e1–7e3). The differences in peak ICa accounted for most of the changes in Itotal across the different groups. The Zmax values for intermediate models reproduced the small shift seen in experiments because ICa was at the correct level at high voltages (-30 mV) when Vlow was at -70mV (Fig 7e3). The other two groups did not produce appropriate Zmax for Vlow = -70mV because either ICa was too small (and hence Itotal too large), resulting in a smaller Zmax (Fig 7e1) or vice versa (Fig 7e2). It was also clear that the more negative voltages allowed for an increase in IH levels and therefore larger contribution to the total current. With Vlow at -70mV, not only was there a larger peak amplitude of IH at the lower voltages, but the current at positive voltages also increased because of the very slow deactivation rate. Consequently, IH did not fully turn off when ICa peaked, so that it also contributes to shaping the upper envelope of the total current. IH kinetics were different across the groups (Fig 7e1–7e3). Taken together with the fact that when IH was removed produced only parameter values with very high and very low (S1 Fig), these data suggest that IH could extend the range of ICa parameters over which MPR could be generated through compensation for variable levels of IH.

The ICa in low models was too small when Vlow was -70 mV, because the low conductance did not allow for a significant contribution from the additional de-inactivation (considering the higher in this group) and therefore the peak current did not increase enough. Consequently, the contribution of IH at low voltages was greater than that of ICa at higher voltages (Fig 7e2). Conversely, in the high group, was more negative and so many more channels were available for de-inactivation and the contribution of ICa at higher voltages was much larger than that of IH 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 fres without influencing Zmax significantly.

The intermediate models were strongly correlated in (R2 = 0.89; p < 0.001 Fig 7f1, and had a stronger correlation in the parameters compared to all models (R2 = 0.65; p < 0.001; Fig 7g). Limiting the optimal models to the intermediate group also revealed a correlation in the parameters (R2 = 0.79; p < 0.001; Fig 7h). This new correlation may be produced by the balance of the amplitudes of IH and ICa at the lower and higher voltages, respectively.

fres and Qz 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 . In particular, we tested for significant sensitivity of fres and QZ when parameters were co-varied in directions parallel (L‖) or perpendicular (L┴) to their respective population correlation lines.

We first examined whether fres and QZ were sensitive to for both high (Fig 8a1), low (Fig 8a2), and intermediate (Fig 8a3) when parameters were moved along L‖ and L┴ (blue and green line; Fig 8a1–8a3). For high and intermediate models, fres sensitivities in the L‖ group were negative and not significantly different (3-way RM ANOVA; N = 50, p > 0.05), but both groups were significantly different from the low group (3-way RM 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 fres when the value of is intermediate or high. For all 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‖ and L┴ for fres (t(93.57) = 28.251, p<0.001). Similarly, for all groups, there was a significant difference in sensitivity between L‖ and L┴ for QZ (t(93.57) = -8.294, p<0.001). Because the difference between L‖ and L┴ for QZ was negative, these results suggest that the correlation determines fres and not QZ (Fig 8b).

thumbnail
Download:
Fig 8. Assessing the dependence of fres and QZ on the linear correlation.

a. Parameter values for each model were changed along a line parallel (‖, blue) to the correlation line (black) or along a perpendicular line (┴, grey). This was done for models with high (cyan; a1), low (red; a2) and intermediate (grey; a3) models. For each model and each line, ‖ or ┴, we fit a line to the relative change in either fres or QZ as a function of the relative change in . b. The sensitivity values of fres or QZ to ‖ or ┴ are shown for the three groups.

https://doi.org/10.1371/journal.pcbi.1005565.g008

We next examined whether fres and QZ were sensitive to the correlation for the three model groups (Fig 9a1–9a3). For all 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‖ and L┴ for fres (t(95.18) = 10.10, p<0.001). Similarly, for all groups, we found a significant difference in sensitivity between L‖ and L┴ for QZ (t(95.18) = -35.62, p<0.001). Therefore, these results suggest that the correlation determines QZ and not fres (Fig 9b).

thumbnail
Download:
Fig 9. Assessing the dependence of fres and QZ on the linear correlation.

a. Parameter values for each model were changed along a line parallel (‖, blue) to the correlation line (black) or along a perpendicular line (┴, grey). This was done for models with high (cyan; a1), low (red; a2) and intermediate (grey; a3) models. For each model and each line, ‖ or ┴, we fit a line to the relative change in either fres or QZ as a function of the relative change in . b. The sensitivity values of fres or QZ to ‖ or ┴ are shown for the three groups.

https://doi.org/10.1371/journal.pcbi.1005565.g009

Finally, we tested the sensitivity of fres and QZ to the correlation in the intermediate 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 fres sensitivity between L‖ and L┴ (t(93.93) = 2.65, p = 0.0095; Fig 10). Although the sensitivity of QZ was not 0 for L‖, the difference in sensitivity values between L‖ and L┴ was also significantly different (t(93.93) = 62.157, p < 0.0001; Fig 10b). These results suggest that, when Vlow is at -70 mV, for this subset of models to shift fres with only small shifts in Zmax, and values must be balanced. It may be possible that the QZ sensitivity is not closer to zero along L‖ because , which is also negatively correlated with , should decrease too to compensate for changes in QZ.

thumbnail
Download:
Fig 10. Assessing the dependence of fres and QZ of the intermediate models on the linear correlation.

a. Parameter values for each model were in the intermediate group (see Fig 7) were changed along a line parallel (‖, blue) to the correlation line (black) or along a perpendicular line (┴, grey). For each model and each line, ‖ or ┴, we fit a line to the relative change in either fres or QZ as a function of the relative change in . b. The sensitivity values of fres or QZ to ‖ or ┴ are shown for the three groups. c. Impedance profiles showing how QZ changes when the parameters vary along a line parallel (blue) or perpendicular (grey) to the correlation line in one optimal model. Arrows show the direction of the movement of Zmax and fres for the change in parameters along ‖ or ┴.

https://doi.org/10.1371/journal.pcbi.1005565.g010

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 IH and ICa [12–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 ICa 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–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 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 ICa can generate MPR in the absence of additional ionic currents [21], that the activation variable diminishes the propensity for MPR and the interaction with IH enhances the dynamic range of parameters producing ICa-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 ICa 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 IH and ICa 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 ICa inactivation or IH 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 IH [4]. It is possible that the kinetic parameters of IH and ICa 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 fres and fφ = 0 are driven by the peak phase of ICa, and that fres and fφ = 0 are nearly equal because of the phase matching of ICa with IL. 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–20]. In fact, as we mentioned above, this is not the case for the ICa model used in [16], although our results on the ICa inactivation time constant are consistent with that study. In these models phase advance for low input frequencies required the presence of IH. 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 fres on the ICa 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 fres 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 fres produced a correlation between the values of time constants of ICa 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 fres constant.

The optimal model parameter sets showed a nonlinear co-variation relationship between the and half-inactivation voltage of ICa. However, the models could be divided into two groups, low and high in each of which this co-variation was close to linear. Interestingly, although ICa alone was the primary current underlying MPR, in the absence of IH (with ) the models were restricted to the high group. A 2D sensitivity analysis showed that co-varying parameters in each groups along their respective correlation lines preserved QZ without affecting fres, 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 had opposite effects on fres between high and low groups. Increasing decreased fres in high models but increased it in low models. A previous modeling study has found that changes in 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 fres [13]. Here we show that decreasing, Vlow, the lower bound of the oscillation voltage of the PD neuron, from -60 to -70 mV, significantly shifted fres to smaller values without affecting Zmax. Within our optimal model parameter sets, we obtained a set of ~1000 models in the intermediate range that produced a similar shift in fres but no change in Zmax. Because Vlow greatly affects both ICa inactivation and IH activation, this indicated a potential interaction between these two currents. In fact, we found that because IH and ICa are activated preferentially in different voltage ranges, their amplitudes needed to be balanced to keep Zmax unchanged when Vlow was decreased. If the ratio of IH to ICa amplitudes is incorrect, then Z will amplify (for high models) or attenuate (for low models). The intermediate models also showed a stronger correlation, which may be important in matching the phase of ICa with that of IL. This group also showed a strong correlation, which may provide a mechanism for controlling the changes in IH amplitude at more negative voltage with similar changes in ICa amplitude at more positive voltages.

In contrast to the findings of Rathour and Narayanan [16], in our optimal models the IH amplitude was not different across the groups with different ICa properties. However, since ICa and IH 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, IH 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 IH, 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 K2SO4 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 (v0 = -45 mV, v1 = 15 mV) and the waveform ranged through frequencies of flo = 0.1 to fhi = 4 Hz over a total duration T = 100 s. Each ZAP waveform was preceded by three cycles of sinusoidal input at flo 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:

thumbnail
Download:
Fig 11. Characterization of impedance amplitude Z(f) and phase φ(f) into target objective functions was performed to constrain the model parameters.

The individual objective functions which collectively measure goodness-of-fit were taken as the distance away from characteristic points along the Z(f) and φ(f) profiles (green circles). a. The attributes used along Z(f) were Z0 = Z(f0) at f0 = 0.1 Hz, Z(f1) at f1 = 4 Hz, maximum impedance Zmax = Z(fres) and the two points of the profile at Z0+QZ/2. QZ = Zmax-Z0. Λ½ is the width of the profile at Z0+QZ/2. b. The attributes used along φ(f) were φ(f0), maximum advance φmax, zero-phase frequency fφ = 0, φf = 2 at 2 Hz and maximum delay φmin.

https://doi.org/10.1371/journal.pcbi.1005565.g011

  • (f0, Z0), where Z0 = Z(f0) and f0 = 0.1 Hz,
  • (fres, Zmax), thereby capturing Qz = Zmax—Z0,
  • (f1, Z(f1)) where f1 = 4 Hz,
  • The two frequencies at which Z = Z0 + QZ / 2. Pinning the profile to these points captures the frequency bandwidth Λ½ which is the frequency range for which f > Z0 + QZ / 2 (Fig 1a).

We also constructed five objective functions to capture the attributes of φ(f) at five points (Fig 11b):

  • (f0, φ(f0)),
  • (fφ = 0, 0), where fφ = 0, is the phasonant frequency
  • (fφmax, φmax) where φmax is the maximum phase advance,
  • (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 ICm is the capacitive current ( in nA), Cm is set to 1 nF and IL is the voltage-independent leak current in nA. The voltage-dependent currents Icurr (ICa or IH) in nA are given by where V is the ZAP voltage input (see below), mcurr is the activation gating variable, hcurr is the inactivation gating variable, is the maximal conductance in μS, Ecurr is the reversal potential in mV, and p and q are non-negative integers. For ICa, p = 3, q = 1 and, for IH, p = 1 and q = 0. The generic equation that governs the dynamics of the gating variables is: where x = mcurr or hcurr, and

The sign of the slope factor (kx) determines whether the sigmoid is an increasing (negative) or decreasing (positive) function of V, and Vx 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 kx of the sigmoid functions , , and were fixed at -8 mV, 6 mV, and -7 mV, respectively. was fixed at -70 mV, using data from experimental measurements in crab [44]. The voltage-dependent time constant for IH was also taken from [44] to be where the range of is given in Table 2.

thumbnail
Download:
Table 2. Limits of parameter values allowed for the PD neuron models.

was fixed at -70 mV since there is little variability in the reporting of this experimental measurement [45, 61]. Voltages are in mV, maximal conductances in μS and time constants in ms.

https://doi.org/10.1371/journal.pcbi.1005565.t002

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 fres but not in Qz 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 x1 is said to dominate solution x2 if it is closer to the target Z(f) and φ(f) profiles in at least one attribute (e.g., fres) and is no worse in any other attributes (e.g., QZ, Z0, 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 fres and QZ 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 pn+1 = exp(±Δpn) p0 with Δpn = 0.001*1.15n and the sign indicating whether the parameter was increased or decreased. To ensure approximate linearity, we added points to the fit until the R2 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.

thumbnail
Download:
Fig 12. Linear fits used to assess the sensitivity of impedance attributes on changes in parameters.

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 (R2<0.98).

https://doi.org/10.1371/journal.pcbi.1005565.g012

We assessed the sensitivity of fres and Qz to parameter pairs (p1 and p2) that were correlated. We first fit a line through the correlated values in the p1-p2 space. We then shifted this line to pass through a subset of 50 random points in p1-p2 space, resulting in a family of parallel lines, L‖. 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‖ or L┴ to calculate value of p2. We fit the relative change in the Z(f) attribute as a linear function of the correlated change in p1 and p2. 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 gCa, each direction, L‖ and L┴, and between each Z attribute, fres and QZ.

For each model, we solved a system of three differential equations for mH, mCa and hCa (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.

Supporting information

S1 Fig. Changing the value of does not change the correlations observed among the model parameters.

a. Correlations shown in Fig 8b with at -70 mV. b. Correlations obtained with set to -96 mV (red dots). MOEA was run only once in this case, compared to 5 times in panel a (hence the difference in the number of points). Black dots are the same as panel a. Note that the values of in this case are about 10 times larger than those in panel a, but the correlations (green boxes) remain intact. More importantly, the range of parameters other than is exactly the same in both cases.

https://doi.org/10.1371/journal.pcbi.1005565.s001

(TIF)

S2 Fig. IH extends the dynamic range of ICa parameters over which ICa-mediated MPR occurs.

Parameter values for the optimal models in space shown for all models (grey dots) and those without IH (blue dots). We removed IH by setting = 0, and ran the MOEA multiple times using the same Z- and φ-profiles to constrain the ICa parameters. A linear fit (green) shows that, when = 0, the relationship between is linear and matches a narrow range of the high values in Fig 6c.

https://doi.org/10.1371/journal.pcbi.1005565.s002

(TIF)

Author Contributions

  1. Conceptualization: DMF HGR FN.
  2. Data curation: DMF.
  3. Formal analysis: DMF.
  4. Funding acquisition: HGR FN TGS.
  5. Investigation: DMF.
  6. Methodology: DMF HGR FN.
  7. Project administration: HGR FN.
  8. Resources: FN.
  9. Software: DMF.
  10. Supervision: HGR FN.
  11. Validation: DMF HGR FN.
  12. Visualization: DMF HGR FN.
  13. Writing – original draft: DMF HGR FN.
  14. Writing – review & editing: DMF HGR FN TGS HT.

References

  1. 1. Hutcheon B, Yarom Y. Resonance, oscillation and the intrinsic frequency preferences of neurons. Trends Neurosci. 2000;23(5):216–22. pmid:10782127.
  2. 2. Hu H, Vervaeke K, Storm JF. Two forms of electrical resonance at theta frequencies, generated by M-current, h-current and persistent Na+ current in rat hippocampal pyramidal cells. J Physiol. 2002;545(Pt 3):783–805. pmid:12482886.
  3. 3. Pike FG, Goddard RS, Suckling JM, Ganter P, Kasthuri N, Paulsen O. Distinct frequency preferences of different types of rat hippocampal neurones in response to oscillatory input currents. J Physiol. 2000;529 Pt 1:205–13. pmid:11080262.
  4. 4. Zemankovics R, Kali S, Paulsen O, Freund TF, Hajos N. Differences in subthreshold resonance of hippocampal pyramidal cells and interneurons: the role of h-current and passive membrane characteristics. J Physiol. 2010;588(Pt 12):2109–32. pmid:20421280.
  5. 5. Erchova I, Kreck G, Heinemann U, Herz AV. Dynamics of rat entorhinal cortex layer II and III cells: characteristics of membrane potential resonance at rest predict oscillation properties near threshold. J Physiol. 2004;560(Pt 1):89–110. pmid:15272028.
  6. 6. Schreiber S, Erchova I, Heinemann U, Herz AV. Subthreshold resonance explains the frequency-dependent integration of periodic as well as random stimuli in the entorhinal cortex. J Neurophysiol. 2004;92(1):408–15. pmid:15014100.
  7. 7. Lampl I, Yarom Y. Subthreshold oscillations of the membrane potential: a functional synchronizing and timing device. J Neurophysiol. 1993;70(5):2181–6. pmid:8294979.
  8. 8. Lampl I, Yarom Y. Subthreshold oscillations and resonant behavior: two manifestations of the same mechanism. Neuroscience. 1997;78(2):325–41. pmid:9145790.
  9. 9. Puil E, Meiri H, Yarom Y. Resonant behavior and frequency preferences of thalamic neurons. J Neurophysiol. 1994;71(2):575–82. pmid:8176426.
  10. 10. Beatty JA, Song SC, Wilson CJ. Cell-type-specific resonances shape the responses of striatal neurons to synaptic input. J Neurophysiol. 2015;113(3):688–700. pmid:25411465.
  11. 11. Song SC, Beatty JA, Wilson CJ. The Ionic Mechanism of Membrane Potential Oscillations and Membrane Resonance in Striatal LTS Interneurons. J Neurophysiol. 2016:jn 00511 2016. pmid:27440246.
  12. 12. Tohidi V, Nadim F. Membrane resonance in bursting pacemaker neurons of an oscillatory network is correlated with network frequency. J Neurosci. 2009;29(20):6427–35. pmid:19458214.
  13. 13. Tseng HA, Martinez D, Nadim F. The frequency preference of neurons and synapses in a recurrent oscillatory network. J Neurosci. 2014;34(38):12933–45. pmid:25232127.
  14. 14. Tseng HA, Nadim F. The membrane potential waveform of bursting pacemaker neurons is a predictor of their preferred frequency and the network cycle frequency. J Neurosci. 2010;30(32):10809–19. pmid:20702710.
  15. 15. Richardson MJ, Brunel N, Hakim V. From subthreshold to firing-rate resonance. J Neurophysiol. 2003;89(5):2538–54. pmid:12611957.
  16. 16. Rathour RK, Narayanan R. Inactivating ion channels augment robustness of subthreshold intrinsic response dynamics to parametric variability in hippocampal model neurons. J Physiol. 2012;590(22):5629–52. pmid:22930270.
  17. 17. Rathour RK, Narayanan R. Homeostasis of functional maps in active dendrites emerges in the absence of individual channelostasis. Proceedings of the National Academy of Sciences of the United States of America. 2014;111(17):E1787–96. pmid:24711394.
  18. 18. Rotstein HG, Nadim F. Frequency preference in two-dimensional neural models: a linear analysis of the interaction between resonant and amplifying currents. J Comput Neurosci. 2014;37(1):9–28. pmid:24254440.
  19. 19. Rotstein HG. Subthreshold amplitude and phase resonance in models of quadratic type: nonlinear effects generated by the interplay of resonant and amplifying currents. J Comput Neurosci. 2015;38(2):325–54. pmid:25586875.
  20. 20. Rotstein HG. Frequency Preference Response to Oscillatory Inputs in Two-dimensional Neural Models: A Geometric Approach to Subthreshold Amplitude and Phase Resonance. J Math Neurosci. 2014;4:11. pmid:24872925.
  21. 21. Hutcheon B, Miura RM, Yarom Y, Puil E. Low-threshold calcium current and resonance in thalamic neurons: a model of frequency preference. J Neurophysiol. 1994;71(2):583–94. pmid:8176427.
  22. 22. Chen Y, Li X, Rotstein HG, Nadim F. Membrane potential resonance frequency directly influences network frequency through electrical coupling. J Neurophysiol. 2016:jn 00361 2016. pmid:27385799.
  23. 23. Buzsaki G. Theta oscillations in the hippocampus. Neuron. 2002;33(3):325–40. pmid:11832222.
  24. 24. Hu H, Vervaeke K, Graham LJ, Storm JF. Complementary theta resonance filtering by two spatially segregated mechanisms in CA1 hippocampal pyramidal neurons. J Neurosci. 2009;29(46):14472–83. pmid:19923281.
  25. 25. Stark E, Eichler R, Roux L, Fujisawa S, Rotstein HG, Buzsaki G. Inhibition-induced theta resonance in cortical circuits. Neuron. 2013;80(5):1263–76. pmid:24314731.
  26. 26. Narayanan R, Johnston D. Long-term potentiation in rat hippocampal neurons is accompanied by spatially widespread changes in intrinsic oscillatory dynamics and excitability. Neuron. 2007;56(6):1061–75. Epub 2007/12/21. pmid:18093527.
  27. 27. Moca VV, Nikolic D, Singer W, Muresan RC. Membrane resonance enables stable and robust gamma oscillations. Cereb Cortex. 2014;24(1):119–42. pmid:23042733.
  28. 28. Tikidji-Hamburyan RA, Martinez JJ, White JA, Canavier CC. Resonant Interneurons Can Increase Robustness of Gamma Oscillations. J Neurosci. 2015;35(47):15682–95. pmid:26609160.
  29. 29. Marder E, Bucher D. Understanding circuit dynamics using the stomatogastric nervous system of lobsters and crabs. Annu Rev Physiol. 2007;69:291–316. pmid:17009928.
  30. 30. Marder E. Variability, compensation, and modulation in neurons and circuits. Proceedings of the National Academy of Sciences of the United States of America. 2011;108 Suppl 3:15542–8. pmid:21383190.
  31. 31. Puil E, Gimbarzevsky B, Miura RM. Quantification of membrane properties of trigeminal root ganglion neurons in guinea pigs. J Neurophysiol. 1986;55(5):995–1016. pmid:3711977.
  32. 32. Gutfreund Y, Yarom Y, Segev I. Subthreshold oscillations and resonant frequency in guinea-pig cortical neurons: physiology and modelling. J Physiol. 1995;483 (Pt 3):621–40. pmid:7776248.
  33. 33. Hutcheon B, Miura RM, Puil E. Subthreshold membrane resonance in neocortical neurons. J Neurophysiol. 1996;76(2):683–97. pmid:8871191.
  34. 34. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol. 1952;117(4):500–44. pmid:12991237.
  35. 35. Golowasch J, Abbott LF, Marder E. Activity-dependent regulation of potassium currents in an identified neuron of the stomatogastric ganglion of the crab Cancer borealis. J Neurosci. 1999;19(20):RC33. Epub 1999/10/12. pmid:10516335.
  36. 36. Khorkova O, Golowasch J. Neuromodulators, not activity, control coordinated expression of ionic currents. J Neurosci. 2007;27(32):8709–18. Epub 2007/08/10. 27/32/8709 [pii] pmid:17687048.
  37. 37. Schulz DJ, Goaillard JM, Marder E. Variable channel expression in identified single and electrically coupled neurons in different animals. Nat Neurosci. 2006;9(3):356–62. pmid:16444270.
  38. 38. Deb K. Multi-objective optimization using evolutionary algorithms. 1st ed. Chichester; New York: John Wiley & Sons; 2001. xix, 497 p. p.
    • 39. Rathour RK, Malik R, Narayanan R. Transient potassium channels augment degeneracy in hippocampal active dendritic spectral tuning. Sci Rep. 2016;6:24678. pmid:27094086.
    • 40. Prinz AA, Bucher D, Marder E. Similar network activity from disparate circuit parameters. Nat Neurosci. 2004;7(12):1345–52. pmid:15558066.
    • 41. Zhao S, Golowasch J. Ionic current correlations underlie the global tuning of large numbers of neuronal activity attributes. J Neurosci. 2012;32(39):13380–8. pmid:23015428.
    • 42. Schulz DJ, Goaillard JM, Marder EE. Quantitative expression profiling of identified neurons reveals cell-specific constraints on highly variable levels of gene expression. Proceedings of the National Academy of Sciences of the United States of America. 2007;104(32):13187–91. pmid:17652510.
    • 43. Rotstein HG, Olarinre M, Golowasch J. Dynamic compensation mechanism gives rise to period and duty-cycle level sets in oscillatory neuronal models. J Neurophysiol. 2016;116(5):2431–52. pmid:27559141.
    • 44. Buchholtz F, Golowasch J, Epstein IR, Marder E. Mathematical model of an identified stomatogastric ganglion neuron. J Neurophysiol. 1992;67(2):332–40. pmid:1373763.
    • 45. Peck JH, Gaier E, Stevens E, Repicky S, Harris-Warrick RM. Amine modulation of Ih in a small neural network. J Neurophysiol. 2006;96(6):2931–40. pmid:16943317.
    • 46. Wang XJ. Neurophysiological and computational principles of cortical rhythms in cognition. Physiol Rev. 2010;90(3):1195–268. pmid:20664082; PMCPMC2923921.
    • 47. Engel TA, Schimansky-Geier L, Herz AV, Schreiber S, Erchova I. Subthreshold membrane-potential resonances shape spike-train patterns in the entorhinal cortex. J Neurophysiol. 2008;100(3):1576–89. pmid:18450582.
    • 48. Rotstein HG. Resonance modulation, annihilation and generation of antiresonance and antiphasonance in 3D neuronal systems: interplay of resonant and amplifying currents with slow dynamics. bioaRxiv 2016;091207. 091207.
    • 49. Prinz AA, Billimoria CP, Marder E. Alternative to hand-tuning conductance-based models: construction and analysis of databases of model neurons. J Neurophysiol. 2003;90(6):3998–4015. Epub 2003/08/29. 00641.2003 [pii]. pmid:12944532.
    • 50. Haas JS, White JA. Frequency selectivity of layer II stellate cells in the medial entorhinal cortex. J Neurophysiol. 2002;88(5):2422–9. pmid:12424283.
    • 51. Taylor AL, Goaillard JM, Marder E. How multiple conductances determine electrophysiological properties in a multicompartment model. J Neurosci. 2009;29(17):5573–86. pmid:19403824.
    • 52. Johnson BR, Kloppenburg P, Harris-Warrick RM. Dopamine modulation of calcium currents in pyloric neurons of the lobster stomatogastric ganglion. J Neurophysiol. 2003;90(2):631–43. Epub 2003/08/09. 90/2/631 [pii]. pmid:12904487.
    • 53. Bucher D, Prinz AA, Marder E. Animal-to-animal variability in motor pattern production in adults and during growth. J Neurosci. 2005;25(7):1611–9. pmid:15716396.
    • 54. Golowasch J. Ionic Current Variability and Functional Stability in the Nervous System. Bioscience. 2014;64(7):570–80. pmid:26069342.
    • 55. Gimbarzevsky B, Miura RM, Puil E. Impedance profiles of peripheral and central neurons. Can J Physiol Pharmacol. 1984;62(4):460–2. pmid:6733594.
    • 56. Mitchell M. An introduction to genetic algorithms. Cambridge, Mass.: MIT Press; 1996. viii, 205 p. p.
      • 57. Malik A, Shim K, Prinz A, Smolinski T G. Multi-objective evolutionary algorithms for analysis of conductance correlations involved in recovery of bursting after neuromodulator deprivation in lobster stomatogastric neuron models. BMC Neurosci. 2014;15(Suppl 1):P370.
      • 58. Deb K, Pratap A, Agarwal S, Meyarivan T. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation. 2002;6(2):182–97. WOS:000175082800006.
      • 59. Nowotny T, Levi R, Selverston AI. Probing the dynamics of identified neurons with a data-driven modeling approach. PLoS One. 2008;3(7):e2627. pmid:18612435.
      • 60. Burden RL, Faires JD. Numerical Analysis: PWS Publishing Company—Boston; 1980.
        • 61. Ballo AW, Nadim F, Bucher D. Dopamine modulation of Ih improves temporal fidelity of spike propagation in an unmyelinated axon. J Neurosci. 2012;32(15):5106–19. pmid:22496556.