Bifurcation analysis of motoneuronal excitability mechanisms under normal and ALS conditions

Abstract

Introduction:

Bifurcation analysis allows the examination of steady-state, non-linear dynamics of neurons and their effects on cell firing, yet its usage in neuroscience is limited to single-compartment models of highly reduced states. This is primarily due to the difficulty in developing high-fidelity neuronal models with 3D anatomy and multiple ion channels in XPPAUT, the primary bifurcation analysis software in neuroscience.

Methods:

To facilitate bifurcation analysis of high-fidelity neuronal models under normal and disease conditions, we developed a multi-compartment model of a spinal motoneuron (MN) in XPPAUT and verified its firing accuracy against its original experimental data and against an anatomically detailed cell model that incorporates known MN non-linear firing mechanisms. We used the new model in XPPAUT to study the effects of somatic and dendritic ion channels on the MN bifurcation diagram under normal conditions and after amyotrophic lateral sclerosis (ALS) cellular changes.

Results:

Our results show that somatic small-conductance Ca²⁺-activated K (SK) channels and dendritic L-type Ca²⁺ channels have the strongest effects on the bifurcation diagram of MNs under normal conditions. Specifically, somatic SK channels extend the limit cycles and generate a subcritical Hopf bifurcation node in the V-I bifurcation diagram of the MN to replace a supercritical node Hopf node, whereas L-type Ca²⁺ channels shift the limit cycles to negative currents. In ALS, our results show that dendritic enlargement has opposing effects on MN excitability, has a greater overall impact than somatic enlargement, and dendritic overbranching offsets the dendritic enlargement hyperexcitability effects.

Discussion:

Together, the new multi-compartment model developed in XPPAUT facilitates studying neuronal excitability in health and disease using bifurcation analysis.

1. Introduction

Bifurcation analysis of neural models is a useful tool for studying the steady-state, non-linear dynamics of neurons and the characteristics of their firing output throughout the stimulus input range. While bifurcation analysis has greatly contributed to our understanding of the neuronal non-linear dynamics underlying membrane oscillations and cell firing under normal (Jalics et al., 2010; Zhou et al., 2018), pathological (Glass and Mackey, 1979; Kaper et al., 2013), and pharmacological (V-Ghaffari et al., 2017) conditions, their use has been limited to single-compartment models (Jalics et al., 2010; Wang et al., 2019; Xing et al., 2022). However, such models lack the non-linear behaviors that arise from dendritic channels (Hounsgaard and Kiehn, 1993; Lee and Heckman, 1999; Li and Bennett, 2007; Carlin et al., 2009). While AUTO (Doedel et al., 2007) and MatCont (Dhooge et al., 2003) are the primary bifurcation tools used in the dynamical systems literature, XPPAUT and its integrated module AUTO (Ermentrout and Mahajan, 2003) are the more common bifurcation analysis software used in the neuroscience literature. Because the number of ordinary differential equations of a high-fidelity neuronal model with multi-compartments and multiple ion channels could be large and XPPAUT—unlike NEURON—lacks tools that abstract the mathematical detail to facilitate model development, the process of developing a multi-compartment model in XPPAUT is cumbersome (N.B., no single multi-compartment model of a neuron has been developed in XPPAUT to date). Thus, extending bifurcation analysis to high-fidelity models with multiple compartments and dendritic channels, which mediate many non-linear firing behaviors in MNs, has been limited. As multi-compartment models are more accurate in simulating the firing behaviors of spinal MNs than reduced models (Hendrickson et al., 2011; Elbasiouny, 2014), bifurcation analysis of more complex models and their behaviors is, therefore, of great importance.

The goals of this study are to (1) develop a multi-compartment computer model that simulates in XPPAUT the non-linear behaviors as empirically measured in spinal MNs, and (2) use the model in XPPAUT to assess the role of somatic and dendritic ion channels in regulating the cell’s repetitive firing using bifurcation analysis under normal and disease conditions. As XPPAUT does not support large cell models with 3D anatomical detail, our first step was to reduce the cell model published by Mousa and Elbasiouny (2020) while preserving as much accuracy in simulating this MN’s firing behaviors as possible. The outcome of this step was a six-compartment (6C) model developed in XPPAUT, which we assessed to verify its electrical properties versus experimental data. To the authors’ best knowledge, this model is the first multi-compartment neuron model with somatic and dendritic channels developed in XPPAUT.

Our results showed that, under normal conditions, somatic small-conductance calcium-activated potassium (SK) channels, which mediate the afterhyperpolarization (AHP) phase of the action potential (AP), support the MN in generating healthy, large-amplitude AP spikes and stable rhythmic cell firing over an extended input current range. Importantly, the presence of somatic SK channels resulted in the emergence of a subcritical Hopf bifurcation node to replace the supercritical Hopf node at the end of the voltage-current (V-I) bifurcation relationship. This bifurcation point results from the negative feedback control between the somatic SK and N-type Ca²⁺ (CaN) channels. Our results also showed that dendritic L-type Ca²⁺ channels increase the cell excitability substantially and shift stable rhythmic firing on the V-I bifurcation relationship to negative currents, allowing the cell to fire continuously in absence of excitatory stimuli. However, dendritic SK channels regulated L-type Ca²⁺ channels activation toward normal levels, yet continued to enable rhythmic cell firing starting at low input currents. Using the new MN model in XPPAUT to study the impact that somatic hypertrophy and dendritic enlargement and overbranching have on cell firing in ALS, our bifurcation analysis showed that the dendritic enlargement has more notable net hypoexcitability effects on cell firing than does somatic enlargement. Interestingly, some hyperexcitability effects also arose from dendritic enlargement, however dendritic overbranching appeared to offset those effects. Together, this work reports the first multi-compartment MN model developed in XPPAUT and makes it available for the scientific community to use. Additionally, our results report the first bifurcation analysis conducted on a neuron model with this higher level of anatomical and ion channel detail under normal and disease conditions. Further, our bifurcation analysis describes novel non-linear behaviors, as well as the roles different ion channels play in regulating cell firing.

