{'fraction_non_alphanumeric': 0.03331426937374893, 'fraction_numerical': 0.009837003145553332, 'mean_word_length': 4.786068828590338, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neuron Signal Propagation Analysis of Cytokine-Storm induced Demyelination

The COVID-19 pandemic has shaken the world unprecedentedly, where it has affected the vast global population both socially and economically. The pandemic has also opened our eyes to the many threats that novel virus infections can pose for humanity. While numerous unknowns are being investigated in terms of the distributed damage that the virus can do to the human body, recent studies have also shown that the infection can lead to lifelong sequelae that could affect other parts of the body, and one example is the brain. As part of this work, we investigate how viral infection can affect the brain by modelling and simulating a neuron's behaviour under demyelination that is affected by the cytokine storm. We quantify the effects of cytokine-induced demyelination on the propagation of action potential signals within a neuron. We used information and communication theory analysis on the signal propagated through the axonal pathway under different intensity levels of demyelination to analyse these effects. Our simulations demonstrate that virus-induced degeneration can play a role in the signal power and spiking rate and the probability of releasing neurotransmitters and compromising the propagation and processing of information between the neurons. We also propose a transfer function that models these attenuation effects that degenerates the action potential, where this model has the potential to be used as a framework for the analysis of virus-induced neurodegeneration that can pave the way to improved understanding of virus-induced demyelination.


INTRODUCTION
The recent outbreak of the coronavirus disease 2019  pandemic caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been impacting society as a whole, leaving a long-lasting impact on people's health. As a result, the scientific community have been joining efforts and resources towards not only the epidemiologic characteristics and transmission dynamics of the virus, but also 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 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)(2)(3)(4) and elicit neurodegeneration (5,6), in which 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 trigger immunotherapy-induced neuroinflammation and, consequently, lead to neurodegeneration (9). Viruses are known for causing dramatic structural and biochemical changes to the host cell, by hĳacking 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 characterised by numerous toxic-metabolic derangements, such as cytokine storms (Figure 1(a)). 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 the infectious agent (10).
Other types of coronavirus have been known to cause demyelination, such as the murine coronavirus (M-CoV). It 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 (10). This behaviour also matches findings on SARS-CoV-1, which reports a decrease in viral titers as clinical disease worsens (11). M-CoV is a type coronavirus of the same genus (betacoronavirus) as SARS-CoV-2, and it is believed to be 43.8%-48% similar to this novel coronavirus (12). 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 virus in human tissues (13).
The virus replication process produces host cells exhaustion that leads to the activation of the immune system, which calls for macrophages' work. 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 tumor necrosis factor alpha (TNF-) aiming to fight the infection. These cytokines exert cytotoxic effects on neurons and glial cells, e.g. oligodendrocytes, where this damages the neurons' myelin sheath providing support and insulation to axons (14).
Then the dynamics of cytokine storms and their efficiency in fighting infections without further cellular damage becomes of increase attention to measure the neuroinflammatory effects of COVID-19 infection. Ludwig et al (15) found out that as the neuropathy is more severe, the higher the serum concentration of TNF-and IL-6. This matches more recent results found by Chen et al (16), which analysed plasma 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. 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 (17) in which they match recordings of 13 different cytokines with a model of non-linear ordinary differential equations. Likewise, Yiu et al (18) 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 structures or non-neuronal cells. For instance, the work of Bitsch et al (19) 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 (20) 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, Molecular Communications (MC) has been improving biological models by accounting for the communication of cells using their signalling mechanisms as information carriers. MC 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. A recent survey (21) 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 synthetic biology-based treatment. Even further targeted drug delivery technology can alleviate the cytokine storms effects on neurons and provoke the restoration of the myelin sheath (22). This paper presents a systems theory-based analysis of the demyelination that is, either directly or indirectly, caused by viral infections from the perspective of MC, more specifically, a neuro-spike MC. The goals of this paper are (1) to propose a mathematical model that represents the cytokine storms stemming from acute infection-derived inflammation; (2) to provide insights on how the demyelination will affect the neuronal information along the axonal pathway; and (3) to provide a transfer function (TF) model that describes the effects of the membrane action potential 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 example, myelin sheaths) are traditionally well-approximated by first-order transfer functions with time delay and the underlying physical rationale, our model joins the family of established, applicable biophysical transfer function models. We expect that these analyses could pave the way for more in-depth studies that can support in vitro and in vivo experimental work on the neurological effects induced by neurotropic viral infections.
The contributions of this work are as follows: • A mathematical model that describes the process of virus-induced demyelination: We present a model starting from the well-documented fact that viral infections trigger cytokine storms. These 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 action potential 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 transfer function 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 neurons' intra-signalling.

