Electromagnetic characteristics of in vivo nerve fibers at the terahertz-far-infrared band

How terahertz signals perform in the neural system has attracted widespread interest in the life sciences community. Relevant experimental reveals that in animal nerve cells, the myelin sheath of the nerve axon has a higher refractive index than the intracellular and extracellular fluids in the Terahertz-far-infrared (THz-FIR) frequency band. This makes THz-FIR wave transmission possible in nerve fibers. Based on this premise, this article carries out the following work from the theoretical level to investigate the electromagnetic (EM) characteristics of in vivo nerve fibers at the THz-FIR band. First, the EM transmission model of the nerve fibers is established and studied theoretically. The dispersion curves of THz-FIR wave modals transmission in nerve fibers are calculated, which predict that nerve fibers can act as dielectric waveguides for transmitting THz-FIR waves and the THz-FIR waves can transmit at speeds up to 10⁸ m/s. Second, a mode matching algorithm is proposed, which is named RNMMA, to calculate the transmission characteristics of THz-FIR waves at the nodes of Ranvier. The scattering matrix obtained from the proposed algorithm is in good agreement with the results from EM simulation software, which reveals how THz-FIR signals are transmitted forward through the nodes of Ranvier with low loss.

Introduction

As a complex physical system, the human brain has often been a popular topic in life science research, especially in the study of neural information generation and transmission in neurobiology (Adolphs, 2015; Liu, 2018). Electrochemical theory played a leading role in the early studies, although it could not perfectly explain many underlying problems, such as the generation and origin of consciousness (Tononi, 2004; Tegmark, 2015). Meanwhile, action potential signaling mainly relies on the transmembrane transport of ions to cause the potential difference between the excitatory site and the non-excited site, thereby forming a local current. However, traditional electrochemical theories cannot explain the transmission mechanism of high-frequency signals in living organisms, such as how THz and FIR signals that are generated by various redox reactions in biological cells are transmitted. In this context, research of the physical mechanisms that coexist with existing electrochemical theories may reveal new prospects for explaining high-level signal transduction mechanisms, including whether there are more efficient physical ways to transmit and store information in the nervous system. Neuro-electromagnetics provides a new way to study the mechanism of the high-level communication and rapid response of the nervous system. However, it also has many problems, such as how to generate and transmit signals efficiently.

Recently, studies on THz-FIR spectroscopy have confirmed that many oxidative metabolic processes in biological cells produce a large number of biological photons, covering the spectrum range of THz-FIR. Many of the biological macromolecules that play an important role in life activities, such as proteins and phospholipids, also have vibration frequencies in the THz-FIR range (Wetzel and LeVine, 1999; Laman et al., 2008; Miller and Dumas, 2010). The ion channels on the surface of the nerve membrane are not only related to action potential activity but also show that the process of ion transport across the membrane generates THz-FIR radiation (Li et al., 2021; Liu et al., 2021; Guo et al., 2022; Hu et al., 2022; Wang et al., 2022). These are several potential sources for the THz-FIR signals transmitting on the nerve fibers. The THz-FIR signals can also affect the activity of the proteins in the cells (Chao et al., 2021), DNA molecular unhelix (Wu et al., 2020), and the selective transmission of ion channels (Li et al., 2021). As information carrier candidates, THz biophotons may become high-speed signal carriers that connect intracellular organelles and cells independently of electrical and chemical transmissions (Xiang et al., 2020).

Some nerve fibers in the vertebrate nervous system are a kind of EM waveguide structure (i.e., wrapped by high resistance and the low capacitance of nerve myelin sheath). This structure can accelerate the conduction of action potential. Meanwhile, in the FIR range, it supports biological photon signal transmission between cells (Kumar et al., 2016; Zangari et al., 2018; Liu et al., 2019). This natural waveguide structure enables signal transmission at a speed of up to 10⁸ m/s, which is far higher than the conduction speed of the action potential. The nerve myelin sheath with a dense membrane structure can isolate the polar water environment on both sides of the membrane and greatly reduce the attenuation of the THz-FIR wave due to polarization loss in the transmission process.

Our current work is focused on how THz-FIR waves behave in nerve fibers, and how the THz-FIR waves hop and exchange energy with periodic nodes of Ranvier structures. In view of this hypothesis, we build a transmission model of nerve fibers and nodes of Ranvier (as shown in Figure 1). The myelinated nerve axon can be regarded as a fiber transmission line that is loaded with layered annular media (i.e., nerve cell fluid, myelin sheath, and extracellular fluid). In contrast to optical fiber, the energies of THz-FIR waves concentrate in the myelin sheath area rather than the fiber core area because of its high dielectric constant; while the physical environment of the nerve fibers are polarity solutions, and there is no obvious boundary. Meanwhile, the node of Ranvier is a periodic structure that is formed by the discontinuity of the myelin sheath, which will inevitably produce energy changes. The main work and innovations of this article are as follows: 1) the field equations and dispersion curves of THz-FIR wave transmission in both myelinated and unmyelinated bare nerve fibers are obtained; 2) the transmission efficiencies of THz-FIR electromagnetic waves in different modes are quantitatively analyzed; and 3) a new mode matching algorithm is proposed that can quickly calculate the scattering parameters and energy exchange of THz-FIR wave transmission at the nodes of Ranvier.

FIGURE 1

FIGURE 1. Schematic diagram of the nerve fiber structure. (A) Cross section of the nerve fiber and (B) 3D perspective graphics of the nerve fiber.

Methods

THz-FIR transmission modes in myelin-sheathed nerve fibers

As shown in Figure 1, the nerve fiber is structurally composed of three regions, and the dielectric constants are represented by ${\epsilon }_{\mathit{1}}$, ${\epsilon }_2$, and ${\epsilon }_3$, respectively. Compared with the intracellular fluid and extracellular fluid of nerve cells (i.e., water as the solvent and have strong absorption characteristics in the THz-FIR band), the myelin sheath is a tight membrane structure with a high dielectric constant and low loss in the THz-FIR band (Liu et al., 2019). The relationship of the dielectric constant of the three-layer media can be expressed as ${\epsilon }_3<{\epsilon }_{\mathit{1}}<{\epsilon }_2$ (Friede and Bischhausen, 1982; Kumar et al., 2016; Liu et al., 2019).

From the completeness of the EM field, it can be known that the THz-FIR transmission modes on the nerve fiber are composed of the guided modes of the finite mode and the integral of the radiation modes over their continuum spectra. Hence, we decompose the fields on the nerve fiber into two parts,

$\mathit{E}\left(r,\phi ,z,t\right)={\displaystyle \sum }_{n=\mathit{1}}^N{C}_n{\mathit{E}}_n\left(r,\phi \right)e^{j\left({\omega t-\beta }_{z,n}z\right)}+\int C_{\xi }{\mathit{E}}_{\xi }\left(r,\phi \right)e^{j\left({\omega t-\beta }_{z,\xi }z\right)}d\xi$(1–1)

$\mathit{H}\left(r,\phi ,z,t\right)={\displaystyle \sum }_{n=\mathit{1}}^N{D}_n{\mathit{H}}_n\left(r,\phi \right)e^{j\left({\omega t-\beta }_{z,n}z\right)}+\int D_{\xi }{\mathit{H}}_{\xi }\left(r,\phi \right)e^{j\left({\omega t-\beta }_{z,\xi }z\right)}d\xi$(1–2)

