¹¹institutetext: Romina Travaglini ²²institutetext: Istituto Nazionale di Alta Matematica, Università degli Studi di Parma
²²email: romina.travaglini@unipr.it

A reaction-diffusion model for relapsing-remitting multiple sclerosis with a treatment term

Romina Travaglini\orcidID0000-0003-4107-1764

Abstract

We present a mathematical study for the development of multiple sclerosis based on a reaction-diffusion system. The model describes interactions among different populations of human cells, motion of immune cells stimulated by cytokines, consumption of myelin sheath due to anomalously activated lymphocytes and its restoration by oligodendrocytes. Successively, we introduce a therapy term representing injection of low-dose IL-2 interleukine. A natural step is then to study the system, investigating the formation of spatial patterns by means of a Turing instability analysis of the problem. In particular, we get spatial patterns oscillating in time that may reproduce brain lesions characteristic of the early stage of the pathology, in both non-treatment and treatment scenarios.

1 Introduction

Multiple Sclerosis (MS) is one of the most severe autoimmune disorders. Its primary pathological characteristic is injury of the myelin sheath, that surrounds axons in the nervous system (CNS) and favors the transmission of cerebral impulses. Such lesions can be observed, through MRI, as development of focal plaques in the white matter.

It is widely accepted that the biological dynamics behind MS involves immune cells as T-cells, B-cells, macrophages and microglia. These cells may be dysfunctionally activated against cells producing myelin or myelin itself. This activation may induce an autoimmune cascade, promoted by production of proinflammatory cytokines, i.e. molecules that attract and stimulate clonal expansion of immune cells.

Medical literature on MS reports a wide range of cases in terms of clinical course, features of lesions and associated irreversible neurological symptoms. Usually, the early phase of the disease is characterized by acute inflammation and formation of active lesions. At the same time, also a restoration process, referred as remyelination, takes place, resulting in the formation of ”shadow” plaques. This demyelination-remyelination process may continue to alternate for months or even years and this phase is indicated as relapsing-remitting MS (RRMS). This stage is usually followed by a second phase, called secondary progressive MS (SPMS), when restoration of active myelin lesions is less frequent, resulting in a progressive neurodegeneration. Further details about medical description of MS may be found in lassmann2005multiple ; lassmann2007immunopathology ; lassmann2012progressive ; mahad2015pathological and references therein.

For our purposes, it must be underlined the fact that self-reactive immune cells can be found in non-pathological conditions as well danke2004autoreactive . But in this case, the action of specific cells called immunosuppressors can inhibit or kill cells presenting the antigen and activated immune cells. In autoimmune conditions as MS, though, both the number and the efficiency of these natural killers are lacking hoglund2013one ; mimpen2020natural ; zozulya2008role .

In a recent work travaglini2023reaction a mathematical model describing peculiar dynamics of MS has been proposed. In particular, starting from kinetic models of active particles applied to autoimmune diseases as delitala2013mathematical ; kolev2018mathematical , followed by ramos2019kinetic ; della2022mathematical , authors derive a system of partial-differential equations of reaction-diffusion type with a chemotaxis term obtained by a diffusive limit in aproper time scaling. Populations of cells and biological substances there considered are: antigen presenting cells (APCs), self-reactive leukocytes (SRLs), immunosuppressive cells (ISCs), cytokines and consumed myelin. Through a Turing instability analysis, they reproduce formation of one-dimensional patterns whose dynamics is comparable to RRMS or SPMS.

The aim of the present work is to focus on the RRMS dynamics, introducing a treatment term representing injection of low-dose interleukin L-2 (IL-2). As reported in literature, in fact, in autoimmune conditions furtado2002interleukin , low-dose IL-2 may stimulate expansion of ISCs klatzmann2015promise and enhance their suppressive function. The manner in which these effects can be achieved in the specific case of MS, coupled with the stimulation of ISC capacity to promote remyelination, is currently under investigation louapre2023randomized .

