On the effect of remyelination in a mathematical model of multiple sclerosis

: In the present work, we consider a mathematical model of multiple sclerosis, extending a model, known in the literature, so that it can account for the process of remyelination. Our model comprises of a reaction–diffusion–chemotaxis system of partial differential equations with a time delay. As a ﬁrst approximation, we consider the model under the assumption of radial symmetry, which is consistent, e.g., with Bal´o’s concentric disease. We conduct numerical experiments in order to study the effect of the remyelination strength on the disease progression. Furthermore, we show that the modiﬁed model has greatly enriched dynamics, which is capable of describing qualitatively different types of multiple sclerosis (according to classical classiﬁcations of the disease progression) as well as giving a better agreement with experimental data.


I. INTRODUCTION
The central nervous system (CNS) consists of the brain and spinal cord, which are responsible for processing and coordinating the flow of information through an organism's body.Neurons are the primary cells of the CNS, and they consist of three main parts: the cell body, dendrites, and axons.Myelin is a fatty substance that surrounds and insulates axons, allowing for faster transmission of signals as well as providing structural support for the axons.A schematic representation of a myelinated nerve cell is given in Fig. 1.
Multiple sclerosis (MS) is a disease of the CNS that disrupts the flow of information within the brain, and between the brain and body.It involves an immunemediated process in which an abnormal response of the body's immune system is directed against the CNS [1].In individuals with MS, the immune system causes inflammation that damages myelin as well as the nerve fibers themselves, and the specialized cells (oligodendrocytes) that produce myelin.When myelin or nerve fibers are damaged or destroyed in MS, messages within the CNS are altered or stopped completely.The damaged areas develop scar tissue which gives the disease its name-multiple areas of scarring or multiple sclerosis.Some of the most common symptoms in patients experiencing multiple sclerosis include, but are not limited to, vision impairment, mobility problems, numbness and tingling, muscle spasms, stiffness and weakness, bladder and bowel problems, speech and swallowing difficulties [3].
The cause of MS is not known, but it is believed to involve genetic susceptibility, abnormalities in the immune system and environmental factors that combine to trigger the disease [4].Diagnosing MS can also become complicated due to the variations in the temporal development of the disease-some patients Fig. 1: Scheme of a myelinated neuron [2].
experience immediate acute effects, while others develop the symptoms insidiously over extended periods of time.Furthermore, while there exists no known cure for MS, disease-modifying therapies, such as immune activity modulation or suppression have been implemented in order to alleviate the symptoms and treat acute attacks [5].However, one should note that such treatments remain in their experimental stage and could in some situation prove to hide certain significant risks [6].
Mathematical modelling could, thus, be very useful for shedding new light on various aspects of the disease, simulating different hypotheses, etc.Unfortunately, besides purely empirical/statistical studies, it is very hard to find many papers, dealing with the topic of mathematical modelling of MS.The works [7][8][9] are among the few, which try to give a (very basic) mathematical description of the underlying mechanisms of MS progression.Those works consider models, described in terms of ODE systems, of the time evolution of most broadly defined cell populations, thought to be important for the disease progression.The other possible view concerns the spatio-temporal dynamics of the lesions generation.To the best of the authors' knowledge, the main works in this direction comprise of the work of Lombardo et al. [10], which steps on the works [11,12].
In the present work, we are interested in the description of the spatio-temporal dynamics regarding most key cells involved in MS pathogenesisoligodendrocytes, macrophages and cytokines.We, thus, step on the model, proposed in [10], and modify it, in order to enrich its dynamics (by accounting for a very important aspect, missing in the original modelremyelination) such that more realistic scenarios for a disease course can be simulated in a broader range of MS types.
The original model from [10] incorporates the pro-cess of demyelination, a crucial component of the disease progression in all MS cases.Based on the types of demyelination and the targets affected by it, four different patterns of lesions have been proposed [13].In this work the accent falls on modelling the progression of the third type of lesions-cases in which the oligodendrocytes are the primary "victims", whereas in type II cases, for example, the main target of the inflammatory response is the myelin.The third type patterns also can be regarded as being formed in the initial stages of disease cases (or to the "pre-demyelinatig" phase) in the scenario of intraindividual homogeneity.
Let us first summarize the main agents, known to be connected with the disease, which are of interest for the present work: • Oligodendrocytes-myelin producing cells which are destroyed by the macrophages; • Macrophages-a type of white blood cells responsible for the destruction of the oligodendrocytes in MS cases; • Cytokines-signaling proteins, the presence of which attracts macrophages.They are produced by both macrophages and destroyed oligodendrocytes.
For type III lesions, macrophages and activated microglia are known to trigger the inflammatory process by producing pro-inflammatory cytokines.Thus, the model in [10] considers the time evolution of those species along with the destroyed oligodendrocytes, since significant oligodendrocyte apoptosis has been observed in this type of lesions.
As noted, we step on the aforementioned model and modify it by adding a remyelination term.It is very important for a variety of cases in some of the main types of MS, as we further explain.As a first approximation, in order to study the dynamics of the extended model, we consider the assumption of radial symmetry.In particular, we explore the possibilities of applying the model to a particularly aggressive demyelinating disease called Baló's concentric disease (BCD), which is often considered to be a variant of multiple sclerosis [14].Baló's disease is most notably characterized by concentric alternating rings of CNS lesions (demyleinated areas) and myelinated white matter [15].This makes it convenient to simplify the model by assuming radial symmetry, since demyelination levels appear nearly identical at any fixed distance from a given center.This process of ring formation has been studied to occur in a wave-like fashion starting outward from a focal point [14].Conclusive information regarding the pathogenesis of such a disease progression remains elusive, however the rings of relatively unharmed white matter have sometimes been hypothesized to be the result of remyelinating processes [16].The present work is structured as follows.Section II deals with the mathematical formulation of the considered model.In particular, we first formulate the original model from [10], which comprises of a reactiondiffusion-chemotaxis system of 3 PDEs and summarize its basic underlying assumptions.Then, we go on to extend the model by including a remyelination term with a time-delay and formulate the newly obtained model in terms of a partial delay differential equation system in operator form.Finally, a 1D coordinate form of the model is derived in polar coordinates under the assumption of radial symmetry.In Section III, we construct an explicit finite difference approximation of the radially-symmetric model, which is further used in section IV to study the effect of the modification we make to the model.To that end, the model solutions are studied, while varying the rate of remyelination in subsection IV-A and the time delay in subsection IV-B.Subsection IV-C emphasizes the enriched dynamics of the model in terms of more types of MS, which can be accounted for and subsection IV-D shows the better descriptive capabilities in terms of a comparison with an MRI scan.Section V summarizes the main findings and outlines some directions for further development of the subject.