where the first terms with the subscript n represent the guided modes at a certain frequency and the second terms with continuous subscript $\xi$ denote the integral over the radiation modes. Here, the integration variable $\xi$ represents the radial wavenumbers of the radiation modes. N is the total number of the guided modes on the nerve fiber at a certain frequency. $C_n$ is the normalized amplitude, $\omega$ is the angular frequency, and ${\beta }_{z,n}$ and ${\beta }_{z,\xi }$ are the propagation constants of the nth or $\xi$th mode along the nerve fiber and have the same value in the three regions. ${\mathit{E}}_n\left(r,\phi \right)$, ${\mathit{H}}_n\left(r,\phi \right)$ describe the electric and magnetic fields of the nth mode in the cross-sectional area of the nerve fiber, which can be expressed as the superpositions of transverse and longitudinal components,

${\mathit{E}}_n\left(r,\phi \right)={\mathit{E}}_{t,n}\left(r,\phi \right)+{\mathit{E}}_{z,n}\left(r,\phi \right),{\mathit{H}}_n\left(r,\phi \right)={\mathit{H}}_{t,n}\left(r,\phi \right)+{\mathit{H}}_{z,n}\left(r,\phi \right)(2)$

The longitudinal components of the electric field ${\mathit{E}}_{z,n}\left(r,\phi \right)$ and magnetic field ${\mathit{H}}_{z,n}\left(r,\phi \right)$ satisfy the Helmholtz equation. Thus, we have

$\begin{array}{c}{\nabla }^{\mathit{2}}{\mathit{E}}_z+{\omega }^{\mathit{2}}\mu \epsilon {\mathit{E}}_{z,n}-j\omega \mu \sigma {\mathit{E}}_{z,n}=0,\\ {\nabla }^{\mathit{2}}{\mathit{H}}_z+{\omega }^{\mathit{2}}\mu \epsilon {\mathit{H}}_{z,n}-j\omega \mu \sigma {\mathit{H}}_{z,n}=0\end{array}(3)$

which can be further expressed as,

$\begin{array}{c}{\mathit{E}}_{\text{iz},n}=C_{\text{in}}{\psi }_{\text{in}}\left({\beta }_{\text{in}}r\right)\left(\begin{array}{c}\sin n\phi \\ \cos n\phi \end{array}\right),i=\mathrm{1,2,3},\\ {\mathit{H}}_{\text{iz},n}=D_{\text{in}}{\psi }_{\text{in}}\left({\beta }_{\text{in}}r\right)\left(\begin{array}{c}\sin n\phi \\ \cos n\phi \end{array}\right),i=\mathrm{1,2,3}\end{array}(4)$

where ${\psi }_{\mathrm{i}\mathrm{n}}$ represents the Bessel function of order n and has different combinations and representations in different regions. The subscript i represents ith region in Figure 1A. ${\beta }_{\mathrm{i}\mathrm{n}}$ represents the transverse propagation constant in the ith region.

Here, we make the approximation ${\epsilon }_{\mathit{1}}\approx {\epsilon }_3$, due to the influence of the high dielectric constant myelin sheaths. However, for the bare nerve axon at the node of Ranvier, the approximation is no longer used because the bare axon is similar to a weakly conducting fiber (which we will discuss in detail later on). For the guided modes in the myelinated axons, the propagation constants ${\beta }_z$ satisfy the relations of $k_0\sqrt{e_{\mathit{1}}}<{\beta }_z<k_0\sqrt{e_{\mathit{2}}}$. The transverse propagation constants of the three regions can be expressed as ${\beta }_{\mathit{1}n}^{\mathit{2}}={\beta }_{z,n}^{\mathit{2}}-k_0^{\mathit{2}}{\epsilon }_{\mathit{1}},{\beta }_{\mathit{2}n}^{\mathit{2}}=k_0^{\mathit{2}}{\epsilon }_{\mathit{2}}-{\beta }_{z,n}^{\mathit{2}},$ and ${\beta }_{\mathit{3}n}^{\mathit{2}}={\beta }_{z,n}^{\mathit{2}}-k_0^{\mathit{2}}{\epsilon }_{\mathit{3}}$, where $k_0=\omega /c_0$. Accordingly, the Bessel function distribution of the corresponding region satisfies ${\psi }_{\mathit{1}n}=I_n\left({\beta }_{\mathit{1}n}r\right)$, ${\psi }_{\mathit{2}n}=J_n\left({\beta }_{\mathit{2}n}r\right)+Y_n\left({\beta }_{\mathit{2}n}r\right)$, and ${\psi }_{\mathit{3}n}=K_n\left({\beta }_{\mathit{1}n}r\right)$, respectively, which means that the field decays exponentially at both sides of the myelin sheath region. The transverse field components are expressed in terms of the longitudinal field components,

${\mathit{E}}_t=\frac{j}{{\omega }^{\mathit{2}}{\mu }_0{\epsilon }_0{\epsilon }_i{-\beta }_z^{\mathit{2}}}\left({-\beta }_z{\nabla }_t{\mathit{H}}_z+{\omega \mu }_0{\mathit{e}}_z\times {\nabla }_t{\mathit{H}}_z\right),{\mathit{H}}_t=\frac{j}{{\omega }^{\mathit{2}}{\mu }_0{\epsilon }_0{\epsilon }_i{-\beta }_z^{\mathit{2}}}\left({-\beta }_z{\nabla }_t{\mathit{E}}_z+{\omega \epsilon }_0{\epsilon }_i{\mathit{e}}_z\times {\nabla }_t{\mathit{E}}_z\right)(5)$

Furthermore, we can obtain the guided mode expressions of 12 field components in the three regions, as follows:

$\left[\begin{array}{l}{\mathit{E}}_{1\mathit{r}}\\ {\mathit{E}}_{1\phi }\\ {\mathit{H}}_{1\mathit{r}}\\ {\mathit{H}}_{1\phi }\end{array}\right]=j\left[\begin{array}{cc}\begin{array}{l}\frac{{\beta }_{z,n}}{{\beta }_{1n}}{I}_n^{\prime }\left({\beta }_{1n}r\right)\sin n\phi \\ \frac{{\beta }_{z,n}n}{{\beta }_{1n}^{2}r}{I}_n\left({\beta }_{1n}r\right)\cos n\phi \\ -\frac{{\omega \epsilon }_0{\epsilon }_{1}n}{{\beta }_{1n}^{2}r}{I}_n\left({\beta }_{1n}r\right)\cos n\phi \\ \frac{{\omega \epsilon }_0{\epsilon }_1}{{\beta }_{1n}}{I}_n^{\prime }\left({\beta }_{1}r\right)\sin n\phi \end{array}& \begin{array}{l}-\frac{\omega \mu n}{{\beta }_1^{2}r}{I}_n\left({\beta }_{1n}r\right)\sin n\phi \\ -\frac{\omega \mu }{{\beta }_{1n}}{I}_n^{\prime }\left({\beta }_{1n}r\right)\cos n\phi \\ \frac{{\beta }_{z,n}}{{\beta }_{1n}}{I}_n^{\prime }\left({\beta }_{1n}r\right)\cos n\phi \\ -\frac{{\beta }_{z}n}{{\beta }_{1n}^{2}r}{I}_n\left({\beta }_{1n}r\right)\sin n\phi \end{array}\end{array}\right]\left[\begin{array}{l}{C}_1\\ D_1\end{array}\right](\mathrm{6a})$

