A Unified Frequency Domain Model to Study the Effect of Demyelination on Axonal Conduction

Introduction

Axons are an integral part of communication in the nervous system. Conduction and the effects of myelination are an integral part of understanding this communication. Extensive research has been going on in the field of axonal communication to model various electrical responses.^1,2 Most of the current models focus on using the traditional time domain approach. Review of the present day classical cable models indicates that simulation, in terms of analog circuit transient responses, results in attenuation losses between the Node of Ranvier (NR) regions. This loss is often modeled as passive resistive networks. This approach does not provide a method to determine the influences of demyelination based on the simulation results. In order to carry out the high-level system analysis of the neural cortex, we need to have intuitive understanding of these models, and they should be easy to understand. The use of frequency domain modeling approach can help meet these goals. We propose using the frequency domain approach to design a model with a system identification approach on the MATLAB/Simulink platform to analyze the myelinated neuronal axon (with N node). The approach differs from purely analog circuit characterization in a number of ways. These include the evaluation and implementation of transfer function blocks, which can be used to represent a segment of the axon and proper signal conditioning to minimize the attenuation losses along the length of the axon, as would normally be observed in passive lumped circuit models.

The formation of myelin is important for the normal function of the mammalian nervous system. While the effects of demyelination can lead to action potential (AP) attenuation and negative transient response along the length of the axon, current models utilizing the cable model have been unable to reproduce such phenomena. In fact, simulations of the current cable models have shown a failure to eliminate the attenuation losses, and thus, modeling the AP propagation along the length of a myelinated axon model has proven difficult. In this article, we use the frequency domain approach to solve these problems and, thus, provide a better and more realistic biological model. We would be able to better understand the effects of demyelination and in turn the model diseases, such as multiple sclerosis (MS), described in the “Multiple Sclerosis” section, by studying the effects of loss in the number of myelin layers (within a particular intermodal region) on AP magnitude and transient response.

We begin by reviewing the present and popular time-based neuron models in the “Multiple Sclerosis” section and compare with the models of frequency domain. In the “Time Domain Vs Frequency Domain” section, we propose the frequency domain modeling and system identification-based transfer function derivation of the various constituents of neuron. Both open loop and closed loop responses of various kinds of stimuli are explored in this section. Next, in the “Proposed Modeling—Architecture” section, we present the simulation results to verify our modeling approach and justify its validity by comparing the results. The “Simulations and Results Comparison” section describes the results of the simulation, and then we conclude our article in the “Conclusion” section.

Multiple Sclerosis

MS is a disease with the loss of myelin across the nerve fibers. To date, there is no research on trying to model the effects of loss of myelin. Using our model, we can start looking into the effects of conduction velocity on various disease characteristics. Eighty-five percent of people with MS have the relapsing-remitting form of MS (as the name suggests that this is the form of MS that relapses and then manifests itself after some time), and it will be very helpful to study the effects of loss of myelin over time using our model. Lesions in MS can occur in various brain regions. Other potential applications of this model could be to study the effects of demyelination on various brain regions. Using the diffusion tensor images of a person with MS, we can model real data and look into predictive capabilities for disease manifestation in the future. This has a huge potential to look into the steps to help mitigate the long-term effects of MS.

Time Domain vs Frequency Domain

Traditional models have been very useful to predict both the static and dynamic properties of the neuron. Time-based ordinary differential equations remain the heart of neuron modeling. It has been shown that the Hodgkin-Huxley equation can be solved by the numerical solutions to get an AP profile for a single nerve fiber.^1,3 We use the popular system identification technique^4,5 to derive the nominal transfer function for the single chamber of a neuron (Fig. 1). Ideally, the pathway of time domain should give us the identical final outcome as the frequency domain method. In this article, we use frequency-based modeling. Figure 2 demonstrates a more direct comparison of the AP found by the proposed frequency domain modeling and two of the most reliable time domain numerical solutions (ODE15s and Runge-Kutta). It should be noted that the transfer function of the system can only be possible if the system can be modeled as linear or at last as piece-wise linear. We list all the assumptions required for the derivation of the transfer function. We observe that our approximate linear time invariant (LTI) system response is in very much agreement when compared to the traditional time-based models.