MATERIALS AND METHODS
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 (18). This model was extended by coupling together with a neuronal model that implements the behaviour of a myelinated axon (23). By extending these models, we developed simulations and the default parameters of the cell (such as the length and diameter of each of their compartments) were 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 action potentials on a Layer 5 (L5) pyramidal cell with its myelin sheath affected. All simulations were performed with extensions to the NEURON Simulator (24).

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, ( ), and its rate of change, Δ ( ), in a vector-matrix form, can be represented as follows 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 (see (18)). Also, and are positive constants that express the sensitivity of the cytokine's acceleration to concentration and rate of change.
The cytokine's response modes are characterised by the eigenvalues, 1 and 2 (rad/day), of the stability matrix of the system, and this is as follows given.
The parameters are chosen to minimize the error between the cytokine concentration and the clinical trial measurements performed in (18). 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.

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 action potentials (AP) or spikes, and avoid significant attenuation due to other synaptic processes. According to Cohen et al (23), the myelin sheath dynamics can be described using circuit theory. This circuit is then coupled with the Hodgkin-Huxley (HH) circuit model, which describes the AP dynamics in neurons ( Figure 2).
Axial resistance, , of the axon core, can be defined as the ratio between axial resistivity, , and the cross-sectional area, and this is represented as where and is the axon core diameter. Let be the periaxonal cross-sectional area, then where the axial resistance in the periaxonal space, , is calculated as the ratio between periaxonal resistivity, , 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 + 2 , thus and the same calculation can be performed for , and . Recognising that a myelin sheath is in-series compaction of layers, the radial resistance of the sheath, , is the sum of the resistance of each myelin membrane, represented as and, the radial capacitance of the myelin sheath, , may vary inversely to the sum of the capacitances of each of its composing membranes, , thus where the resistance and capacitance of a single myelin membrane are and , respectively. Based on this the expression of the function can be represented in terms of the number of myelin lamellae, , as follows For a more detailed analysis of the myelin sheath's modelling, the reader is referred to the work of Cohen et al (23).

Cytokine-induced Demyelination
As mentioned in the Introduction, a cytokine storm is released to fight infections and, in the process, can damage healthy tissues.
There are numerous pieces of evidence linking cytokine storms to neurodegeneration (15,(17)(18)(19). 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 sheaths. Applying linear regression ( = 0.508, = 0.001) on the data provided by Ludwig et al (15), we establish a proportional relationship between the severity of the neuropathy ( ) and the serum concentration of TNF-, ( ), as and in this relationship, the stronger the cytokine storm is, the more severe the degeneration. The severity of the neuropathy was defined according to (15) as a scoring function, where scoring nerve functions of '0' is 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. The severity that initially ranges between 0 and 8 in total should be re-scaled to match a normal as reference.

A Linear Model of Demyelination
Based on a given healthy neuron, we want to predict the signal effects on spike trains that are expected from a demyelinated neuron. Let , the number of the myelin sheaths, be a tunable parameter of the model (in this section, we refer to as for clarity of index notation). This mechanism is depicted in the block diagram in Figure 1(b).
Observing faulty signals resulting from the propagation of standard signals through a transfer function superimposed on the original system is a well-established concept; examples of its application include modelling of cables (25). It makes intuitive sense to observe such a transfer function as taking the health system's output and delivering the faulty, unhealthy signal as the faulty output, which results in signal degradation for a worsening channel it traverses. The opposite process, in which an "adaptor" transfer function 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. Now that we have established the processing chain, the choice for the transfer function is made by observing the general trends in the output signals for various values of . Namely, as seen in Figure 3 the transfer function emulating the effect of myelin deficiency needs to allow for the attenuation of the signal, widening of the spikes, and overall propagation delay. An obvious candidate is the traditional First Order Plus Time Delay (FOPTD) function, given by where we assume that the parameters , and have different values for different values of . While those values could be tabulated and looked up for specific values of , 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 = 0 · , = · and = 0 · ), 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 sheathes. For example, if the healthy baseline signal is achieved with = 13, emulation of = 10 can be interpreted as undoing the effect of Δ = 3 sheathes, i.e. the effect of a cascade of three blocks representing a "sheath undo function". This approximation is good for larger Δ , as it allows for the approximations like (1 +¯) Δ ≈ 1 + Δ , i.e. giving rise to a ∼ exp( ) law. In the Results and Discussion section, 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 transfer function gives us the final form of our transfer function, which is represented as The final form of the equation 11 suggests the rationale behind suggesting exponential behaviour of = log instead of itself: in the -domain, we obtain a function of the form exp( + ) + , which in turn in Fourier domain corresponds to exp( + ) + .
Knowing that the behaviour of the myelin circuit originates from connections of (purely passive, real impedance) and (purely active, imaginary impedance), it is expected to observe this symmetric real-imaginary coupling of terms.
FOPTD is not an uncommon choice in biophysics, where it has been used to model glucose control (26), because they are simple for identification (27) and for quick and accurate tuning of the controllers that can regulate their behaviour (28). 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): Here, RMSE , stands for the RMSE of the output of our transfer function compared to the actual output signal for sheaths, i.e. if the is the number of samples, which is represented as Analogously, we can find the value for RMSE , which corresponds to the RMSE of the output signal produced by another model compared to the th actual output signal. The quantity is positive where our model is more accurate (has lower RMSE) than the model we are comparing to it.

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 , where = {1, 2, 3, ..., } identifies the order of each spike and, as the peak times at the output. Thus, we define the relative mean time shift, , as analogously, we can define the relative mean amplitude shift, , with and as the peak amplitudes of spike at the input and output, 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 . 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 in the soma of the neuron. Later, analogous to the relationship between and the latency, we decided to look into a potential attenuation on the power of the signal, , as we have already discussed a relatively heavy mean attenuation in the membrane potential, . Let us define as where in a set of samples, [ ] corresponds to the potential of the membrane at the -th sample.
As there are changes in specific time points where spikes peak, it is most likely that there will be changes in the spiking rate. For that reason, we calculated the approximated spiking rate from the mean inter-spike interval (ISI). Consider a spike train with spikes and, between each spike there is an interval of value where = {1, 2, 3, ..., − 1}. We can define an equation for calculating the spiking rate ( ) from the ISI, thus which, according to Yu et al (29), there is a linear relationship ( 2 = 0.9304) between the spiking rate, , and the probability of releasing glutamate into the synaptic cleft, ( ). Glutamate is an excitatory neurotransmitter that is present in most synaptic connections. Its deficit can dysregulate the depolarisation of the postsynaptic neuron and, as a consequence, compromise the transfer of information in a neuronal network. This relationship is expressed as We also quantified the attenuation of the signal, 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, where and 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.      Lastly, we also analysed the relation between the reference spike train ( = 13) and all other demyelinating scenarios. The intention is to understand how much power is being transferred between each pair of signals. With that in mind, we applied a coherence ( ) metric, which is described as where ( ) is the cross-spectral density between the two signals and, ( ) and ( ) are the power spectrum densities of input and output, respectively.

