
==== Front
Open Biol
Open Biol
RSOB
royopenbio
Open Biology
2046-2441
The Royal Society

rsob240138
10.1098/rsob.240138
18130Methods and Techniques
Methods and Techniques
In silico modelling of neuron signal impact of cytokine storm-induced demyelination
In silico modelling of neuron signal impact of cytokine storm-induced demyelination
Adonias Geoflly L. 1 Formal analysis Methodology Software Validation Visualization Writing – original draft Writing – review and editing mail@gadonias.com

Siljak Harun 2 Formal analysis Methodology Resources Software Validation Writing – original draft Writing – review and editing harun.siljak@tcd.ie

Balasubramaniam Sasitharan 3 Conceptualization Funding acquisition Investigation Project administration Supervision Writing – original draft Writing – review and editing sasi@unl.edu

https://orcid.org/0000-0002-9765-7660
Barros Michael Taynnan 4 Conceptualization Investigation Methodology Supervision Writing – original draft Writing – review and editing m.barros@essex.ac.uk

1 Walton Institute for Information and Communication Systems Science, South East Technology University , Waterford, Ireland
2 Department of Electronic and Electrical Engineering, Trinity College Dublin , Dublin, Ireland
3 School of Computing, University of Nebraska-Lincoln , Lincoln, NE, USA
4 School of Computer Science and Electronic Engineering, University of Essex , Colchester, UK
2024
04 9 2024 September 4, 2024
04 9 2024 September 4, 2024
14 9 24013822 5 2024 May 22, 2024
16 7 2024 July 16, 2024
17 7 2024 July 17, 2024
© 2024 The Author(s).
2024
https://creativecommons.org/licenses/by/4.0/ Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

In this study, we develop an in silico model of a neuron’s behaviour under demyelination caused by a cytokine storm to investigate the effects of viral infections in the brain. We use a comprehensive model to measure how cytokine-induced demyelination affects the propagation of action potential (AP) signals within a neuron. We analysed the effects of neuron-neuron communications by applying information and communication theory at different levels of demyelination. Our simulations demonstrate that virus-induced degeneration can play a role in the signal power and spiking rate, which compromise the propagation and processing of information between neurons. We propose a transfer function to model the weakening effects on the AP. Our results show that demyelination induced by a cytokine storm not only degrades the signal but also impairs its propagation within the axon. Our proposed in silico model can analyse virus-induced neurodegeneration and enhance our understanding of virus-induced demyelination.

neuron
; action potential
; demyelination
; molecular communications
; cytokine storm
; Hodgkin–Huxley
National Science Foundation (NSF) 2316960 Balasubramaniam Sasitharan
==== Body
pmc1. Introduction

The recent outbreak of the coronavirus disease 2019 (COVID-19) pandemic caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has shaken society as a whole, leaving a long-lasting impact on people’s health. As a result, scientists from multidisciplinary fields have been joining efforts and resources towards the study of not only the epidemiological characteristics and transmission dynamics of the virus, but also of the physiological damage to the human body as the infection is prolonged. SARS-CoV-2 is well known for affecting primarily the respiratory system, and it can potentially leave lifelong sequelae in tissues and organs. Furthermore, there has been an increase in works that suggest SARS-CoV-2 may be able to invade the nervous system [1–4] and elicit neurodegeneration [5,6] where sequelae may have other detrimental effects on patients’ lives post-infection.

Viruses that present the ability to infect nerve cells are known to exhibit neurotropic properties and can also be called neuroinvasive. By infecting cells in the nervous system and replicating themselves within it, these viruses can negatively impact neurological functions and even cause severe nerve damage by triggering a pro-inflammatory immune response [6]. Unfortunately, SARS-CoV-2 is not the only virus that exhibits this kind of behaviour; for example, it has been shown that the Zika virus (ZKV) [7] can infect the peripheral nervous system (PNS) and sometimes spread to the central nervous system (CNS). Furthermore, viruses such as the human immunodeficiency virus (HIV) [8] can infect the CNS and cause neuroinflammation, induced by an immune response of the body, consequently leading to neurodegeneration [9]. A growing amount of data from the past few years support the hypothesis that chronic damage caused by different infectious agents can lead to neurodegeneration [10], as early experimental models of virus-induced demyelination have indicated [11].

Viruses are known for causing dramatic structural and biochemical changes to the host cell, by hijacking and exhausting its machinery for replication until, eventually, the cell is killed. This viral manipulation of a host cell provokes neuroinflammatory defence mechanisms that can be characterized by numerous toxic-metabolic derangements, such as cytokine storms (figure 1a ). Such inflammation could potentially lead to several types of neurodegeneration, including demyelination. As cytokines are released to fight the infection, healthy tissues could be affected as a ‘collateral damage’ of the fight against infectious agents [12]. Other types of coronavirus have been known to cause demyelination, such as the murine coronavirus (M-CoV), which has been identified to cause demyelinating disease and, even after the virus is cleared from the CNS, the demyelination can continue for a few months [12]. This behaviour also matches findings on SARS-CoV-1, which reports a decrease in viral titers as clinical disease worsens [13]. M-CoV is a type of coronavirus of the same genus (Betacoronavirus) as SARS-CoV-2, and it is believed to be 43.8–48% similar to this novel coronavirus [14]. SARS-CoV-2, when compared to SARS-CoV-1, triggers lower levels of interferons and pro-inflammatory cytokines and chemokines. However, it is capable of infecting and replicating a significantly higher amount of viruses in human tissues [15].

Figure 1. Pathogenesis of a virus-induced demyelination. (a) Depiction of the viral entry into the respiratory system, passage through the BBB and effects of demyelination to the neurons. (b) The system model captured by our model in this paper.

Pathogenesis of a virus-induced demyelination.

This article presents a systems theory-based analysis of the demyelination (either directly or indirectly) caused by viral infections from the perspective of molecular communication (MC), and more specifically a neuro-spike MC. The goals of this article are: (i) to propose a model extension coupled with computational analysis on the Hodgkin–Huxley (HH) formalism to account for the effects of cytokine storms indirectly caused by viral infections; (ii) to provide insights on how the demyelination will affect the neuronal information along the axonal pathway; and (iii) to provide a transfer function (TF) model that describes the effects of the membrane action potential (AP) caused by the degeneration of the myelin sheath. This TF can be considered fundamentally as a model of the demyelination process itself. It is intrinsically tied with the behaviour of equivalent resistance–capacitance (RC) circuits, but at the same time has a reduced complexity as it opts for exponential asymptotics. Given that complex cascades of equivalent simple blocks (in our scenario, myelin sheaths) are traditionally well approximated by first-order TFs with time delay and the underlying physical rationale, our model joins the family of established, applicable biophysical TF models. We expect that these analyses not only could pave the way for more in-depth studies that can support in vitro and in vivo experimental works on the neurological effects induced by neurotropic viral infections but also serve as a prediction tool that can help to guide future experimental works. It may even be necessary that, as the demyelinating scenarios get more and more biologically plausible, future works may have to fine-tune some of the parameters (e.g. regarding the type of cell or the intensity of cytokine storm-induced demyelination) of the proposed model in order to have as accurate an analysis as possible.