In this paper, we perform a reduction of the reaction-diffusion system derived in travaglini2023reaction and introduce the new term representing effects of therapy. Then we propose a Turing instability analysis of the system, in order to find conditions on parameters leading to formation and restoration of plaques. Results will be validated and compared by means of numerical simulations in a two-dimensional domain.

2 The model

We recall here the macroscopic reaction-diffusion model derived in travaglini2023reaction that describes formation and restoration of myelin plaques in multiple sclerosis, in a non-dimensional form. The system of reaction-diffusion equations reads as

$\displaystyle\frac{\partial A}{\partial t}$	$\displaystyle=\,$	$\displaystyle 1+\beta\,A\,R-A\,S-\zeta A,$
$\displaystyle\frac{\partial S}{\partial t}$	$\displaystyle=\,$	$\displaystyle\mu\,A\,S-\,S,$
$\displaystyle\frac{\partial R}{\partial t}$	$\displaystyle=\,$	$\displaystyle\nabla_{\mathbf{x}}\cdot\left(\Phi_{0}(R)\,\nabla_{\mathbf{x}}\,R-{\chi}\,\Phi_{1}(R)R\,\nabla_{\mathbf{x}}\,C\right)+\eta\,A\,R-\psi\,R\,S-\theta\,R,$
$\displaystyle\frac{\partial C}{\partial t}$	$\displaystyle=\,$	$\displaystyle\delta\Delta_{\mathbf{x}}\,C+A\,R-\tau\,C,$
$\displaystyle\frac{\partial E}{\partial t}$	$\displaystyle=\,$	$\displaystyle\frac{\gamma\,R}{\xi+R}\,R\,(1-E)-\lambda\,E,$

where involved quantities are macroscopic densities for

$A$ – Self-antigen presenting cells (SAPCs).
$S$ – Immunosuppressive cells (ISCs).
$R$ – Self-reactive leukocytes (SRLs).
$C$ – Cytokines.
$E$ – Destroyed myelin,

depending on $t\in\mathbb{R}_{0}^{+}$ and $\mathbf{x}\in\Omega$ , with $\Omega$ bounded domain in $\mathbb{R}^{2}$ . Moreover, we choose

\Phi_{0}(y)=\cos\left(\frac{\pi}{2}\,y\right)+\frac{\pi}{2}\,y\,\sin\left(\frac{\pi}{2}\,y\right),\quad\Phi_{1}(y)=\cos\left(\frac{\pi}{2}\,y\right).

A detailed discussion about more general choice for spatial domain and functions can be found in travaglini2023reaction .

As first, since the aim of the present work is to focus on the effects of low-dose IL-2 on ISCs, SRLs and myelin sheath, we assume that the population of SAPCs constitutes a background which SRLs and ISCs interact with. Henceforth, we suppose that SAPCs density is constant at equilibrium, i.e.

\displaystyle A

\displaystyle=\,

\displaystyle\frac{1}{\zeta-\beta\,R+S}.

Moreover, we introduce the parameter $\alpha\in[0,1)$ representing low-dose IL-2 delivery, obtaining the system

$\displaystyle\frac{\partial S}{\partial t}$	$\displaystyle=$	$\displaystyle\,S\,\left(\frac{\mu}{\zeta-\beta\,R+S}-\Gamma\right),$	(1)
$\displaystyle\frac{\partial R}{\partial t}$	$\displaystyle=\,$	$\displaystyle\nabla_{\mathbf{x}}\cdot\left(\Phi_{0}(R)\,\nabla_{\mathbf{x}}\,R-{\xi}\,\Phi_{1}(R)R\,\nabla_{\mathbf{x}}\,C\right)$	(2)
		$\displaystyle\,+R\left(\frac{\eta}{\zeta-\beta\,R+S}-\Psi\,S-\theta\,\right),$	(2)