II. MATHEMATICAL FORMULATION
The original model, proposed in [10], consists of three partial differential equations, describing the spatio-temporal dynamics of the following three species: • m(X,t)-density of activated macrophages; • c(X,t)-concentration of chemoattractants (cytokines); • d(X,t)-density of the destroyed oligodendrocytes, where (X,t)∈ Ω × R + , Ω is a bounded region of R n , n = 1, 2. X and t signify the respective position in space and moment in time.
The basic assumptions that underlie the model are as follows: • For the evolution of the macrophages: their flux is caused by random movement as well as directed movement (called chemotaxis) in the direction of cytokines; the natural growth of the population of macrophages is logistic; a Monod-type model [17] that models the prevention of overcrowding is considered.
• For the cytokines: their flux is caused only by random motion; their production depends on both destroyed oligodendrocytes and macrophages; their natural decay is proportional to their concentration; the evolution of cytokines may be considered to be on a different time scale than the other two, i.e., compared to the other two considered populations, cytokine concentration evolves more quickly.• For the destroyed oligodendrocytes: they are static; their destruction is caused by coming in contact with macrophages; the damage that macrophages cause is modelled by an increasing function with saturation for high values of the macrophages' density.After a proper nondimensionalization, as done in [10], the model takes the form where and χ ∞ is the chemoattraction coefficient, τ is a timescaling parameter for the cytokine-related processes, is the cytokine diffusion coefficient, β is the rate of cytokine production by macrophages, δ is the rate of cytokine production by oligodendrocytes, and r is the intensity, with which the macrophages damage the oligodendrocytes.All the mentioned parameters are positive constants.
This model was discussed extensively in [10] and [18], here we shall highlight some of the most important results obtained in these works.The system (1) has two uniform steady states: the disease-free equilibrium P 0 = (0, 0, 0) and the non-trivial point P * = (1, β + δ, 1).The disease-free equilibrium is always unstable and it can be shown that in order for Turing instability to be possible (i.e. to allow the possibility of pattern formation) the inequality χ ∞ > χ c must be satisfied, where Upon extensive search in specialized literature and estimation of the unavailable parameters, the intervals for the parameters in the nondimensional form of the model (1), proposed by the authors of [10], are presented in Table I.Upon inspection of the model, proposed above, it can easily be deduced that the third equation implies a time derivative of the destroyed oligodendrocyte function that is (for appropriate initial conditions) non-negative at all times.This, in turn, guarantees that the function d(t) is monotonously increasing.Such a scenario proves to be too unrealistic in many cases, however, since remyelination is regularly examined to be a concomitant process with multiple sclerosis in many patients [19].Hence, this renders the model in its current form unsuitable for those cases and we now modify the model, in order to account for remyelination.
Remyelination in MS patients is known to sometimes begin taking place several weeks after lesion genesis [20], which indicates a certain delay between pathological effects and CNS reaction.One factor that could explain such a delay is the theory that Oligodendrocyte progenitor cells (OPCs) are mainly responsible for the generation of new myelin [21].OPCs are a group of slowly dividing non-neuronal cells which, after reaching maturity, are hypothesized to lose their division capabilities and become fully functioning oligodendrocytes creating new myelin.This scenario would explain the delay in the remyelinating processes, since the following two steps would need to take place before new myelin can be produced: 1) recruitment phase-OPCs proliferate in order to populate demyelinated areas; 2) differentiation phase-first OPCs differentiate into premyelinating oligodendrocytes which in turn contact demyelinated axons and differentiate into mature, myelinating oligodendrocytes that form functional myelin sheaths.
Having the aforementioned in mind, we modify the third equation in system (1) by adding a remyelination term.To obtain the desired effect, we turn to a delay partial differential equation.
The levels of remyelination in affected areas have been shown to correlate significantly to the concentration of macrophages and oligodendrocytes in the lesion [22].For simplicity, we shall propose a connection only between the part of destroyed oligodendrocytes and the rate of myelin restoration.We shall assume that the term responsible for remyelination depends on a previous moment in time, with ∆ 0 being the delay.We postulate the rate of remyelination to be proportional to the part of the destroyed oligodendrocytes at the delayed moment with a scale factor of µ.
Therefore, we consider the model where D ∆0 is a time-delay operator, defined with We allow that r = r(t) > 0 and µ = µ(t) ≥ 0, ∀t be time-dependent as we discuss in the later sections with all other parameters being positive constants.
In order to write the model in coordinate form, we use polar coordinates: The form of the differential operators ∇ and ∆ in polar coordinates is as follows [23]: ( Now we shall study the special case, in which we impose an assumption of radial symmetry, i.e. m, c and d do not depend on the angle ϕ (this is, e.g., consistent with the pattern formation in Baló's concentric disease).Thus, m, c, and d are only functions of ρ and t and we obtain the following coordinate form of the model, which is 1D with respect to space: where χ and F are defined with (2), ρ ∈ [0, R], t ∈ [0, T ], and all time-dependent quantities are evaluated at time t, except for the last term in the third equation (thus, we state explicitly the point of time only in this term and omit the arguments of all other functions for simplicity of notation).We close the system by imposing homogeneous Neumann boundary conditions for m and c at the right boundary: symmetry boundary conditions at the left boundary: and initial conditions: For d, we need to impose an initial history function instead of an initial condition.We do this in the following manner: Note that under these conditions, the modified model's solution will be identical to the solution of (1) in the time interval [0, ∆ 0 ).This lets us directly compare the solutions of the two models.

III. FINITE DIFFERENCE APPROXIMATION
In this section, we construct an explicit finite difference approximation of the model ( 6)- (10).First we discretize time and space by introducing the uniform mesh ωs × ωt , where and h s and h t are the respective space and time discretization steps.Next we introduce the approximate solutions of the model which are defined for the nodes of the mesh:

A. Approximation of the main equations
Now we are ready to approximate equations ( 6) by substituting finite difference formulae for the time and spatial derivatives in such a way that second order of local approximation is obtained with respect to the spatial discretization step.
For the time derivatives, we utilise a forward difference formula of the form [24,25]: which for sufficiently smooth functions f has an approximation error O(∆t).We approximate the second order spatial derivatives, resulting from the modelling of the diffusive and chemotactic processes, using the following formula: which has been shown to have an approximation error O(∆ρ 2 ) for sufficiently smooth functions α and f [26,27].Thus, we obtain: where θ := ∆ 0 h t (in the numerical experiments we shall use such a value of ∆ 0 , so that ∆ 0 /h t is a whole number), and r j := r(t j ), µ j := µ(t j ).

B. Approximation of the initial and boundary conditions
We approximate the initial conditions exactly through the functions given in ( 9): The boundary conditions can be approximated with the following finite difference formulae: • for the left boundary we apply the following threepoint right-sided finite difference formula: • for the right boundary we apply the following three-point left-sided finite difference formula: Both formulae have local approximation error O(∆ρ 2 ) for sufficiently smooth functions f [28].Thus, we obtain the following finite difference approximations of the boundary conditions: Order of Convergence Fig. 2: Estimated order of convergence for the numerical scheme ( 11)-( 14).

C. Practical estimate of the order of convergence
We study the convergence of the finite difference scheme numerically, by using Runge's method over three nested meshes [24,27].We show second order of convergence with respect to h s .Let us note that h t is taken in the numerical experiments of the order O(h 2 t ) from considerations of stability and, thus, we study the convergence with respect to h s .
The core concept of the method is the following-we solve the system numerically over three nested meshes with step sizes h s , h s /2, and h s /4, respectively.If we denote the approximate solution at the common points of the three meshes (i.e., the nodes of the coarse mesh) with y j i,hs , y j i,hs/2 , and y j i,hs/4 , respectively, one can derive the following estimate of the order of convergence p at a given point: .
Here, we present the results from one such experiment, which is representative for all parameter values that are further used in the paper (thus, we postpone the assignment of specific values to the model parameters until the next section, even though, formally speaking, we present the results from a single numerical experiment).To estimate the order of convergence we conduct an experiment, which illustrates the convergence order of the presented numerical scheme, given sufficiently smooth initial conditions.Using Runge's method, we are able to show that in almost all nodes of the space discretization at least second order of convergence is observed (see Fig. 2).

IV. NUMERICAL EXPERIMENTS
In this section, we present results from the implementation of the finite difference scheme, derived in Section III.We set the following initial conditions: corresponding to the presence of macrophages near the centre of the lesion at time t = 0 with all oligodendrocytes being still healthy.Let us further fix the following parameter values, which are used in what follows, unless stated otherwise: The numerical domain is defined with T = 20, R = 20 and the discretization steps are chosen to be h s = 0.2, h t = 0.0039.For all the figures, illustrating the results from the numerical experiments, we have reverted back to dimensional variables in order to obtain more informative images of disease progression.Thus, let us note that in dimensional quantities T = 20 corresponds approximately to 13.2 days and R = 20 is equal to roughly 5 mm.Since the authors of [10] give intervals as estimations for the dimensional values of b/ν and α/ν, we assume, for definiteness, these values to be respectively b/ν = 10 −6 and α/ν = 10 −2 .
A. Numerical study of the effect of remyelination strength on the model solutions First, we conduct several experiments to determine the influence of the remyelination strength on the quantity of destroyed oligodendrocytes.We assume for those experiments that µ is a constant.We compare the results for three possible values of µ: µ = 0 (i.e.no remyelination), and µ = 1, µ = 10.
We present the results for the density of macrophages at time t = T in Fig. 3. Here, there is no distinctive difference in the quantitative behaviour (in terms of number of clusters of macrophages and values of the local maxima) for varying remyelination intensity.However heterogeneities in the spatial distribution of the macrophage clusters are observed.
For the cytokine distributions at time t = T , see Fig. 4. Here, in addition to the varying spatial distribution, we can also clearly document significant quantitative variation between the three experiments, specifically that higher remyelination intensity leads to lower cytokine levels.For the third cell species-the destroyed oligodendrocytes, a more detailed analysis is in order, since the effects of the varying remyelination coefficent has a direct relation to their concentration in contrast to the two previous cell types.Furthermore, we aim to ultimately recreate different MS types, and the most direct route to that goal would be to better understand the most direct cause (lack of myelin produced by functioning oligodendrocytes) of the clinical symptoms observed.The graphs for the destroyed oligodendrocytes at time points t = 7.5, t = 12, t = 20 can be seen in Figures 5a-5c, respectively.
From these graphs, several effects of higher remyelination rates can be documented, namely:  intensive remyelination; • Lower demyelination levels in the most negatively affected areas; • Larger myelinated areas separating the more severely demyelinated ones.The most obvious consequence of the introduction of remyelination into the model is the overall significantly lower levels of destroyed oligodendrocytes.A somewhat more indirect effect of the remyelinating term is observed in the density of the macrophages and in the concentration of cytokines.These are phenomena that can hopefully be observed, once more clarity is present on the mechanisms and quantification of remyelinating processes.
Since the model assumes that cytokines are being produced by destroyed oligodendrocytes, the simulations with non-zero values for µ leads to lower cytokine levels as a consequence of the fewer destroyed oligodendrocytes observed due to remyelination.This, in turn, would appear to also influence the diffusionbased spatial dynamics of the cytokines, as visible from Fig. 4.
The leftward shift of the cytokine levels' extrema appear to consequently influence the macrophages' chemotactic processes by forming the clusters in different positions as compared to the non-remyelinating case (note how the maxima of the cytokine concentration and macrophage density appear in what seem to be identical positions for the different remyelination levels in Figures 3 and 4).
Also worth pointing out is that higher remyelination levels appear to result in overall slower processes of disease progression.
Upon further inspection including longer periods of time (i.e. more than 100 days) no equilibrium spatial profile of the destroyed oligodendrocytes could be observed for µ > 0 in contrast to the result proven for the original model in [10].In fact, the destroyed oligodendrocytes appear to move in clusters towards the focal point of the inflammation (i.e.ρ = 0) with some pairs of clusters merging into a single one while new clusters spontaneously arise from healthy oligodendrocytes.Analogous behaviour is displayed by the macrophages and cytokines.The results of these numerical experiments lead us to believe that, although still lacking an analytical proof, the modified model tends to no steady state for this type of choice for the parameters-including positive constant values for µ and r (see Fig. 6).The latter observation could possibly be connected to the relapses and remissions in the active forms of MS as we further discuss in subsection IV-C.

B. Numerical study of the effect of remyelination delay on the model solutions
Now we shall study how varying the values of ∆ 0 affects the model solutions.For the following experiments the same parameter values as before were retained with µ = 10.The large value for µ is chosen with the prospect of stronger remyelination potentially reacting more sensitively to delay variation and thus illustrating the effects of ∆ 0 more clearly.
In the experiments we focus on the destroyed oligodendrocytes due to their importance in disability development and since the macrophages and cytokines display analogous reaction to varying delay rates.We compare three cases with ∆ 0 = 10h t , ∆ 0 = 20h t and ∆ 0 = 40h t respectively, which translates to 37, 74 and 148 minutes in dimensional quantities.The main consequence of larger ∆ 0 values is best seen when two clusters are in the process of merging.These phenomena are observed at an earlier moment in time when ∆ 0 is large (see Fig. 7).This indicates that disease progression for longer remyelination delays is faster when compared to a faster recuperative reaction from the CNS.Although the relationships between the remyelination delay with disease progression are still in- sufficiently well understood, it would seem biologically plausible that bigger delay in recuperative processes would lead to more rapid disease progression, e.g due to the longer initial period of uninterrupted degenerative processes until the first effects of remyelination take place, among other reasons.The varying of the values µ and ∆ 0 is also of practical interest since older patients typically display weaker regenerative processes [29] which would translate to large values of ∆ 0 and smaller values for µ in terms of our model.

C. Ability to describe different types of MS
A tool for measuring the extent of disease progression was developed in 1983 by neurologist John Kurtzke [30].Named the Extended Disability Status Scale (EDSS), this scale takes into account examinations of various functional systems in the brain and yields a score ranging from 0-lack of any diagnosed disability, to 10-death caused by the disease.We shall most prominently distinguish between the following types of MS: • Relapsing-remitting MS (RRMS): This is the most common form.Around 85% of people with MS are initially diagnosed with RRMS.This type of MS involves episodes of new or increasing symptoms, followed by periods of remission, during which symptoms go away partially or totally.In a more recent 2014 work [31], a further modification to MS classification was presented without coming in conflict or abolishing the established in 1996 disease subtypes.Essentially, the existing MS types were further subdivided in accordance to their status regarding the so called "disease activity" and "disease progression".Thus, the authors speak about active/nonactive and progressive/non-progressive disease, meaning as follows: • Active disease-relapses, acute or subacute episodes of new or increasing neurologic dysfunction followed by full or partial recovery, in the absence of fever or infection.• Progressive disease-steadily increasing objectively documented neurologic dysfunction/disability without unequivocal recovery (fluctuations and phases of stability may occur).In this work we shall use the term "worsening" instead of "progressive" as to avoid confusion with the other progressive disease types.
For our next simulations, we shall vary the model parameters in order to show the highly enriched dynamics of the model (with respect to the model without remyelination) and its ability to mirror more adequately the different types of multiple sclerosis.In order to do so, we monitor the temporal evolution of the destroyed oligodendrocytes at various fixed distances from the focal point, since it has been well established in neuroscience that different parts of the brain relate to different functionality and thus studying a localized part of the CNS can be expected to be directly connected to concrete symptoms which can then be directly measured.1) PPMS (non-active worsening MS): The primaryprogressive MS type is the only one well described by the original model (i.e without remyelination), since it obviously leads to monotonously increasing values of d.For this reason we set the remyelinating coefficient µ = 0, rendering the oligodendrocytes incapable of recuperation once after being destroyed.All the other parameters are the same as in the previous examples.The results of the experiment are shown in Fig. 8 2) Benign RRMS (active and non-worsening MS): For this experiment, we set µ = 10 and reuse all other parameter values from previous simulations.In the case of constant coefficients with a non-negative value for µ, an interesting pattern is observed, and namely reminiscent of the relapses and remissions usually documented on a macro-level (see Fig. 9).It should be noted that the lack of significant residual damage after each relapse can be linked to the benign MS type, a somewhat loosely defined variant of the RRMS, characterized by extremely slow to no disability accumulation and asymptomatic periods between relapses [32].
3) SPMS (MS transitioning from active and nonworsening to non-active and worsening): Our next goal is to recreate a scenario, which emulates a transition from a relapsing-remitting start with no significant worsening after the first relapses to a progressive disease course, uninterrupted by relapses (namely, the

SPMS subtype).
For the purpose of recreating the CNS's decreasing remyelinating capacity, we utilise a monotonously decreasing function of time for the remyelinating coefficient and more precisely µ(t) = max 10 T − t T , 0 .
The results indeed show a qualitative behaviour close to the desired one, as shown in Fig. 10.Note, that the subtle increase in destroyed oligodendrocytes between the first and second relapse could indicate that the reported "stability" between relapses in the initial stages of SPMS might in fact prove to simply be a misinterpretation of the much more slowly progressing disability between these periods.4) PRMS (active and worsening MS): As we aim to recreate the final of the main MS types-PRMS, we modify the damaging coefficient r to also be a time-dependent function as opposed to being a constant.More precisely we assume that r follows an increasing  linear law of the form r(t) = 60 t T .In this manner, in retaining the other coefficient values from the previous experiment, and more specifically the monotonously decreasing remyelinating coefficient µ(t), we achieve a behaviour of the destroyed oligodendrocytes, which mimics the PRMS disease type, characterized by ongoing continuous worsening, intertwined with concomitant relapses, see Fig. 11.A comparison of the spatial distribution of the destroyed oligodendrocytes at two different time points for the SPMS and PRMS are given in Fig. 12. Near disease onset the two types of MS progress in a similar manner (see Fig. 12a), but after a certain point the differences between them become more apparent (see Fig. 12b).We can see that PRMS disease progression is more rapid than in SPMS, which is also a well documented clinical fact [33].