The contributions of this work are as follows:

— A mathematical model that serves as a building block for a TF of the signal propagation. We present a model inspired by the well-documented fact that viral infections trigger cytokine storms and lead us towards a general TF that accounts for the propagation of the signal. The cytokine storms are counter-measures of the immune system against the infection. We investigate the signal power, the signal attenuation and magnitude-squared coherence (MSC) analysis of the axonal pathway as a communication channel.

— An analysis of the effects of demyelination on the neuronal AP propagation. We conduct a variety of analyses, such as latency, attenuation and spiking rate, on the spike trains that pass through a demyelinated pathway. Furthermore, we also analyse how the intensity of cytokines storms correlates with the amount of attenuation present in the signal at the output of the neuron.

— A TF model that accounts for the effects of the demyelination-induced attenuation. We propose a transfer function that describes the transition of a healthy to an unhealthy neuron. The model accounts for sheath-by-sheath demyelination which affects membrane potential, peak times and spike width. This should lead to more in-depth analysis and open a new view on demyelination modelling, especially those triggered by a neurotropic viral infection and how it affects the transmembrane molecular exchange (i.e. exchange of ions and release of neurotransmitters) of the neurons.

2. Background and literature review

The virus replication process exhausts host cells leading to the activation of the immune system, which calls for the work of macrophages. Those are types of cells of the immune system, and their resident in the brain is the microglia [1,4]. Macrophages secret pro-inflammatory cytokines such as interleukin-1 (IL-1), interleukin-6 (IL-6) and tumour necrosis factor-alpha (TNF- α ) aiming to fight the infection. These cytokines exert cytotoxic effects on neurons and glial cells (e.g. oligodendrocytes), damaging the myelin sheaths which are responsible for providing support and insulation to axons [16]. The dynamics of cytokine storms and their efficiency in fighting infections without further cellular damage become of increased attention to measure the neuroinflammatory effects of COVID-19. Ludwig et al. [17] found that as the neuropathy becomes more severe, the serum concentration of TNF- α and IL-6 becomes higher. This matches more recent results found by Chen et al. [18] and Song et al. [19] with respect to cytokine levels concerning the severity of COVID-19 cases. They identified a positive correlation between the levels of cytokines and the severity of the COVID-19 cases, meaning that the more severe the COVID-19 cases were, the higher were the levels of cytokines, which can lead to demyelinating lesions [5,20].

On the other hand, the dynamics of a cytokine storm can be similar for different inflammatory scenarios, yielding the need for a more general approach that evaluates cytokine dynamics. An example is the work of Waito et al. [21] in which they match recordings of 13 different cytokines with a model of non-linear ordinary differential equations. Likewise, Yiu et al. [22] analysed the dynamics of cytokine storms and provided evidence for how cytokines induce or inhibit other cytokines. Additionally, some pieces of literature report the effects of cytokine storms on specific neuronal and non-neuronal structures. For instance, the work of Bitsch et al. [23] shows a negatively correlated relationship between the amount of microglia-produced TNF- α and the concentration of myelin oligodendrocyte glycoprotein (MOG). It shows that the more TNF- α there is, the less MOG oligodendrocytes will produce, compromising the myelin sheath structure. On the other hand, Redford et al. [24] showed how the number of axons found in sciatic nerves is affected. They found that the more the concentration of TNF- α is increased, the more axons are found to be damaged. However, the complete biophysical models of these various effects to neurons caused by infection, particularly infections from COVID-19, require urgent attention since correct treatment procedures for acute infection damage can benefit from mathematical modelling.

In the past decade, MC has been improving biological models by accounting for the communication of cells using their signalling mechanisms as information carriers [25–27]. MC is a new field that is looking to characterize and engineer biological cells using concepts from communication engineering and networking [28–30], it bridges electrical and communications engineering, molecular biology and biomedical engineering, and provides complete end–end models of biophysical transmission of molecules, their propagation, and their reception [31–35]. A recent survey [36] reports numerous works concerning the use of MC for the analysis and modelling of infectious diseases. However, accounts for the effects of infections are still missing from a biophysical approach even within MC models, as the COVID-19 effects are many, and the molecular interactions with the body can be used to predict the behaviour of a population of cells, even tissues and organs. The emergence of novel MC models for biophysical processes, such as the demyelination induced by COVID-19, delivers an in-depth analysis of tissue behaviour that is needed for treatments based on synthetic biology. Even further targeted drug delivery technology can alleviate the cytokine storms effects on neurons and provoke the restoration of the myelin sheath [37].

3. End-to-end computational model for a cytokine storm-induced demyelination

Models that describe the dynamics and evolution of cytokine concentrations, without regard to the cells that secrete or are affected by them, have been proposed by Yiu et al. [22]. This model was extended by coupling together with a neuronal model that implements the behaviour of a myelinated axon [38]. By extending these models, we performed simulations with the default parameters of the cell (such as the length and diameter of each of their compartments) based on the original model. Therefore, the myelin sheath properties concerning their morphological characteristics were not modified or changed. We studied and quantified the propagation of APs on a Layer 5 (L5) pyramidal neuron mimicking damage to its myelin sheath. All simulations were performed with extensions to the NEURON Simulator [39].

3.1. Cytokine signalling in microglia

The growth and decay of an individual cytokine’s response to its given initial state are first represented by a second-order, linear, time-invariant ordinary differential equation. Denoting the serum concentration ρ(t) in pg ml−1, and its rate of change Δρ(t) , it can be represented in a vector-matrix form as follows:

(3.1) [ρ˙(t)Δρ˙(t)]=[01−a−b][ρ(t)Δρ(t)],[ρ(0)Δρ(0)]given,

where the initial concentration, ρ(0)=0 , is referenced to the cytokine’s basal level, and the initial rate of change, Δρ(0) , is stimulated by the TGN1412 infusion (refer Yiu et al. [22] for further details). Also, a and b are positive constants that express the sensitivity of the cytokine’s acceleration to concentration and rate of change.

The cytokine’s response modes are characterized by the eigenvalues, λ1 and λ2 (rad day−1), of the stability matrix of the system, and this is as follows:

(3.2) [ρ˙(t)Δρ˙(t)]=[01−λ1λ2(λ1+λ2)][ρ(t)Δρ(t)],[ρ(0)Δρ(0)]given.

The parameters are chosen to minimize the error between the cytokine concentration and the clinical trial measurements performed by Yiu et al. [22]. For TNF- α , λ1=λ2=−2.63 and Δρ(0)=32821 . In this work, we will be investigating the inflammatory effects of TNF- α , as this cytokine is well-known for its pro-inflammatory properties.

3.2. Conduction through a myelinated axon

It is well known that some neurons contain a myelin sheath wrapped around sections of their axons. Myelin sheath helps propagate electrical impulses, known as APs or spikes, and avoid significant attenuation due to parallel synaptic processes. According to Cohen et al. [38], the myelin sheath dynamics can be described using circuit theory (figure 2). This circuit is then coupled with the HH circuit model [40], which describes the membrane potential dynamics in neurons.