$\displaystyle\frac{\partial C}{\partial t}$	$\displaystyle=\,$	$\displaystyle\delta\Delta_{\mathbf{x}}\,C+\frac{R}{\zeta-\beta\,R+S}-\tau\,C,$	(3)
$\displaystyle\frac{\partial E}{\partial t}$	$\displaystyle=\,$	$\displaystyle\frac{\gamma\,R}{\xi+R}\,R\,(1-E)-\Lambda\,E,$	(4)

where $\Gamma=(1-\alpha)$ , $\Psi=\psi\,(1+\alpha)$ , $\Lambda=\lambda\,(1+\alpha)$ .

2.1 Turing instability

We study macroscopic system (1)-(4) by means of Turing instability analysis. More in detail, we aim at finding a particular range for parameters leading to the appearance of oscillating spatial patterns that may reproduce the appearance (and possible reconstruction) of myelin plaques.

We consider system (1)-(4) taking as initial data

\mathbf{U}(0,\mathbf{x})=\mathbf{U}_{0}(\mathbf{x})\geq 0,\mbox{ with }R_{0}(\mathbf{x})\leq 1,\,E_{0}<1,

and imposing zero-flux at the boundary

\Big{(}\Phi_{0}(R)\,\nabla_{\mathbf{x}}\,R-{\xi}\,\Phi_{1}(R)R\,\nabla_{\mathbf{x}}\,C\Big{)}\cdot{\bf\widehat{n}}=0,\quad\nabla_{\mathbf{x}}C\cdot{\bf\widehat{n}}=0,

with ${\bf\widehat{n}}$ being the external unit normal to the boundary $\partial\Omega$ . Turing instability turing52 occurs when a spatially homogeneous steady state for the system without spatial gradients turns to be unstable when diffusive and chemotactic terms are added. To this aim, we individuate an equilibrium $U_{1}=(S_{1},\,R_{1},\,C_{1},\,E_{1})$ for the system (1)-(4), in spatially homogeneous conditions, that is biologically relevant, i.e. that belongs to the set $\mathcal{E}=\left\{0<R(t,\mathbf{x})\leq 1,\,C(t,\mathbf{x})>0,\,0<E(t,\mathbf{x})<1\right\}.$ Thus we find

U_{1}=\,\left(\frac{\Gamma\,\eta-\theta\,\mu}{\mu\,\Psi},\,\frac{\Gamma\,\zeta\,\mu\,\Psi+\Gamma^{2}\,\eta-\Gamma\,\theta\,\mu-\mu^{2}\,\Psi}{\beta\,\mu\,\Psi},\,\frac{1}{\mu\tau}\,R_{1},\,\frac{R_{1}^{2}\gamma}{R_{1}^{2}\,\gamma+R_{1}\,\Lambda+\Lambda\,\xi}\right),

that belongs to the set $\mathcal{E}$ when the following conditions on parameters are satisfied

\theta<\frac{\Gamma\,\eta}{\mu},\quad\bar{\theta}>\theta>\bar{\theta}-{\beta\Psi},\quad\mbox{with}\quad\bar{\theta}:=\frac{\Gamma\,\eta}{\mu}+\frac{\Psi(\Gamma\,\zeta-\mu)}{\Gamma}.

(5)

We find it convenient to take $\Gamma\,\zeta>\mu$ , in such a way conditions in (5) become

\frac{\Gamma\,\eta}{\mu}>\theta>\bar{\theta}-{\beta\Psi}.

Linearizing system (1)-(4) around equilibrium $U_{1}$ , we have the Jacobian matrix defined as