D. Better agreement with experimental data
Next, we shall compare the results from the original model ( 3) and the modified one (6) from another point of view-their agreement with experimental data.In Fig. 13 we compare the solutions of the models with and without remyelination with parameter values taken from the numerical experiments in section IV and remyelination coefficients respectively µ = 0 and µ = 1.Since we are mainly interested with the progression of the disease, we depict only the graphs of the destroyed oligodendrocytes.The result is visualised for time t = 20.
As seen from the illustrations, the effects of remyelination lead to a much steeper slope between the remyelinated and demyelinated zones, as well as areas of the nervous system that remain almost practically unscathed.This is much more strongly reminiscent of the images of the nervous system observed in patients with Baló's concentric sclerosis (see Fig. 14a).Fig. 14b compares the results obtained from the model respectively with and without remyelination, with the former bearing a significantly larger resemblance to the experimentally obtained images, mainly due to the presence of completely myelinated areas as is the case in Fig. 14a, while the original models displays only partially or totally demyelinated areas.

V. CONCLUSION AND DISCUSSION
In the current work, we have proposed a modification of an existing PDE model of multiple sclerosis in order to enable the observance of a wider variety of MS subtypes by the means of a single mathematical model.We have achieved this mainly through the addition of a remyelination term, which in turn relies upon a time delay operator for the purpose of modelling the organism's latency in reacting to the accrual of CNS damage.Also a key assumption of our modification to the model was to allow for the respective coefficients of demyelination and remyelination intensity to vary over  time, a property of MS which has been documented e.g.regarding the decreasing intensity of the inflammatory reactions and recuperation in SPMS [35].
The model, in its new form, was capable of displaying behaviour, qualitatively reminiscent of that found in RRMS, PPMS, SPMS and PRMS which encompasses all of the most commonly described types of MS.However, numerous challenges in regards to "perfecting" the modelling of MS still remain unresolved and require intensive future study and consideration.
From a medical point of view, the underlying processes need to be more thoroughly studied and described for the purpose of further modifying the model's Once more light has been shed on the mechanisms behind the pathogenesis of MS and a more suitable model has been formulated, one of the next steps would be to examine each of the parameter's influence on the results of the simulation, so that the model will not only be capable of displaying satisfactory qualitative behaviour but also more realistic quantitative results can be obtained, making it possible to predict the disease progression in individual-specific cases.
From a mathematical point of view, further development of the theory behind delay differential equations is necessary in order to better evaluate the practical properties of such models, such as conducting stability analysis and determining the conditions for convergence or monotonicity for difference schemes, applied to the model.