Figure 2. Schematic of the equivalent myelinated axon circuit, with added periaxonal ( Pa ) and paranodal ( Pn ) axial resistances [38].

Schematic of the equivalent myelinated axon circuit.

Axial resistance, Ri ( Ω⋅cm−1 ), of the axon core, can be defined as the ratio between axial resistivity, ri ( Ω⋅cm ), and the cross-sectional area of the axon core, and this is described as

(3.3) Ri=4riπd2,

where d (nm) is the axon core diameter. Let δpa (nm) be the radius of the periaxonal space, then

(3.4) δpa=12[−d+d2+(4rpaπRpa)],

where the axial resistance in the periaxonal space, Rpa (G Ω⋅cm−1 ), is calculated as the ratio between periaxonal resistivity, rpa ( Ω⋅cm ), and periaxonal cross-sectional area. In this case, the axon core cylinder is surrounded by both the periaxonal and perinodal spaces forming a ‘larger’ axon cylinder of diameter d+2δpa , thus

(3.5) Rpa=rpaπδpa(d+δpa),

and an analogous calculation can be performed for δpn (nm), rpn ( Ω⋅cm ) and Rpn (T Ω⋅cm−1 ).

Recognizing that a myelin sheath is an in-series compaction of n layers, the radial resistance of the sheath, Rmy (k Ω⋅cm2 ) is the sum of the resistance of each myelin membrane, Rmm (k Ω⋅cm2 ), formulated as

(3.6) Rmy=∑i=1nRmmi,

and the radial capacitance of the myelin sheath, Cmy ( μF⋅cm−2 ), may vary inversely to the sum of the capacitances of each of its composing membranes, Cmm ( μF⋅cm−2 ). Thus,

(3.7) 1Cmy=∑i=1n1Cmmi,

where the resistance and capacitance of a single myelin membrane are Rmm and Cmm , respectively. Based on this, equations (3.6) and (3.7) can be represented in terms of the number of myelin lamellae, nmy , as follows:

(3.8) nmy=Rmy2Rmm=Cmm2Cmy.

For a more detailed analysis of the myelin sheath’s modelling, the reader is referred to the work of Cohen et al. [38]. Each of the lamellae that compose the whole myelin sheath is modelled with the same set of parameters. This is basically an assumption that each myelin lamellae is identical to the other and that the degeneration caused by the demyelination affects the myelin sheath proportionally.

3.3. Cytokine-induced demyelination

As previously discussed in figure 1, the literature indicates that, as microglia cells release pro- and anti-inflammatory cytokines to fight infection, demyelination may also occur as a side effect and, consequently, compromise the neuronal signal propagation. There are numerous pieces of evidence linking cytokine storms to neurodegeneration [17,21–23]; however, to the best of our knowledge, none of them goes as far as linking the storm’s intensity to an approximate number of myelin lamellae, nmy . A linear regression ( R=0.508 , p=0.001 ) applied by Ludwig et al. [17] to their own data reveals a proportional relationship between the severity of the neuropathy, ζ (in our scenario, the neuropathy is the demyelination) and the serum concentration of TNF- α , ρ(t) , as

(3.9) ζ(t)=20.727−0.9228⋅ρ(t).

Even though the cytokine dynamics by Yiu et al. [22] are triggered by a specific antibody, the proportionality presented in equation (3.9) is also suggested by the works of Hartung [41] and Empl et al. [42]. In this relationship, the stronger the cytokine storm is, the more severe the degeneration of the myelin sheaths. Furthermore, in this work, we are using data from Cohen et al. [38] as a reference to indicate healthy myelin with 13 lamellae. In the study by Ludwig et al. [17], the severity of the neuropathy was defined by scoring four nerve functions as ‘ 0 ’ for typical values of the nerve, ‘ 1 ’ for affected nerves (either decreased amplitude and/or decreased nerve conduction velocity) and ‘ 2 ’ for no stimulation possible. As we are looking at it from the perspective of nmy , we had to re-scale these scores by applying linear interpolation to indicate approximately how many sheaths are being degraded. Therefore, building on top of their findings we hypothesize that the worst-case scenario would have a score of 8 (four nerve functions multiplied by two—no stimulation possible) which indicates nmy=0 and the best-case scenario would have a score of 0 indicating nmy=13 . Similar experiments conducted by the scientific community potentially on top of this work can re-scale ζ(t) depending on their reference of normal myelination as the number of myelin sheaths may vary among different types of neurons from different parts of either the CNS or the PNS.

The model proposed by Ludwig et al. [17] addresses neuropathies in the PNS in which peripheral neurons are myelinated by Schwann cells. On the other hand, this work examines neuropathies in the CNS where myelination is mainly controlled by oligodendrocytes. Even though oligodendrocytes and Schwann cells can express different types of myelin proteins, a major difference would be that oligodendrocytes can myelinate several axons at the same time while Schwann cells can wrap around only a single axon at a time [43,44]. In terms of myelination processes, both types of cells are quite similar and the model itself already leaves room for improvement as the parameters may be finely tuned as soon as new evidence is presented in the literature.

3.4. A linear model of demyelination

The linearization process aims to find a linear approximation of a nonlinear system (i.e. the HH model) at an equilibrium point. This means that for small variations around said point, the linear system should behave similarly to the nonlinear system [45]. The linear model is the result of the process of linearization applied to the conventional (non-linear) HH model. In this process, the dynamics of the ionic channels are ‘lost’, the components that represent each channel (e.g. variable conductance and a voltage source) are ‘reduced’ to a static conductance and an inductance. As the demyelination only affects the RC circuits coupled to the internodal compartments, the electronic elements for the axon should remain in their default parameters and the myelin sheath RC circuits be affected by the inflammatory effects of a cytokine storm. Based on a given healthy neuron, we want to predict the signal effects on spike trains that are expected from a demyelinated neuron. Let n , the number of the myelin sheaths, be a tuneable parameter of the model (in this section, we refer to nmy as n for clarity of index notation). This mechanism is depicted in the block diagram in figure 1b .

In telecommunication systems, the study of damage and degradation often asks for a model of damage effects which could explain how healthy (i.e. signals characteristic to the system without damage) transition into faulty signals (i.e. signals characteristic to the damaged system). Observing faulty signals resulting from the propagation of standard signals through a TF superimposed on the original system is a well-established concept; examples of its application include modelling of cables [46]. It makes intuitive sense to observe such a TF as taking the healthy system’s output and delivering the faulty, unhealthy signal as the faulty output, which results in signal degradation for a worsening channel it traverses. In our case, this means that we will construct a TF which, for an input that represents a healthy neuronal signal, delivers an output equivalent to that of a demyelinated neuron. The opposite process, in which an ‘adaptor’ TF would convert a signal from a deteriorated neuron into one of a healthy neuron, is anti-causal as it would have to introduce negative time delays in the signal.