\mathbb{A}=\left(\begin{array}[]{cccc}-\frac{\Gamma^{2}\,S_{1}}{\mu}&\frac{\Gamma^{2}\,\beta\,\,S_{1}}{\mu}&0&0\\ \\ -\frac{R_{1}\left(\Gamma^{2}\eta+\mu^{2}\Psi\right)}{\mu^{2}}&\frac{\Gamma^{2}\,\beta\,\eta\,R_{1}}{\mu^{2}}&0&0\\ \\ -\frac{\Gamma^{2}\,R_{1}}{\mu^{2}}&\frac{\Gamma}{\mu}+\frac{\Gamma^{2}\,\beta\,{R_{1}}}{\mu^{2}}&-\tau&0\\ \\ 0&\frac{\gamma\,\Lambda\,\ R_{1}(2\,\xi+R_{1})}{(\xi+R_{1})(\xi\,\Lambda+R_{1}(\Lambda+\gamma\,R_{1}))}&0&-\frac{\gamma\,R_{1}^{2}}{\xi+R_{1}}-\Lambda\end{array}\right).

The eigenvalues of the matrix $\mathbb{A}$ are $-\tau$ and $-\Lambda-\frac{\gamma R_{1}^{2}}{\xi+R_{1}}$ , and the eigenvalues $\sigma_{1},\,\sigma_{2}$ of the 2x2 square submatrix in the top left corner. It holds that

\sigma_{1}+\sigma_{2}=\frac{\Gamma^{2}}{\mu}\left(\frac{R_{1}\,\eta\,\beta}{\mu}-S_{1}\right),\quad\sigma_{1}\,\sigma_{2}=\frac{\Gamma^{2}\,\beta\,\Psi}{\mu}\,R_{1}\,S_{1}>0.

Thus, the equilibrium $U_{1}$ is locally asymptotically stable holding the condition

\theta<\tilde{\theta},\quad\tilde{\theta}=\frac{\Gamma\,\eta}{\mu}+\frac{\eta\,\Psi\,(\Gamma\,\zeta-\mu)}{\Gamma\,(\eta-\mu)}.

Now, we consider the complete system with the spatial derivative terms. When linearizing, we get the diffusion matrix $\mathbb{D}$ having all entries null except for elements in positions (2,2), (2,3) and (3,3), that are equal to $\Phi_{0}(R_{1})$ , $-\Lambda\Phi_{1}(R_{1})R_{1}$ and $\delta$ , respectively. Formation of spatial pattern may arise when it is possible to find an interval $(k_{1},k_{2})$ such that for any $k_{1}\leq k\leq k_{2}$ , $Re(\lambda_{k})>0$ , being $\lambda_{k}$ an eigenvalue for the matrix $\mathbb{A}-k^{2}\mathbb{D}$ . We can individuate only a sufficient condition to have eigenvalues with positive real part, that is $det(\mathbb{A}-k^{2}\mathbb{D})<0$ , where

\det({\mathbb{A}}-k^{2}{{\mathbb{D}}})=\frac{\Gamma^{2}\,S_{1}}{\mu}\left(\Lambda+\frac{\gamma R_{1}^{2}}{\xi+R_{1}}\right)\,h(k^{2}),

h(k^{2}):=k^{4}\delta\,\Phi_{0}(R_{1})+k^{2}q+\tau\,R_{1}\,\beta\,\Psi,

q=R_{1}\,\delta\,\beta\,\Psi-\frac{\Gamma\,\chi}{\mu}\,\Phi_{1}(R_{1})\,R_{1}+\Phi_{0}(R_{1})\,\tau.

It is straightforward to observe that the condition $\det({\mathbb{A}}-k^{2}{{\mathbb{D}}})<0$ is satisfied if and only if we have the following two

q<0\mbox{ and }q^{2}-4\,R_{1}\,\delta\,\beta\,\Psi\,\Phi_{0}(R_{1})\tau>0,

that, after computations, lead to the a threshold value for $\chi$ that reads as

\chi>\mu\,\frac{2\sqrt{\Phi_{0}(R_{1})\,R_{1}\,\delta\,\tau\,\beta\,\Psi}+\delta\,R_{1}\,\beta\,\Psi+\Phi_{0}(R_{1})\,\tau}{\Gamma\,\Phi_{1}(R_{1})\,R_{1}}.