OPEN IN VIEWER

Figure 1.

Frequency-domain open loop model for neuron for N Nodes of Ranvier.

OPEN IN VIEWER

Figure 2.

Comparison of the action potential in both conventional time domain vs proposed frequency domain for a single chamber.

In numerical analysis, the Runge-Kutta methods are an important family of implicit and explicit iterative methods, which are used in temporal discretization for the approximation of solutions of ordinary differential equations, while ODE15s method is a numerical MATLAB-based, nonlinear, time-based solution method.

The fundamental advantage of this LTI-based frequency domain analysis is that if we know the transfer function of a block or chamber, we can get the overall system response simply by cascading the individual transfer function blocks. Also, this analysis enables us to have direct control over each localized portion of axonal conduction as opposed to overall time-based models. Some of the key differences between these two approaches are as follows. In a time-based system, the response to myelinated portions is only modeled as a basic static repeater, while in the frequency domain analysis, the most essential phase response is modeled. This is very important for the stability analysis of closed loop neuron system. Using the frequency-based model, we can get a better logical intuition for various parameters. Similarly, the portions of dendrites can be made much more featured when compared to a passive conduction model in the time-based systems. ⁶

Proposed Modeling—Architecture

The basic foundation of this study is to assume the neural conduction model as a LTI system. Because the transfer function is only valid for an LTI system, if the system is nonlinear, then the task becomes to convert the nonlinear system into a linear response first. Suitable assumptions have to be made in order to derive the LTI-based models. Proposed assumptions are as follows:

I = C V m d t + g K n 4 ( V m − V k ) + g N a m 3 h ( V m − V N a ) + g l ( V m − V l )

(1)

To check the validity and universal convergence of these assumptions, we model the transfer function by minimizing Equation 1 (by ignoring sodium and potassium activation and inactivation) as follows:

I = C V m d t + g l ( V m − V l )

(2)

Thus, this system is reduced to a first-order system with finite output response (Fig. 3). By the theory of control analysis, ⁷ this result proves the existence of valid and stable transfer function (following our set of heuristic assumptions).

OPEN IN VIEWER

Figure 3.

First-order response of a reduced form of Hodgkin-Huxley equation (as described in Equation 2). This shows the effect of neglecting the ion activation and inactivation channels. (A) Reduced circuit diagram of the single neuron chamber. (B) First order response with only leakage component into consideration.

Before we move to the actual model development, it should be noted that there are some direct correspondences between the time and frequency domains. A delay in time domain is a phase delay in frequency domain, and the degree of freedoms in the time-dependent variable is directly correlated with the number of poles of the system.

V ( s ) I ( s ) = 5453 s − 3258 s 3 + 4.29 s 2 + 32.63 s + 41.1

(3)

In Equation 3, we see the overall single chamber transfer function of the neuron. The transfer function has been derived by the system identification method in MATLAB. We can see that the transfer function has an right hand plane (RHP) zero and two left hand plane (LHP) poles (dominant). This creates an approximate transfer function relating to the input injection current with the membrane potential with the appropriate phase responses for any localized region of axonal conduction. The overall neuron can be divided into three important regions: soma (the cell body), the periodic unmyelinated portions of axon, and the periodic myelinated portions of the axons. These three portions are illustrated in Figure 4. This figure shows the frequency domain circuit model for the unmyelinated section. The inputs and outputs are buffered from the external (chemical/electrical) environment using buffers, and the central response is governed by the LTI transfer function. The parameters of this transfer function depend on the radius and the length of the periodic myelinated portions. Figure 5 shows the proposed frequency domain model for the myelinated section of the axon. There is a transport delay caused by this section due to the repetition of the AP from its receiving end to its propagating end. The delay of this section depends on the radii of unmyelinated and myelinated portions and the length of the capsule as shown in the figure. The overall transfer function of the entire length will be the result of the cascade of such sections in alternate fashion (cascade of the sections of Figs. 4 and 5).

