Analysis of Family Structures Reveals Robustness or Sensitivity of Bursting Activity to Parameter Variations in a Half-Center Oscillator (HCO) Model

Visual Abstract


Visual Abstract
The underlying mechanisms that support robustness in neuronal networks are as yet unknown. However, recent studies provide evidence that neuronal networks are robust to natural variations, modulation, and environmental perturbations of parameters, such as maximal conductances of intrinsic membrane and synaptic currents. Here we sought a method for assessing robustness, which might easily be applied to large brute-force databases of model instances. Starting with groups of instances with appropriate activity (e.g., tonic spiking), our method classifies instances into much smaller subgroups, called families, in which all members vary only by the one parameter that defines the family. By analyzing the structures of families, we developed measures of robustness for activity type. Then, we applied these measures to our previously developed model database, HCO-db, of a two-neuron half-center oscillator (HCO), a neuronal microcircuit from the leech heartbeat central pattern generator where the appropriate activity type is alternating bursting. In HCO-db, the maximal conductances of five intrinsic and two synaptic currents were varied over eight values (leak reversal potential also varied, five values). We focused on how variations of particular conductance parameters maintain normal alternating bursting activity while still allowing for functional modulation of period and spike frequency. We explored the trade-off between robustness of activity type and desirable change in activity

Significance Statement
Robustness is an attribute of living systems and mathematical models that describe them. We developed a method for assessing the robustness of activity types (e.g., bursting), which can be applied to brute-force databases of neuronal model instances in which biologically relevant parameters are varied, and where sensitivity analyses are conceptually and practically difficult to apply. By organizing all instances with appropriate activity into families, in which all members vary only by the one parameter defining the family, we developed measures of robustness for activity type based on family structure and address a fundamental challenge to robustness, modulation, which, by changing parameters, may alter activity type. The method determines which parameters predictably alter activity characteristics, (e.g., burst period), without changing activity type.