3 Numerical simulations

Analytical results obtained in the previous section are now applied to the study of quantities involved in system (1)-(4) in space and time. We perform numerical simulation using VisualPDE walker2023visualpde in a square domain of size $L=10$ . In particular, we focus on the dynamics of destroyed myelin $E$ . We fix the following parameters

\begin{array}[]{c}\beta=0.8,\quad\zeta=2.05,\quad\mu=2,\quad\tau=0.5,\quad\eta=1,\quad\delta=0.1,\quad\theta=0.4,\\ \\ \psi=0.3,\quad\gamma=2,\quad\xi=1,\quad\lambda=0.3.\end{array}

(6)

Vales above individuate an interval for the treatment parameter $\alpha$ that is $0\leq\alpha<0.095$ . We start by simulating the non-treatment, setting $\alpha=0$ and as initial conditions a perturbation of equilibrium $U_{1}$ for quantities $S,R,C$ while $E$ is taken uniformly null at the beginning. In this case we have a critical $\chi_{c}\approx 4.59$ , then we fix $\chi=7$ and we obtain results shown in Figure 1. It can be observed that the myelin plaques (areas where the density of $E$ is higher) have an oscillatory behavior in time, as happens in the case of RRMS.

Refer to caption — Figure 1: Evolution of quantity $E$ (destroyed myelin), taking values as in (6), $\alpha=0$ , $\chi=7$ .

For the second case, instead, we take $\alpha=0$ and initial conditions as in the previous case, then, when lesions appear, we increase $\alpha=0.05$ . The critical value for $\chi$ now is $\chi_{c}\approx 6.37$ , hence we keep $\chi=7$ , obtaining dynamics reported in Figure 2. It is possible to observe that, at time $t=85$ plaques completely disappear, differently from the first case. Moreover, lesions take a longer time to form again, we underline that this effect is indirectly induced by a lower decreasing and a major suppressive efficiency of ISCs, given by parameters $\Lambda$ and $\Psi$ , respectively. This fact affects negatively the growth of SRLs and, consequently the consumption of sane myelin.

4 Conclusions

In this work, we considered the model proposed in travaglini2023reaction , where a reaction-diffusion model for the study of multiple sclerosis is derived from kinetic level and analyzed. In particular, we focused on results that may describe the first phase of multiple sclerosis, called relapsing-remitting phase, characterized by formation and restoration of plaques in the white matter of the central nervous system. Our aim was, on the one hand to reproduce the behavior of the model performing simulations on a two-dimensional domain, on the other hand to study the effects of possible treatment terms. For this reason, we performed a reduction of the original system, introducing a term representing injection of low-dose interleukin L-2 and we went through a Turing instability analysis. We found conditions on parameters allowing for formation of patterns and we validated results numerically. We reproduced the periodic formation and restoration of myelin plaques, confirming that the treatment term helps the remyelination and delays demyelination.

Acknowledgements.

Tha author is a post-doc fellow of the National Institute of Advanced Mathematics (INdAM), Italy. The research was carried out in the frame of activities sponsored by the Cost Action CA18232.

\ethics

Competing Interests The research work was supported by INdAM (National Institute of Advanced Mathematics), by the Portuguese national funds (OE), through FCT/MCTES FCT/MCTES (Fundação para a Ciência e a Tecnologia) Projects UIDB/00013/2020, UIDP/00013/2020,
PTDC/03091/2022 (“Mathematical Modelling of Multi-scale Control Systems: applications to human diseases (CoSysM3)”), and by University of Parma through the action Bando di Ateneo 2022 per la ricerca co-funded by MUR-Italian Ministry of Universities and Research - D.M. 737/2021 - PNR - PNRR - NextGenerationEU” (project: ”Collective and self-organised dynamics: kinetic and network approaches”)