RESULTS AND DISCUSSION
Let us consider the proportionality between the severity of the neuropathy and the number of myelin lamellae, , as described in the Cytokine-induced Demyelination section. Our analysis consisted of evaluating the neuronal behaviour and spike propagation under normal circumstances ( = 13), followed by an analysis on the demyelination by decreasing . The objective is to understand what happens to the neuronal information when travelling through a demyelinated axonal pathway. The neuron receives an external current of 3 nA for 15 ms at time = 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 (Figure 1(b)).

Analysis of a Demyelination-induced Channel Attenuation
The results for and from Equation (14) are shown in Figure 4(a). It shows that, although time has quite slight delays as we decrease , the amplitude of the membrane potential is a bit more drastically attenuated. In other words, is almost eight times more degraded than in the worst scenario ( = 1). 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 (30,31). Each curve's shape resembles a capacitor charging (amplitude) or discharging  (time) in RC circuits. This comes as no surprise since both the axonal compartment and the myelin sheath are roughly expressed as RC circuits ( Figure 2). Furthermore, as latency is an operation on time of the peak, it is expected the curve behaves similarly as it does for , as shown in Figure 4(b). This time, we notice that there is a subtle "lag" for the signal to travel across the axon even under normal circumstances matching findings in the literature (32) 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. On the other hand, Figure 4(c) 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 action potentials. 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 (31). 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. Figure 4(d) depicts how the spiking rate is affected by the demyelination. As we expected, as the spikes start to get wider and far from each other, the spiking rate gets lower. This corroborates findings on demyelination-induced effects on spiking rate (33). The rate at which a neuron fires action potentials is significant for modulating and encoding neuronal information in cognitive, sensory and motor functions.
Not surprisingly, when we calculate the probability for releasing glutamate in Figure 4(e), it shows that the probability curve has the same shape as the spiking rate when increasing . This is due to the linear relationship expressed in Equation (17). There is a subtle decrease in ( ), which may be aggravated depending on how a specific neuron encodes the information. This matches findings in the literature that state axonal coding is very sensitive to subtle changes in the pattern of the stimuli (34).
Furthermore, we decided to investigate the attenuation caused by the demyelination from a more generic communication systems point of view as expressed in Equation 18. The results presented in Figure 4(f) 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 (31).
Lastly, Figure 5 shows the results for the signal coherence. In Figure 5(a), we show the coherence with regards to our worst scenario, = 1, in which there clearly are more oscillations in lower frequency bands when compared to Figure 5(b). Figure 5(b) shows the coherence for a light demyelination, = 12, where only 1 sheath has been removed with regards to a healthy scenario, = 13. Figure 5(c), describes the mean and standard deviation of how the coherence values (1 ≤ ≤ 12) fluctuate in the 0 − 50kHz spectrum. Neuronal coherence has been known to serve neuronal communication as an indicator of the efficiency of the exchange of information (35). 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.