$\left[\begin{array}{l}{\mathit{E}}_{2r}\\ {\mathit{E}}_{2\phi }\\ {\mathit{H}}_{2r}\\ {\mathit{H}}_{2\phi }\end{array}\right]=j\left[\begin{array}{cccc}-\frac{{\beta }_{z,n}}{{\beta }_{2n}}{J}_n^{\prime }\left({\beta }_{2n}r\right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{n}\phi & -\frac{{\beta }_{z,n}}{{\beta }_{2n}}{Y}_n^{\prime }\left({\beta }_{2n}r\right)\sin n\phi & \frac{\omega \mu n}{{\beta }_{2n}^{2}r}{J}_n\left({\beta }_{2n}r\right)\sin n\phi & \frac{\omega \mu n}{{\beta }_{2n}^{2}r}{Y}_n\left({\beta }_{2n}r\right)\sin n\phi \\ -\frac{{\beta }_{z,n}n}{{\beta }_{2n}^{2}r}{J}_n\left({\beta }_{2n}r\right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{n}\phi & -\frac{{\beta }_{z,n}n}{{\beta }_{2n}^{2}r}{Y}_n\left({\beta }_{2n}r\right)\cos n\phi & \frac{\omega \mu }{{\beta }_{2n}}{J}_n^{\prime }\left({\beta }_{2n}r\right)\cos n\phi & \frac{\omega \mu }{{\beta }_{2n}}{Y}_n^{\prime }\left({\beta }_{2n}r\right)\cos n\phi \\ \frac{{\omega \epsilon }_0{\epsilon }_{2}n}{{\beta }_{2n}^{2}r}{J}_n\left({\beta }_{2n}r\right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{n}\phi & \frac{{\omega \epsilon }_0{\epsilon }_{2}n}{{\beta }_{2n}^{2}r}{Y}_n\left({\beta }_{2n}r\right)\cos n\phi & -\frac{{\beta }_{z,n}}{{\beta }_{2n}}{J}_n^{\prime }\left({\beta }_{2n}r\right)\cos n\phi & -\frac{{\beta }_{z,n}}{{\beta }_{2n}}{Y}_n^{\prime }\left({\beta }_{2n}r\right)\cos n\phi \\ -\frac{{\omega \epsilon }_0{\epsilon }_2}{{\beta }_{2n}}{J}_n^{\prime }\left({\beta }_{2n}r\right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{n}\phi & -\frac{{\omega \epsilon }_0{\epsilon }_2}{{\beta }_{2n}}{Y}_n^{\prime }\left({\beta }_{2n}r\right)\sin n\phi & \frac{{\beta }_{z,n}n}{{\beta }_{2n}^{2}r}{J}_n\left({\beta }_{2n}r\right)\sin n\phi & \frac{{\beta }_{z,n}n}{{\beta }_{2n}^{2}r}{Y}_n\left({\beta }_{2n}r\right)\sin n\phi \end{array}\right]\left[\begin{array}{l}{C}_{21}\\ C_{22}\\ D_{21}\\ D_{22}\end{array}\right](\mathrm{6b})$

$\left[\begin{array}{l}{\mathit{E}}_{3r}\\ {\mathit{E}}_{3\phi }\\ {\mathit{H}}_{3r}\\ {\mathit{H}}_{3\phi }\end{array}\right]=j\left[\begin{array}{cc}\begin{array}{l}\frac{{\beta }_{z,n}}{{\beta }_{3n}}{K}_n^{\prime }\left({\beta }_{3n}r\right)\sin n\phi \\ \frac{{\beta }_{z,n}n}{{\beta }_{3n}^{2}r}{K}_n\left({\beta }_{3n}r\right)\cos n\phi \\ -\frac{{\omega \epsilon }_0{\epsilon }_{3}n}{{\beta }_{3n}^{2}r}{K}_n\left({\beta }_{3n}r\right)\cos n\phi \\ \frac{{\omega \epsilon }_0{\epsilon }_3}{{\beta }_{3n}}{K}_n^{\prime }\left({\beta }_{1}r\right)\sin n\phi \end{array}& \begin{array}{l}-\frac{\omega \mu n}{{\beta }_{3n}^{2}r}{K}_n\left({\beta }_{3n}r\right)\sin n\phi \\ -\frac{\omega \mu }{{\beta }_{3n}}{K}_n^{\prime }\left({\beta }_{3n}r\right)\cos n\phi \\ \frac{{\beta }_{z,n}}{{\beta }_{3n}}{K}_n^{\prime }\left({\beta }_{3n}r\right)\cos n\phi \\ -\frac{{\beta }_{z,n}n}{{\beta }_{3n}^{2}r}{K}_n\left({\beta }_{3n}r\right)\sin n\phi \end{array}\end{array}\right]\left[\begin{array}{l}{C}_3\\ D_3\end{array}\right](\mathrm{6c})$

The dispersion curves of THz-FIR wave transmitting in nerve fibers can be obtained by using the tangential continuous boundary conditions for the electric and magnetic fields on $r=R_{\mathit{1}}$ and $r=R_{\mathit{2}}$, and the dispersion equation of nth HEM mixing mode is as follows:

$\left|\begin{array}{cccccccc}{I}_n\left({\beta }_{\mathit{1}n}{R}_{\mathit{1}}\right)& 0& {-J}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& {-Y}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& 0& 0& 0& 0\\ 0& I_n\left({\beta }_{\mathit{1}n}{R}_{\mathit{1}}\right)& 0& 0& {-J}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& {-Y}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& 0& 0\\ \frac{{\beta }_{z}n}{{\beta }_{\mathit{1}n}^{\mathit{2}}{R}_{\mathit{1}}}{I}_n\left({\beta }_{\mathit{1}n}{R}_{\mathit{1}}\right)& -\frac{\omega \mu }{{\beta }_{\mathit{1}n}}{I}_n^{\prime }\left({\beta }_{\mathit{1}n}{R}_{\mathit{1}}\right)& \frac{{\beta }_{z}n}{{\beta }_{\mathit{2}n}^{\mathit{2}}{R}_{\mathit{1}}}{J}_n\left({\beta }_{\mathit{1}n}{R}_{\mathit{1}}\right)& \frac{{\beta }_{z}n}{{\beta }_{\mathit{2}n}^{\mathit{2}}{R}_{\mathit{1}}}{Y}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& -\frac{\omega \mu }{{\beta }_{\mathit{2}n}}{J}_n^{\prime }\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& -\frac{\omega \mu }{{\beta }_{\mathit{2}n}}{Y}_n^{\prime }\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& 0& 0\\ \frac{{\omega \epsilon }_0{\epsilon }_{\mathit{1}}}{{\beta }_{\mathit{1}n}}{I}_n^{\prime }\left({\beta }_{\mathit{1}}{R}_{\mathit{1}}\right)& -\frac{{\beta }_{z}n}{{\beta }_{\mathit{1}n}^{\mathit{2}}{R}_{\mathit{1}}}{I}_n\left({\beta }_{\mathit{1}n}{R}_{\mathit{1}}\right)& \frac{{\omega \epsilon }_0{\epsilon }_{\mathit{2}}}{{\beta }_{\mathit{2}n}}{J}_n^{\prime }\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& \frac{{\omega \epsilon }_0{\epsilon }_{\mathit{2}}}{{\beta }_{\mathit{2}n}}{Y}_n^{\prime }\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& -\frac{{\beta }_{z}n}{{\beta }_{\mathit{2}n}^{\mathit{2}}{R}_{\mathit{1}}}{J}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& -\frac{{\beta }_{z}n}{{\beta }_{\mathit{2}n}^{\mathit{2}}{R}_{\mathit{1}}}{Y}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{1}}\right)& 0& 0\\ 0& 0& J_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& Y_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& 0& 0& {-K}_n\left({\beta }_{\mathit{3}n}{R}_{\mathit{2}}\right)& 0\\ 0& 0& 0& 0& J_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& Y_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& 0& {-K}_n\left({\beta }_{\mathit{3}n}{R}_{\mathit{2}}\right)\\ 0& 0& -\frac{{\beta }_{z}n}{{\beta }_{\mathit{2}n}^{\mathit{2}}{R}_{\mathit{2}}}{J}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& -\frac{{\beta }_{z}n}{{\beta }_{\mathit{2}n}^{\mathit{2}}{R}_{\mathit{2}}}{Y}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& \frac{\omega \mu }{{\beta }_{\mathit{2}n}}{J}_n^{\prime }\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& \frac{\omega \mu }{{\beta }_{\mathit{2}n}}{Y}_n^{\prime }\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& -\frac{{\beta }_{z}n}{{\beta }_{\mathit{3}n}^{\mathit{2}}{R}_{\mathit{2}}}{K}_n\left({\beta }_{\mathit{3}n}{R}_{\mathit{2}}\right)& \frac{\omega \mu }{{\beta }_{\mathit{3}n}}{K}_n^{\prime }\left({\beta }_{\mathit{3}n}{R}_{\mathit{2}}\right)\\ 0& 0& -\frac{{\omega \epsilon }_0{\epsilon }_{\mathit{2}}}{{\beta }_{\mathit{2}n}}{J}_n^{\prime }\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& -\frac{{\omega \epsilon }_0{\epsilon }_{\mathit{2}}}{{\beta }_{\mathit{2}n}}{Y}_n^{\prime }\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& \frac{{\beta }_{z}n}{{\beta }_{\mathit{2}n}^{\mathit{2}}{R}_{\mathit{2}}}{J}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& \frac{{\beta }_{z}n}{{\beta }_{\mathit{2}n}^{\mathit{2}}{R}_{\mathit{2}}}{Y}_n\left({\beta }_{\mathit{2}n}{R}_{\mathit{2}}\right)& -\frac{{\omega \epsilon }_0{\epsilon }_{\mathit{3}}}{{\beta }_{\mathit{3}n}}{K}_n^{\prime }\left({\beta }_{\mathit{1}}{R}_{\mathit{2}}\right)& \frac{{\beta }_{z}n}{{\beta }_{\mathit{3}n}^{\mathit{2}}{R}_{\mathit{2}}}{K}_n\left({\beta }_{\mathit{3}n}{R}_{\mathit{2}}\right)\end{array}\right|=0(7)$

It is worth noting that when n≠0, the aforementioned equation has non-zero solutions, indicating that the THz-FIR waves on the nerve fibers are transmitting in mixed modes. When n=0, the THz-FIR waves will no longer transmit in the mixed modes but are decomposed into independent TE and TM modes, and the field components satisfy ${\mathit{E}}_z={\mathit{E}}_{\phi }={\mathit{H}}_r=0$ or ${\mathit{H}}_z={\mathit{H}}_{\phi }={\mathit{E}}_r=0$. Therefore, in the case of n = 0, the left-hand determinant can be rewritten as a 4 × 4 order determinant.

For the radiation modes, it is used to describe the radiation fields outside and inside of the myelin sheath, which are necessary supplements to the guided modes to provide a completely orthogonal set of modes. The fields of the radiation modes will not decay exponentially for increasing values of r outside of the nerve fiber (i.e., extend to infinity). Because there is no need to limit the function describing the radiation modes to those decaying exponentially for large r, a combination of Bessel and Neumann functions is used to describe the radiation modes. It is worth noting that the radiation modes also satisfy the boundary conditions and the fields remain finite on-axis at r = 0. The propagation constants of radiation modes satisfy ${\beta }_z<k_0\sqrt{e_{\mathit{1}}}$, and the transverse propagation constants can be expressed as ${\beta }_{\mathit{1}n}^{\mathit{2}}=k_0^{\mathit{2}}{\epsilon }_{\mathit{1}}-{\beta }_{z,n}^{\mathit{2}},{\beta }_{\mathit{2}n}^{\mathit{2}}=k_0^{\mathit{2}}{\epsilon }_{\mathit{2}}-{\beta }_{z,n}^{\mathit{2}}$, and ${\beta }_{\mathit{3}n}^{\mathit{2}}=k_0^{\mathit{2}}{\epsilon }_{\mathit{3}}-{\beta }_{z,n}^{\mathit{2}}$ in three regions. Accordingly, the Bessel function ${\psi }_{\mathrm{i}\mathrm{n}}$ in Eq. 4 can be written as ${\psi }_{\mathit{1}n}=J_n\left({\beta }_{\mathit{1}n}r\right)$, ${\psi }_{\mathit{2}n}=J_n\left({\beta }_{\mathit{2}n}r\right)+Y_n\left({\beta }_{\mathit{2}n}r\right)$, and ${\psi }_{\mathit{3}n}=J_n\left({\beta }_{\mathit{3}n}r\right)+Y_n\left({\beta }_{\mathit{3}n}r\right)$. Unfortunately, in the general field solution of the radiation mode, there are 10 unknowns and only eight boundary conditions. Therefore, the propagation constants of the radiation modes can theoretically take continuous values (Dietrich and Marcuse, 1970; Marcuse, 1973).

Here, we ignore the evanescent radiation modes of the nerve fiber, which are the modes with imaginary values of the propagation constants ${\beta }_z$ and do not propagate in the z-direction. However, omitting evanescent radiation modes keeps the set of orthogonal modes from being complete.

Nodes of Ranvier mode matching algorithm

In the vertebrate nervous system, nerve fibers are not completely enclosed by myelin sheaths, which periodically disintegrate to form structures called the node of Ranvier. The membranes of the axon that are exposed at the narrow node of Ranvier contain large numbers of ion channels, which enable the action potential to jump over the nodes. A schematic diagram of the structure is shown in Figure 2. The calculations of the electric and magnetic fields at the node of Ranvier are very complicated because of the discontinuity of the longitudinal section of the THz-FIR wave in the transmission direction. To analyze the THz-FIR wave transmission and field mode transformation at the nodes of Ranvier, we propose the nodes of Ranvier pattern matching algorithm (RNMMA), which theoretically consists of three steps: 1) the RNMMA begins by expanding the tangential components of the electric and magnetic fields at the node in terms of the normal modes in each region; 2) the scattering matrices on the boundaries of III and IV are obtained by using mode orthogonality and power normalization; and 3) the total scattering matrix of one node of Ranvier can be obtained from the matrix cascade, and then the distributions and mode characteristics of the electric and magnetic fields at the nodes of Ranvier can be solved. To simplify the analysis, we made the boundaries of III and IV perpendicular to the z-direction, while the boundaries are not strictly perpendicular physiologically.