References

[1] N. A. Danke, D. M. Koelle, C. Yee, S. Beheray, and W. W. Kwok. Autoreactive t cells in healthy individuals. J. Immunol., 172(10):5967–5972, 2004.
[2] M. Delitala, U. Dianzani, T. Lorenzi, and M. Melensi. A mathematical model for immune and autoimmune response mediated by T-cells. Comput. Math. Appl., 66(6):1010–1023, 2013.
[3] R. Della Marca, M. d. P. Machado Ramos, C. Ribeiro, and A. J. Soares. Mathematical modelling of oscillating patterns for chronic autoimmune diseases. Math. Meth. Appl. Sci., 45(11):7144–7161, 2022.
[4] G. C. Furtado, M. A. C. de Lafaille, N. Kutchukhidze, and J. J. Lafaille. Interleukin 2 signaling is required for CD4+ regulatory T cell function. J. Exp. Med., 196(6):851–857, 2002.
[5] R. A. Høglund, T. Holmøy, H. F. Harbo, and A. A. Maghazachi. A one year follow-up study of natural killer and dendritic cells activities in multiple sclerosis patients receiving glatiramer acetate (GA). PLoS One, 8(4):e62237, 2013.
[6] D. Klatzmann and A. K. Abbas. The promise of low-dose interleukin-2 therapy for autoimmune and inflammatory diseases. Nat. Rev. Immunol., 15(5):283–294, 2015.
[7] M. Kolev and I. Nikolova. A mathematical model of some viral-induced autoimmune diseases. Math. Applicanda, 46(1), 2018.
[8] H. Lassmann. Multiple sclerosis pathology: evolution of pathogenetic concepts. Brain Pathol., 15(3):217–222, 2005.
[9] H. Lassmann, W. Brück, and C. F. Lucchinetti. The immunopathology of multiple sclerosis: an overview. Brain Pathol., 17(2):210–218, 2007.
[10] H. Lassmann, J. Van Horssen, and D. Mahad. Progressive multiple sclerosis: pathology and pathogenesis. Nat. Rev. Neurol., 8(11):647–656, 2012.
[11] C. Louapre, M. Rosenzwajg, M. Golse, A. Roux, F. Pitoiset, L. Adda, N. Tchitchek, C. Papeix, E. Maillart, A. Ungureanu, et al. A randomized double-blind placebo-controlled trial of low-dose interleukin-2 in relapsing–remitting multiple sclerosis. J. Neurol., pages 1–12, 2023.
[12] D. H. Mahad, B. D. Trapp, and H. Lassmann. Pathological mechanisms in progressive multiple sclerosis. Lancet Neurol., 14(2):183–193, 2015.
[13] M. Mimpen, J. Smolders, R. Hupperts, and J. Damoiseaux. Natural killer cells in multiple sclerosis: a review. Immunol. Lett., 222:1–11, 2020.
[14] J. M. Oliveira and R. Travaglini. Reaction–diffusion systems derived from kinetic theory for Multiple Sclerosis. Math. Models Methods Appl. Sci., 34(07):1279–1308, 2024.
[15] M. M. Ramos, C. Ribeiro, and A. Soares. A kinetic model of t cell autoreactivity in autoimmune diseases. J. Math. Biol., 79(6-7):2005–2031, 2019.
[16] A. M. Turing. The chemical basis of morphogenesis. Phil. Trans. R. Soc. Lond. B, 237:37–72, 1952.
[17] B. J. Walker, A. K. Townsend, A. K. Chudasama, and A. L. Krause. VisualPDE: rapid interactive simulations of partial differential equations. Bull. Math. Biol., 85(11):113, 2023.
[18] A. L. Zozulya and H. Wiendl. The role of regulatory T cells in multiple sclerosis. Nat. Clin. Pract. Neurol., 4(7):384–398, 2008.