2. Materials and methods

In the present study, we used a fatigue-resistant (FR) cat MN model of the medial gastrocnemius (MG) muscle, which was modeled in high detail in Mousa and Elbasiouny (2020). This high-fidelity model incorporates known non-linear firing behaviors of MNs. This model was developed in the NEURON software (Hines and Carnevale, 1997), which does not feature bifurcation analysis. As XPPAUT (Ermentrout and Mahajan, 2003), and its bifurcation analysis AUTO module, do not have tools to facilitate the development of high-fidelity neuronal models with detailed anatomy and multiple ion channels, bifurcation analysis of the full 3D model of Mousa and Elbasiouny (2020) is infeasible, due to the hundreds of balance and ordinary differential equations to be written in XPPAUT. To overcome that, the first goal of the present study was to reduce the 3D Mousa and Elbasiouny (2020) model into a simpler multi-compartment model with all non-linear MN behaviors retained and implemented in XPPAUT, where bifurcation analysis is possible. As Auto-07p usage in neuroscience is still limited and XPPAUT’s time domain simulations were needed for the model verification process and comparison against the NEURON simulations, we used XPPAUT (version 8.0) and AUTO 2000 in our simulations. All figures were constructed using Python3 (Van Rossum and Drake, 2010), Matplotlib (Hunter, 2007), Seaborn (Waskom, 2021), and Pandas (The Pandas Development Team, 2022).

2.1. Model description

Following the methodology of Elbasiouny (2014) in reducing 3D models into reduced models that still demonstrate highly accurate firing behaviors, we reduced the high-fidelity computer model of Mousa and Elbasiouny (2020) into a six-compartment (6C) model (Figure 1A, 6C model morphology, and electrical parameters are shown in Tables 1, 2, respectively). This 6C model was able to simulate experimental data and the 3D model behaviors with acceptable accuracy (Table 3). In the reduction process, the simplified 6C model was designed to resemble the 3D model in area distribution (Figure 1C), electrical distance distribution (Figure 1B), and electrical properties, which were optimized to fall within the 95% confidence interval of experimental data [passive and active membrane properties, rheobase, action potential (AP) and afterhyperpolarization (AHP) properties, frequency-current (FI) relationship properties] and to match the 3D model properties as much as possible.

FIGURE 1

Mousa and Elbasiouny (2020) (blue). (C) A comparison between the cumulative area and physical path distance of the model compartments in the reduced 6C model in XPPAUT (red) and the 3D model of Mousa and Elbasiouny (2020) (blue). In panels (B,C), the x-axes are normalized to the longest physical path distance of each model.

Similar to the 3D model of Mousa and Elbasiouny (2020), the 6C model was formed of initial segment/axon hillock (IS/AH), soma, and four dendrites (Figure 1A). All compartments of the 6C model had passive leak channels. The soma had the following additional channels: fast sodium (Naf), delayed rectifier potassium (Kdr), N-type calcium (CaN), and small conductance calcium-activated potassium (SK_AHP) channels. The IS/AH had the following channels (in addition to leak channels): Naf, persistent sodium (Nap), and Kdr. Dendrites 0 and 1 had only leak channels, whereas dendrites 2 and 3 also contained L-type Ca²⁺ (CaL) and small conductance calcium-activated potassium (SK_L) channels. Importantly, the location and conductance of dendritic CaL and SK_L channels were retained in the 6C model, as in the 3D model of Mousa and Elbasiouny (2020), to generate comparable Ca persistent inward currents (Ca²⁺ PIC). Specifically, the 3D model had dendritic regions of high (26% of the dendritic area located between 0.44λ and 0.6λ from the soma) and low (10% of the dendritic area located between 0.6λ and 1.1λ from the soma) Ca²⁺ PIC and SK conductances. To mimic that, the distal dendrites of the 6C model were simulated with two compartments, one with high Ca²⁺ PIC and SK conductances (dendrite 2, with 26% of the dendritic surface area and distance between 0.45λ and 0.55λ from the soma) and another with low Ca²⁺ PIC and SK conductances (dendrite 3, with 10% of the dendritic surface area and distance between 0.55λ and 0.93λ from the soma). In this way, the reduced 6C model had dendritic conductances very similar to those of the 3D model of Mousa and Elbasiouny (2020).

All ion channels were conductance-based, following the Hodgkin and Huxley formalism (Hodgkin and Huxley, 1952), such that activation and inactivation states are voltage-dependent, except for SK channels which use saturation function that depends on the calcium concentration in the compartment (Mousa and Elbasiouny, 2020). The cell membrane specific capacitance was 1 μF/cm². All model equations are listed in the Supplementary material.

2.2. Bifurcation diagrams

In this paper, the V-I bifurcation curves have four colors: (1) red traces indicate stable equilibrium points, (2) black traces indicate unstable equilibrium points, (3) green traces indicate stable periodic solutions; the maximal and minimal value of the green curve at a single current value refers to the maximal and minimal amplitude of the stable oscillation at this current value, and (4) blue traces indicate unstable periodic solutions. In the bifurcation diagram, the firing range is the range where stable periodic solutions exist.

2.3. XPPAUT simulation settings