FIGURE 2

FIGURE 2. Schematic diagram of one node of Ranvier structure, and electromagnetic wave transmissions and reflections.

As shown in Figure 2, the structure of the node of Ranvier (④ and ⑤ regions) can be regarded as an unmyelinated neural axon, whose structure is similar to an optical fiber. The permittivity of the extracellular fluid in ⑤ region is the same as that in ③ region, and the permittivity in ④ region is the same as that in ① region. The field components can also be described by Eqs 1–5. The Bessel function in ④ and ⑤ regions can be expressed as ${\psi }_{4n}=J_n\left({\beta }_{4n}r\right)$ and ${\psi }_{5n}=K_n\left({\beta }_{5n}r\right)$ for guided modes, which can ensure that the field propagates in region ④ and decays exponentially in the r-direction in region ⑤. The propagation constants satisfy $k_0\sqrt{{\epsilon }_5}<{\beta }_{z,n}<k_0\sqrt{{\epsilon }_{\mathit{3}}}$, and the transverse propagation constants can be expressed ${\beta }_{4n}^{\mathit{2}}=k_0^{\mathit{2}}{\epsilon }_4-{\beta }_{z,n}^{\mathit{2}}$ and ${\beta }_{5n}^{\mathit{2}}={\beta }_{z,n}^{\mathit{2}}-k_0^{\mathit{2}}{\epsilon }_5$. When ${\beta }_{z,n}<k_0\sqrt{{\epsilon }_5}$, the transmission modes of the THz-FIR wave generated by the node of Ranvier are radiation modes. The transverse propagation constants can be expressed as ${\beta }_n^{\mathit{2}}=k_0^{\mathit{2}}{\epsilon }_4-{\beta }_{z,n}^{\mathit{2}}$ and ${\beta }_{5n}^{\mathit{2}}=k_0^{\mathit{2}}{\epsilon }_5-{\beta }_{z,n}^{\mathit{2}}$. The field components in the r-direction can be described by a combination of Bessel and Neumann functions, ${\psi }_{4n}=J_n\left({\beta }_{4n}r\right)$ and ${\psi }_{5n}=J_n\left({\beta }_{5n}r\right)+Y_n\left({\beta }_{5n}r\right)$.