Introduction
Robustness is a fundamental feature of complex systems (Kitano, 2004) like neuronal networks, yet remains only operationally defined, and the underlying mechanisms that support it are largely unknown (Gutierrez et al., 2013;Marder et al., 2015). If robustness for a neuronal network is defined as the maintenance of a desirable activity state in the face of parameter variation (e.g., maximal conductance of intrinsic membrane and synaptic currents), then abundant experimental evidence (Swensen and Bean, 2005;Schultheiss et al., 2012;Talebi and Baker, 2012;Tang et al., 2012;Langen et al., 2013;Caplan et al., 2014;Dethier et al., 2015;) suggests that neuronal networks are robust to animal-to-animal variations, modulation, and environmental perturbations of these parameters.
Modulation poses a particular challenge to network robustness, because it often must modify network activity without changing activity type (e.g., changing the period of a rhythmically active network without disrupting its rhythmicity or phase). Some modeling studies have addressed the mechanisms underlying such functional robustness (Goldman et al., 2001(Goldman et al., , 2003. Also, sensitivity analysis is commonly used to assess the influence of a parameter on activity characteristics and type in neuronal models (Olypher and Calabrese, 2007), but is difficult to apply to multidimensional parameter spaces that neurons occupy, though there has been some notable progress Drion et al., 2015).
The purpose of this study was twofold, with a focus on assessing robustness. First, we developed a strategy for assessing the robustness of model neurons applicable to large databases of model instances based on a grid structure (brute-force). We classified instances in the database into small groups called families in which members vary only by the value of the parameter that defines the family, and we used the family structures to develop new measures for assessing the robustness (or sensitivity) of elec-trical activity to changes in model parameters. Second, we applied our measures to an existing model of a twoneuron half-center oscillator (HCO).
Central pattern generators (CPGs; Marder, 2011; pace adaptable rhythmic behaviors, such as walking and breathing. Their inherent rhythmicity and adaptability to perturbations results from the interaction of their intrinsic and synaptic properties with neuromodulatory and sensory inputs. CPGs exhibit remarkably robust activity types (e.g., rhythmic alternating bursting), yet extensive modulation of activity characteristics (e.g., burst period) occurs. Thus, they represent an excellent test bed for exploring the interplay of robustness and modifiability.
Reciprocally inhibitory neurons (often autonomous bursters) called HCOs are prevalent circuit building blocks of CPGs that assure robust alternating bursting (Selverston et al., 2000;Cymbalyuk et al., 2002). HCO models display a wide range of bursting activity when intrinsic and synaptic conductances of the neurons are varied (Prinz et al., 2004). We used the HCO model of Hill et al., 2001', which reproduces the electrical activity observed in the leech heartbeat CPG and consists of a pair of reciprocally inhibitory model neurons represented as single isopotential compartments with intrinsic and synaptic conductances of the type defined by Hodgkin and Huxley (1952). Each model neuron contains eight voltage-dependent currents, five inward currents (I Na , a fast Na ϩ current; I P , a persistent Na ϩ current; I CaF , rapidly inactivating low-threshold Ca current; I CaS , slowly inactivating low-threshold Ca current; and I h , a hyperpolarization-activated cation current), three outward currents (I K1 , a delayed rectifier-like K current; I K2 , a persistent K current; and I KA , a fast transient K current), and two types of inhibitory synaptic transmission: graded (I SynG ) and spike mediated (I SynS ).
To explore the HCO parameter space [maximal conductances (g values) of intrinsic and synaptic currents], we simulated ϳ10.5 million model instances, whose characteristics we recorded into a database named HCO-db (Doloc-Mihu and Calabrese, 2011). Here, we systematically explored the parameter space of two identified groups from the database: realistic HCOs or (rHCOs), which show normal physiological activity with 99,066 instances; and functional HCOs (fHCOs), which show nonphysiological but functional alternating bursting activity with 1.1 million instances. We subdivided instances from both groups into families of instances that vary by only one parameter that defines the family. By examining family structures and patterns, we developed new measures of robustness of bursting activity to alterations of model parameters. Using this analysis, we showed that the HCO model is robust to variations of I h but is highly sensitive to the I P . Moreover, the burst period is reliably and predictably regulated by modulating g h , suggesting that it is ideal for neuromodulation. This new analysis also allowed us to iden-tify ensembles of instances that show typical robust physiological activity for future analysis of parameter variations.

Half-center oscillator (HCO) model
We used a half-center oscillator model (Hill et al., 2001) that produces rhythmic alternating bursting activity resembling the electrical activity in the heartbeat CPG of the leech. This model is publicly available on ModelDB repository (https://senselab.med.yale.edu/ModelDB; accession #19698). The model consists of a two reciprocally inhibitory model interneurons, represented as single isopotential electrical compartments with intrinsic and synaptic membrane conductances of the Hodgkin and Huxley (1952) type. Each model neuron contains, in addition to a leakage current, eight voltage-dependent currents, five inward currents (I Na , I P , I CaF , I CaS , and I h ), and three outward currents (I K1 , I K2 , and I KA ). In what follows, this half-center oscillator model is simply referred to as the model or the model neurons, and the currents are referred to by their letters. The model has two types of inhibitory synaptic transmission between the two interneurons: SynG and SynS. The graded transmission SynG was modeled as a postsynaptic conductance controlled by presynaptic Ca 2ϩ concentration and the spike-mediated transmission SynS was modeled as a postsynaptic conductance triggered by presynaptic spikes. The values for the maximal conductances and the leak reversal potential are the free parameters in the model. For our canonical model, these values are g CaS ϭ 3.2 nS, g h ϭ 4 nS, g P ϭ 7 nS, g K2 ϭ 80 nS, g Leak ϭ 8 nS, g SynS ϭ 60 nS, g SynG ϭ 30 nS, g Na ϭ 200 nS, g CaF ϭ 5 nS, g K1 ϭ 100 nS, g KA ϭ 80 nS, and E leak ϭ Ϫ60 mV (Doloc-Mihu and Calabrese, 2011). The kinetics, voltage dependencies, reversal potentials of the intrinsic currents, and the synaptic connections of the HCO model interneurons have all been verified and previously adjusted to fit the biological data of leech interneurons (Hill et al., 2001). The differential equations of the model are given in the study by Doloc-Mihu and Calabrese, 2014. In previous work, Doloc-Mihu and Calabrese, 2011 performed extensive simulations of this model by systematically varying eight key parameters (a brute-force approach). All model simulations were started from the same initial conditions, which were different for each of the two neurons and were obtained by running the canonical HCO model (Hill et al., 2001) for 200 s, such that one of the two neurons was in its bursting state and the other one was being inhibited. The same parameter values were used in each of the paired model neurons. The eight parameters varied were as follows: seven maximal conductances (g SynS , g SynG , g Leak , g P , g CaS , g h , and g K2 ), across 0%, 25%, 50%, 75%, 100%, 125%, 150%, and 175% of their canonical values, and E leak across Ϫ70, Ϫ65, Ϫ60, Ϫ55, and Ϫ50 mV. After changing a parameter, we ran each model instance for 100 s to allow the system to establish stable activity, and then we ran it for another 100 s, from which we recorded the voltage traces of the electrical activity corresponding to its paired neurons and the corresponding spike times. The firing char-acteristics were analyzed and recorded in a database. By systematically varying the eight key parameters (a bruteforce approach) in all possible combinations, we developed a database of 10,485,760 simulated model instances named "HCO-db" (available upon request; see description in Doloc-Mihu and Calabrese, 2011). The resulting parameter space includes 10,321,920 instances, which have at least one synaptic component present and thus are potential working HCOs, and 163,840 isolated neuron instances, which contain twin neurons without any synaptic interaction.

Definitions
In voltage traces, we recognized a spike only if the potential waveform crossed a threshold of Ϫ20 mV. We defined a burst as having at least three spikes and a minimum interburst interval of 1 s. We defined the burst period as being the interval between the middle spikes of two consecutive bursts. Phase was calculated on a per burst basis, as being the delay from the middle spike of a burst of neuron B to the middle spike of the preceding burst of neuron A divided by the interval from this middle spike of the next burst of neuron A to the middle spike of the preceding burst of neuron A. The duty cycle was defined as the percentage of the period occupied by a burst. Each burst spike frequency is defined as the number of spikes within the burst divided by the burst duration. The spike frequency of a neuron is defined as the mean of all burst spike frequencies divided by the mean of all burst durations.
We defined a half-center oscillator instance (here referred to simply as HCO) as having two model interneurons each showing bursting activity with at least two bursts in a 40 s time interval, with each burst having normal spikes (coefficient of variation of the amplitudes of the spikes within any burst, Ͻ0.07); a small variation of period (coefficient of variation of period, Ͻ 0.05); relative phase in the range of 0.45-0.55; and at least one synaptic component present (either g SynS 0, or g SynG 0, or both g SynS 0 and g SynG 0). We considered an rHCO instance as being an HCO that showed physiological bursting corresponding to that observed in leech oscillator heart interneurons. Precisely, it was an HCO with period between 5 and 15 s, average spike frequency between 8 and 25 Hz, and duty cycle between 50% and 70%. All HCO instances that did not meet the criteria to be designated as realistic (i.e., rHCOs) were designated as fHCOs because they still maintained functional alternating bursting. Note that rHCOs and fHCOs are nonoverlapping subsets of HCOs.
We defined an isolated neuron instance (isolated neuron) as having two identical interneurons (though started with different initial conditions, but otherwise identical), and no synaptic interaction (i.e., g SynS ϭ 0 and g SynG ϭ 0). We defined a burster instance as being an isolated neuron instance for which both neurons had at least two bursts, each with normal spikes, and regular periods (as defined above for the HCOs). We defined a realistic burster as being a burster that showed realistic bursting corresponding to isolated leech oscillator heart interneurons. Pre-cisely, it was a burster with a period between 5 and 15 s, and an average spike frequency between 8 and 25 Hz. Note that realistic bursters can be thought of as being HCOs with no synaptic connections.
We define a family as being a subset (of a group) of instances that all have the same parameter values except the one that defines the family (e.g., all realistic bursters that vary only by g P constitute a family of P). We define a family sequence as the order of values (e.g., increasing) of its defining parameter. Note that each family member has a unique g value of the defining parameter, and all the g values within the same family form the family sequence. A family sequence is delimited by a beginning and an ending parameter value (permissible range), not by the grid limits of the database. If the family sequence includes all grid points in the database within the permissible range, this family is called "noninterrupted"; if not all grid points within the permissible range are members of the family sequence, then the family is called "interrupted." In the latter case, we say that the family sequence has broken down.
By using our definitions from above as criteria, we identified from our HCO-db database those four groups that include simulated instances showing appropriate burst characteristics (period, spike frequency, duty cycle): 99,066 rHCOs (realistic); 1,103,073 fHCOs (functional not realistic); 307 realistic bursters; and 117 bursters that are not realistic. By querying our HCO-db database, we analyzed (with our own Java and Matlab scripts) the sensitivity of the leech burst characteristics to changes in maximal conductances for the realistic groups of instances.

Sensitivity classification
For each family, plotting the burst characteristics (period, spike frequency) versus the g value of each family member and then connecting the points obtained corresponding to the adjacent family members via lines forms a curve. We built and used our own Matlab scripts to automatically separate the families according to their curve steepness calculated as the maximum slope. As a family curve is not smooth but is made of several connected line segments, it was important first to assess its monotonicity. Then, the algorithm can be applied easily. The algorithm used the polyfit and polyval functions in Matlab (code inspired from Matlab Central, Jean-Luc Dellis "getthetangent": http://www.mathworks.com/ matlabcentral/fileexchange/23799-getthetangent/content/get thetangent.m) to fit a smooth curve (polynomial of degree 5 was enough for most h families) to our family curves. The slope of a curve at a certain point is the slope of the tangent at that point. So, the algorithm finds the tangent to the curve with the maximum slope (calculated by using the derivative of the curve). The tangent to a family curve gave the slope or the angle of the decay (with respect to the x-axis). Our algorithm first found the families that decline steeply and classified them as the high-sensitivity group; then it found the families that had a gradual decline and classified them as the low-sensitivity group; and finally the rest of the families were classified as the medium-sensitivity group. We set up a threshold of Ϫ2.5 for the high-sensitivity group (angle range, ϳ68 -90°clockwise), and a maximum of Ϫ0.4 for the low-sensitivity group (angle range, ϳ0 -22°clockwise).

Results
We started our study by first searching for mechanisms involving correlated conductance parameters (g values) that influence the robustness of activity type, here realistic HCO bursting activity, since parameter correlation is one mechanism that produces and maintains robustness. Recent work (Doloc-Mihu and Calabrese, 2014) investigated the potential relationships between parameters that maintain bursting activity in isolated neurons from the HCO model (bursters, see Definitions, in Materials and Methods). The authors found a linearly correlated set of three maximal conductances (of leak current, g Leak ; of a persistent K current, g K2 ; and of a persistent Na ϩ current, g P ) that maintains bursting activity in burster (including realistic burster) model instances, therefore underlying the robustness of bursting activity. In addition, they found that bursting activity was very sensitive to individual variation of these parameters; only correlated changes could maintain the activity type. Now we ask whether or not there is a similar linear correlation between maximal conductances of the g Leak of the g K2 and of the g P that maintains realistic HCO bursting activity.

Nonlinear correlation between g P, g Leak, g K2, and E leak for realistic HCOs
To address this question, we developed a Matlab script to visualize the following five characteristics of a dataset at once: three parameters, which form a 3D parameter space of the data (here, rHCOs), the number of instances projected onto each point in this space given by the size of each point, and pie charts of a fourth parameter showing all instances projected onto each point in the 3D space. Here, the three parameters forming the 3D space were g P , g K2 , and g Leak , and the fourth parameter was E leak . The plot obtained is shown in Figure 1. We selected these parameters inspired by the research in the study by Doloc-Mihu and Calabrese, 2014. The pie chart was split into five slices according to the number of values possible for E leak , with each slice having a different color. If there was no instance projected into the 3D space for a particular value of the fourth parameter, then its corresponding slice was not shown in the pie chart. For a better visualization of the points projected onto the 3D plot ( Fig. 1), we used the natural logarithm to size the pies (as radius). For example, the largest point from the plot included 2,440 projected realistic HCO instances, with 488 per each E leak slice, and its pie size was ln(2440) ϭ 7.7998. The plot in Figure 1 shows that there is no linear correlation among g P , g K2 , and g Leak that is similar to the one observed for bursters in the study by Doloc-Mihu and Calabrese, 2014, which showed a clustering of the points along the main diagonal. The 3D shape from Figure 1 has a complicated contour, similar to a wedge, and reveals potential nonlinear correlation among these three maximal conductances. A similar wedge-like shape correlation was observed by Goldman et al. (2001, their Fig. 2) among maximal conduc-tances of A, Ca, Na, and KCa currents in a theoretical study on lobster superior temporal gyrus neurons.
Expansions of all the E leak pie charts revealed that the main diagonal is a border between E leak values: above the main diagonal were predominantly more positive values of E leak (Ϫ50 and Ϫ55 mV); along the diagonal were all E leak values; and below the diagonal were mostly more negative values (Ϫ70 and Ϫ65 mV). The number of instances seemed to be biased toward the main diagonal as the biggest pies can be seen along it. That is, realistic HCO activity is influenced (produced) by almost equal (including Ϯ25% and Ϯ50% variations) values of g P , g K2 , and g Leak (range, 0 -175%), along with middle values of E leak (Ϫ55, Ϫ60, and Ϫ65 mV).
However, while a linear correlation among g P , g K2 , and g Leak was not observed for the realistic HCO instances, these three parameters (along with E leak ) defined a limited parameter space in which realistic activity was observed. Now we wanted to ask how other parameters (working individually) influence the alternating bursting activity in the HCOs, given the constraints imposed by E leak , g Leak , g K2 , and g P .

Increasing g h promotes bursting
Here, we answer this question by first looking at how varying the value of g h affects bursting activity in HCOs and realistic HCOs. From now on, we separate the HCOs into two nonoverlapping groups, the rHCOs and Plot of all rHCO instances projected onto the 3D space given by the maximal conductances of I P , I K2 , and I Leak. Each point displays a pie chart of the E leak of all instances from the group having the same values of g P , g K2 , and g Leak as the respective projected point. The number of instances projected onto each point in the 3D space is shown by the size of the E leak pie (as radius). For a better visualization of the points, we applied the natural logarithm to the sizes of the pies.

New Research
fHCOs, which include all HCOs that are not realistic HCOs.
Many fHCOs (461,724 instances) had only the duty cycle outside the physiological range [0.5-0.7] for the three criteria used to determine realistic status. Also, there were many (155,556 instances) fHCOs that had only spike frequencies that were too high (Ͼ25 Hz) or too low (Ͻ8 Hz). There were relatively fewer fHCOs (23,523 instances) that had burst periods that were too low (Ͻ5 s) or too high (Ͼ15 s). A good number of fHCOs (84,586 instances) did not satisfy any of the physiological criteria. Also, there were fHCOs that did not meet two of these three realistic requirements, in every combination: 30,751 fHCOs failed period and spike frequency; 87,025 fHCOs failed period and duty cycle; and 259,908 fHCOs failed spike frequency and duty cycle. Figure 2 illustrates the number of instances within each of the rHCO and fHCO groups versus with their period values for each value of g h . Figure 2 reveals right-skewed distributions of the fHCO and rHCO instances for each h value. For fHCOs, increasing the value of g h has the following two effects: it monotonically increases the number of instances within the group (from 94,639 instances for g h ϭ 0 to 181,083 for g h ϭ 175%), and it increases the number of instances with faster bursting activity (periods between 5 and 10 s). For the latter, we used the skewness Matlab function to show that the fHCOs distributions have positive monotonically increasing skewness values (from 0.95 for g h ϭ 0 to 1.93 for g h ϭ 175%) as the value of g h increases. For rHCOs, Figure 2 shows a peak in the number of these instances for g h ϭ 75%. Similar to the distributions of the fHCOs, the distributions of the rHCOs show positive monotonically increasing skewness values (from 1.38 for g h ϭ 0 to 2.9 for g h ϭ 175%) as the value of g h increases. We conclude that a larger g h value promotes functional and realistic HCO bursting.
The amount of h current thus appears to be important in regulating alternating bursting in mutually inhibitory heart interneurons. Moreover, increasing g h values increases the number of instances with shorter burst periods (between 5 and 10 s). Next, we investigated in more detail the effect of each parameter on the realistic bursting activity.

Robustness as defined by the family size
For each parameter, we queried the HCO-db database to build up the families (of instances) existing for that particular parameter within the rHCO group (we also build up the families within fHCOs; data not shown). Then, we partitioned the rHCOs into families according to the number of members in the family (see Definitions). Table 1 summarizes, for each parameter, the number of families of each size (number of members) within the rHCOs. From Table 1, we can see that P has many families (on the order of thousands) with a small number of members (one and two members) and no large families (with three to eight members). K2 shows a family structure similar to that of P, with thousands of families with a small number of members (one to three members), but also has a small number of families with four and five members. Leak, CaS, and SynS share a similar family structure, with thousands of families with a small number of members (one to three members), and a good number of families with four to six members. Finally Figure 2. The number of instances within the rHCO group and fHCO group vs their period values for each value of g h . fHCOs are in blue and rHCOs are in red. Period values were between 1 and 40 s for the fHCOs and between 5 and 15 s for the rHCOs. Plots were obtained using the bar Matlab function (R2014a). Each bar shows the period values within one unit or between [n, n ϩ 1). The bars corresponding to the realistic instances were shifted with 0.3 on the x-axis for better visualization.
g h values are given as the percentage of the canonical (100%). structure, including a large number of families with seven and eight members, and thousands of the rest.

Family sequence
In this work, we used the g values within the same family in increasing order of their values to form the family sequence (see Definitions). By analyzing the family sequences, we revealed that many families have internal gaps or breaks. A family sequence is delimited by a beginning and an end parameter value (permissible range), not by the grid limits of the database. A family is noninterrupted if its sequence includes all grid points in the database between these limits; otherwise, the family is called interrupted (i.e., if there is at least one missing grid point from the family sequence between these limits). In Table 1, we present in parentheses the number of noninterrupted families. The percentage of noninterrupted families within the total number of families of a given number of members varied for different parameters. For g K2 , it showed a steady increase with increasing number of family members from 89.9% to 100%. For g P , we had all families as noninterrupted families (100%). For the rest of the parameters, the percentage of noninterrupted families fluctuated (i.e., changed nonmonotonically), as follows: between 58.7% and 86.7% for g Leak (corresponding to families with five and two members, respectively); between 73.2% and 100% for g CaS (for families with three and six members, respectively); between 51.2% and 70% for g SynS (for families with five and six members, respectively); between 32.3% and 100% for g h (for families with six and eight members, respectively); and between 17.5% and 100% for g SynG (for families with five and eight members, respectively). Interestingly, both g h and g SynG showed a similar tendency: both had a large number of families with any given number of members, and, with the exception of their families with eight members, both displayed very low numbers of noninterrupted families (Ͻ56% for g h and Ͻ42.2% for g SynG ).
Here, we used a large number of families with many members (best if all eight members are present) as a measure of the robustness of realistic bursting activity over the defining parameter changing within these families. Moreover, noninterrupted families indicated strong robustness to the variation of the defining family parameter. Using this criterion, from above we can conclude that realistic HCO activity is robust to changes in h and SynG, and that it is sensitive to changes in P and maybe K2. Next, we looked at the family beginning and end values of the defining g , for each g value, along with their potential sequence interruptions to assess potential robustness or sensitivity.

Robustness as defined by family sequence breakdown
For each parameter, we plotted the distributions of its families that were interrupted. The presence of interruptions in a family sequence indicates that realistic HCO activity is not robust over the entire permissible range of the family's defining parameter, with all other parameters remaining constant. Conversely, noninterrupted (or con-tinuous) families indicate robustness to the variation of the defining family parameter. Figure 3 illustrates the distributions of interrupted and noninterrupted h families for the rHCOs. There were a total of 25,057 families with two to seven members, of which 11,066 were noninterrupted and 13,991 were interrupted. As the permissible range of the family sequence increased-move along any given row (beginning point) toward the right-most column (largest end point)-the number of interrupted families increases, while the number of noninterrupted families tends to decrease then increase for the largest sequence end point (g h ϭ 175%). For large permissible ranges, there are a large number of families with multiple interruptions in their sequences. Conversely, the smaller the permissible range, the larger the proportion of noninterrupted families. This analysis also points out that with a large permissible range (Ͼ4) of g h , regardless of the beginning or end point permissible, there are a considerable number of families with a sizeable number of members (more than four), albeit some with interruptions, and there are many families that have such permissible ranges (4,669). What is somewhat surprising about these results is that there are so many interruptions in family sequences, and for large permissible ranges there can be several interruptions. It is thus important to explore whether these interruptions lead to a change in activity type or reflect a more subtle change like a movement beyond the physiological range criteria (e.g., spike frequency) for an rHCO, while maintaining fHCO activity. Figure 4 illustrates the distributions of interrupted and noninterrupted CaS families for the rHCOs. There were a total of 24,110 CaS families with two to six members, of which 18,157 were noninterrupted and 5,953 were interrupted. No family is thus able to span the entire range of g CaS values tested. Figure 4 shows that the number of interrupted families generally decreases as the permissible range increases (look across rows), which is in stark contrast to h families (Fig. 3), except for sequences that start at g CaS ϭ 0, where there is a peak for sequences that end at g CaS ϭ 75%. Similarly, the number of noninterrupted families decreases as the permissible range increases (look across rows). This analysis points out that there are no families with seven or eight members and few with six members (4 families) with a modest number of families with five members (67 families), most of which are noninterrupted (63 families). Still, as with h families, there can be several interruptions in family sequences except for families with permissible ranges of approximately four. As with h families, it is important to explore whether these interruptions lead to a change in activity type or reflect a more subtle change like a movement beyond the physiological range criteria (e.g., spike frequency) for an rHCO, while maintaining fHCO activity.

Sequence breakdown versus changes in g P
For rHCOs, we obtained only families with two members for P, and all these families are noninterrupted (Table 1).
There were no families of P that started at 0. The peak in the number of P families was obtained for families starting at 50% and ending at 75% (1,614 families), followed by families starting at 75% and ending at 100% (861 families).

Sequence breakdown versus changes in g K2
For rHCOs, there were a total of 14,497 K2 families with two to five members, of which 13,144 families were noninterrupted. The number of K2 families, both noninterrupted and interrupted, increased monotonically as the end point of the sequence increased, with the largest number (4,600 and 741, respectively) of these families for g K2 ϭ 175%. There was an increase in the number of families as the g K2 end point increased, which suggests that more g K2 promotes realistic HCO bursting.

Sequence breakdown versus changes in g Leak
For rHCOs, there were a total of 14,557 leak families with two to six members, of which 12,273 were noninterrupted. The maximum number of noninterrupted families of any size (number of members) was when the sequence started with no leak present (g Leak ϭ 0), and as the amount of g Leak increased the number of noninterrupted families decreased monotonically. Lower g Leak values (0 -100%) seemed to promote more noninterrupted families and, thus, resulted in more robust realistic HCO activity.

Sequence breakdown versus changes in g SynS
There were a total of 25,128 SynS families with two to six members, of which 15,432 were noninterrupted. For each number of members, the number of noninterrupted families increased monotonically and steeply as the sequence beginning point increased and also as the sequence end point increased (peak at 175%). These observations sug-gest that stronger spike-mediated synaptic transmission promoted robust realistic HCO activity.

Sequence breakdown versus changes in g SynG
There were a total of 22,335 SynG families with two to seven members, of which 6,437 were noninterrupted. Like h, SynG had many families (1,023 families) with all eight  possible members (noninterrupted sequence). As the sequence beginning point increased (peak at 0), both the number of noninterrupted families and interrupted families decreased monotonically and steeply (from 1,904 and 8,267 for g SynG ϭ 0 respectively, to 384 and 218, respectively, for g SynG ϭ 150% and g SynG ϭ 125%). So, it appears that more g SynG did not promote robustness.

Robustness as shown by missing family members
In the following, we analyze the missing members from all rHCO families (noninterrupted and interrupted); the missing members outside the family's permissible range (beginning and end points); and the missing members in the family sequence to determine the proportion of missing members that maintain functional HCO activity (regular alternating bursting but not realistic). We use this as a measure of the robustness of alternating bursting activity state.

Bursting activity shows robustness to changes in g h
The plots in Figure 5 illustrate the distributions of the missing members from h families of the rHCOs for each g h value. In these plots, the bars at g h ϭ 0% and g h ϭ 175% values contain only missing members outside the family's permissible range, and all bars but these extremes are mixed. As family size decreases, more of the missing values illustrated are from beyond the family's permissible range. We grouped together all families with the same number of members (and thus the same number of missing members) and showed their data on the same subplot. We color coded the activity type of the missing members. It is easy to see that yellow, the color corresponding to the fHCOs, is predominant, which reveals that most (from 85% for families with seven members to 68% for families with one member) of the missing members from the interrupted h families of the rHCOs are in fact fHCO instances. Table 2 shows the number of missing (h) family members of rHCOs that have functional HCO bursting grouped by the physiological criterion or criteria that they fail. For all families, regardless of size (number of members), the biggest number of missing instances from rHCO interrupted families did not meet the duty cycle criterion. A very small number of these missing instances from rHCO families did not satisfy any of the three realistic criteria. In addition, Table 1 shows that there was a large number of multimember h families of the rHCOs (with   more than a third of them being noninterrupted), including 144 perfect h families with all eight members present (circuit stability/strength). These observations show that variations in g h might move HCO bursting outside the physiological range but still keep it in the functional HCO range, which shows the robustness of the functional alternating bursting activity to these variations.

Bursting activity shows robustness to changes in g CaS , if it is sufficiently high
Distributions (data not shown) of the missing members from CaS families of the rHCOs for each g CaS value revealed that, if the model neurons had low g CaS values (0, 25%, and 50%), most of the missing members from the CaS families (ϳ61%) had continuous spiking activity, and that if the model neurons had sufficiently high g CaS (at least 75%), then most (between 79% for families with one member and 95.6% for families with five members) of the missing members from these families had fHCO bursting activity. That is, a minimum of g CaS ϭ 75% was the cutoff for system robustness; changes in g CaS values Ͻ75% will result in changes in the activity type, and changes in g CaS values starting from 75% will result in keeping the functional HCO bursting activity.
Bursting activity shows sensitivity (no robustness) to changes in g P Table 1 shows that the rHCOs had P families with only one or two members, and none of these are interrupted. The plots in Figure 6 show the distributions of the missing members from interrupted P families of the rHCOs for each g P value; thus, the plots in Figure 6 had seven and six missing members, respectively. These plots reveal that the activity of the rHCOs was quite sensitive to changes in the g P value: any change in g P value was very likely (Ͼ85.6%) to move the activity type outside functional HCO bursting. For g P values Ͻ100%, the activity type became predominantly silent (60%) or spiking (12%), and for g P values Ն100% the bursting activity became (between 28.8% for families with one member and 66.7% for families with two members) asymmetric (e.g., one cell is silent while the other cell is spiking; Doloc-Mihu and Calabrese, 2011), continuous spiking (between 25.8% for families with one member and 17.4% for families with two members), functional (13.6%), or silent (8%, but only for families with one member). In summary, the activity became predominantly nonfunctional (i.e., no alternating bursting).
Bursting activity shows more sensitivity (than robustness) to changes in g K2 Table 1 shows that K2 had mainly noninterrupted families (13,094), with some interrupted families (1,403). The distributions of the missing members from K2 families for the rHCOs (data not shown) reveal that rHCO activity was quite sensitive to changes in g K2 value: a change in g K2 value was very likely (Ͼ81%) to move the activity type outside functional HCO bursting. For low g K2 values (0%  or 25%), the activity type became mainly spiking or plateau, and for the rest of the g K2 values the activity became mostly asymmetric, to a lesser extent spiking, or in a significant number of cases fHCO (varying from 16% to 26% across g K2 values), which was very similar to the sensitivity to changes in P current. Increasing the amount of g K2 changed the activity type from spiking to bursting, but mostly asymmetric bursting and not fHCO activity.
Bursting activity shows moderate sensitivity to changes in g Leak Table 1 shows that Leak had mostly noninterrupted families (12,235) with a fairly large proportion of interrupted families (2,322). Distributions of the missing members from Leak families for the rHCOs (data not shown) reveal that the activity of the rHCOs was quite sensitive to changes in g Leak values, as the proportion of missing members that showed functional HCO activity was low, varying between 21% and 51% across g Leak values. Increasing g Leak values Ͼ100% increased the proportion of missing members with spiking activity and decreased the proportion of missing members with asymmetric activity (i.e., outside the fHCO activity). It seems that g Leak ϭ 100% was the cutoff for model sensitivity: any g Leak value above it decreased the proportion of missing family members with fHCO bursting. This observation indicates that high g Leak values move the model neurons outside the fHCO bursting range, but lower g Leak values can maintain activity inside the fHCO bursting. Table 1 shows that SynS had many noninterrupted families (15,432) but also many interrupted families (9,696). Distributions of the missing members from SynS families for the rHCOs (data not shown) revealed that low values of g SynS (Ͻ100%) produced mostly (between 57.7% and 93.7% across g SynS values) missing members with spiking activity, and that for sufficient values of g SynS (Ն100%) most (between 64.5% and 96.8%) missing members of the SynS families had fHCO bursting activity. As the amount of g SynS increased, there were substantially more missing members with fHCO activity and fewer missing members with spiking activity within each group of families with the same number of missing members. This means that a strong spike-mediated synapse was necessary (at least g SynS ϭ 100%) to maintain fHCO bursting activity. In addition, the relatively large number of multimember SynS families of the rHCOs shows that rHCO bursting activity was quite robust to variations of the spike-mediated synapse (g SynS ); variations of g SynS , especially toward large values (at least g SynS ϭ 100%) maintained either realistic or functional HCO bursting activity. Table 1 shows that SynG had mainly interrupted families (15,898) with a significant proportion of noninterrupted families (6,437). Distributions of the missing members from SynG families for the rHCOs (data not shown) revealed that most (between 75% and 81%) of the missing members from the SynG families of the rHCOs were in fact fHCO instances. This observation indicates that vary-ing g SynG does not disrupt fHCO bursting activity. Based also on the large number of multimember SynG families of the rHCOs, including 1,023 families with all eight members present, it seems that in the HCO model realistic bursting activity was robust to variations of the g SynG ; variations in g SynG values maintained either realistic or functional HCO bursting activity of the neurons.

Sensitivity of burst characteristics to variations of g
To assess the sensitivity of a characteristic (period, spike frequency) to the variation of a parameter, we plotted the burst characteristic (period or spike frequency) versus the corresponding maximal conductance values of all family members for each of the families (of the parameter considered) that have the same number of members. Then, we connected via a line the burst characteristic values corresponding to two consecutive maximal conductance values (adjacent family members); thus, for each family a curve was plotted. We developed Matlab scripts to automatically analyze each such family curve. First, our scripts checked whether a family curve was monotonic or not, by simply comparing the burst characteristic values of each two adjacent family members (i.e., by calculating their difference). If, for a family, all these differences have the same sign (either positive or negative), the family curve is considered monotonic; otherwise, it is considered nonmonotonic. We separated the families into subgroups according to their monotony (e.g., all four member h families showing monotonically decreasing curves corresponding to period were grouped together). Then, our scripts calculated the average change in burst characteristic for each such subgroup of families (sum over all families of the difference of the last and first members' burst characteristic values in each family divided by the number of families). Figure 7, A and B, shows the period plots obtained for h families with eight members and with four members, respectively, for the rHCO group. However, we obtained similar plots for all possible h family sizes of the rHCOs. For all families of h, such plots (Fig. 7A,B) showed that increasing the maximal conductance (g h ) of the hyperpolarization-activated current monotonically decreases the burst period of the rHCOs, which confirms and extends the results in the study by Hill et al. (2001) stated for the canonical model when one parameter was varied at a time (more restricted parametric space). This decrease occurs because of the ability of the h current to depolarize the inhibited neuron and advance the transition to the burst phase. Figure 7, C and D, shows the spike-frequency plots obtained for h families with eight members and with four members, respectively, for the rHCO group. For most families of h (except for 8 out of 2,703 families with five members, for 22 of 4,322 families with four members, for 63 of 6,325 families with three members, and for 353 of 9,895 families with two members), spike-frequency plots (Fig. 7C,D) revealed that increasing the maximal conductance (g h ) of the 560 hyperpolarization-activated current monotonically increased the spike frequency of the rH-COs. As the family size increased, the average increase in spike frequency for each family increased, from ϳ0.29 Hz for h families with two members to 0.71 Hz for h families with eight members (amount calculated as the ratio between the average increase in frequency for all families having same number of members and the number of families). However, as stated in the previous section, increasing g h did not increase the burst period but decreased it, due to the ability of I h to promote the escape of the inhibited neuron. The increase in spike frequency happens because by increasing g h the inhibited neuron escapes earlier and starts bursting earlier, which makes the burst duration shorter, leading to higher average spike frequencies.

Increasing g CaS increases both period and spike frequency
For most families of CaS, such period plots (data not shown) revealed that increasing the maximal conductance of the slow Ca current (g CaS ) monotonically increased the cycle period of the rHCOs. For many of these families, the increase was almost linear. However, for few CaS families-2 families (of 929) of four members, 35 families (of 5,263) of three members, and 137 families (of 17,851) of two members-increasing g CaS decreased (not monotonically) the period. We conclude that these results confirm and extend the previous results in the study by Hill et al. (2001).
For all families of CaS, spike-frequency plots (data not shown) revealed that increasing the maximal conductance of the slow Ca current (g CaS ) monotonically increased the spike frequency of the rHCOs. For most families, the increase was big (with an average increase from 5 Hz for families with three members up to Ͼ11 Hz for families with six members). This result confirmed the results in the study by Hill et al. (2001) regarding the influence of CaS on spiking activity of the canonical model and explains the effect on period because increased spike frequency leads to greater inhibition of the opposite cell of the HCO.

Increasing g P increases both period and spike frequency
For most (except for 34 of 3,048) families (of two members) of P, plots (data not shown) revealed that increasing the maximal conductance of the persistent Na ϩ current (g P ) increased the period of the rHCOs. For many of these families, the increase was large (Ͼ9.77 s); the average (over all families with two members) increase of the period was 4.255 s.
For the rHCOs, spike-frequency plots (data not shown) of P families revealed that increasing the maximal con- ductance of the persistent Na ϩ current (g P ) increased the spike frequency. For many P families, the increase was noteworthy, with an average increase of 5.37 Hz. The effect on spike frequency accounts for the effect on period through increased inhibition of the opposite cell of the HCO.
Increasing g K2 decreases the period For the rHCOs, plots (data not shown) of period versus g K2 values for all families of K2 revealed that for most (except for 19 of 13,163 families with two members) families of K2, increasing the maximal conductance of the persistent K current (g K2 ) monotonically decreased the period of the rHCOs. For many of these families, the decrease was large (average of 3.32 s for families with two members, of 4.87 s for families with three members, and of 6.59 s for families with four members). Once again, our results confirm and extend the previous results of the study by Hill et al. (2001).
For many families of K2 (except for 424 of 13,163 families with two members, and for 4 of 1,299 families with three members), spike-frequency plots (data not shown) showed that an increase in the maximal conductance of the persistent K current (g K2 ) led to a decrease in the spike frequency. The decrease was Ͼ2 Hz on average, and it occurs due to the ability of I K2 to decrease the peak of the (slow-wave) oscillation during the burst. Again, the change in spike frequency with g K2 accounts for the change in period due to a change in inhibition of the opposite cell in the HCO.

Increasing g Leak decreases both period and spike frequency
For the rHCOs, plots (data not shown) of period versus g Leak values for all Leak families revealed that for most (except for 345 of 11,873 families with two members, for 1 of 2,133 families with three members, and for 1 of 452 families with four members) families increasing g Leak decreased the period. The average decrease of the period was 2.46 s for families with two members and went up to 5.78 s for families with six members.
For most Leak families (except for 13 of 11,873 families with two members), spike-frequency plots (data not shown) of spike frequency versus g Leak values revealed that increasing g Leak decreased the spike frequency. The average decrease of the spike frequency was 2.72 Hz for families with two members, and it increased up to 7.55 Hz for families with six members. Again, the change in spike frequency with g Leak accounts for the change in period due to a change in inhibition of the opposite cell in the HCO.

Increasing g SynG results in negligible changes in period and spike frequency
Plots (data not shown) of period versus g SynG values for all rHCO families of SynG revealed that for many families (except for 1,241 of 7,047 families with two members, for 1,099 of 5,086 families with three members, for 1,114 of 4,124 families with four members, for 1,038 of 3,046 families with five members, for 799 of 1,971 families with six members, for 572 of 1,061 families with seven members, and for 761 of 1,023 families with eight members) increasing g SynG increased the period. As the number of members within families increased, there was an increase in the number of families for which increasing g SynG decreased the period. However, the change (either increase or decrease) in period was small (Ͻ0.4 s on average), and therefore we conclude that SynG does not play an important role in controlling period.
For most families of SynG, spike-frequency plots (data not shown) revealed that increasing g SynG changed the spike frequency of the rHCOs negligibly (average change between 0.18 and 0.31 Hz for families with two and seven members, respectively). For most families (except for 2,291 of 7,047 families with two members, for 1,522 of 5,086 families with three members, for 999 of 4,124 families with four members, for 704 of 3,046 families with five members, for 409 of 1,971 families with six members, for 214 of 1,061 families with seven members, and for 170 of 1,023 families with eight members), increasing g SynG slightly increased the spike frequency. We conclude from the small changes in burst characteristics (period and spike frequency) when g SynG is varied that, while some amount of graded transmission may contribute to alternating bursting (maybe 25% is good enough), there is no added benefit to having higher g SynG values.

Increasing g SynS increases the period, but decreases the spike frequency
For almost all rHCO families of SynS, plots (data not shown) revealed that increasing the maximal conductance of the spike-mediated synapse (g SynS ) monotonically increased the burst period of the rHCOs (except for 74 of 16,950 families with two members, and for 3 of 6,390 families with three members). The average increase of period was from 2 s for families with two members to 5 s for families with six members.
For many SynS families (Ͼ55%), spike-frequency plots (data not shown) showed that increasing the maximal conductance of the spike-mediated synapse (g SynS ) decreased the spike frequency of the rHCOs. The decrease started from an average of 0.35 Hz for families with two members, and it increased to an average of 2.5 Hz for families with six members. The larger the family size, the larger was the number of the families that show this decrease in spike frequency (from 55% of families with two members to 90% of families with six members).

Analysis of family structure shows how g h influences period of realistic HCO instances
To analyze how burst characteristics in rHCOs are influenced by varying a g value, we took advantage of family structure and made graphs like those in Figure 7, in which we plotted a burst characteristic [e.g., period (Fig. 7A,B) or spike frequency (Fig. 7C,D)] vs g values for families of different sizes. By looking across families, we could spot trends in these relationships. We will focus here on period versus g h values for h families of rHCO instances. All these plots show monotonically decreasing curves of period versus g h values for each h family, regardless of family size (we checked the monotonies using our Matlab scripts). We observed that these monotonically decreasing curves show different rates of decline, indicating dif-ferent sensitivities of period to changes in g h . Some seemed to decline very steeply with increasing g h , and some seemed to decline quite gradually.
Using our own Matlab scripts (see Sensitivity classification), we automatically separated the h families according to their curve steepness (period sensitivity) into the following three groups: high-sensitivity (steep slope), medium-sensitivity, and low-sensitivity (shallow slope). Figure 8A shows the split between the curves of the h families with eight members whose curves were illustrated in Figure 7A. Our algorithm found 35 families with high-sensitivity, 57 with medium-sensitivity, and 52 with low-sensitivity. Table 3 provides the number of curves for each sensitivity group for all h families of rHCOs. It is easy to see that the group with steep slopes had the smallest number of families (curves). Most families with many (six to eight) members were in the medium-sensitivity group, otherwise the low-sensitivity group prevailed. For each sensitivity group (high, medium, and low), the plots in Figure 8B show the ranges of each parameter for h families with eight members. Plots in Figure 8C show the parameter ranges for the h families with four members (Fig. 7B, family curves). In the plots of Figure 8, B and C, the size of each point shows how many distinct families of the respective family group have a certain parameter value. For better visualization, the points in Figure 8B have been scaled up five times because there was a small number of families for many parameter values, and black diamond shapes show the median values for each parameter.
For the high-sensitivity group (Fig. 8 B1,C1), the reversal potential E leak showed a full range. For h families with many members (six to eight members), E leak (Ϫ55 or Ϫ60 mV) showed the greatest number of families; for the rest of h families (with three to five members), a hyperpolarized E leak (Ϫ65 or Ϫ70 mV) showed a slight increase in the number of families. g Leak showed also a full range of its values, with an increase in the number of families as g Leak decreases. g K2 had values between 25% and 175%, which meant that K2 must be present for rHCO h families to have high h sensitivity. An increased amount of g K2 seemed to increase the number of families, with a maximum reached at 150% (median, 125%). g P had a range similar to g K2 , with values between 25% and 175% (median, 75%). g CaS showed a full range of its values for h families with three members, and a range between 25% and 150% for the other h families. More interestingly, for g CaS values between 25% and 75%, we observed the most h families, while for higher values we observed just a few families (peak at 50%, which was also the median value). g SynS seemed to be needed (at least 25% present) to have full stability (eight member families); otherwise, it showed the whole range of values (e.g., families with four members; Fig. 8C1). g SynG and g h showed the full range, with an almost equal number of families at each value.
For the medium-sensitivity group (Fig. 8B2,C2), the reversal potential E leak showed a full range for all h families except the ones with eight members, for which the range was from Ϫ60 to Ϫ50 mV. There was a monotonic decrease in the number of families as the reversal potential becomes hyperpolarized. g Leak showed a full range for its values for all h families except those families with eight members, for which the range was from 25% to 175%. g K2 had values between 25% and 175%, except for one case of a family with five members, which meant that K2 must be present to have rHCOs h families with mediumsensitivity. An increased amount of g K2 seemed to increase the number of families, with a maximum reached at 150%; the median value was at 125%. g P had a range similar to that of g K2 , with values between 25% and 175% (median, 75%). g CaS showed a range between 25% and 175% for the h families with three to seven members, and a range between 25% and 100% for the families with eight members. g SynS seemed to be needed (at least 75% present) to have full stability (eight-member families), otherwise it showed a whole range of values. g SynG and g h showed the whole range of values for all families of different sizes.
For the low-sensitivity group (Fig. 8B3,C3), the reversal potential E leak showed a full range for families with fewer than seven members, and a range between Ϫ60 and Ϫ50 mV for families with many members (seven and eight members). More interestingly, for this sensitivity group there was a monotonic decrease in the number of families as the reversal potential becomes hyperpolarized. g Leak showed a full range for families with fewer than seven members. g K2 had values between 25% and 175% for h families with six to eight members, and a full range for the remaining h families (three to five members), although analysis showed a very small number of families where there is no K2 (Fig 8C3, this point in the plot is too small to see and represents one family). An increased amount of g K2 seemed to increase the number of families, and the median value was at 125%. g P had values between 25% and 175% for families with five to eight members, and full range for the remaining h families; the median value was 75%. g CaS showed a full range for families with three members, and a range between 25% and 175% for the remaining h families. The median values were at 100% for families with three to five members and at 75% for the others (six to eight members). g SynS seemed to be needed (at least 75% present) to have full stability (eight-member families), otherwise it showed the whole range of values (data not shown). g SynG and g h showed the whole range for all families of different sizes.
There is no set or range of parameters that we can detect that characterizes the period sensitivity (period vs g h ). Moreover, we observed no linear correlations of parameters in any of the period sensitivity groups. For example, we applied principal component analysis (PCA) to each of the three sensitivity subgroups of eight-member h families of rHCOs in Figure 8, A and B (data not shown), and found no linear correlations.
There are, however, a few observations that we can make from the plots in Figure 8 about how parameters contribute to period sensitivity. The graded synapse g SynG seems to not have any influence on period sensitivity as varying it does not change the number of families (constant). g K2 , g P , and g CaS seem to be needed to produce rHCOs (values of at least 25%). As expected, E leak and . Lines connect adjacent members of an h family forming a curve. Curves were split into three groups according to their period sensitivity (slopes), as follows: high (steep slopes; 1); medium (medium slopes; 2); and low (shallow slopes; 3). B, Parameter values vs the number of distinct families for all h families with eight members separated into these three groups of sensitivity. The size of each point quantifies how many different families have that parameter value (here, we scale up the sizes by five times for a better display). Black diamond shapes show the median values for each parameter. C, Parameter values vs the number of distinct families from all h families with four members separated into the three groups of sensitivity. For simplicity, plots do not show data for g SynG since the graded transmission did not have a significant influence on burst characteristics (see Sensitivity of period to variations of g's, and Sensitivity of spike frequency to variations of g's).
g Leak seem to affect activity in the same way; the number of rHCOs increases for high sensitivity or decreases for low sensitivity for both parameters at the same time. g K2 and g P also affect activity but differently, as it seems that one compensates for the effect of the other; this result (correlated triplet of g Leak , g K2 , and g P ) was observed (Doloc-Mihu and Calabrese, 2014) in regular bursting neurons (isolated HCO cells). The amount of g CaS does seem to have an effect on sensitivity type. There are few highsensitivity h families with three members with zero g CaS and few low-sensitivity h families with three members with zero g K2 and g P . Except for h families with three members where it shows the full range from 0% to 175%, it seems that the amount of g CaS shows the smallest range for the high-sensitivity group (25-75%) for families with eight members, and 25-150% for families with four to six members, with most families having values between 50% and 100%. Except for families with eight members, the lowand medium-sensitivity families show a range of 25-175% for g CaS . For the high-sensitivity families with different numbers of members, the median is consistently at 50% of g CaS ; for medium-sensitivity families, it is at 50% for h families with many members (seven and eight), and at 75% for the remaining h families (3-6 members); and for the low-sensitivity families, it is either at 75% or at 100% of g CaS .

Analysis of family structure helps identify ensembles of model instances for modeling studies
Our family analysis can help to identify parameter sets that approximate the physiological activity of any neuronal or network model more closely than hand-tuning methods and more efficiently than search methods, such as genetic algorithms (Holland, 1992). We propose selecting these robust parameter sets from the noninterrupted families with a large number of members (more than four members); each family identifies one such potential set of parameters. Our first step for identifying a good ensemble of model instances to best approximate the HCO physiological bursting was to select from the set of noninterrupted large member families a subset that have nonzero values for each parameter (i.e., g 0) other than the defining parameter. For example, there were 119 such families from the 144 h families of realistic HCOs with eight members. Then, for our next step we selected from the above subset only those families that include members whose corresponding isolated neurons show a specific activity type.
For example, in Figure 9, A1 and B1, we show (in red) two h families of realistic HCOs selected from the families with eight and five members, respectively. We selected these families such that the isolated neurons corresponding to their members have spiking and realistic bursting, respectively. In Figure 9, A2 and B2, we show the voltage traces of both neurons for some members of these two h families; 22 s was represented in each trace. One can see how varying g h influences the bursting activity. The first family (Fig. 9A) shows a moderate sensitivity to variations of g h , and the second one (Fig. 9B) shows a high sensitivity (busting speeds up fast by modulating h).
The voltage traces shown in Figure 9, A2 and B2, illustrate activity closely corresponding to leech physiological HCO activity (period between 5 and 15 s; average spike frequency between 8 and 25 Hz; duty cycle between 50% and 70%), as indeed each member of a rHCO family must. Model instances from the family in Figure  9A2 are of particular interest to us since they have the strong synaptic inhibition, observed experimentally. By having a duty cycle Ͻ50% (49.38%), the canonical model of Hill et al. (2001) would not be included in our realistic HCO group but in the functional HCO group (it belongs to the set of h families of fHCOs with six members). Like the canonical model of Hill et al. (2001), our model instances from Figure 9A2 have isolated neurons with spiking activity type.
All the h families of rHCO model instances from Figure  9A1 thus better approximate the physiological activity of the leech HCO system than the canonical model of Hill et al. (2001), because every member of every family conforms to the physiological activity. Similarly, other large h families of rHCOs (Fig. 9B1, four-member families) conform to the physiological activity. These families thus comprise a useful ensemble of model instances for exploring the mechanisms of the physiological activity. We conclude that our family-based method is a good tool for identifying such ensembles.

Discussion
We used an HCO model that replicates the rhythmic alternating bursting of mutually inhibitory interneurons of the leech heartbeat CPG to investigate the robustness (or sensitivity) and modulability of bursting activity to maximal conductance (g ) variations. We systematically explored the parameter space of two groups from the 10.5 million instances comprising the entire HCO model space, whose characteristics were previously recorded into a database, as follows: rHCOs with 99,066 instances; and fHCOs with 1,103,073 instances (Doloc-Mihu and Calabrese, 2011). To analyze such a large number of instances, we developed a new method (family based) that allowed us to classify instances from these two large groups (rHCOs and fHCOs) into much smaller subgroups of instances called families. Families are groups of model instances that differ only by the value of the one parameter that defines the family. By examining family structures and patterns, we developed new measures of robustness (or sensitivity) of bursting activity to changes in model parameters and investigated how such variations can be harnessed to modulate activity characteristics without compromising robustness. These measures can be easily adapted to other forms of electrical activity and other parameter types, and applied to brute-force databases. Our decision to limit our parameter variations to variations of g values while holding constant other kinetic parameters of the ionic currents was based on our available voltage-clamp analyses, which used average data to determine activation/inactivation and temporal characteristics of the currents. This decision limits the impact of our analysis (Amendola et al., 2012) but not the ability to apply our method to other systems where varying these parameters seems appropriate.

Robustness measures
We investigated independently the role of each conductance in the robust maintenance of functional bursting activity. Families organize model instances (rHCOs and fHCOs) into small subsets of instances sharing the same seven parameter values. Using families, we defined ro-bustness with a number of different measures and analyzed all families of each parameter with respect to these measures to find potential patterns. Our first measure was based on the family size: the greater number of these large families, the greater robustness. rHCO families with many members (seven or eight members of possible eight members) all maintain realistic bursting activity over a large range of values of the parameter. By this measure, realistic HCO activity was robust to changes in g h and g SynG , but it was very sensitive to changes in g P , and to a lesser extent to changes in g K2 . Our previous results (Doloc-Mihu and Calabrese, 2014) show that isolated neurons that are realistic bursters had only CaS families with one to three members, meaning that the activity of the realistic bursters was quite sensitive to changes in the value of g CaS . By adding inhibition (HCO configuration), the number of CaS families with many members (one to six) increased, making the system less sensitive to variations in the value of g CaS . However, the small number (four) of CaS families with six members for rHCOs suggested that the system is not as robust to variations in g CaS values as it is to variations in g h values.
Our second measure of robustness was based on family sequence (the sequence of values of the varied parameter within the family), as follows: a noninterrupted family sequence shows robustness (not broken by changing the defining family parameter value over its permissible range) and an interrupted family sequence shows sensitivity (susceptible to a change in activity type or nonphysiological burst characteristics by changing the defining family parameter values). We found that the higher the values of g h within the family, the more robust the sequence; thus, the higher likelihood of finding realistic HCO activity. Conversely, the lower the values of g h within the family, the more sensitive the sequence. Lower g Leak values (0 -100%) resulted in more robust realistic HCO activity (more noninterrupted families than for higher g Leak values). Stronger (Ն60 nS) spike-mediated synaptic transmission promoted robust realistic HCO activity. The system is quite robust to variations of graded transmission (SynG).
Our third measure of robustness was based on the activity type shown by the missing family members in both interrupted and noninterrupted families (i.e., both within and outside the permissible range): best if missing members have functional HCO activity, because by varying the parameter the HCO bursting is still kept functional. Plots of distributions of missing points from h families of rHCOs revealed that HCO bursting activity was robust to variations of g h , as variations of g h maintained fHCO bursting activity, though they might interrupt or terminate a realistic HCO family sequence. We conclude that g h is a potential target for modulation, since the HCO system shows robustness to its variations. The h families of the realistic bursters group (isolated neurons) were in general small (64 families with two and three members) with very few larger families (6 families with four members and 1 family with five members), with most missing members silent (data not shown; Doloc-Mihu and Calabrese, 2014). A change in g h has a significant chance of disrupting regular bursting activity in synaptically isolated neurons transforming their activity into silence, but by adding inhibition and forming an HCO, such variations of g h maintain functional HCO bursting activity. Based on these observations, it seems that in the heartbeat HCO inhibition makes the activity very robust to variations of g h .
Plots of distributions of missing points from CaS families of rHCOs revealed that HCO bursting activity was robust to variations of g CaS , if g CaS was Ն75% of the canonical value. Realistic HCO activity is robust with a low g Leak (0 -100%) present in the system (missing members have mostly fHCO activity); g Leak Ͼ100% favors spiking and silent activity. The distribution of missing points from P families of rHCOs revealed that changes in the amount of g P moved the HCO bursting activity outside the fHCO bursting range (to spiking, silent, or asymmetric activity). We conclude that g P is not good as a target for modulation since the activity type of the system is very sensitive to its variations. Changes in the amount of g K2 also disrupted HCO bursting activity, but not quite as severely as g P . Finally, distributions of the missing members from SynS families for the rHCOs show that at least 100% of g SynS was necessary to maintain fHCO bursting. Strong spike-mediated synaptic transmission promotes functional HCO activity.
For h families of rHCOs, most of the missing members (increasing percentage, from 62.3% for families with two members to 85.6% for families with seven members) have fHCO bursting activity for which only their duty cycle is outside the permissible physiological range (0.5-0.7). This result pertains to all h families of rHCOs (Table 3). The percentage of these missing members of h families that had burst periods either too long or too short to have rHCO activity (only their period outside the permissible range of 5-15 s) varied between 6.8% for families with seven members and 10.3% for families with five members. The number of missing members that spike too slowly (spike frequency Ͻ8 Hz) or too quickly (spike frequency Ͼ25 Hz) during the burst diminishes (from 8.7% for families with two members to 2% for families with seven members). For all h families, Ͻ0.3% of missing members from h families of rHCOs went outside the permissible realistic HCO range for all three criteria.
A corollary of our measures of robustness, which directly indicates the suitability of a parameter for functional modulation of the bursting activity, was based on a reliable and predictable (monotonic) variation in burst characteristics (period and spike frequency) within the realistic and functional range when a parameter was individually varied. Several clear patterns have emerged for rHCOs. Increasing g h increased the spike frequency moderately and decreased the period, speeding up the realistic HCO bursting activity. Increasing g CaS increased the spike frequency more strongly and increased the period. Increasing g P strongly increased the spike frequency and the period. For many K2 families, increasing g K2 decreased the spike frequency and the period. For most Leak families, increasing g Leak decreased the spike frequency and the period. Increasing g SynG showed negligible changes in spike frequency and period, and increasing g SynS decreased the spike frequency and increased the period. All these results on the influence of parameter variation on burst characteristics confirmed and extended the previous results in the study by Hill et al. (2001) obtained in a more restricted parametric space (variation of a single parameter over a range similar to the one presented here but on a background of canonical values for all other parameters). So, to decrease period, for example, we can increase the amount of g h , g K2 , and g Leak or decrease the amount of g CaS , g P , and g SynS . If we consider the other measures of robustness and this measure of modulatory effectiveness, then varying g h appears ideal for modulating period because the HCO system is very robust to the variation in g h and g h consistently and substantially modulated period with minimal adverse effects on spike frequency or duty cycle.
For a specific (chosen) parameter, t, we propose the following simple formula as a general measure of robustness for large databases: where X t, n is the number of families with more than a chosen number of members (n, best when n is large); Y t, n is the number of noninterrupted families with more than n members; Z t, n is the number of interruptions (missing families members) that do not change the activity type; and w X , w Y , w Z are chosen weights that show the importance of each individual robustness measure in assessing the final robustness of the parameter t, with w X ϩ w Y ϩ w Z ϭ 1. The formula is a simple weighted cumulative sum of the three measures of robustness that we proposed above. The weights indicate the importance of each measure of robustness within the system and depend on what the user values as robust in their system.
In our study here, we chose to assess the robustness for the realistic HCOs. We selected as a parameter t, the maximal conductance of the hyperpolarization-activated cation (h) current, g h . For n ϭ 4, our numbers for calculating total robustness follow. First, we selected from our rHCOs all h families with more than four members and obtained the first measure of robustness, X h, 4 ϭ 4,669. The total number of noninterrupted families is Y h ϭ 11,210, and the number of interruptions that keep the functional bursting activity, fHCO is Z h ϭ 98,882. Next, Y h, 4 ϭ 1,729, and Z h, 4 ϭ 8,945. Finally, we chose to assess robustness by giving equal importance to our first two measures and by not considering the third one (i.e., w X ϭ w Y ϭ 0.5, w Z ϭ 0, and R h, 4 ϭ 3,199). To put this number in perspective, it can be normalized to the user's goal. For example, normalizing to the total number of families with more than four members shows that ϳ70% (0.685) of the h families with more than four members maintain realistic alternating bursting activity when g h varies (i.e., no interruptions).

Robustness and parameter correlations
We emphasize that we previously found no linear correlations among parameters in the rHCO group (Doloc-Mihu and Calabrese, 2014). Figure 1 reveals a nonlinear correlation among E leak , g Leak , g K2 , and g P for the rHCOs (c.f. Goldman et al., 2001), but no other nonlinear correlations were observed. Nor were there any linear correlations (via the PCA method) or other nonlinear correlations in any subgroup tested. For example, we applied PCA to each of the three sensitivity subgroups of eight-member h families of rHCOs of Figure 8, A and B (data not shown), and found no linear correlations.
How then would parameter correlations affect robustness as determined by our cumulative measure? Our previous work showed that for the component neurons of the HCO to be endogenous bursters, there must be a strict linear correlation among g Leak , g K2 , and g P . We enforced this correlation in our database by requiring endogenous bursting of the component neurons in all families considered. We then used the same measures as above and obtained X h, 4 enf ϭ 53, Y h, 4 enf ϭ 51, Z h, 4 enf ϭ 2, and R h, 4 enf ϭ 52. These numbers indicate that this correlation, should it be biologically enforced, would limit the robustness. This finding is consistent with the observation that endogenous bursting in heart interneurons is very sensitive to changes in leak (e.g., as caused by sharp micro-electrode penetration) and that it is not necessary for robust alternating bursting activity (Sorensen et al., 2004) or when h current is modulated (Tobin and Calabrese, 2006). Moreover, the vast majority [94,487 (95.37%)] of rHCOs in our database are made up of component neurons that are spiking. It is interesting to note that such biologically enforced correlations have been observed in the stomatogastric nervous system (Goaillard et al., 2009;Tobin et al., 2009) and cardiac ganglion of crustaceans (Ball et al., 2010). In the cardiac ganglion at least, such correlations appear to increase robustness.

Period sensitivity
We separated our h families with eight members into three groups according to their period sensitivity to increasing g h (high, medium, and low sensitivity). For high-sensitivity families, increasing g h speeds up bursting strongly; a large decrease of period with increasing g h occurs (typically before 50% g h ), then this decrease moderates at higher g h values. For low-sensitivity families, increasing g h speeds up bursting more uniformly; period decreases moderately but almost linearly for all g h values. The medium-sensitivity families are intermediate; period decreases steadily with increasing g h , but there is neither a sudden drop nor a range of weak period decrease. This splitting seems functional to us and shows how parameters interact within specific ranges to produce these different types of sensitivities.
We put in the medium-sensitivity group all those families that were in neither the high-sensitivity nor the lowsensitivity group. One could argue that defining this group of medium sensitivity might not be germane to the analysis, because the set of h families illustrates period curves whose slopes (with respect to the horizontal axis) occupy the entire spectrum of angles. The decision to keep three groups versus two groups (high and low-sensitivity) was based on visual inspection of 100 randomly selected slopes. This process helped us to set the criteria for slope angles of the two important cases of low and highsensitivity. The results show that the region with the steepest slope (high-sensitivity region) can appear anywhere, but in ϳ98% of the cases it appears at the very beginning of the curve [i.e. at low values of g h (Ͻ100%)]. For the medium-sensitivity group, a steep slope region occurs at higher g h (100 -150%) or not at all. For the low-sensitivity group, there is no steep slope region, only a nearly constant slope. Note that all curves show monotonically decreasing periods with the increase in g h .
Splitting the families into three sensitivity groups, moreover, helped us to define how other parameters define period sensitivity to the variation of g h by analysis of the ends of a spectrum. We used this analysis of h families to investigate the effects of background conductances on the period of realistic HCO bursting. To sum up this period sensitivity analysis, several parameters influence the realistic HCO burst period. Shorter periods seem to be influenced by K2 and P working together against Leak, with CaS also having some effect. Longer periods seem to be affected by K2 and E leak working together against P, by the amount of CaS, which has the same influence as K2, and also by the spike-mediated synaptic transmission SynS. From this analysis, we conclude that period is influenced by groups of parameters but it is easily and predictably controllable by modulating g h .
It appears that in the HCO system there are several different mechanisms that influence robust realistic bursting. One mechanism that we found involves E leak , I Leak , I P , and I K2 working together to compensate for the variations in each to keep bursting functional. In isolated bursting neurons, these parameters interact in a linearly correlated way (Doloc-Mihu and Calabrese, 2014). Linking these isolated bursting neurons by mutually inhibitory synapses into HCOs increases the system robustness and causes these parameters to interact in a nonlinear way.
Another mechanism we found involves specific parameters working individually to change the burst characteristics (here, we focused on burst period and spike frequency). Our results show that realistic, physiological bursting activity is robust to changes in the amount of h current, and that h is a great target for modulating period (because varying the amount of g h changes the burst period in a significant and predictable way). This result confirms and supplements the conclusions of previous theoretical studies (Nadim et al. 1995;Hill et al., 2001) and extends them to combinations of key parameters that are varied together. Moreover, the peptide myomodulin upmodulates the h current in heart interneurons that comprise HCOs and consistently speeds the burst period Tobin and Calabrese, 2006). Recent experiments (Tsuno et al., 2013) show that I h is a good modulation target for the phase of neurons in rat cortex, emphasizing the general importance of this current for modulation.

Applicability of the robustness measures
Last, we emphasize that the measures of robustness we developed here are easily adaptable to other neuronal and network models. Brute-force databases, such as HCOdb, lend themselves to the family analysis we performed here. By identifying families and applying our measuresthe number of large families, the number of noninterrupted families, and the missing family members that show functional albeit not physiological activity (i.e., no change in activity type)-one can identify parameters that can be safely modulated. Then, using large families one can identify those parameters that consistently modulate an adaptable and desirable activity characteristic. Family analysis can also allow the identification of robust parameter sets that more closely approximate physiological activity of any neuronal or network model, and thus allow the construction of an ensemble of physiological model instances for further mechanistic studies. For example, the ensemble of physiological model instances with large noninterrupted families we have identified will allow us to ask whether for effective, functional modulation there are advantages to covarying currents over a single particularly efficacious current like h in the leech heartbeat HCO.