We used the following parameters in XPPAUT: dt is 0.025 ms, meth = backeul, parmax = 1, parmin = −1, NPR = 100,000, NMAX = 10,000,000, NTST = 100. Regarding parmax and parmin, the maximal injected current during all simulations is 50 nA, so a scaling factor is used to map the current range to the parameter range. NTST can be higher in some simulations to avoid mx (failed to converge) error in XPPAUT. The simulation ran using a personal computer with 8 Gb RAM, Intel Core i7-5500 U CPU, and Ubuntu 20.04 OS. Long pulses are generated using unit step pulses, and FI curves are generated using ramp pulses with a slope of 4 nA/s and a maximal current injection of 25 nA.

3. Results

3.1. Development and verification of the MN model in XPPAUT

To study the role of different ion channels in regulating MN repetitive firing, bifurcation analysis of a MN model with sufficient anatomical and ion channel details in XPPAUT was needed. As XPPAUT does not facilitate simulations of large cell models with 3D anatomical detail, the first goal of this study was to reduce the MN model published in Mousa and Elbasiouny (2020) into an equivalent one that could be implemented in XPPAUT for bifurcation analysis. The Mousa and Elbasiouny (2020) model was chosen as the reference 3D model because it is currently the cell model that simulates the most non-linear electrical and firing behaviors of spinal MNs (e.g., firing bistability, voltage hysteresis, synaptic amplification, PIC warm-up, dendritic staircase currents, dendritic persistent inward and outward currents). We reduced the 3D, 731-compartment anatomy of the Mousa and Elbasiouny (2020) model into an equivalent 6-compartment (6C) model that could be implemented in XPPAUT (Figure 1A). Six were the fewest compartments that could match most of the 3D model’s electrical properties within the 95% confidence interval of experimental data (see Table 3). The 6C model’s morphological parameters were set to replicate the 3D model dendritic area spatial distribution, such that (1) distal dendrites contain CaL and SK_L channels that mediate the Ca²⁺ PIC (see section “Materials and methods” for detail), (2) Electrotonic distance profile matches the original 3D model (Figure 1B), and (3) the cumulative distribution of the dendritic surface area matches the experimental 3D morphology of Cullheim et al. (1987) model (Figure 1C). Also, all channel kinetic parameters were preserved and translated into the XPPAUT environment to develop a 25-state-variable model. To confirm that the 6C MN model was implemented successfully in XPPAUT, a NEURON version of the 6C reduced model was also implemented and the results of both were compared and confirmed to be identical. The 6C model in XPPAUT was then used for bifurcation analysis using the Auto module. With somatic/axon hillock/dendritic compartments, 25 state variables and multiple ion channels, this cell model is the most comprehensive model implemented in XPPAUT to date. This enabled us to perform bifurcation analysis of one neuron with unprecedented details. The following sections present the bifurcation analyses of the 6C model in XPPAUT, in which we added one compartment and one ion channel at a time to examine the effect of each channel on cell firing.

3.2. Somatic SK channels support stable rhythmic cell firing over extended input range

To examine the effect of somatic ion channels on cell firing, we conducted the first bifurcation analysis on a single compartment model composed of only the soma of the 6C model in XPPAUT, with only its somatic leak, Na, and K channels included (Figure 2A). These are the minimal channels needed for generating an AP. We started with this single compartment soma model to serve as a comparison of the following bifurcation graphs, in which additional channels are added to the model one at a time. When bifurcation analysis of the single compartment soma model was conducted, the typical V-I bifurcation graph of a Hodgkin–Huxley model with a firing behavior starting and ending between subcritical and supercritical Hopf bifurcation points, respectively, was obtained (Figure 2C). In this bifurcation diagram, the subthreshold membrane depolarization with no cell firing at low input current is shown by the increasing red trace on the left until rhythmic cell firing is evoked at the cell rheobase (shown by the green traces). As Na channels inactivate gradually with increasing input current, the height of AP decreases gradually (shown by the decreasing limit cycles voltage amplitude in Figure 2C) until cell firing dampens completely (shown by the red trace on the right in Figure 2C).

FIGURE 2

When CaN channels were added to the soma model (Figure 2B), a similar voltage amplitude (but with a narrower limit cycles’ range) and a bit smaller current range were seen in the bifurcation graph (Figure 2C). The decrease in limit cycles voltage amplitude and current range is due to the additional depolarization provided by the CaN channels, which reduced the AP amplitude more. However, when the somatic SK channels were next added to the soma model (Figure 3A), which mediate the spike AHP (Figure 3B, blue traces), the limit cycles voltage amplitude and current range increased substantially, and the supercritical Hopf node at the end of the firing range was replaced with a subcritical node (Figure 3E, the “with SK” trace). Given the negative feedback control exerted on the membrane potential by the somatic SK and N-type Ca (CaN) channels, the subcritical bifurcation node in the V-I bifurcation graph reflects this membrane potential non-linearity.

FIGURE 3

In the time domain, the negative feedback control provided by somatic SK channels on the membrane potential was evident. For instance, without somatic SK channels included, the frequency-current (F-I) relationship in response to a triangular current ramp reached very high, non-physiological firing rates, and had a steep slope (Figure 3C, red trace). Also, long current pulses evoked dwarf APs firing at a high rate (Figure 3D, red trace). However, the addition of SK channels regulated the firing of the neuron model and brought the F-I relationship slope and long pulse firing down to the physiological firing range (blue traces in Figures 3C, D). Interestingly, at 40 nA, the model without SK channels had fully inactivated Na channels with no cell firing (Figure 3D, lower panel). However, the AHP mediated by SK channels helped relieve some Na channels from inactivation and evoked cell firing, which is also reflected in the bifurcation diagram at 40 nA (Figure 3E, there are no limit cycles for the model without SK at 40 nA). Collectively, somatic SK channels support the generation of full (large amplitude) and stable rhythmic cell firing over an extended input current range and generate a subcritical Hopf bifurcation point in the V-I bifurcation graph.