The tangential components of electric and magnetic fields at the boundary can be expressed in terms of each normal mode:

${\mathit{E}}_t={\displaystyle \sum }_{N,p}{C}_N^p{\mathit{E}}_{t,N}^p\left(r,\phi \right)e^{{-\mathit{jp}\beta }_{\text{zN}}z}(8)$

${\mathit{H}}_t={\displaystyle \sum }_{N,p}{C}_N^p{\mathit{H}}_{t,N}^p\left(r,\phi \right)e^{{-\mathit{jp}\beta }_{\text{zN}}z}(9)$

where p (p = +, −) represents the propagation direction of the wave, N indicates the mode number, $A_N^p$ is the amplitude coefficient of each mode, and ${\mathit{E}}_{t,N}^p\left(r,\phi \right)$ and ${\mathit{H}}_{t,N}^p\left(r,\phi \right)$ are the transverse electric field and magnetic field of each mode in different regions. The pattern functions of different regions all satisfy the orthogonality theorem of patterns,

$\underset{S}{\oint }{\mathit{e}}_z\cdot \left({\mathit{E}}_{t,m}^q\times {\mathit{H}}_{t,n}^{p^{*}}+{\mathit{E}}_{t,n}^{p^{*}}\times {\mathit{H}}_{t,m}^q\right)\text{ds}={4\delta }_{\mathit{mn}}{S}_0\left(p,q\right)P_n(10)$

where ${\delta }_{\mathrm{m}\mathrm{n}}$ is the Kronecker delta. When $\mathrm{m}=\mathrm{n}$, ${\delta }_{\mathrm{m}\mathrm{n}}=\mathit{1}$, otherwise ${\delta }_{\mathrm{m}\mathrm{n}}=0$, which means the different modes are orthogonal to each other. Meanwhile, p and q represent the propagation directions of the modes. When p and q are both positive, $S=\mathit{1}$; when p and q are both negative, $S=-\mathit{1}$; otherwise, $S=0$. This means the same modes of forwarding and backward transmissions are also orthogonal to each other. The scattering matrix (SM) is the premise of analyzing the THz-FIR wave propagation mode at the node of Ranvier, the core of which is establishing the coupling matrix (CM). N and M modes are, respectively, considered in the axon and node of Ranvier regions on both sides of the III boundary. By applying the continuity conditions (${\mathit{E}}_{t,\mathrm{a}\mathrm{x}\mathrm{o}\mathrm{n}}={\mathit{E}}_{t,\mathrm{N}\mathrm{R}}$ and ${\mathit{H}}_{t,\mathrm{a}\mathrm{x}\mathrm{o}\mathrm{n}}={\mathit{H}}_{t,\mathrm{N}\mathrm{R}}$) at the junction (z = 0), we obtain the equation for the electric field as