In this regard, when we speak of our model, we have the TF in mind. Our model’s input is the series of spikes produced by a healthy pre-synaptic neuron and the model’s output is the series of spikes a demyelinated neuron would fire, as illustrated in figure 1b . The model parameters are factors in the TF; their value depends on the number of myelin sheaths we desire the demyelinated neuron to have. Hence, the number of myelin sheaths is a tuneable scalar value characterizing the model. Now that we have established the processing chain, the choice for the TF is made by observing the general trends in the output signals for various values of n . Namely, as shown in figure 3, the TF emulating the effect of myelin deficiency needs to allow for the attenuation of the signal, widening of the spikes and overall propagation delay. This delay is caused by the demyelination process and the proposed TF does not account for random delays [47]. An obvious candidate is the traditional first order plus time delay (FOPTD) function [48], given by

Figure 3. Mapping of the model components and the effects on the propagated signal, and illustration of signal behaviour as it transitions from a healthy state to an unhealthy state.

Mapping of the model components and the effects on the propagated signal.

(3.10) Wn(s)=kn1+Tnse−τns,

where we assume that the parameters kn , Tn and τn have different values for different values of n . While those values could be tabulated and looked up for specific values of n , we set a more ambitious goal of determining the laws according to which they change. Here, we make a hypothesis that they follow exponential law (namely, log⁡kn=a0⋅arn , Tn=To⋅Trn and τn=τ0⋅τrn ), based on the following reasoning.

Turning a healthy signal into a deteriorated one in our model is a process of cancelling the effect of passing through some of the sheaths. For example, if the healthy baseline signal is achieved with n=13 , emulation of n=10 can be interpreted as undoing the effect of Δn=3 sheaths (i.e. the effect of a cascade of three blocks representing a ‘sheath undo function’). This approximation is good for larger Δn , as it allows for the approximations like (1+Ts‾)Δn≈1+TΔns , i.e. giving rise to a Tn∼en law. In §4, we revisit this hypothesis, both in terms of the exponential law’s existence and in terms of the domain of accuracy. Introducing exponential laws for the coefficients in this TF gives us the final form of our TF, which is represented as

(3.11) Wn(s)=ea0⋅arn1+T0⋅Trnse−τ0⋅τrns=ea0⋅arn−τ0⋅τrns1+T0⋅Trns.

The final form of equation (3.11) suggests the rationale behind suggesting exponential behaviour of a=log⁡k instead of k itself: in the s -domain, we obtain a function of the form e(α+βs)γ+δs , which in turn in the Fourier domain corresponds to e(α+jβω)γ+jδω . Knowing that the behaviour of the myelin circuit originates from connections of R (purely passive, real impedance) and C (purely active, imaginary impedance), it is expected to observe this symmetric real-imaginary coupling of terms.

The knowledge about the equivalent model of the myelin circuit, a ladder of resistors and capacitors supports the choice of FOPTD, and this is known as a convenient model for large RC circuits [49]. Namely, the circuit consists of linear components and as such can be represented by a linear model. Furthermore, with the increasing order of such a linear model, there is a necessity to replace it with a simpler first-order model such as FOPTD as it grows and will keep the complexity of the model low while retaining accuracy.

FOPTD is not an uncommon choice in biophysics, where it has been used to model glucose control [50], because they are simple for identification [48] and for quick and accurate tuning of the controllers that can regulate their behaviour [51]. Nonetheless, our model may help in designing chemical control loops for myelin reinforcement in a similar manner.

To verify the quality of the model, we introduce a metric based on root mean square error (RMSE):

(3.12) Mn=20log10⁡(RMSEN,nRMSEW,n).

Here, 𝑅𝑀𝑆𝐸W,n stands for the RMSE of the output of our TF Wn compared to the actual output signal for n sheaths if m is the number of samples, which is represented as

(3.13) RMSEW,n=1m∑i≥1(xn,i−xW,i)2.

Analogously, we can find the value for 𝑅𝑀𝑆𝐸N,n which corresponds to the RMSE of the output signal produced by another model N compared to the n th actual output signal. The quantity Mn is positive where our model Wn is more accurate (has lower RMSE) than the model N we are comparing to it.

Certainly, the parameters might eventually change; however, qualitatively speaking, changes in axon structure from different neurons should exhibit similar behaviour to each other. For example, some delay on signal propagation may be introduced due to different axonal lengths [52], but it should not affect the shape of the APs as it is the case with axons that have been demyelinated. Furthermore, the TF works on sub-threshold neuronal signalling and this stimulation technique usually does not lead to the firing of APs even though it can still fire depending on the intensity of the stimuli. Our aim is to apply a similar approach to Khodaei & Pierobon [53] to study the impact on the signalling caused by the demyelination that may have been indirectly caused a viral infection.

3.5. Signal analysis

Visually, we first noticed subtle shifts in amplitude (peak potential reached by the membrane) and in time (spikes were taking longer to reach their peak values). We then decided to quantify these shifts, both in amplitude and in time, on average, by proposing a metric that we called 'relative mean shift'. For the analysis on time shift, we consider the points in time where each spike peaked at the input as Tink , where k={1,2,3,...,K} identifies the order of each spike and, Toutk as the peak times at the output. Thus, we define the relative mean time shift, δt¯ (ms), as

(3.14) δt¯=1K∑k=1K(Toutk−Tink).

Analogously, we can define the relative mean amplitude shift, δv¯ (mV), with Vink and Voutk as the peak amplitudes of spike k at the input (spike peak observed in the soma) and output (spike peak observed in the axon), respectively.

As there is an average shift in time inside the channel itself due to demyelination, we also expect an increase in latency as we decrease nmy . In order words, latency is the time interval between the input and the output, and it often occurs due to the channel or network’s intrinsic characteristics. In our scenario, we are looking for the time the first spike peaked at the output, concerning the point when this same spike peaked at the soma of the neuron. Later, analogous to the relationship between δt¯ and the latency, we decided to look into a potential attenuation on the power of the signal, P (mW), as we have already discussed a relatively heavy mean attenuation in the membrane potential, δv¯ . Let us define P as

(3.15) P=limJ→∞(12J+1∑j=−JJ|x[j]|2),

where in a set of J samples, x[j] corresponds to the potential of the membrane at the jth sample.

We also quantified the attenuation of the signal, A in dB/ 100 m, which is the gradual loss of power of a signal over its propagation through the channel. Depending on the attenuation coefficient, one can calculate more accurately the attenuation in a specific material. For our analysis, we used the generic form of attenuation for RF cables. This decision is based on the fact that the axonal pathway is modelled using cable theory as a leaky cable. Thus,

(3.16) A=10⋅log10⁡(PiPo),

where Pi and Po , in W , are the input and output power of the signal. In this analysis, the input is the spike train through a healthy neuron and the output would be the spike train through a demyelinated one.

Finally, we also analysed the relation between the reference spike train ( nmy=13 ) and all other demyelinating scenarios ( 1≤nmy≤12 ). The intention is to understand how much power is being transferred between each pair of signals. With that in mind, we applied a coherence ( Cxy ) metric, which is described as

(3.17) Cxy(ω)=Sxy(ω)2Sxx(ω)Syy(ω),