3.3. Dendritic L-type Ca²⁺ channels enable MN firing at lower currents

According to Li and Bennett (2007), L-type Ca²⁺ and SK channels are present on the dendrites of MNs, which were simulated in Mousa and Elbasiouny (2020). The currents of both channels summate to produce Ca²⁺ persistent inward currents (Ca²⁺ PICs). In this section, we first investigate the effects of the L-type Ca²⁺ channels alone. To that end, we added four dendritic compartments to the somatic section to form the complete 6C model (Figure 4A), then added L-type Ca²⁺ channels to the distal dendrites (dendrites 2 and 3 in Figure 4A) to replicate the 3D model of Mousa and Elbasiouny (2020). The somatic spike and AHP of the complete 6C model before and after adding the dendritic L-type Ca2⁺ channels are shown in Figure 4B. Without L-type Ca²⁺ channels, the model with passive dendrites showed a similar bifurcation diagram to the one with only a somatic compartment (compare the “with SK” trace in Figure 3E to the “Without CaL” trace in Figure 4E). However, the limit cycles were a little shortened in both voltage and current ranges, as the somatic current now activates a larger cell area with larger capacitance, which in turn causes a smaller AP height. The firing behavior during the current ramp (F-I relationship) was similar to that of the somatic compartment model, but with lower firing rates due to the larger cell surface area and capacitance.

FIGURE 4

When L-type Ca²⁺ channels were added to the passive dendrites (Figure 4E), the limit cycles were shifted in the negative current direction (by > 12 nA), allowing the cell to fire repetitively at zero current (i.e., without input). This behavior was also seen in the F-I relationship, where the cell continued to fire on the descending current ramp until zero current was injected (Figure 4C). Long pulse simulations showed similar results, where the cell continued to fire repetitively after the current pulse was terminated (Figure 4D, see the arrow showing self-sustained firing). The addition of L-type Ca²⁺ channels increased the MN maximal firing rate from 37.8 to 82.4 Hz on the current ramp (Figure 4C) and increased the average firing rate during long pulses (Figure 4D).

The bifurcation analysis shows that, at low current amplitudes (< 10 nA), two stable states exist in the MN model with dendritic Ca channels (Figure 4E, the “with CaL” trace): one non-oscillatory state representing the subthreshold membrane depolarization (the red line) and another oscillating state with limit cycles representing the cell firing (between the green lines). The presence of these two stable states indicates that a MN with dendritic Ca²⁺ PIC may fire or not depending on its initial conditions, which is reflected in the F-I relationship (Figure 4C, blue trace) with no cell firing in the range < 10 nA on the ascending firing segment but with repetitive cell firing in the same range during the descending firing segment. The temporal activation of Ca²⁺ PIC is seen during the long pulse simulations, in which a 30 nA current pulse evoked repetitive cell firing at the beginning of the pulse; then it changed to a non-oscillatory steady state membrane depolarization at the end of the pulse (due to full Na channels inactivation). In sum, L-type Ca²⁺ channels shift in the MN bifurcation diagram to lower currents, allowing the cell to fire at negative currents (i.e., in absence of excitation) and at higher firing rates.

3.4. Dendritic SK channels regulate L-type Ca²⁺ channels activation

Dendritic L-type Ca²⁺ channels co-localize with SK channels on the dendrites to regulate the amplitude of Ca²⁺ PIC (Li and Bennett, 2007). When SK channels were added to the model’s distal dendrites (Figure 5A), the bifurcation diagram shifted to the right, as it partially deactivated L-type Ca²⁺ channels (i.e., reduced the Ca²⁺ PIC) and stopped cell firing at less negative current (Figure 5E). This was observed in the F-I relationship (Figure 5C, blue trace), in which dendritic SK channels allowed graded activation of the Ca²⁺ PIC, causing the cell to fire at lower firing rates and for bistability to emerge (i.e., non-overlapping ascending and descending F-I segments in Figure 5C, blue trace). In absence of dendritic SK channels, L-type Ca²⁺ channels fully activate immediately, leading to Ca²⁺ PIC saturation upon cell recruitment and lack of firing bistability (i.e., overlapping ascending and descending F-I segments) (Figure 5C, red trace). The effect of dendritic SK channels is primarily on repetitive firing as they had negligible effect on single spike properties (Figure 5B).

FIGURE 5

The effects of graded Ca²⁺ PIC activation by dendritic SK channels were seen in the time domain long pulse simulations (Figure 5D, blue traces) which show no self-sustained firing after the termination of a 20 nA current pulse (as L-type Ca²⁺ channels were not fully activated), but with self-sustained firing present after the termination of a 30 nA current pulse (as L-type Ca²⁺ channels were fully activated by the very high amplitude pulse). In absence of dendritic SK channels, self-sustained firing was always present regardless of pulse amplitude (Figure 5D, red traces). Interestingly, and contrary to their somatic counterpart, dendritic SK channels did not extend the limit cycles or the cell firing range (as in Figure 3E). They only shifted the bifurcation diagram rightward, offsetting some of the Ca²⁺ PIC effects.

3.5. Examining mechanisms of MN excitability dysfunction in ALS