$\begin{array}{cc}0\le r<R_{\mathit{1}}& {\displaystyle \sum }_{n=\mathit{1}}^N\left(C_{\mathit{1},n}^{+}\pm C_{\mathit{1},n}^{-}\right){\psi }_{\mathit{1}t,n}\left(r,\phi \right)\\ R_{\mathit{1}}\le r<R_{\mathit{2}}& {\displaystyle \sum }_{n=\mathit{1}}^N\left(C_{\mathit{2},n}^{+}\pm C_{\mathit{2},n}^{-}\right){\psi }_{\mathit{2}t,n}\left(r,\phi \right)\\ R_{\mathit{2}}\le r& {\displaystyle \sum }_{n=\mathit{1}}^N\left(C_{\mathit{3},n}^{+}\pm C_{\mathit{3},n}^{-}\right){\psi }_{\mathit{3}t,n}\left(r,\phi \right)\end{array}\}=\{\begin{array}{cc}{\displaystyle \sum }_{m=\mathit{1}}^M\left(D_{4,m}^{+}\pm D_{4,m}^{+}\right){\psi }_{t,m}\left(r,\phi \right)& 0\le r<R_{\mathit{1}}\\ {\displaystyle \sum }_{m=\mathit{1}}^M\left(D_{5,m}^{+}\pm D_{5,m}^{+}\right){\psi }_{5t,m}\left(r,\phi \right)& R_{\mathit{1}}\le r\end{array}(11)$

where ${\psi }_{\mathrm{i}\mathrm{t},n}\left(r,\phi \right)$ represents tangential electric field (E) and tangential magnetic field (H). By making use of the property of mode orthogonality, we can derive a set of an equation involving the unknown coefficients. To this end, the aforementioned equation is multiplied ${\psi }_{\mathrm{i}\mathrm{t},n}\left(r,\phi \right)$ in the different regions and integrated with the whole interface. We can get

${\displaystyle \sum }_{i=\mathrm{1,2,3}}{\displaystyle \sum }_{n=\mathit{1}}^N{F}_{i,n}\left(C_{i,n}^{+}{+C}_{i,n}^{-}\right)={\displaystyle \sum }_{j=\mathrm{1,2,3}}{\displaystyle \sum }_{i=\mathrm{4,5}}{\displaystyle \sum }_{m=\mathit{1}}^M{f}_{\mathit{ij},\mathit{nm}}\left(D_{i,m}^{+}{+D}_{i,m}^{+}\right)(12)$

${\displaystyle \sum }_{i=\mathrm{1,2,3}}{\displaystyle \sum }_{n=\mathit{1}}^N{F}_{i,n}^h\left(C_{i,n}^{+}{-C}_{i,n}^{-}\right)={\displaystyle \sum }_{j=\mathrm{1,2,3}}{\displaystyle \sum }_{i=\mathrm{4,5}}{\displaystyle \sum }_{m=\mathit{1}}^M{f}_{\mathit{i}\mathit{j},\mathit{n}\mathit{m}}^h\left(D_{i,m}^{+}{-D}_{i,m}^{+}\right)(13)$

where $F_{i,n}$ and $f_{\mathrm{i}\mathrm{j},\mathrm{n}\mathrm{m}}$ are the CM, and their representations are as follows:

$F_{i,n}={\int }_0^{\mathit{2}\pi }{\int }_{R_{i-\mathit{1}}}^{R_i}{\psi }_{\mathit{it},n}^e{\psi }_{\mathit{it},n}^{e}r\mathit{drd}\theta ,\mathrm{i}=\mathrm{1,2,3}$(14–1)

$f_{\text{ji},\text{nm}}={\int }_0^{\mathit{2}\pi }{\int }_{R_{j-\mathit{1}}}^{R_j}{\psi }_{\mathit{jt},m}^e{\psi }_{\mathit{it},n}^{e}r\mathit{drd}\theta ,j=\mathrm{4,5}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{i}=\mathrm{1,2,3}$(14–2)

$F_{i,n}^h={\int }_0^{\mathit{2}\pi }{\int }_{R_{i-\mathit{1}}}^{R_i}{\psi }_{\mathit{it},n}^h{\psi }_{\mathit{it},n}^{h}r\mathit{drd}\theta ,\mathrm{i}=\mathrm{1,2,3}$(14–3)

$f_{\text{ji},\text{nm}}^h={\int }_0^{\mathit{2}\pi }{\int }_{R_{j-\mathit{1}}}^{R_j}{\psi }_{\mathit{jt},n}^h{\psi }_{\mathit{it},m}^{h}r\mathit{drd}\theta ,j=\mathrm{4,5}\mathrm{a}\mathrm{n}\mathrm{d}i=\mathrm{1,2,3}$(14–4)

where ${\psi }_{\mathrm{i}\mathrm{t},n}^e$ represents the transverse electric field component of the nth mode in the ith region of the axon region. ${\psi }_{\mathrm{i}\mathrm{t},n}^h$ represents the transverse magnetic field component of the nth mode in the ith region of the axon region. Meanwhile, ${\psi }_{\mathrm{j}\mathrm{t},m}^e$ represents the transverse electric field component of the mth mode in the jth region of the NR region. Finally, ${\psi }_{\mathrm{j}\mathrm{t},m}^h$ represents the transverse magnetic field component of the mth mode in the jth region of the NR region.

Therefore, the SM at the Node of Ranvier can be obtained by using the method in Ref. (Itoh, 1989), as follows:

${\mathit{S}}_{22\_\text{III}}={\left(\mathit{I}+{\mathit{M}}_{y\_\text{III}}^t{\mathit{M}}_{x\_\text{III}}\right)}^{-\mathit{1}}\left(\mathit{I}-{\mathit{M}}_{y\_\text{III}}^t{\mathit{M}}_{x\_\text{III}}\right)$(15-1)

${\mathit{S}}_{21\_\text{III}}=\mathit{2}{\left(\mathit{I}+{\mathit{M}}_{y\_\text{III}}^t{\mathit{M}}_{x\_\text{III}}\right)}^{-\mathit{1}}{\mathit{M}}_{y\_\text{III}}^t{\mathit{M}}_{x\_\text{III}}$(15-2)

${\mathit{S}}_{12\_\text{III}}={\mathit{M}}_{x\_\text{III}}\left(\mathit{I}+{\mathit{S}}_{22\_\text{III}}\right)$(15-3)

${\mathit{S}}_{11\_\text{III}}={\mathit{M}}_{x\_\text{III}}{\mathit{S}}_{21\_\text{III}}-\mathit{I}$(15-4)

where

${\mathit{M}}_{x\_\mathrm{I}\mathrm{I}\mathrm{I}}={\left({\mathit{Z}}_{\mathit{a}\mathit{x}\mathit{o}\mathit{n}}^{\mathit{2}}{\mathit{F}}_{\mathit{axon}}\right)}^{-\mathit{1}}\left({\mathit{f}}_{\mathit{N}\mathit{R}}{\mathit{Z}}_{\mathit{axon}}{\mathit{Z}}_{\mathit{NR}}\right)$(16-1)

${\mathit{M}}_{y\_\mathrm{I}\mathrm{I}\mathrm{I}}={\left({\mathit{Z}}_{\mathit{axon}}^{\mathit{2}}{\mathit{F}}_{\mathit{axon}}^h\right)}^{-\mathit{1}}\left(f_{\mathit{NR}}^h{\mathit{Z}}_{\mathit{axon}}{\mathit{Z}}_{\mathit{NR}}\right)$(16-2)

and I is the unit matrix. Each term in the CM is given by Eq. 14, and its matrix expression is as follows:

${\mathit{F}}_{\text{axon}}={\left[\begin{array}{ccc}{F}_{\mathrm{1,1}}{+F}_{\mathrm{2,1}}{+F}_{\mathrm{3,1}}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & F_{\mathit{1},N}{+F}_{\mathit{2},N}{+F}_{\mathit{3},N}\end{array}\right]}_{\left(N\times N\right)}$(17-1)

${\mathit{F}}_{\text{axon}}^h={\left[\begin{array}{ccc}{F}_{\mathrm{1,1}}^h{+F}_{\mathrm{2,1}}^h{+F}_{\mathrm{3,1}}^h& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & F_{\mathit{1},N}^h{+F}_{\mathit{2},N}^h{+F}_{\mathit{3},N}^h\end{array}\right]}_{\left(N\times N\right)}$(17-2)

${\mathit{f}}_{\text{NR}}={\left[\begin{array}{ccc}{f}_{\mathrm{41,11}}{+f}_{\mathrm{52,11}}{+f}_{\mathrm{53,11}}& \cdots & f_{\mathrm{41,1}M}{+f}_{\mathrm{52,1}M}{+f}_{\mathrm{53,1}M}\\ ⋮& \ddots & ⋮\\ f_{41,N\mathit{1}}{+f}_{52,N\mathit{1}}{+f}_{53,N\mathit{1}}& \cdots & f_{41,\mathit{NM}}{+f}_{52,\mathit{NM}}{+f}_{53,\mathit{NM}}\end{array}\right]}_{\left(N\times M\right)}$(17-3)

${\mathit{f}}_{\text{NR}}^h={\left[\begin{array}{ccc}{f}_{\mathrm{41,11}}^h{+f}_{\mathrm{52,11}}^h{+f}_{\mathrm{53,11}}^h& \cdots & f_{\mathrm{41,1}M}^h{+f}_{\mathrm{52,1}M}^h{+f}_{\mathrm{53,1}M}^h\\ ⋮& \ddots & ⋮\\ f_{41,N\mathit{1}}^h{+f}_{52,N\mathit{1}}^h{+f}_{53,N\mathit{1}}^h& \cdots & f_{41,\mathit{NM}}^h{+f}_{52,\mathit{NM}}^h{+f}_{53,\mathit{NM}}^h\end{array}\right]}_{\left(N\times M\right)}$(17-4)

${\mathit{Z}}_{\text{axon}}=\left[\begin{array}{ccc}{P}_{\mathit{1},\mathit{axon}}^{-0.5}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & P_{N,\mathit{axon}}^{-0.5}\end{array}\right]$(17-5)

${\mathit{Z}}_{\text{NR}}=\left[\begin{array}{ccc}{P}_{\mathit{1},\mathit{NR}}^{-0.5}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & P_{M,\mathit{NR}}^{-0.5}\end{array}\right]$(17-6)

where ${\mathit{Z}}_{\mathrm{a}\mathrm{x}\mathrm{o}\mathrm{n}}$ and ${\mathit{Z}}_{\mathrm{N}\mathrm{R}}$ are the normalized matrices, which can be determined by Eq. 10. The SM of the boundary IV can be solved by the same methods. Then, the SM of the entire NR can be obtained by the cascade formula of two-port networks (Chu and Itoh, 1986).

Results and discussion

Dispersion curves and field distributions of different modes

The dispersion curves of the guided modes in the nerve fiber can be obtained by solving Eq. 7. The results are shown in Figure 3. Here, the radii of the three regions of the model are set at $R_{\mathit{1}}=8\mu \mathrm{m}$, $R_{\mathit{2}}=10\mu \mathrm{m}$, and $R_{\mathit{3}}=20\mu \mathrm{m}$ (Kwon et al., 2017). The refractive indexes of the medium in the three regions are set as $n_{\mathit{1}}=2.1$, $n_{\mathit{2}}=\mathit{3}$, and $n_{\mathit{3}}=\mathit{2}$ (Liu et al., 2019), and the dielectric constants can be expressed as ${\epsilon }_i=n_i^{\mathit{2}}$. Meanwhile, we will further investigate the effect of the difference in the dielectric constant between the myelin sheath, and the internal and external fluids on the transmission of the THz-FIR wave.

FIGURE 3

FIGURE 3. (A) Dispersion curves of different guided modes in the myelinated nerve fibers within the frequency range of 0–30 THz, where m is the number of angular variations; (B) phase velocity of different guided modes in the myelinated nerve fibers within the frequency range of 0–30 THz; (C) dispersion curves of different guided modes in the unmyelinated nerve fibers within the frequency range of 0–30 THz; and (D) phase velocity of different guided modes in the unmyelinated nerve fibers within the frequency range of 0–30 THz.

B is the relative propagation constant, which can be represented as $B={\beta }_z/k_0$. $k_0$ is the vacuum propagation constant and is expressed as $k_{=0}=\omega \sqrt{{\epsilon }_0{\mu }_0}$. As shown in Figure 3, we only show the case of m = 0, 1, 2, and 3, but in practice the mode order n can take any value. Because the outer boundary of the nerve fiber is a metal tube wall with incomplete conductance, the tangential electric field on the outer boundary surface is not zero and the transmitted electromagnetic waves will no longer be pure TE mode or TM mode. Similar to the structure of optical fibers, the main transmission mode of terahertz waves in myelinated nerve fibers is mainly HEM mixed mode. According to the dominance of electric and magnetic fields, the modes obtained under the nerve fiber are denoted as the HEmn and EHmn hybrid modes when m > 1 (or TE0n or TM0n modes).

As can be seen from the dispersion curves in Figures 3A, C, both myelinated and unmyelinated nerve fibers transmit independent patterns, and each mode has a different propagation constant and corresponding transverse cutoff wavenumber at a given frequency. Hence, for any given excitation signal, the mode that satisfies the transmission characteristics can be found by the dispersion curve. Meanwhile, axons encased in high-refractive myelin sheaths greatly increase the number of modes that can be transmitted in nerve fibers compared to bare axons. By using Eq. 18, we further investigated the phase velocities of different modes of propagation in myelin-wrapped and bare axons. The calculated results are shown in Figures 3B, D. The transmission speed of THz waves in nerve fibers can reach 108 m/s, which is much higher than the propagation speed of action potential in myeloid axons (order of magnitude 102 m/s). Therefore, this unique biological structure is more conducive to supporting long-distance and high-speed communication between nerve cells.

${\mathit{v}}_p=\frac{\omega }{{\beta }_z}(18)$

To further analyze the number of guided modes that nerve fibers can transmit, we investigate the effect of the difference between the dielectric constant of the myelin sheath and internal fluid on the number of guided modes when f = 30 THz. The results are shown in Figure 4.

FIGURE 4

FIGURE 4. Difference between the dielectric constant of the myelin sheath and internal fluid on the number of guided modes at f = 30 THz.

The results show that the number of guided modes increases with the increase of the difference between the dielectric constants of the myelin sheath and internal fluid, which indicates that the nerve fiber can transmit THz signals to a greater extent in the frequency band with a higher permittivity of the myelin sheath. In addition, the lateral field distributions of different modes can be obtained by taking the results of Figure 3 into Eq. 6, the results are shown in Figure 5. Here, Figures 5A1–E1 show the electric field intensities and their vector diagrams of the first five modes of THz-FIR wave transmission at the nerve fiber port. Figures 5A2–E2 show the magnetic field intensities and their vector diagrams.

FIGURE 5

FIGURE 5. The lateral field distributions of different modes. (A1-E1) Electric field intensities and their vector diagrams of the first five modes; (A2-E2) magnetic field intensities and their vector diagrams of the first five modes.

The results show that the energies of electric and magnetic fields are mainly concentrated in the myelin region in both fundamental and higher modes. To compare the influence of different myelin dielectric constants and modes on the energy distribution, we conducted the following supplementary studies. The results are shown in Figure 6.

FIGURE 6

FIGURE 6. Effects of changes in the permittivity of myelin sheaths on energy distribution in axons, myelin sheaths, and fluid regions.

It can be seen from Figure 6 that the energy in the myelinated region increases as the difference in dielectric constants increases. At the same frequency, lower-order modes with high propagation constants concentrate energy more efficiently in the myelin sheath region.

Transmission characteristics of THz-FIR waves in nerve fibers

The complex permittivity of a polar solution can be expressed as:

$\epsilon ={\epsilon }^{\prime }-j\left({\epsilon }^{″}+\frac{\sigma }{\omega }\right)(19)$

where ${\epsilon }^{″}$ is the imaginary part of the permittivity, representing the polarization loss in the nerve fiber. Meanwhile, $\sigma$ is the electrical conductivity, representing the ohmic loss of the nerve fiber. During the simulation, we take ${\epsilon }^{″}=0.1$, $\sigma =30.9\mathrm{S}/\mathrm{m}$. The normalized Poynting vector of the longitudinal transmission THz-FIR waves in the nerve fibers are shown in Figure 7. It can be seen that myelinated nerve fibers can greatly reduce the loss caused by the absorption of polar liquids, such as water.

FIGURE 7

FIGURE 7. Transmission efficiency of different modes of THz-FIR waves in myelinated and unmyelinated nerve fibers at f = 30 THz.

We further investigate the effects of different myelin sheath thickness and dielectric constant on THz-FIR wave transmission. The theoretical calculation results of HE11 mode are shown in Figure 8.

FIGURE 8

FIGURE 8. (A) Changes in the permittivity of the myelin sheath, and $d_{myelin}=2\mu \mathrm{m}$; (B) changes in the thickness of myelin sheath, and $\Delta \epsilon =5$.

Figure 8 shows the relationship between the maximum attenuation of THz-FIR waves traveling along the 100 μm long axon with the difference between the dielectric constants of the myelin sheath and the internal and external fluids, and the thickness of the myelin sheath when the frequency is 30 THz. It can be seen that with the increase of the permittivity and thickness of the myelin sheath, the attenuation of THz-FIR waves is greatly reduced. This indicates that the bearing capacity of the nerve fibers in different positions in the organism for the transmission of the THz-FIR wave is different. Thicker nerve fibers (e.g., squid’s spinal nerve) facilitate long-distance transmission of THz-FIR wave.

Transmission characteristics at nodes of Ranvier

As an important physical parameter, the SM can be used to calculate the energy transfer process at the nodes of Ranvier. The calculation process of the SM is shown in Eqs 16 and 17. Here, we set the gap width of the nodes of Ranvier to 2 μm, the thickness of myelin sheath to 2 μm, the length of the axon to 100 μm, and the permittivity of three regions to ${\epsilon }_{\mathit{1}}=4$, ${\epsilon }_{\mathit{2}}=9$, and ${\epsilon }_{\mathit{3}}=4.41$, respectively. We then calculated the S11 curve and S21 curve of the node of Ranvier and compared them with the calculation results of the CST simulation software. The results are shown in Figure 9.

FIGURE 9

FIGURE 9. S-parameters of THz-FIR wave propagation at the nodes of Ranvier calculated by RNMMA and CST.

Compared with the result calculated by CST software, the results of the theoretical calculation are consistent with those from CST simulation software, which verifies the accuracy of our method. At f = 12.5 THz, S11 has the minimum value. This indicates that the reflection of THz-FIR transmission is the smallest at this time, which is mainly related to the width of the nodes of Ranvier. At the same time, S21 is greater than −5dB, which indicates that the THz-FIR wave has good transmission performance through the Ranvier junction. Meanwhile, the results also show that the S21 curve maintains a relatively stable value in the range of 10–30 THz and its value is greater than −5dB. This indicates that the node of Ranvier gap width (2 μm) has a weak influence on the transmission of THz waves. This happens because when the periodic gap width is much smaller than the wavelength of the THz-FIR wave, the THz-FIR wave can pass through the structure with minimal loss. This method provides theoretical support for further research on the interaction between THz-FIR waves and ion motion at the nodes of Ranvier.

Conclusion

In this article, we investigated the propagation characteristics of THz-FIR waves at the nerve fibers and the nodes of Ranvier by theoretical analysis. We aimed to study the new physical mechanism of high-frequency signal transmission in the nervous system. First, the exact expressions of the guiding mode and radiation mode in the nerve fiber have been obtained by theoretical analysis and the dispersion curves of the transmission of the THz-FIR wave in the nerve fiber have been obtained by applying proper boundary conditions. Myelinated nerve fibers can support more modes of transmission than unmyelinated bare axons, which greatly enhances the possibility that THz-FIR waves can travel through nerve fibers. It was found from the calculation results of field distributions that the myelinated nerve fibers could concentrate the field energy within the myelin region to a great extent, thus avoiding the strong absorption of THz-FIR waves by the surrounding polar liquids. Meanwhile, the thickness and permittivity of myelin sheath had a considerable influence on the propagation of THz-FIR wave along the nerve fiber. With the increase in the thickness and permittivity of the myelin sheath, the attenuation of the terahertz wave in longitudinal transmission gradually decreases.

We established a pattern-matching algorithm for the first time, named RNMMA, which was used to investigate the electromagnetic wave propagation and reflection characteristics of longitudinal discontinuities at the axon and the nodes of Ranvier. The accuracy of our algorithm has been verified by comparing it with the calculation results of CST.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

Author contributions

YG and SW designed the research and revised the manuscript. LG and DX conducted the study and wrote the manuscript. KW, YS, QZ, HN, and CL validated the data and helped revise the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (T224100002, 61921002, and 61988102).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Adolphs, R. (2015). The unsolved problems of neuroscience. Trends cognitive Sci. 19, 173–175. doi:10.1016/j.tics.2015.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Chao, T., Guang-Xu, Z., Jun, H., Jun-Hong, L., Wu, K., Ji, W., et al. (2021). Nonthermal and reversible control of neuronal signaling and behavior by midinfrared stimulation. Proc. Natl. Acad. Sci. U. S. A. 70, e2015685118. doi:10.1073/pnas.2015685118

CrossRef Full Text | Google Scholar

Chu, T. S., and Itoh, T. (1986). Generalized scattering matrix method for analysis of cascaded and offset microstrip step discontinuities. IEEE Trans. Microw. Theory Tech. 34, 280–284. doi:10.1109/tmtt.1986.1133323

CrossRef Full Text | Google Scholar

Dietrich and Marcuse (1970). Radiation losses of the dominant mode in round dielectric waveguides. Bell Syst. Tech. J. 49, 1665–1693. doi:10.1002/j.1538-7305.1970.tb04284.x

CrossRef Full Text | Google Scholar

Friede, R. L., and Bischhausen, R. (1982). How are sheath dimensions affected by axon caliber and internode length? Brain Res. 235, 335–350. doi:10.1016/0006-8993(82)91012-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, L., Bo, W., Wang, K., Wang, S., and Gong, Y. (2022). Theoretical investigation on the effect of terahertz wave on Ca2+ transport in the calcium channel. Iscience 25, 103561–110042. doi:10.1016/j.isci.2021.103561

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, Z.-H., Lv, W.-P., Hui, D.-X., Wang, X.-J., and Wang, Y.-N. (2022). Permeability enhancement of the KcsA channel under radiation of a terahertz wave. Phys. Rev. E 105, 024104. doi:10.1103/physreve.105.024104

PubMed Abstract | CrossRef Full Text | Google Scholar

Itoh, T. (1989). Numerical techniques for microwave and millimeter wave passive structures. New Jersey, U.S.: John Wiley & Sons.

Google Scholar

Kumar, S., Boone, K., Tuszynski, J., Barclay, P., and Simon, C. (2016). Possible existence of optical communication channels in the brain. Sci. Rep. 6, 36508. doi:10.1038/srep36508

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwon, J., Kim, M., Park, H., Kang, B.-M., Jo, Y., Kim, J.-H., et al. (2017). Label-free nanoscale optical metrology on myelinated axons in vivo. Nat. Commun. 8, 1832–1839. doi:10.1038/s41467-017-01979-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Laman, N., Harsha, S. S., Grischkowsky, D., and Melinger, J. S. (2008). High-resolution waveguide THz spectroscopy of biological molecules. Biophysical J. 94, 1010–1020. doi:10.1529/biophysj.107.113647

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Chang, C., Zhu, Z., Sun, L., and Fan, C. (2021). Terahertz wave enhances permeability of the voltage-gated calcium channel. J. Am. Chem. Soc. 143, 4311–4318. doi:10.1021/jacs.0c09401

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, G. (2018). The conjectures on physical mechanism of vertebrate nervous system. Chin. Sci. Bull. 63, 3864–3865. doi:10.1360/n972018-01143

CrossRef Full Text | Google Scholar

Liu, G., Chang, C., Qiao, Z., Wu, K., Zhu, Z., Cui, G., et al. (2019). Myelin sheath as a dielectric waveguide for signal propagation in the mid‐infrared to terahertz spectral range. Adv. Funct. Mat. 29, 1807862. doi:10.1002/adfm.201807862

CrossRef Full Text | Google Scholar

Liu, X., Qiao, Z., Chai, Y., Zhu, Z., Wu, K., Ji, W., et al. (2021). Nonthermal and reversible control of neuronal signaling and behavior by midinfrared stimulation. Proc. Natl. Acad. Sci. U. S. A. 118, e2015685118. doi:10.1073/pnas.2015685118

PubMed Abstract | CrossRef Full Text | Google Scholar

Marcuse, D. (1973). Coupled mode theory of round optical fibers. Bell Syst. Tech. J. 52, 817–842. doi:10.1002/j.1538-7305.1973.tb01992.x

CrossRef Full Text | Google Scholar

Miller, L. M., and Dumas, P. (2010). From structure to cellular mechanism with infrared microspectroscopy. Curr. Opin. Struct. Biol. 20, 649–656. doi:10.1016/j.sbi.2010.07.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Tegmark, M. (2015). Consciousness as a state of matter. Chaos, Solit. Fractals 76, 238–270. doi:10.1016/j.chaos.2015.03.014

CrossRef Full Text | Google Scholar

Tononi, G. (2004). An information integration theory of consciousness. BMC Neurosci. 5, 42–22. doi:10.1186/1471-2202-5-42

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, K., Wang, S., Yang, L., Wu, Z., Zeng, B., and Gong, Y. (2022). THz trapped ion model and THz spectroscopy detection of potassium channels. Nano Res. 15, 3825–3833. doi:10.1007/s12274-021-3965-z

CrossRef Full Text | Google Scholar

Wetzel, D. L., and Levine, S. M. (1999). Imaging molecular chemistry with infrared microscopy. Science 285, 1224–1225. doi:10.1126/science.285.5431.1224

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, K., Qi, C., Zhu, Z., Wang, C., Song, B., and Chang, C. (2020). Terahertz wave accelerates DNA unwinding: A molecular dynamics simulation study. J. Phys. Chem. Lett. 11, 7002–7008. doi:10.1021/acs.jpclett.0c01850

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiang, Z., Tang, C., Yan, L., Chang, C., and Liu, G. (2020). A new viewpoint and model of neural signal generation and transmission: Signal transmission on myelinated neuron. Chin. Phys. B 29, 108701. doi:10.1088/1674-1056/abaee9

CrossRef Full Text | Google Scholar

Zangari, A., Micheli, D., Galeazzi, R., and Tozzi, A. (2018). Node of ranvier as an array of bio-nanoantennas for infrared communication in nerve tissue. Sci. Rep. 8, 539. doi:10.1038/s41598-017-18866-x

PubMed Abstract | CrossRef Full Text | Google Scholar