The Linear Model Identification and Verification
Let us observe the changes that the output of a neuron goes through when varies, and understand these dynamics as depicted in Figure 3. Reduction in the number of myelin sheaths ( , 1 ≤ ≤ 13) causes an increasing delay in the signal, i.e., spikes start    later in neurons with less myelin, and they take longer to reach the peak value. 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, FWHM) also in Figure 3. This solution was adopted by identifying a suitable transfer function. In Figures 6(a) and 6(b) we verify that time intervals represented here correspond to Δ = − 13 for 1 ≤ ≤ 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 < 13).
At this point, we can say that (1) for 1 ≤ ≤ 10, we see an exponential decay in the "lag" as grows; (2) spike onsets reach the values observed in = 13 one by one (1 st spike for = 10, 2 nd for = 11, 3 rd for = 12) while all the peaks reach the = 13 time values in the same case of = 12. The first conclusion suggests that we will have a transport delay term − in the transfer function, in which the delay will be an exponential function of . The second conclusion suggests that this term cannot explain all of the dynamics: some of the lag is contributed by a real pole −1/ , i.e. a term (1 + ) −1 in the transfer function. Furthermore, after = 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 = 10). Again, given the linearity of the log plot, i.e. exponential nature of the curve, it is expected that is an exponential function of .
In Figures 6(c) and 6(d) we observe the behaviour of logarithms of spike amplitudes and pulse widths in the region of interest, which suggests (1) time-invariance of the system, as all three pulses collapse in the same amplitude curve, and (2) that the change in the pulse width requires the real pole −1/ . This reasoning, graphically presented in Figure 3, confirms our hypothesis about the applicability of the FOPTD transfer function (10) and its exponential coefficients form (11).
The identified parameters of the model 11 are 0 = 0.35, = 0.7, 0 = 54.42, = 0.66, 0 = 20.27 and = 0.8. Those values were found with the Levenberg-Marquardt (36) numerical optimisation. The iterative procedure was conducted by determining the values of , , for = 6, then using those values as initial guesses for = 5 and = 7, and subsequent values. The relationship between exponential approximation, or linearisation in log domain, and the best choice of coefficients without exponential law assumption is shown in Figure 7(a).
It is expected that this would be a good approximation for the observed signals in the "exponential domain", 1 ≤ ≤ 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 (12) introduced earlier. Figure 7(b), representing for 1 ≤ ≤ 12, gives the answer to the following question: if one ignores the variability in and replaces every output signal ( ) with (a) one of a completely deteriorated neuron = 1, (b) averagely damaged neuron = 6, or (c) a healthy neuron = 13, how high is the amplitude of error, compared to that of our model. As expected, for = this ratio goes down to −∞ dB as the approximation with exact signals is perfect. However, for any other value, even ± 1, our approximation is superior (i.e. above 0 dB).

CONCLUSION
In this work, we proposed an end-to-end model that takes advantage of the fact that viruses can invade the nervous system and affect the brain. This model expresses the dynamics of a cytokine storm and its relation to the amount of demyelination that may subject to a neuron. From evidence found in the literature, we developed a model capable of mimicking virus-induced demyelination's degenerative effects. We also proposed a transfer function aiming towards a linear model of the demyelination. 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. The results show that the demyelination induced by the cytokine storm not only degrades the signal but also compromises the communications among neurons. The signal is attenuated and shifted in time and influences the release of excitatory neurotransmitters into the synaptic cleft. This whole analysis led to the development of a transfer function 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 transfer function models. We believe that the proposed models will contribute to bioengineering approaches for neurodegeneration, especially demyelinating disease. For future work, we plan to couple our computational modelling with wet-lab experiments to assess and improve our model's accuracy. Such a technique could pave the way for more sophisticated and precise approaches for the treatment of neurodegeneration.

AUTHOR CONTRIBUTIONS
GA performed the simulations, data analysis and wrote the first draft of the manuscript. HS performed the linear model analysis. MB and SB led work development. All authors contributed to manuscript writing and revision. All authors have also read and approved the submitted version.