To illustrate the utility of the 6C model in XPPAUT in studying neuronal excitability dysfunction in neurodegenerative diseases, we examined the effect of pathological MN anatomical changes in ALS by conducting bifurcation analysis on the 6C model in XPPAUT under ALS conditions. Specifically, in the early stages of ALS, spinal MNs undergo anatomical changes, such as somatic (Shoenfeld et al., 2014; Dukkipati et al., 2018) and dendritic enlargement (Amendola and Durand, 2008), as well as dendritic overbranching (Filipchuk and Durand, 2012) in multiple ALS mouse models. While these individual cellular changes are expected to cause cell hypoexcitability, it is unknown whether these changes are of similar magnitude and impact on MN excitability. To examine the effects of somatic vs. dendritic size enlargement, we increased the soma size by 10% [based on experimental data of P10 SOD MN increases of 6–20%, (Dukkipati et al., 2018; Figure 6A)] and compared that to a 30% increase in the dendritic size [based on experimental data of P10 SOD MN increases of 25–90%, (Amendola and Durand, 2008)] via the addition of one dendritic branch (Figure 6B). In these simulations, all other model parameters, including the ionic conductances, were unchanged.

FIGURE 6

Interestingly, relative to the normal model, the increase in soma size shifted the bifurcation diagram to the right, causing the cell recruitment current to increase and the limit cycles’ height to also increase at any given current (Figure 7A, blue label). These hypoexcitability effects were confirmed in the time domain, as the AP height was increased and self-sustained firing was absent (Figure 7B, blue labels), recruitment current was increased, and firing bistability (measured as the difference between the ascending and descending FI relationship segments) was decreased (Figure 7C, blue labels), and cell firing rates were generally decreased (Figures 7C, D, blue labels).

FIGURE 7

Conversely, the increase in dendritic area had mixed effects on the cell excitability. For instance, relative to the normal model, the cell recruitment current was increased, the limit cycles’ height and range were reduced, and the subcritical Hopf bifurcation node in the bifurcation diagram was exacerbated, effectively shrinking the operational firing range of the cell significantly (Figure 7A, purple label). These are hypoexcitability changes whose effects were observed in the time domain as a reduction in the AP height (Figure 7B, purple label), lack of cell ability to fire repetitively due to Na channel inactivation (Figure 7B, purple label), and a significant increase in cell recruitment current (Figure 7C, purple label). Paradoxically, the increase in dendritic area increased the cell’s repetitive firing on the ascending F-I segment relative to the normal model (compare the purple to black segments in the left part of Figure 7D), which is a hyperexcitability change. This increase in ascending firing rates was due to the significant increase in cell capacitance due to the dendritic enlargement, leading to an increase in the cell time constant (i.e., slower activation of the cell); thereby causing Ca²⁺ PIC full activation before cell recruitment, which inactivated Na channels and reduced the cell firing range. Collectively, these results show that somatic and dendritic enlargements have different effects on MN excitability in ALS, with the dendritic enlargement having much stronger hypoexcitability effects on MN excitability than those from soma size increase, primarily because of the significant increase in cell capacitance by the dendritic enlargement.

To examine the effects of dendritic overbranching on cell excitability, we compared the ALS model when the dendritic area was enlarged by 30% via one (Figure 6B) vs. three dendritic branches (Figure 6C). Interestingly, although the dendritic area was increased by the same percentage in both ALS models and both models experienced hypoexcitability effects relative to the normal condition, the hypoexcitability effects were less in the model with dendritic overbranching (Figure 8). Specifically, relative to the one-branch ALS model, the three-branch ALS model had larger operational firing range (as the subcritical Hopf node in the bifurcation diagram shifted rightward, Figure 8A—orange label) and repetitive cell firing on long pulses and ramps was maintained longer (Figures 8B, D—orange label) and at less firing rate (Figures 8C, D—orange label). Together, these results show that dendritic overbranching attenuates the hyperexcitability effects of the dendritic enlargement, leading to less firing rates.

FIGURE 8

4. Discussion

The present study has several highly novel aspects. First, it presents the development of the first (to our knowledge) multi-compartment neuron model in XPPAUT (model code is publicly available for researchers to download and use¹). Second, using this novel model, we were able, for the first time, to study the effects of somatic and dendritic ion channels on the bifurcation diagrams of spinal MNs. Third, our bifurcation analysis reports novel bifurcation behaviors, such as the generation of a new subcritical Hopf bifurcation node, resulting from somatic SK channels. To demonstrate the utility of the novel model in XPPAUT in investigating neurodegenerative diseases, we used the model to study and contrast the separate effects that somatic and dendritic enlargements have on MN excitability in ALS. Our results show that, under normal conditions, somatic SK and dendritic L-type Ca²⁺ channels have the strongest effects on the bifurcation diagram of spinal MNs. Specifically, somatic SK channels extend the limit cycles range and generate a new subcritical Hopf bifurcation node in the V-I bifurcation diagram of the MN. L-type Ca²⁺ channels do not influence the limit cycles’ shape but shift them to negative currents (reflecting self-sustained cell firing in absence of input). Dendritic SK channels—a less powerful player than their somatic counterpart—lessen the Ca²⁺ PIC effects and shift the limit cycles back toward positive currents. Under ALS conditions, our bifurcation analysis shows that, while somatic and dendritic enlargements have net hypoexcitability effects on the cell, those from the dendritic enlargement are more drastic, reduce the cell operational firing range, and induce some hyperexcitability effects. Importantly, dendritic overbranching appeared to offset the dendritic enlargement hyperexcitability effects. On the other hand, somatic enlargement lowers the cell firing without compromising the cell firing range. Accordingly, the novel XPPAUT 6C model has demonstrated the potential to improve our understanding of how pathological changes impact neuronal excitability in neurodegenerative diseases.

4.1. Negative feedback loops are different