where Sxy(ω) is the cross-spectral density between the two signals, and Sxx(ω) and Syy(ω) are the power spectrum densities of input and output, respectively.

4. Computational results and discussion

Let us consider the proportionality between the severity of the neuropathy, ζ(t) , and the number of myelin lamellae, nmy , as described in §3.3. Our analysis consisted of an evaluation of the neuronal behaviour and spike propagation under normal circumstances ( nmy=13 ), followed by an analysis on the demyelination by decreasing nmy . The objective is to understand what happens to the neuronal information when travelling through a demyelinated axonal pathway. In other words, we are looking at the signalling within a single neuron that has been impacted by the demyelination on its own, rather than signals that have changed from pre-synaptic neurons that may be affected by cytokine storms. The neuron receives an external current of 3 nA for 15 ms starting at time t=2.5 ms in a 20 ms simulation. The spikes evoked under normal circumstances are considered our input of the channel, and the spikes at the far end of the axon are the output. The axon itself is our communication channel, while the demyelination is acting as an attenuator for the channel, as illustrated in figure 1.

4.1. Analysis of a demyelination-induced channel attenuation

The results for δv¯ and δt¯ from equation (3.14) are shown in figure 4a , d , respectively. The former shows that membrane potential is on average affected by an eightfold decrease. This matches fundamental computational neuroscience theory on neuronal modelling which states that a neuron gradually leaks a small amount of the input signal as it travels through it, and this ‘leak’ worsens on unmyelinated cables [54,55]. On the other hand, from the latency data (figure 4e ), we note that there is a subtle ‘lag’ for the signal to travel across the axon even under normal circumstances matching findings in the literature [56] for models of CNS demyelination. This is most likely due to a maximum conduction speed inherent in the axonal membrane itself. As we go from our worst scenario towards a regular healthy myelin sheath, there is a massive decrease of about 73% in the latency.

Figure 4. Measurements of attenuation and delay. All results except for (e) use the signal nmy=13 as their input and 1≤nmy≤12 as the output; (e) uses the recordings from the soma of the neuron and is compared with 1≤nmy≤13 . (a) Relative mean amplitude shifts, (b) spiking rate, (c) signal power, (d) relative mean time shift, (e) relative mean latency and (f) attenuation.

Measurements of attenuationMeasurements of attenuation and delay.

Figure 4b depicts how the spiking rate is affected by the demyelination. As we expected, as the spikes start to get wider and further from each other, the spiking rate gets lower. This corroborates findings on demyelination-induced effects on spiking rate [57]. The rate at which a neuron fires APs is significant for modulating and encoding neuronal information in cognitive, sensory and motor functions.

Figure 4c shows that the decrease in signal power is quite subtle, and from the worst-case scenario to the best, there is a difference of less than 20μW . This indicates how demyelination affects the energy consumption per unit time used to propagate the axon’s APs. As the signal starts to get degraded, it is less and less likely a spike would be evoked at the post-synaptic neurons connected to a demyelinated cell [55]. Demyelination does not affect only the post-synaptic neuron by reducing the chances of evoking post-synaptic potential, but it can compromise the spiking rate of the demyelinated neuron itself. Furthermore, we decided to investigate the attenuation caused by the demyelination from a more generic communication systems point of view as expressed in equation (3.16). The results presented in figure 4f show how the signal is more attenuated as we remove each myelin sheath. It not only shows consistency between our results and validates our hypothesis but also supports findings on failing pre-synaptic AP due to demyelinating diseases [55]. We decided to sweep through the entire range of the number of myelin sheaths instead of emphasizing the cytokine storm time-dependency because it made more sense as we wanted to identify every potential demyelination effect to the axonal pathway.

Finally, figure 5 shows the results for the signal coherence. In figure 5a , we show the coherence with regard to our worst scenario, nmy=1 , in which there clearly are more oscillations in lower frequency bands when compared to figure 5b . Figure 5b shows the coherence for a light demyelination, nmy=12 , where only one sheath has been removed with regard to a healthy scenario, nmy=13 . Figure 5c describes the mean and standard deviation of how the coherence values ( 1≤nmy≤12 ) fluctuate in the 0–50 kHz spectrum. Neuronal coherence has been known to serve neuronal communication as an indicator of the efficiency of the exchange of information [58]. In this work, it is noticeable how coherence measurements show higher instability for severely demyelinated neurons in comparison with light demyelination for low- and mid-frequency ranges. As demyelination worsens, so does the reliability of the information going down the axonal channel. All coherence plots showed some fluctuations for the higher end of the frequency range. We believe the fluctuations in the coherence plots is a finite window effect. In other words, the fast Fourier transform (FFT) implicitly filters the data with a rectangular time-domain filter, which is a sinc-shaped filter in the frequency domain. For this reason, we are bound to get those lobes.

Figure 5. Signal coherence between nmy=13 and (a) nmy=1 ; (b) nmy=12 ; and (c) mean and standard deviation of the signal coherence with 1≤nmy≤12 . (a) Signal coherence with n my = 1 (b) Signal coherence with n my = 12. (c) Mean and standard deviation of the signal coherence for all n my.

Signal coherence.

4.2. The linear model identification and verification