(cells per mm 2 )Fig. 5 :
Fig. 5: Progression of the density of destroyed oligodendrocytes for various values of µ
(a) Lesions in Baló's concentric sclerosis [34].(b) Baló's disease, as represented by the destroyed oligodendrocytes d at time t = 20 (13.2 days), obtained in the numerical experiments without (left, µ = 0) and with (right, µ = 1) remyelination.The grey scale is defined as follows.White corresponds to no destroyed oligodendrocytes, while black corresponds to completely demielynated areas.

Fig. 14 :
Fig. 14: Comparison between the observed CNS damage in Baló's concentric disease and the results yielded by the numerical simulation (white areas represent more severe oligodendrocyte damage).

Table I :
[10]mated possible intervals for the nondimensional version of the model[10].
Secondary progressive MS (SPMS): At first, people will experience episodes of relapse and remission, but then the disease will start to progress steadily.
• Primary progressive MS (PPMS): Symptoms worsen progressively, without early relapses or remissions.Some people may experience times of stability and periods when symptoms worsen and Biomath 12 (2023), 2309267, https://doi.org/10.55630/j.biomath.2023.09.267 9/13 then get better.Around 15% of people with MS have PPMS.• Biomath 12 (2023), 2309267, https://doi.org/10.55630/j.biomath.2023.09.267 12/13 G. M. Bazlyankov, T. B. Ivanov, On the effect of remyelination in a mathematical model of multiple sclerosis equations to more truthfully describe the real interactions between the main agents involved in MS.This would in all likelihood contribute to the elimination of inaccuracies from the model's solution such as the observed minima with close to zero values preceding the relapses in e.g.Figures 10 and 11.