Neuronal firing is the outcome of many interacting non-linear membrane properties and ionic mechanisms, and bifurcation analysis provides a means to study and quantify the effect(s) each membrane property has on cell firing. As XPPAUT does not have tools to facilitate the development of high-fidelity neuronal models with detailed anatomy and multiple ion channels, bifurcation analysis has been traditionally conducted on reduced neuronal models of a single compartment with only a few ion channels (White et al., 1995; Wang et al., 2019, 2021; Xing et al., 2022). While bifurcation analysis of reduced neuronal models has been instrumental in shaping our understanding of the neuronal non-linear dynamics underlying membrane oscillations and cell firing under normal (Jalics et al., 2010; Zhou et al., 2018), pathological (Glass and Mackey, 1979; Kaper et al., 2013), and pharmacological (V-Ghaffari et al., 2017) conditions, extending this analysis to higher-fidelity models with multiple compartments and dendritic channels, which mediate many non-linear firing behaviors in MNs (Hounsgaard and Kiehn, 1993; Li and Bennett, 2007; Carlin et al., 2009) has been limited. As multi-compartment models are more accurate in simulating the firing behaviors of spinal MNs than reduced models (Hendrickson et al., 2011; Elbasiouny, 2014), bifurcation analysis of their behaviors is, therefore, of a great importance. Accordingly, the multi-compartment model developed in the present study in XPPAUT is a novel tool that allows the scientific community to study the effects that somatic and dendritic ion channels have on the bifurcation diagrams of neurons. Applied on a spinal MN, our bifurcation analysis showed a subcritical Hopf node in the cell V-I bifurcation diagram, resulting from somatic SK channels. A subcritical Hopf node is a bifurcation node in which unstable limit cycles surround a stable equilibrium point. Our analysis shows that the subcritical Hopf node is mediated by somatic SK channels due to its negative feedback effects on the membrane potential. Interestingly, while dendritic SK and L-type Ca²⁺channels have similar negative feedback effects between them, their interaction did not generate a subcritical Hopf node in the bifurcation diagram (no subcritical Hopf bifurcation node emerged when somatic SK channels were deactivated but dendritic SK channels were present, data not shown). However, because CaN channels inactivate whereas L-type Ca²⁺channels do not, the negative feedback loop between the somatic SK and CaN channels is dynamic (i.e., changes in magnitude with time), whereas that between the dendritic SK and L-type Ca²⁺channels is static. This difference in negative feedback loop dynamics is probably responsible for the appearance of the subcritical Hopf node with the somatic, but not dendritic, SK channels. Together, these results indicate that not all negative feedback loops have similar effects on cell firing.

4.2. Neuronal input dynamics

Neurons (including MNs and interneurons) receive inputs from presynaptic sources and these synaptic inputs could be constant at times or dynamic at other times. Because bifurcation analysis examines steady state responses of neurons, we tested the model with long current pulses (to evoke steady-state cell firing) to relate the model firing behavior to the bifurcation diagrams (see long pulse simulations in Figures 3D, 4D, 5D, 7B, 8B). To also study dynamic firing of the cell, we tested the MN model with increasing/decreasing triangular current ramps (see current ramp simulations in Figures 3C, 4C, 5C, 7C, D, 8C, D). Therefore, the analysis in the present study covered both constant and dynamic responses of the MN model. While conducted on a MN model, this work could be extended to other types of neurons or interneurons (e.g., central pattern generator interneurons) in the nervous system.

4.3. Bifurcation analysis of diseased neurons

In many neurological conditions, neurons typically experience excitability dysfunction to which somatic and dendritic anatomical changes contribute, such as in ALS (Amendola and Durand, 2008; Elbasiouny et al., 2012; Shoenfeld et al., 2014; Dukkipati et al., 2018), spinal cord injury (Kitzman, 2005), and Alzheimer’s disease (Anderton et al., 1998). For bifurcation analysis to examine neuronal excitability dysfunction in these conditions, multi-compartment models that represent dendritic anatomy and ion channels must be developed in XPPAUT. The multi-compartment model in XPPAUT developed in the present study is a first step in that direction. We demonstrated the utility of that model in investigating the roles somatic and dendritic anatomical changes play in MN excitability dysfunction in ALS. One primary result of our work is that the response of the cell to somatic enlargement is very different from that to dendritic enlargement. While both were expected to reduce the excitability of the cell, somatic enlargement increased the cell recruitment current and reduced the cell firing rate without affecting the operational firing range of the cell (i.e., the range of stable limit cycles along the x-axis in Figure 7A). Dendritic enlargement, on the other hand, had more drastic hypoexcitability effects on the cell, e.g., increased cell recruitment current (Figures 7A, C), and lowered descending firing rate (Figures 7C, D). Dendritic enlargement also reduced the MN’s operational firing rate substantially (Figure 7A), yet in some respects induced cell hyperexcitability by increasing the cell’s ascending firing rate (Figures 7C, D). This increase in firing rate was due to the increase in cell capacitance because of the dendritic enlargement, leading to increased cell time constant; thereby causing Ca PIC full activation before cell recruitment, which inactivated Na channels and reduced the cell firing range. The mixed dendritic enlargement excitability effects is important as it illustrates how one cellular change could induce both hypo- and hyperexcitability effects on the MN, which is highly relevant to the conflicting literature about MN excitability dysfunction in ALS [for review, see Elbasiouny (2022)]. In our ALS simulations, we tried to keep the percent changes we tested within the range of what has been reported experimentally. For dendritic enlargement, we only tested 30%, whereas the mean of reported data is ∼58%, because the cell stopped firing when the dendritic size was increased by > 40%. This prevented us from conducting bifurcation analysis of the effects that dendritic enlargement has on cell firing. Therefore, our reported effects on dendritic enlargement are likely to be underestimated. When both the soma and dendrites are enlarged, the cell’s intrinsic excitability dropped significantly such that the cell was no longer able to fire action potential (data not shown). Given that it is experimentally challenging to study the separate effects of neuronal anatomical changes, bifurcation analysis is an important tool that fulfills this role, and can provide deep insight on disease pathogenesis.