Let us observe the changes that the output of a neuron goes through when n varies, and understand these dynamics as depicted in figure 3. Reduction in the number of myelin sheaths ( n , 1≤n≤13 ) causes an increasing delay in the signal, such that spikes start later in neurons with less myelin, and they take longer to reach the peak value (in this section, we also refer to nmy as n for improved consistency with §3.4. On the other hand, in terms of the shape of the spikes, we observe the effect on the spike height and the spike width (full width at half maximum) also in figure 3. This solution was adopted by identifying a suitable TF.

In figure 6a , b , we verify that time intervals represented here correspond to Δt=tn−t13 for 1≤n≤12 , i.e. they are the ‘lag’ observed between 13-sheath neuron, which will be the input to our model and other analysed scenarios that the model needs to approximate well, given the value of n<13 .

Figure 6. Exponential relationships between signals from different demyelinated neurons compared to the healthy n=13 case. (a) Time offset of spike peaks. (b) Time offset of spike onsets.

Exponential relationships between signals from different demyelinated neurons.

At this point, we can say that: (i) for 1≤n≤10 , we see an exponential decay in the ‘lag’ as n grows; and (ii) spike onsets reach the values observed in n=13 one by one ( 1st spike for n=10 , 2nd for n=11 , 3 ⁣rd for n=12 ) while all the peaks reach the n=13 time values in the same case of n=12 . The first conclusion suggests that we will have a transport delay term e−τns in the TF, in which the delay τn will be an exponential function of n . The second conclusion suggests that this term cannot explain all of the dynamics: some of the lag is contributed by a real pole −1/Tn , i.e. a term (1+Tns)−1 in the TF. Furthermore, after n=10 , the delay term vanishes and the only effect seen is the one of the pole (we will ignore this effect, as our approximation focus will be for the interval up to n=10 ). Again, given the linearity of the log plot, i.e. exponential nature of the curve, it is expected that Tn is an exponential function of n .

In figure7a , b , we observe the behaviour of logarithms of spike amplitudes and pulse widths in the region of interest, which suggests (i) time-invariance of the system, as all three pulses collapse in the same amplitude curve, and (ii) that the change in the pulse width requires the real pole −1/Tn . This reasoning, graphically presented in figure 3, confirms our hypothesis about the applicability of the FOPTD TF (3.10) and its exponential coefficients from equation (3.11).

Figure 7. Exponential relationships between signals from different demyelinated neurons compared to the healthy n=13 case. (a) Ratio of spike amplitudes. (b) Ratio of pulse widths.

Exponential relationships between signals from different demyelinated neurons.

The identified parameters of the model equation 3.11 are a0=0.35 , ar=0.7 , τ0=54.42 , τr=0.66 , T0=20.27 and Tr=0.8 . Those values were found with the Levenberg–Marquardt [59] numerical optimization. The iterative procedure was conducted by determining the values of kn,τn,Tn for n=6 , then using those values as initial guesses for n=5 and n=7 , and subsequent values. The relationship between exponential approximation, or linearization in the log-domain, and the best choice of coefficients without exponential law assumption is shown in figure 8a .

Figure 8. Estimation and performance for a linear model of demyelination. (a) Estimated parameters of the TF, and their exponential approximation. (b) Ratio of the RMSE for approximating a demyelinated with the single sheath case, average demyelinated neuron, or a healthy neuron and RMSE for approximating the demyelinated neurons with our TF.

Estimation and performance for a linear model of demyelination.

It is expected that this would be a good approximation for the observed signals in the ‘exponential domain’, 1≤n≤10 . While the approximation can be accurate outside of this domain as well, we focus on applicability within the range, and we verified it using the RMSE metric introduced earlier in equation (3.12). For larger values of n , our numerical estimation was found τ to be zero, hence it could not be shown in a log plot.

Figure 8b , representing Mn for 1≤n≤12 , gives the answer to the following question: if one ignores the variability in n and replaces every output signal (xn) with (i) one of a completely deteriorated neuron N=1 , (ii) averagely damaged neuron N=6 or (iii) a healthy neuron N=13 , how high is the amplitude of error, compared to that of our model. As expected, for n=N this ratio goes down to −∞ dB as the approximation with exact signals is perfect. However, for any other value, even N±1 , our approximation is superior (i.e. above 0 dB). It is important to emphasize that, as observed in figures 4–8, this demyelinating behaviour can be caused by other sources of degeneration. However, there are no claims those results are caused only by virus infections, rather we can say that it is clear from the literature that there is a specific demyelination caused by viruses, and we proposed a phenomenological model that takes into account the potential effects indirectly caused by a viral infection into the nervous system.

5. Conclusion

In this work, we propose an in silico model to quantify the impact of infection on APs. This model describes the dynamics of a cytokine storm and its relation to the extent of demyelination in a neuron. From evidence found in the literature, our phenomenological model is capable of mimicking cytokine-storm-induced demyelination’s degenerative effects. We also proposed a TF aiming towards a linear model of the demyelination process. Using traditional control and systems theory, we propose a FOPTD function that could help the design of chemical control loops for the reinforcement of myelin. It is important to emphasize that, even though we believe the linear approach described in this work can be a stepping stone for more complex systems, the whole system itself is highly nonlinear. Thus, this should be taken into account as it can affect the accuracy of future synthetically engineered therapeutics.

Although the discovery and design of new drugs requires multidisciplinary teams working together, it is important to not only understand what computational tools can offer to improve the accuracy of proposed approaches but also to reproduce cells and molecules interactions as an area of computer-aided drug development. We believe that this model can help biotechnologists as well as pharmacologists to design drugs that will be able to minimize and, hopefully, neutralize the impact of cytokine storms on myelin sheaths. This could be achieved by adding new modules to the model to facilitate the remyelination of axons, and even though the process of remyelination rarely regenerates the myelin sheaths back to their original state [60], this approach could help determine the treatment strategy needed to remyelinate the neuron, aiming to restore it as closely as possible to a healthy, fully myelinated state.

Our results show that demyelination induced by a cytokine storm not only degrades the signal but also impairs its propagation within the axon. This whole analysis led to the development of a TF that fundamentally represents the process of demyelination itself. It not only decreases the level of complexity of the system linking itself with the behaviour of an RC circuit, but also underlies the physical rationale of the system by applying biophysically plausible TF models. We believe that the proposed models will contribute to bioengineering approaches for neurodegeneration, especially demyelinating disease.

For future work, we plan to validate our computational modelling with wet-lab experiments to assess and improve the model’s accuracy. We are interested in analysing and modelling additional experiments, such as whether cytokine diffusion behaves differently in various parts of the brain due to differing diffusion coefficients. This could lead to models that incorporate multimodal data from imaging, genetic or molecular analysis, and enhanced longitudinal data. Also, different viral infections can trigger cytokine storms with different concentration levels. We expect that our proposed technique could lead to more sophisticated and precise approaches for treating neurodegeneration through computational assessment of brain damage from infection.

Acknowledgement

Figures 1 and 2 were created with BioRender.com.

Ethics

This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility

This article has no additional data.

Declaration of AI use

We have not used AI-assisted technologies in creating this article.

Authors’ contributions

G.L.A.: formal analysis, methodology, software, validation, visualization, writing—original draft, writing—review and editing; H.S.: formal analysis, methodology, resources, software, validation, writing—original draft, writing—review and editing; S.B.: conceptualization, funding acquisition, investigation, project administration, supervision, writing—original draft, writing—review and editing; M.T.B.: conceptualization, investigation, methodology, supervision, writing—original draft, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

This publication has emanated from research conducted with the financial support of National Science Foundation (NSF) under grant number 2316960.
==== Refs
References

1. Zubair AS , McAlpine LS , Gardin T , Farhadian S , Kuruvilla DE , Spudich S . 2020 Neuropathogenesis and neurologic manifestations of the coronaviruses in the age of coronavirus disease 2019: A review. JAMA Neurol. 77 , 1018–1027. (10.1001/jamaneurol.2020.2065)32469387
2. Sanclemente-Alaman I , Moreno-Jiménez L , Benito-Martín MS , Canales-Aguirre A , Matías-Guiu JA , Matías-Guiu J , Gómez-Pinedo U . 2020 Experimental models for the study of central nervous system infection by SARS-cov-2. Front. Immunol. 11 , 2163. (10.3389/fimmu.2020.02163)32983181
3. Montalvan V , Lee J , Bueso T , De Toledo J , Rivas K . 2020 Neurological manifestations of COVID-19 and other coronavirus infections: a systematic review. Clin. Neurol. Neurosurg. 194 , 105921. (10.1016/j.clineuro.2020.105921)32422545
4. Song E et al . 2021 Neuroinvasion of SARS-cov-2 in human and mouse brain. J. Exp. Med. 218 , e20202135. (10.1084/jem.20202135)33433624
5. Zanin L , Saraceno G , Panciani PP , Renisi G , Signorini L , Migliorati K , Fontanella MM . 2020 SARS-Cov-2 can induce brain and spine demyelinating lesions. Acta Neurochir. 162 , 1491–1494. (10.1007/s00701-020-04374-x)32367205
6. Wu Y , Xu X , Chen Z , Duan J , Hashimoto K , Yang L , Liu C , Yang C . 2020 Nervous system involvement after infection with COVID-19 and other coronaviruses. Brain Behav. Immun. 87 , 18–22. (10.1016/j.bbi.2020.03.031)32240762
7. Oh Y et al . 2017 Zika virus directly infects peripheral neurons and induces cell death. Nat. Neurosci. 20 , 1209–1212. (10.1038/nn.4612)28758997
8. Valcour V et al . 2012 Central nervous system viral invasion and inflammation during acute HIV infection. J. Infect. Dis. 206 , 275–282. (10.1093/infdis/jis326)22551810
9. Cheng Y , Skinner DD , Lane TE . 2018 Innate immune responses and viral-induced neurologic disease. J. Clin. Med. 8 , 3. (10.3390/jcm8010003)30577473
10. De Chiara G , Marcocci ME , Sgarbanti R , Civitelli L , Ripoli C , Piacentini R , Garaci E , Grassi C , Palamara AT . 2012 Infectious agents and neurodegeneration. Mol. Neurobiol. 46 , 614–638. (10.1007/s12035-012-8320-7)22899188
11. Dal Canto MC , Rabinowitz SG . 1982 Experimental models of virus-induced demyelination of the central nervous system. Ann. Neurol. 11 , 109–127. (10.1002/ana.410110202)6280582
12. Stohlman SA , Hinton DR . 2001 Viral induced demyelination. Brain Pathol. 11 , 92–106. (10.1111/j.1750-3639.2001.tb00384.x)11145206
13. Perlman S , Dandekar AA . 2005 Immunopathogenesis of coronavirus infections: implications for SARS. Nat. Rev. Immunol. 5 , 917–927. (10.1038/nri1732)16322745
14. Geldenhuys M et al . 2018 A metagenomic viral discovery approach identifies potential zoonotic and novel mammalian viruses in Neoromicia bats within South Africa. PLoS One 13 , e0194527. (10.1371/journal.pone.0194527)29579103
15. Chu H et al . 2020 Comparative replication and immune activation profiles of SARS-cov-2 and SARS-cov in human lungs: an ex vivo study with implications for the pathogenesis of COVID-19. Clin. Infect. Dis. 71 , 1400–1409. (10.1093/cid/ciaa410)32270184
16. Merson TD , Binder MD , Kilpatrick TJ . 2010 Role of cytokines as mediators and regulators of microglial activity in inflammatory demyelination of the CNS. Neuromolecular Med. 12 , 99–132. (10.1007/s12017-010-8112-z)20411441
17. Ludwig J , Binder A , Steinmann J , Wasner G , Baron R . 2008 Cytokine expression in serum and cerebrospinal fluid in non-inflammatory polyneuropathies. J. Neurol. Neurosurg. Psychiatry. 79 , 1268–1274. (10.1136/jnnp.2007.134528)18550631
18. Chen G et al . 2020 Clinical and immunological features of severe and moderate coronavirus disease 2019. J. Clin. Invest. 130 , 2620–2629. (10.1172/JCI137244)32217835
19. Song P , Li W , Xie J , Hou Y , You C . 2020 Cytokine storm induced by SARS-cov-2. Clin. Chim. Acta. 509 , 280–287. (10.1016/j.cca.2020.06.017)32531256
20. Zoghi A , Ramezani M , Roozbeh M , Darazam IA , Sahraian MA . 2020 A case of possible atypical demyelinating event of the central nervous system following COVID-19. Mult. Scler. Relat. Disord. 44 , 102324. (10.1016/j.msard.2020.102324)32615528
21. Waito M , Walsh SR , Rasiuk A , Bridle BW , Willms AR . 2016 A mathematical model of cytokine dynamics during a cytokine storm. In Mathematical and computational approaches in advancing modern science and engineering (ed. J Bélair ), pp. 331–339. Cham, Switzerland: Springer International Publishing. (doi:10.1007\%2F978-3-319-30379-6_31)
22. Yiu HH , Graham AL , Stengel RF . 2012 Dynamics of a cytokine storm. PLoS One 7 , e45027. (10.1371/journal.pone.0045027)23049677
23. Bitsch A , Kuhlmann T , Da Costa C , Bunkowski S , Polak T , Brück W . 2000 Tumour necrosis factor alpha mrna expression in early multiple sclerosis lesions: correlation with demyelinating activity and oligodendrocyte pathology. Glia 29 , 366–375. (10.1002/(sici)1098-1136(20000215)29:4<366::aid-glia7>3.0.co;2-y)10652446
24. Redford EJ , Hall SM , Smith KJ . 1995 Vascular changes and demyelination induced by the intraneural injection of tumour necrosis factor. Brain 118 , 869–878. (10.1093/brain/118.4.869)7655885
25. Barros MT , Silva W , Regis CDM . 2018 The multi-scale impact of the alzheimer’s disease on the topology diversity of astrocytes molecular communications nanonetworks. IEEE Access 6 , 78904–78917. (10.1109/ACCESS.2018.2885518)
26. Adonias GL , Yastrebova A , Barros MT , Koucheryavy Y , Cleary F , Balasubramaniam S . 2020 Utilizing neurons for digital logic circuits: a molecular communications analysis. IEEE Trans. Nanobioscience 19 , 224–236. (10.1109/TNB.2020.2975942)32092011
27. Adonias GL , Yastrebova A , Barros MT , Balasubramaniam S , Koucheryavy Y . A logic gate model based on neuronal molecular communication engineering. In Proceedings of the 4th Workshop on Molecular Communications, Linz, Austria, pp. 15–16. https://molecularcommunications.org/wp-content/uploads/2022/01/Technical24.pdf.
28. Barros MT , Doan P , Kandhavelu M , Jennings B , Balasubramaniam S . 2021 Engineering calcium signaling of astrocytes for neural-molecular computing logic gates. Sci. Rep. 11 , 595. (10.1038/s41598-020-79891-x)33436729
29. Akyildiz IF , Pierobon M , Balasubramaniam S . 2019 An information theoretic framework to analyze molecular communication systems based on statistical mechanics. Proc. IEEE 107 , 1230–1255. (10.1109/JPROC.2019.2927926)
30. Adonias GL , Duffy C , Barros MT , McCoy CE , Balasubramaniam S . 2021 Analysis of the information capacity of neuronal molecular communications under demyelination and remyelination. IEEE Trans. Neural Syst. Rehabil. Eng. 29 , 2765–2774. (10.1109/TNSRE.2021.3137350)34932481
31. Kuscu M , Dinc E , Bilgin BA , Ramezani H , Akan OB . 2019 Transmitter and receiver architectures for molecular communications: a survey on physical design with modulation, coding, and detection techniques. Proc. IEEE 107 , 1302–1341. (10.1109/JPROC.2019.2916081)
32. Hou P , Eckford AW , Zhao L . 2020 Analysis and design of two-hop diffusion-based molecular communication with ligand receptors. IEEE Access 8 , 189458–189470. (10.1109/ACCESS.2020.3032009)
33. Egan M , Kuscu M , Barros MT , Booth M , Llopis-Lorente A , Magarini M , Martins DP , Schäfer M , Stano P . 2023 Toward interdisciplinary synergies in molecular communications: perspectives from synthetic biology, nanotechnology, communications engineering and philosophy of science. Life 13 , 208. (10.3390/life13010208)36676156
34. Soldner CA et al . 2020 A survey of biological building blocks for synthetic molecular communication systems. IEEE Commun. Surv. Tutorials 22 , 2765–2800. (10.1109/COMST.2020.3008819)
35. Aghababaiyan K , Shah-Mansouri V , Maham B . 2018 Axonal channel capacity in neuro-spike communication. IEEE Trans. Nanobioscience 17 , 78–87. (10.1109/TNB.2018.2800899)29570078
36. Barros MT , Veletić M , Kanada M , Pierobon M , Vainio S , Balasingham I , Balasubramaniam S . 2021 Molecular communications in viral infections research: modeling, experimental data, and future directions. IEEE Trans. Mol. Biol. Multiscale Commun. 7 , 121–141. (10.1109/TMBMC.2021.3071780)35782714
37. Chahibi Y . 2017 Molecular communication for drug delivery systems: a survey. Nano Commun. Netw. 11 , 90–102. (10.1016/j.nancom.2017.01.003)
38. Cohen CCH , Popovic MA , Klooster J , Weil MT , Möbius W , Nave KA , Kole MHP . 2020 Saltatory conduction along myelinated axons involves a periaxonal nanocircuit. Cell 180 , 311–322.(10.1016/j.cell.2019.11.039)31883793
39. Carnevale NT , Hines ML . 2009 The NEURON book, 1st ed. New York, NY: Cambridge University Press.
40. Hodgkin AL , Huxley AF . 1952 A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117 , 500–544. (10.1113/jphysiol.1952.sp004764)12991237
41. Hartung HP . 1993 Immune-mediated demyelination. Ann. Neurol. 33 , 563–567. (10.1002/ana.410330602)8498835
42. Empl M , Renaud S , Erne B , Fuhr P , Straube A , Schaeren–Wiemers N , Steck AJ . 2001 TNF-alpha expression in painful and nonpainful neuropathies. Neurol. 56 , 1371–1377. (10.1212/WNL.56.10.1371)
43. Bhatheja K , Field J . 2006 Schwann cells: origins and role in axonal maintenance and regeneration. Int. J. Biochem. Cell Biol. 38 , 1995–1999. (10.1016/j.biocel.2006.05.007)16807057
44. Baumann N , Pham-Dinh D . 2001 Biology of oligodendrocyte and myelin in the mammalian central nervous system. Physiol. Rev. 81 , 871–927. (10.1152/physrev.2001.81.2.871)11274346
45. Adonias GL , Siljak H , Barros MT , Marchetti N , White M , Balasubramaniam S . 2020 Reconfigurable filtering of neuro-spike communications using synthetically engineered logic circuits. Front. Comput. Neurosci. 14 , 556628. (10.3389/fncom.2020.556628)33178001
46. Pintelon R , Van Biesen L . 1990 Identification of transfer functions with time delay and its application to cable fault location. IEEE Trans. Instrum. Meas. 39 , 479–484. (10.1109/19.106276)
47. Aghababaiyan K , Shah-Mansouri V , Maham B . 2020 Capacity and error probability analysis of neuro-spike communication exploiting temporal modulation. IEEE Trans. Commun. 68 , 2078–2089. (10.1109/TCOMM.2019.2962805)
48. Sundaresan KR , Krishnaswamy PR . 1978 Estimation of time delay time constant parameters in time, frequency, and laplace domains. Can. J. Chem. Eng. 56 , 257–262. (10.1002/cjce.5450560215)
49. Juchem J , Dekemele K , Chevalier A , Loccufier M , Ionescu CM . 2019 First order plus Frequency dependent delay modeling: new perspective or mathematical curiosity? In 2019 IEEE International Conference on Systems, Man and Cybernetics (SMC), Bari, Italy, pp. 2025–2030. (10.1109/SMC.2019.8914386). https://ieeexplore.ieee.org/xpl/mostRecentIssue.jsp?punumber=8906183.
50. Farmer TG , Edgar TF , Peppas NA . 2009 Effectiveness of intravenous infusion algorithms for glucose control in diabetic patients using different simulation models. Ind. Eng. Chem. Res. 48 , 4402–4414. (10.1021/ie800871t)20161147
51. Padma Sree R , Srinivas MN , Chidambaram M . 2004 A simple method of tuning PID controllers for stable and unstable FOPTD systems. Comput. Chem. Eng. 28 , 2201–2218. (10.1016/j.compchemeng.2004.04.004)
52. Debanne D , Campanac E , Bialowas A , Carlier E , Alcaraz G . 2011 Axon physiology. Physiol. Rev. 91 , 555–602. (10.1152/physrev.00048.2009)21527732
53. Khodaei A , Pierobon M . 2016 An intra-body linear channel model based on neuronal subthreshold stimulation. In ICC 2016–2016 IEEE International Conference on Communications, Kuala Lumpur, Malaysia, pp. 1–7. (10.1109/ICC.2016.7511483)
54. Eliasmith C , Anderson CH . 2003 Neural engineering: computation, representation, and dynamics in neurobiological systems. Cambridge, MA: MIT Press.
55. Hamada MS , Popovic MA , Kole MHP . 2017 Loss of saltation and presynaptic action potential failure in demyelinated axons. Front. Cell. Neurosci. 11 , 45. (10.3389/fncel.2017.00045)28289377
56. Bando Y , Takakusaki K , Ito S , Terayama R , Kashiwayanagi M , Yoshida S . 2008 Differential changes in axonal conduction following CNS demyelination in two mouse models. Eur. J. Neurosci. 28 , 1731–1742. (10.1111/j.1460-9568.2008.06474.x)18973589
57. Coggan JS , Prescott SA , Bartol TM , Sejnowski TJ . 2010 Imbalance of ionic conductances contributes to diverse symptoms of demyelination. Proc. Natl Acad. Sci. USA 107 , 20602–20609. (10.1073/pnas.1013798107)20974975
58. Schoffelen JM , Oostenveld R , Fries P . 2005 Neuronal coherence as a mechanism of effective corticospinal interaction. Science 308 , 111–113. (10.1126/science.1107027)15802603
59. Moré JJ . 1978 The Levenberg–Marquardt algorithm: implementation and theory. In Numerical analysis (ed. GA Watson ), pp. 105–116. Berlin, Germany: Springer. (10.1007/BFb0067700)
60. Franklin RJM , Ffrench-Constant C . 2008 Remyelination in the CNS: from biology to therapy. Nat. Rev. Neurosci. 9 , 839–855. (10.1038/nrn2480)18931697