OPEN IN VIEWER

Figure 4.

Architecture of the implemented model of the Nodes of Ranvier. This shows both side and cross-sectional views illustrating the parameters required for computation.

OPEN IN VIEWER

Figure 5.

Architecture of the implemented model of the myelination sheath between two Nodes of Ranvier. This shows both side and cross-sectional views illustrating the parameters required for computation.

Figure 6 illustrates the open loop characteristics of a single axon/dendrite chamber. As shown in Figure 6B, the transfer function has one zero and three poles. Also from Figure 6A, it can be seen that the phase margin is 33°, which guarantees the stability of the transfer function.

OPEN IN VIEWER

Figure 6.

Simulated frequency response of a single chamber of a neuron variation of action potential velocity vs Node of Ranvier. (A) Magnitude and phase response and (B) transfer function and pole zero plot.

The poles and zeroes shown in this study are dependent parameters (Table 1), which change with respect to time and internal and external chemical concentrations. The actual response to each chamber is described by the overall sum of the cascade of such open loop transfer functions. The Node of Ranvier regions simply duplicate the signals and put the replicated signals on the output. The actual response to each chamber response is described by the overall sum of the cascade of these open loop transfer functions.

OPEN IN VIEWER

Table 1.

Values of parameters.

PARAMETER	Value
g Na	12
g K	14
m, n, h	0.6, 0.02, 0.5
NMax	10
L, Lsh, Ln	150 um, 3 um, 2 um
Rsh	4 um
R n	3 um
r sh	0.01 m-ohm

Simulations and Results Comparison

The simulation for this study has been done using the MAT-LAB and Simulink models following the framework that has been formulated in the “Time Domain Vs Frequency Domain” and “Proposed Modeling—Architecture” sections. Figure 7 shows the simulated variation of AP velocity vs Node of Ranvier (N). This is achieved by the Simulink modeling of a single nerve fiber with a cascade of N nodes (based on the transfer function as derived in the “Time Domain Vs Frequency Domain” section). From Figure 7A, we can see that if we keep the length of the fiber fixed, ⁸ then there is an improvement in the AP conduction velocity. But as N becomes large, there are diminishing returns. Thus, we can always predict an optimum Nodes of Ranvier, N, which will help us to get maximum conduction velocity. From the set of parameter values listed in Table 1, we conclude the optimum N value to be 8.

OPEN IN VIEWER

Figure 7.

Simulated variation of action potential velocity vs Node of Ranvier. (A) The case when the length of axon is fixed. (B) The case when distance between inter nodes is constant.

Similarly, Figure 7B demonstrates the case when the overall nerve length is not fixed, but the internodal distance is constant. ⁷ In this case, the improvement due to the Nodes of Ranvier (myelination) is much more than the case illustrated in Figure 7A. Our proposed modeling approach predicts precise improvements over the conventional methods.^6–10 As our study is focused on AP propagation, we assume that a stimulus of sufficient strength is present to elicit an AP response. Implementation of axonal dimensions allows the realistically relevant propagation velocity, and we furthermore eliminate the attenuation in AP seen in passive circuit models. By treating each internodal region as the locations where demyelination will introduce transient delays, one can see AP propagation delays at each Node of Ranvier once the attenuation effects have been modeled and introduced into this segment as an effort for future work.

Conclusion

The proposed model is intended to minimize the attenuation losses and analytical complexity that comes by using a time-dependent cable model.^5,11 Additionally, it is desired to compare the contributions of myelin toward AP propagation velocity along the length of the axon. All outputs and values to be reported are with respect to the proposed model shown in Figures 4 and 5. As predicted, the current model is capable of eliciting the AP propagation along the length of the axon. Utilization of the transfer function block set, derived from the process described earlier, at the stimulation site NR segment results in the AP propagation along each Node of Ranvier as shown in Figure 7. This study provides an alternate, more intuitive method of analyzing the axonal conduction that will enable the important analysis of myelinated and unmyelinated fibers. ¹²