Another primary result of our work is that dendritic overbranching is a hypoexcitability mechanism in itself that offsets the dendritic enlargement hyperexcitability effects resulting from the cell increased capacitance. As dendritic overbranching leads to the flow of less current through more dendritic branches, this mechanism generates less dendritic membrane depolarization; thereby, less dendritic PIC and firing rate. This result could partially explain why SOD MNs, which experience dendritic enlargement and overbranching, had no increased dendritic Ca transients, relative to WT, despite having larger Ca PIC (Quinlan et al., 2011, 2015). Additionally, this result highlights the importance of simulating dendritic branching in reduced models of MNs, as single-branch models will not simulate the dendritic active properties accurately, as shown in this and earlier studies (Hendrickson et al., 2011; Elbasiouny, 2014). An XPPAUT MN model with multiple dendritic branches, as the one developed in the present study, is therefore an important tool for more accurate bifurcation analysis of MN behaviors.

5. Conclusion

In conclusion, the novel multi-compartment MN model developed in XPPAUT in the present study, which is publicly available to the scientific community, is expected to greatly expand our capabilities to study neuronal function under normal and disease conditions.

References

• 1

  AmendolaJ.DurandJ. (2008). Morphological differences between wild-type and transgenic superoxide dismutase 1 lumbar motoneurons in postnatal mice.J. Comp. Neurol.511329–341. 10.1002/cne.21818

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AndertonB. H.CallahanL.ColemanP.DaviesP.FloodD.JichaG. A.et al (1998). Dendritic changes in Alzheimer’s disease and factors that may underlie these changes.Prog. Neurobiol.55595–609. 10.1016/S0301-0082(98)00022-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  CarlinK. P.BuiT. V.DaiY.BrownstoneR. M. (2009). Staircase currents in motoneurons: insight into the spatial arrangement of calcium channels in the dendritic tree.J. Neurosci.295343–5353. 10.1523/JNEUROSCI.5458-08.2009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  CullheimS.FleshmanJ. W.GlennL. L.BurkeR. E. (1987). Three-dimensional architecture of dendritic trees in type-identified alpha-motoneurons.J. Comp. Neurol.25582–96. 10.1002/cne.902550107

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  DhoogeA.GovaertsW.KuznetsovY. A. (2003). MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs.ACM Trans. Math. Softw.29141–164. 10.1145/779359.779362

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  DoedelE. J.ChampneysA. R.DercoleF.FairgrieveT. F.KuznetsovY. A.OldemanB.et al (2007). AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Equations.Montreal: Concordia University.

  • Google Scholar
  • View reference in article

• 7

  DukkipatiS. S.GarrettT. L.ElbasiounyS. M. (2018). The vulnerability of spinal motoneurons and soma size plasticity in a mouse model of amyotrophic lateral sclerosis.J. Physiol.5961723–1745. 10.1113/JP275498

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  ElbasiounyS. M. (2014). Development of modified cable models to simulate accurate neuronal active behaviors.J. Appl. Physiol.1171243–1261. 10.1152/japplphysiol.00496.2014

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  ElbasiounyS. M. (2022). Motoneuron excitability dysfunction in ALS: pseudo-mystery or authentic conundrum?J. Physiol.6004815–4825. 10.1113/JP283630

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  ElbasiounyS. M.QuinlanK. A.EissaT. L.HeckmanC. J. (2012). Electrophysiological Abnormalities in SOD1 Transgenic Models in Amyotrophic Lateral Sclerosis: The Commonalities and Differences.London: IntechOpen.

  • Google Scholar
  • View reference in article

• 11

  ErmentroutB.MahajanA. (2003). Simulating, analyzing, and animating dynamical systems: a guide to XPPAUT for researchers and students.Appl. Mech. Rev.56B53–B53. 10.1115/1.1579454

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  FilipchukA. A.DurandJ. (2012). Postnatal dendritic development in lumbar motoneurons in mutant superoxide dismutase 1 mouse model of amyotrophic lateral sclerosis.Neuroscience209144–154. 10.1016/J.NEUROSCIENCE.2012.01.046

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  FleshmanJ. W.SegevI.BurkeR. B. (1988). Electrotonic architecture of type-identified alpha-motoneurons in the cat spinal cord.J. Neurophysiol.6060–85. 10.1152/jn.1988.60.1.60

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  FoehringR. C.SypertG. W.MunsonJ. B. (1986). Properties of self-reinnervated motor units of medial gastrocnemius of cat. I. Long-term reinnervation.J Neurophysiol.55931–946. 10.1152/JN.1986.55.5.931

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  GlassL.MackeyM. C. (1979). Pathological conditions resulting from instabilities in physiological control systems.Ann. N. Y. Acad. Sci.316214–235. 10.1111/j.1749-6632.1979.tb29471.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  HendricksonE. B.EdgertonJ. R.JaegerD. (2011). The capabilities and limitations of conductance-based compartmental neuron models with reduced branched or unbranched morphologies and active dendrites.J. Comput. Neurosci.30301–321. 10.1007/s10827-010-0258-z

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  HinesM. L.CarnevaleN. T. (1997). The NEURON simulation environment.Neural Comput.91179–1209.

  • Google Scholar
  • View reference in article

• 18

  HochmanS.McCreaD. A. (1994a). Effects of chronic spinalization on ankle extensor motoneurons II. Motoneuron electrical properties.J. Neurophysiol.711468–1479. 10.1152/jn.1994.71.4.1468

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  HochmanS.McCreaD. A. (1994b). Effects of chronic spinalization on ankle extensor motoneurons III. Composite Ia EPSPs in motoneurons separated into motor unit types.J. Neurophysiol.711480–1490. 10.1152/jn.1994.71.4.1480

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  HodgkinA. L.HuxleyA. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve.J. Physiol.117500–544. 10.1113/jphysiol.1952.sp004764

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  HounsgaardJ.KiehnO. (1993). Calcium spikes and calcium plateaux evoked by differential polarization in dendrites of turtle motoneurones in vitro.J. Physiol.468245–259. 10.1113/jphysiol.1993.sp019769

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  HunterJ. D. (2007). Matplotlib: a 2D graphics environment.Comput. Sci. Eng.990–95. 10.1109/MCSE.2007.55

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  JalicsJ.KrupaM.RotsteinH. G. (2010). Mixed-mode oscillations in a three time-scale system of ODEs motivated by a neuronal model.Dyn. Syst.25445–482. 10.1080/14689360903535760

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  KaperT. J.KramerM. A.RotsteinH. G. (2013). Introduction to focus issue: rhythms and dynamic transitions in neurological disease: modeling, computation, and experiment. Chaos An Interdiscip.J. Nonlinear Sci.23:46001. 10.1063/1.4856276

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  KernellD. (1965). High-frequency repetitive firing of cat lumbosacral motoneurones stimulated by long-lasting injected currents.Acta Physiol. Scand.6574–86.

  • Google Scholar
  • View reference in article

• 26

  KitzmanP. (2005). Alteration in axial motoneuronal morphology in the spinal cord injured spastic rat.Exp. Neurol.192100–108. 10.1016/j.expneurol.2004.10.021

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  LeeR. H.HeckmanC. J. (1999). Paradoxical effect of QX-314 on persistent inward currents and bistable behavior in spinal motoneurons in vivo.J. Neurophysiol.822518–2527. 10.1152/jn.1999.82.5.2518

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  LiX.BennettD. J. (2007). Apamin-sensitive calcium-activated potassium currents (SK) are activated by persistent calcium currents in rat motoneurons.J. Neurophysiol.973314–3330. 10.1152/jn.01068.2006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  ManuelM. (2005). How much afterhyperpolarization conductance is recruited by an action potential? A dynamic-clamp study in cat lumbar motoneurons.J. Neurosci.258917–8923. 10.1523/JNEUROSCI.2154-05.2005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  MousaM. H.ElbasiounyS. M. (2020). Dendritic distributions of L-type Ca ²⁺ and SK _L channels in spinal motoneurons: a simulation study.J. Neurophysiol.1241285–1307. 10.1152/jn.00169.2020

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  QuinlanK. A.LamanoJ. B.SamuelsJ.HeckmanC. J. (2015). Comparison of dendritic calcium transients in juvenile wild type and SOD1G93A mouse lumbar motoneurons.Front. Cell. Neurosci.9:139. 10.3389/fncel.2015.00139

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  QuinlanK. A.SchusterJ. E.FuR.SiddiqueT.HeckmanC. J. (2011). Altered postnatal maturation of electrical properties in spinal motoneurons in a mouse model of amyotrophic lateral sclerosis.J. Physiol.5892245–2260. 10.1113/jphysiol.2010.200659

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  ShoenfeldL.WestenbroekR. E.FisherE.QuinlanK. A.TysselingV. M.PowersR. K.et al (2014). Soma size and Cav1.3 channel expression in vulnerable and resistant motoneuron populations of the SOD1G93A mouse model of ALS.Physiol. Rep.21–13. 10.14814/phy2.12113

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  The Pandas Development Team. (2022). Pandas-Dev/Pandas: Pandas. Available online at: https://zenodo.org/record/7344967#.Y9JiuHZBzIU(accessed November 22, 2022). 10.5281/zenodo.3509134

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  Van RossumG.DrakeF. L.Jr. (2010). The Python Language Reference.Amsterdam: Python Software Foundation.

  • Google Scholar
  • View reference in article

• 36

  V-GhaffariB.KouhnavardM.ElbasiounyS. M. (2017). Mixed-mode oscillations in pyramidal neurons under antiepileptic drug conditions.PLoS One12:e0178244. 10.1371/journal.pone.0178244

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  WangT.ChengZ.BuR.MaR. (2019). Stability and Hopf bifurcation analysis of a simplified six-neuron tridiagonal two-layer neural network model with delays.Neurocomputing332203–214. 10.1016/j.neucom.2018.12.005

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  WangT.WangY.ChengZ. (2021). Stability and hopf bifurcation analysis of a general tri-diagonal BAM neural network with delays.Neural Process. Lett.534571–4592. 10.1007/s11063-021-10613-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  WaskomM. (2021). seaborn: statistical data visualization.J. Open Source Softw.6:3021. 10.21105/joss.03021

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  WhiteJ. A.BuddeT.KayA. R. (1995). A bifurcation analysis of neuronal subthreshold oscillations.Biophys. J.691203–1217. 10.1016/S0006-3495(95)79995-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  XingR.XiaoM.ZhangY.QiuJ. (2022). Stability and hopf bifurcation analysis of an (n + m)-neuron double-ring neural network model with multiple time delays.J. Syst. Sci. Complex.35159–178. 10.1007/s11424-021-0108-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  ZengelJ. E.ReidS. A.SypertG. W.MunsonJ. B. (1985). Membrane electrical properties and prediction of motor-unit type of medial gastrocnemius motoneurons in the cat.J. Neurophysiol.531323–1344. 10.1152/jn.1985.53.5.1323

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  ZhouY.VoT.RotsteinH. G.McCarthyM. M.KopellN. (2018). M-Current expands the range of gamma frequency inputs to which a neuronal target entrains.J. Math. Neurosci.81–32. 10.1186/S13408-018-0068-6/FIGURES/16

  • CrossRef
  • Google Scholar
  • View reference in article