1 Introduction

Multiple Sclerosis is a demyelinating disorder affecting the central nervous system and causing severe and progressive physical and neurological impairment. More specifically, it is characterized by inflammation and demyelination, resulting in the formation of focal areas of myelin loss in the white matter of the brain, called plaques or lesions [1]. Different histological patterns of plaques have been identified [2], but it is commonly believed that the onset of any typical new lesion is characterized by the same pathological changes [3].

Some mathematical models able to reproduce many of the typical pathological features of Multiple Sclerosis have been proposed and investigated; they are usually constituted by proper systems of PDEs [4,5,6,7], but also systems of ODEs [8] or stochastic models [9] are able to describe the relapse–remittance dynamics in patients affected by Multiple Sclerosis.

The model proposed in [7] is a PDE system describing the interaction between the density of activated immune cells (macrophages) \(m = m(t, x)\), the concentration of chemical species (cytokine) \(c = c(t, x)\) secreted by the immune cells, and the density of the destroyed oligodendrocytes \(d = d(t, x)\). Proper Turing instability analysis of the proposed reaction–diffusion system is performed, showing the formation of spatial patterns when the chemotactic coefficient overcomes a proper threshold; two-dimensional numerical simulations show the appearance of different kinds of patterns (as circular rings, or small clusters, etc.), similar to those really observed in a damaged brain (like Balo’s plaques). Spatial modulation of the Turing–type structures through Eckhaus and zigzag instability is also addressed in [10]. Rigorous results on this system have been proven in [11, 12], concerning the global existence of weak and strong solutions, uniform bounds in time for such solutions, and nonlinear stability results under proper assumptions on the chemotaxis term. On the other hand, the stability analysis fits the setting of PDEs–ODEs systems of the recent paper [13].

In this paper, we consider a modification of the model proposed in [7], changing the term describing the production and saturation of the activated macrophages. In [7], the rate of macrophages activation is modelled by means of a logistic functional form describing the proliferation and saturation effects taken into account also in previous ODE models for acute inflammation [14, 15]. However, even if it is commonly conjectured that activation of microglia is responsible for the appearance of Multiple Sclerosis lesions [16], the underlying mechanism still remains unknown [7]. It is also plausible that a minimal amount of activated macrophages is needed to start and sustain the inflammatory response. This case was not considered in [7], since they modelled the macrophage activation using a logistic growth function. In the present work, we thus take into account a different model for macrophages activation, by including the so-called Allee effect [17, 18], which is a growth function used in population dynamics to take into account undercrowding effects (see [19,20,21] and references therein). Mainly two types of Allee effects have been considered in the literature [22], namely the weak and the strong (or critical) Allee effect, depending on the fact that the prey growth function is non-negative or negative, respectively, for small population densities. Examples of the use of the Allee effect in biomedical applications can be found for instance in [23], to model in a liver inflammation the fact that a very small virus population can be eliminated by local immune reactions, and not necessarily increases as prescribed by a logistic growth.

The aim of this work is to investigate how the choice of a different model for macrophage activation can change the pattern scenarios presented and discussed in [7]. To this end, we will perform a weakly nonlinear analysis of the modified reaction–diffusion system, together with a bifurcation analysis. The weakly nonlinear approach has been adopted in [6, 7, 24, 25] and it is based on the multiple scales method; suitable expansions in terms of a parameter measuring the dimensionless distance with respect to a suitable critical bifurcation value allow to derive amplitude equations which yield the form of the pattern close to the bifurcation threshold. We perform this asymptotic procedure up to the third-order accuracy, getting a Stuart-Landau equation for the amplitude of spatially periodic solutions. On the other hand, we also perform a bifurcation analysis by computing the bifurcation structure and its deformation under parameter variations. In this regard, we exploit the numerical continuation software pde2path [26,27,28], based on a FEM discretisation of the stationary problem. The software has been recently used beyond its standard setting, for instance in cross-diffusion [29,30,31] and fractional diffusion problems [32].

The paper is organised as follows. In Sect. 2, the mathematical model with the Allee effect is introduced; in Sects. 3 and 4, we perform the Turing instability and the weakly nonlinear analysis, respectively. Numerical results, namely bifurcation diagrams obtained via pde2path and numerical simulations, are shown in Sect. 5. Finally, some concluding remarks can be found in Sect. 6. The Matlab scripts for the numerical continuation are freely accessible in the GitHub folder [33].

2 The mathematical model

We want to describe the initial stage of the disease, which is characterized by the production of pro-inflammatory cytokines \(\tilde{c}\) by macrophages and activated microglia \(\tilde{m}\) and by the emergence of lesions due to apoptosis of oligodendrocyte \(\tilde{d}\). We suppose that all the concentrations depend on the position \({\textbf{X}}\in \Omega \subset {\mathbb {R}}^n\) and time \(T\in {\mathbb {R}}^+\).

We consider the following system of partial differential equations

$$\begin{aligned} \begin{aligned}&\dfrac{\partial \tilde{m}}{\partial T}\,=\,D\Delta _{\textbf{X}}\tilde{m}+\tilde{G}_j(\tilde{m})-\nabla _{\textbf{X}}\cdot (\Psi (\tilde{m})\nabla _{\textbf{X}}\tilde{c}),&\text {on}\,\Omega \times {\mathbb {R}}^+,\\&\dfrac{\partial \tilde{c}}{\partial T}\,=\,\dfrac{1}{\nu }\,[\varepsilon \Delta _{\textbf{X}}\tilde{c}+\mu \tilde{d}-\alpha \tilde{c}+b\tilde{m}],&\text {on}\,\Omega \times {\mathbb {R}}^+,\\&\dfrac{\partial \tilde{d}}{\partial T}\,=\,\kappa \tilde{F}(\tilde{m})\tilde{m}(\bar{d}-\tilde{d}),&\text {on}\,\Omega \times {\mathbb {R}}^+,\\&\dfrac{\partial \tilde{m}}{\partial n}=\dfrac{\partial \tilde{c}}{\partial n}=\dfrac{\partial \tilde{d}}{\partial n}=0,&\text {on}\,\partial \Omega \times {\mathbb {R}}^+,\\&\tilde{m}({\textbf{X}},0)=\tilde{m}_{in}({\textbf{X}}),\, \tilde{c}({\textbf{X}},0)=\tilde{c}_{in}({\textbf{X}}),\, \tilde{d}({\textbf{X}},0)=\tilde{d}_{in}({\textbf{X}}),\,&\text {on}\,\Omega , \end{aligned} \end{aligned}$$

being \(\Omega \) a bounded and connected domain. The evolution of macrophages/microglia density is governed by random movements with constant diffusion coefficient D, a production/decay term \(\tilde{G}_j\) (where \(j=1,2\) denotes two different functions detailed below) and a chemotactic motion towards regions with a high concentration of chemoattractants; the function \(\Psi (\tilde{m})=\psi \tilde{m}/(\bar{m}+\tilde{m})\), where \(\bar{m}\) is the characteristic density of macrophages and \(\psi \) is the maximal chemotactic rate. The cytokines evolution is ruled by random diffusion with diffusion coefficient \(\varepsilon \) and by production/decay terms with constant proliferation coefficients \(\mu ,\, b\) and death coefficient \(\alpha \). For what concerns oligodendrocytes, their evolution is governed by a production term, where \(\kappa \) measures the destructive strength of the macrophages, \(\tilde{F}(\tilde{m})=\tilde{m}/(\bar{m}+\tilde{m})\) and \(\bar{d}\) is the characteristic density of oligodendrocytes.

Regarding the production term for macrophages \(\tilde{G}_j\), we shall consider, contrast and compare the effects of two different growth functions; the first choice is a logistic term

$$\begin{aligned} \tilde{G}_1(\tilde{m})=\lambda \tilde{m}(\bar{m}-\tilde{m}), \end{aligned}$$

where \(\lambda \) is the production rate; alternatively, the production rate can also depend on the concentration of macrophages and we shall consider a cubic function

$$\begin{aligned} \tilde{G}_2(\tilde{m})=\tilde{\lambda }\tilde{m}(\bar{m}-\tilde{m})(\tilde{m}-\hat{m}), \end{aligned}$$

where \(\tilde{\lambda }\) is the rate of increase in presence of an Allee effect, described by the additional parameter \(\hat{m}<\bar{m}\). When \(\hat{m}\le 0\), the Allee effect is weak, namely the growth rate is reduced but still positive; when \(\hat{m}> 0\), we are in presence of a strong Allee effect, namely \(\hat{m}\) describes a threshold in the concentration of macrophages allowing their growth.

We introduce the scaled variables and parameters

$$\begin{aligned} \begin{aligned} m=\tilde{m}/\bar{m},\;d = \tilde{d}/\bar{d},\;c&=\dfrac{\alpha }{b\bar{m}}\,\tilde{c},\;t=\lambda \bar{m}T,\;{\textbf{x}}=\sqrt{\dfrac{\lambda \bar{m}}{D}}\,{\textbf{X}},\\ \chi =\dfrac{\bar{b}}{\alpha D}\,\psi ,\;\tau =\dfrac{\lambda \bar{m}}{\alpha }\,\nu ,\;\epsilon&=\dfrac{\lambda \bar{m}}{\alpha D}\,\varepsilon ,\;\beta =b/\bar{b},\;r=\kappa /\lambda ,\;\delta =\dfrac{\bar{d}}{\bar{m}\bar{b}}\,\mu , \end{aligned} \end{aligned}$$

where \(\bar{b}\) is a typical production rate for macrophages.

The governing equations can be rewritten in the nondimensional form

$$\begin{aligned} \begin{aligned}&\dfrac{\partial m}{\partial t}\,=\,\Delta _{\textbf{x}}m+G_j(m)-\nabla _{\textbf{x}}\cdot (\Phi (m)\nabla _{\textbf{x}}c),&\text {on}\,\Omega \times {\mathbb {R}}^+,\\&\dfrac{\partial c}{\partial t}\,=\,\dfrac{1}{\tau }\,[\epsilon \Delta _{\textbf{x}}c+(\delta d-c+\beta m)],&\text {on}\,\Omega \times {\mathbb {R}}^+,\\&\dfrac{\partial d}{\partial t}\,=\,r F(m)m(1-d),&\text {on}\,\Omega \times {\mathbb {R}}^+,\\&\dfrac{\partial m}{\partial n}=\dfrac{\partial c}{\partial n}=\dfrac{\partial d}{\partial n}=0,&\text {on}\,\partial \Omega \times {\mathbb {R}}^+,\\&m({\textbf{x}},0)=m_{in}({\textbf{x}}),\, c({\textbf{x}},0)=c_{in}({\textbf{x}}),\, d({\textbf{x}},0)=d_{in}({\textbf{x}}),\,&\text {on}\,\Omega , \end{aligned} \end{aligned}$$
(1)

with \(\chi ,\tau ,\epsilon ,r>0\) and \(\beta ,\delta \ge 0\); the production term \(G_j,\, j=1,2\) for logistic and cubic growth function has the following form

$$\begin{aligned} G_1(m)=m(1-m)\;\text { and }\;G_2(m)=\Lambda m(1-m)(m-M), \end{aligned}$$

respectively, where \(\Lambda \) is a proper nondimensional production rate and \(M=\hat{m}/\bar{m}<1\).

The rigorous results proven in [11, 12] for the model in [7] concerning the global existence of weak and strong solutions, uniform bounds in time for such solutions, also apply in this case.

3 Turing instability

In this section, we analyse the linear stability of system (1) by focusing on the formation of stationary patterns due to the destabilisation of a homogeneous steady state. First, we look for equilibria of the homogeneous system

$$\begin{aligned} \begin{aligned} \dfrac{\partial m}{\partial t}\,&=\,G_j(m),\\ \dfrac{\partial c}{\partial t}\,&=\,\dfrac{1}{\tau }(\delta d-c+\beta m),\\ \dfrac{\partial d}{\partial t}\,&=\,r F(m)m(1-d). \end{aligned} \end{aligned}$$

We observe that the line \(\{(0,\delta d,d), \,d\in {\mathbb {R}}^+\}\) of steady states characterized by the absence of macrophages is stable for both logistic growth and a weak Allee effect, and unstable for a strong Allee effect. The coexistence state \(P_*=(m_*,c_*,d_*)=(1,\beta +\delta ,1)\) is a stable equilibrium in all cases (logistic growth, weak and strong Allee effect). When a strong Allee effect is considered, we have an additional unstable steady state \(P_\sharp =(m_\sharp ,c_\sharp ,d_\sharp )=(M,\beta M+\delta ,1)\).

Now we investigate the conditions leading to the destabilisation of \(P_*\) in presence of spatial diffusion and the chemotactic term.

We denote by \(\varvec{\bar{w}}=(w_m,w_c,w_d)=(m-m_*,c-c_*,d-d_*)\) the perturbation with respect to the homogeneous state and we consider for it the following linearised problem

$$\begin{aligned} \dfrac{\partial \varvec{w}}{\partial t}\,=\,{\textbf{J}}^\prime \varvec{\bar{w}}+{\textbf{D}}^\prime \Delta _{\textbf{x}}\varvec{\bar{w}}, \end{aligned}$$

where, for \(j=1,2\) we have

$$\begin{aligned} {\textbf{J}}^\prime \,=\, \begin{pmatrix} -[\Lambda (1-M)]^{j-1} &{} 0 &{} 0\\ \\ \dfrac{\beta }{\tau } &{} -\dfrac{1}{\tau } &{} \dfrac{\delta }{\tau }\\ \\ 0 &{} 0 &{} -\dfrac{r}{2} \end{pmatrix} ,\quad {\textbf{D}}^\prime \,=\, \begin{pmatrix} 1 &{} -\dfrac{\chi }{2} &{} 0\\ \\ 0 &{} \dfrac{\epsilon }{\tau } &{} 0\\ \\ 0 &{} 0 &{} 0 \end{pmatrix}. \end{aligned}$$
(2)

From now on, we assume that \(\Omega \) is a rectangular domain in \({\mathbb {R}}^n\); we look for solutions of the form \(\varvec{\bar{w}}\,\propto \text {exp}(\lambda t+i {\textbf{k}}\cdot {\textbf{x}})\) and compute the eigenvalues \(\lambda \). One is easily given by \(-r/2\); the others are given as proper functions of \(k=|{\textbf{k}}|\), by solving the dispersion relation

$$\begin{aligned} \det (\lambda {\textbf{I}}-{\textbf{J}}+k^2{\textbf{D}})\,=\,0, \end{aligned}$$

where

$$\begin{aligned} {\textbf{J}}= \begin{pmatrix} -[\Lambda (1-M)]^{j-1} &{} 0\\ \\ \dfrac{\beta }{\tau } &{} -\dfrac{1}{\tau } \end{pmatrix}\,\text { and }\, {\textbf{D}}={\textbf{D}}(\chi )= \begin{pmatrix} 1 &{} -\dfrac{\chi }{2}\\ \\ 0 &{} \dfrac{\epsilon }{\tau } \end{pmatrix}. \end{aligned}$$
(3)

We obtain

$$\begin{aligned} \lambda ^2+g(k^2)\lambda +h(k^2)\,=\,0, \end{aligned}$$

where

$$\begin{aligned} g(k^2)\,=\,k^2\text {tr}({\textbf{D}})-\text {tr}({\textbf{J}}),\quad h(k^2)\,=\,\det ({\textbf{D}})k^4+qk^2+\det ({\textbf{J}}), \end{aligned}$$

being

$$\begin{aligned} q\,=\,\dfrac{2(1+\xi ^2)-\chi \beta }{2\tau }, \quad \xi =\sqrt{[\Lambda (1-M)]^{j-1}\epsilon }. \end{aligned}$$

The Turing instability is guaranteed by a positive eigenvalue for at least one \(k>0\). Since for all the values of the parameters

$$\begin{aligned} g(k^2)\,=\,k^2(1+\epsilon /\tau )+[\Lambda (1-M)]^{j-1}+1/\tau \,>\,0, \end{aligned}$$

the presence of a positive eigenvalue is guaranteed by \(h<0\) for some k. Being h a second-order polynomial in \(k^2\) with positive highest and lowest order coefficients, we require the minimum of h to be negative. The minimum is reached at

$$\begin{aligned} k^2_c\,=\,-\dfrac{q}{2\det ({\textbf{D}})}\,=\,-\dfrac{2(1+\xi ^2)-\chi \beta }{4\epsilon }. \end{aligned}$$
(4)

For this to be positive, it is required

$$\begin{aligned} 2(1+\xi ^2)-\chi \beta \,<\,0. \end{aligned}$$

Here we see that it is not diffusion but the chemotactic term the key ingredient destabilising the homogeneous steady state. We consider \(\chi \) as the bifurcation parameter, obtaining the threshold value

$$\begin{aligned} \chi \,>\,\bar{\chi }\,=\,\dfrac{2(1+\xi ^2)}{\beta }. \end{aligned}$$

We look for a critical value for this parameter by imposing

$$\begin{aligned} h(k_c^2)\,=\,0, \end{aligned}$$

and we find the threshold for the bifurcation parameter \(\chi \)

$$\begin{aligned} \chi _c\,=\,\dfrac{2(1+\xi )^2}{\beta }\,>\,\bar{\chi } \end{aligned}$$
(5)

and hence, substituting \(\chi _c\) in (4), we have

$$\begin{aligned} k_c^2\,=\,\dfrac{\xi }{\epsilon }=\sqrt{\dfrac{[\Lambda (1-M)]^{j-1}}{\epsilon }}. \end{aligned}$$
(6)

Note that \(\chi _c\) decreases as M or \(\beta \) increases, and it increases as \(\epsilon \) or \(\Lambda \) increases.

Bifurcation points may occur for \(\chi >\chi _c\). More precisely, the mode associated to the eigenvalues \(k^2\) of the Laplacian destabilises at \(\chi =\chi _{A,k^2}\), where

$$\begin{aligned} \chi _{A,k^2}=\dfrac{\epsilon k^4+\tau \left( 2+\epsilon [\Lambda (1-M)]^{j-1}\right) k^2+[\Lambda (1-M)]^{j-1}}{\beta k^2}, \end{aligned}$$
(7)

obtained by solving \(h(k^2)=0\) with respect to \(\chi \). Note that, in the logistic case, this reduces to

$$\begin{aligned} \chi _{L,k^2}=\dfrac{\epsilon k^4+\tau \left( 2+\epsilon \right) k^2+1}{\beta k^2}. \end{aligned}$$

4 Weakly nonlinear analysis

We are now interested in deriving the Stuart–Landau [7, 24] equation for the amplitude of spatially periodic solutions of system (1) in the one-dimensional case (\(n=1\)) in the domain \(\Omega =[0,\bar{L}]\). We use a multiple-scale method for performing a weakly nonlinear analysis around the stationary state \(P_*\). Since there is no diffusion in the equation for d, we shall consider the reduced system in the variables \(\varvec{w}=(w_m,w_c)\). Anyway, repeating the same procedure as follows for the whole set of variables \(\varvec{w}=(w_m,w_c,w_d)\), it is easy to show that this gives a null solution for \(w_d\). We rewrite the reduced system as

$$\begin{aligned} \dfrac{\partial \varvec{w}}{\partial t}\,=\,{\mathcal {L}}_\chi \varvec{w}+{\mathcal {N}}\varvec{w}, \end{aligned}$$
(8)

where the linear operator is defined as \({\mathcal {L}}_\chi ={\textbf{J}}+{\textbf{D}}\partial _{xx}\) (\({\textbf{J}}\) and \({\textbf{D}}\) are given in (3)) while \({\mathcal {N}}\) keeps track of nonlinear contributions.

We introduce a small parameter \(\eta ^2=(\chi -\chi _c)/\chi _c\) measuring the distance of the bifurcation parameter from its critical value; we consider the following expansions

$$\begin{aligned} \begin{aligned} \varvec{w}&=\eta \varvec{w}_1+\eta ^2\varvec{w}_2+\eta ^3\varvec{w}_3+O(\eta ^4),\\ \chi&=\chi _c+\eta ^2\chi _2+O(\eta ^4), \end{aligned} \end{aligned}$$

where \(\varvec{w}_i=(w_m^{(i)},w_c^{(i)}), \, i=1,2,3\); for what concerns the time dependence of the solutions, we assume a multiple scale dependence \(\varvec{w}=\varvec{w}(T_2,T_4,\ldots )\), where \(T_{2i}=\eta ^{2i} t, \, i=1,2\), so that the time derivative can be expanded as follows

$$\begin{aligned} \dfrac{\partial }{\partial t}=\eta ^2\dfrac{\partial }{\partial T_2}+\eta ^4\dfrac{\partial }{\partial T_4}+O(\eta ^5). \end{aligned}$$

We can substitute the previous expansions in (8) and collect the terms of the same order.

At the first order of accuracy, we have to solve

$$\begin{aligned} {\mathcal {L}}_{\chi _c}\varvec{w}_1=\varvec{0}; \end{aligned}$$
(9)

by imposing no flux boundary conditions, we look for solutions of the form

$$\begin{aligned} w_\ell ^{(1)}\,=\,\rho _\ell A(T_2,\ldots )\cos (k_cx),\qquad \ell =m,c. \end{aligned}$$
(10)

System (9) admits non trivial solutions; in particular, the following condition holds

$$\begin{aligned} \rho _m=\dfrac{1+\xi }{\beta }\,\rho _c. \end{aligned}$$

At the second order of accuracy, we have to solve a non-homogeneous system

$$\begin{aligned} {\mathcal {L}}_{\chi _c}\varvec{w}_2=\varvec{F}(\varvec{w}_1), \end{aligned}$$
(11)

where \(\varvec{F}\) contains only constants and terms proportional to \(\cos (2 k_c x)\). For the Fredholm alternative, the solvability condition for (11) is \(<\varvec{F},\varvec{\psi }>=0\), where \(\varvec{\psi }\in \ker ({\mathcal {L}}^*)\) and \({\mathcal {L}}^*\) is the adjoint operator of \({\mathcal {L}}_{\chi _c}\), and it is automatically satisfied, since \(\varvec{\psi }=(\psi _1\cos (k_cx),\psi _2\cos (k_cx))\) with

$$\begin{aligned} \psi _2\,=\,\dfrac{\tau }{\beta \epsilon }\,\xi (1+\xi )\psi _1. \end{aligned}$$

The solution can be easily written as

$$\begin{aligned} w^{(2)}_\ell \,=\,A^2\left[ \mu _\ell +\theta _\ell \cos (2k_cx)\right] ,\qquad \ell =m,c, \end{aligned}$$
(12)

with

$$\begin{aligned} \begin{aligned} \mu _m\,&=\,-\dfrac{\rho _m^2}{2}\,\left( 1+\omega \right) ,&\qquad \theta _m\,&=\,-\dfrac{\rho _m^2}{18}\left( \omega -\dfrac{1}{\xi }\right) (1+4\xi )&\\ \mu _c\,&=\,-\beta \,\dfrac{\rho _m^2}{2}\,\left( 1+\omega \right) ,&\qquad \theta _c\,&=\,-\beta \,\dfrac{\rho _m^2}{18}\left( \omega -\dfrac{1}{\xi }\right)&\end{aligned} \end{aligned}$$

and \(\omega =(j-1)/(1-M)\). At the third order, we have to solve

$$\begin{aligned} {\mathcal {L}}_{\chi _c}\varvec{w}_3=\varvec{G}(\varvec{w}_1,\varvec{w}_2), \end{aligned}$$
(13)

where the right hand side can be rewritten in the form

$$\begin{aligned} \varvec{G}(\varvec{w}_1,\varvec{w}_2)=\left[ {\textbf{G}}^0\dfrac{dA}{dT_2}+{\textbf{G}}^1A+{\textbf{G}}^2A^3\right] \cos (k_cx)+{\textbf{G}}^*, \end{aligned}$$
(14)

with

$$\begin{aligned} \begin{aligned} {\textbf{G}}^0=(G_m^0,G_c^0)&=\rho _m\left( 1,\dfrac{\beta }{1+\xi }\right) \\ {\textbf{G}}^1=(G_m^1,G_c^1)&=\rho _m\left( -\dfrac{\chi _2}{2}\dfrac{\beta }{\epsilon }\dfrac{\xi }{1+\xi },0\right) \\ {\textbf{G}}^2=(G_m^2,G_c^2)&=\Bigg (2\dfrac{\xi ^2}{\epsilon }\left( 1+\omega \right) \rho _m\left( \mu _m+\dfrac{\theta _m}{2}\right) +\dfrac{3}{4\epsilon }\omega \xi ^2\rho _m^3\\&\quad -\dfrac{1}{4}\chi _ck_c^2\left[ \rho _m\theta _c+\rho _c\left( \mu _m-\dfrac{\theta _m}{2}\right) -\dfrac{1}{8}\rho _m^2\rho _c\right] ,0\Bigg ) \end{aligned} \end{aligned}$$

and \({\textbf{G}}^*\propto \cos (3k_cx)\) satisfies the Fredholm condition. By imposing the solvability condition at this order, we obtain the Stuart–Landau equation for the amplitude

$$\begin{aligned} \dfrac{dA}{dT_2}\,=\,\sigma A-LA^3, \end{aligned}$$

where

$$\begin{aligned} \sigma \,=\,-\dfrac{G^1_m\psi _1+G^1_c\psi _2}{G^0_m\psi _1+G^0_c\psi _2}\,\text { and }\,L\,=\,\dfrac{G^2_m\psi _1+G^2_c\psi _2}{G^0_m\psi _1+G^0_c\psi _2}. \end{aligned}$$

We can observe that

$$\begin{aligned} \sigma =\dfrac{1}{2}\beta \chi _2\dfrac{\xi }{(1+\xi )(\epsilon +\tau \xi )}>0; \end{aligned}$$

in fact, we consider the perturbation \(\chi _2>0\) because Turing instability may occur only for \(\chi >\chi _c\). For what concerns L, it can be rewritten as

$$\begin{aligned} L=-\dfrac{\rho _m^2}{144(\epsilon +\tau \xi )}p(\xi ,\omega ), \end{aligned}$$

and its sign depends on the sign of the polynomial

$$\begin{aligned} p(\xi ,\omega )=4\omega (8\omega +9)\xi ^3+(152\omega ^2+122\omega +63)\xi ^2-(46\omega +55)\xi +2. \end{aligned}$$

In absence of Allee effect, i.e. when we consider the logistic function \(G_1(m)\), we have

$$\begin{aligned} \xi =\sqrt{\epsilon },\;\omega =0, \end{aligned}$$

and the polynomial p reduces to

$$\begin{aligned} p(\sqrt{\epsilon })=63\epsilon -55\sqrt{\epsilon }+2; \end{aligned}$$

therefore, we obtain that

$$\begin{aligned} L>0\,\Longleftrightarrow \, p< 0\,\Longleftrightarrow \,\epsilon \in (0.0014,0.6972). \end{aligned}$$

In presence of Allee effect, namely when we consider the cubic function \(G_2(m)=\Lambda m(1-m)(m-M)\), we have that the variables \(\xi \) and \(\omega \) depend on the choice of parameters M and \(\Lambda \). In order to contrast and compare the outcomes of different growth terms \(G_j, \,j=1,2\) varying these two parameters, we identified two different ways to set the parameter \(\Lambda \). We explore the following cases, whose differences are shown in Fig. 1:

  • Case 1: for a fixed value of \(M<1\), we set the parameter \(\Lambda =\Lambda _1:=(1-M)^{-1}\) in order to reproduce the bifurcation values of the logistic case, namely \(\chi _{A,k^2}=\chi _{L,k^2}\) (see equation (7)).

  • Case 2: for a fixed value of \(M<1\), we set the parameter \(\Lambda \) in order to reproduce the maximal growth rate of the logistic function

    $$\begin{aligned} \max _{0\le m\le 1}G_2(m)=\max _{0\le m\le 1}G_1(m)=\frac{1}{4}. \end{aligned}$$

    This leads to the value

    $$\begin{aligned} \Lambda _2=\left( 4m_2(1-m_2)(m_2-M)\right) ^{-1}, \quad m_2=\dfrac{(M+1)+\sqrt{M^2-M+1}}{3} >0 . \end{aligned}$$

    By simple algebra, one gets that for all \(M<1\), the value \(m_2 \in (M,1)\), and \(m_2 =1\) if and only if \(M=1\); therefore, \(\Lambda _2>0\) for all \(M<1\). Moreover, taking into account that \(4x(1-x)\le 1 \), it follows that

    $$\begin{aligned} \Lambda _2(1-M)>1,\quad \forall M<1, \end{aligned}$$

    and hence, for a fixed \(k^2\), the bifurcation values \(\chi _{A,k^2}\) are greater than the corresponding values for the logistic case, namely \(\chi _{A,k^2}>\chi _{L,k^2}\) (see equation (7)).

The regions corresponding to supercritical (\(L>0\)) and subcritical (\(L<0\)) bifurcation at \(\chi _c\) are reported in Fig. 2  varying \(\epsilon \) for the logistic growth function and in the \((M,\epsilon )\)-plane when the Allee effect is considered, in both cases \(\Lambda _1\) and \(\Lambda _2\).

Fig. 1
figure 1

Comparison of the profile of logistic growth function \(G_1\) (grey line) with cubic growth functions (black lines) with parameter \(\Lambda \) as in Case 1 (\(G_2(m;\Lambda _1)\), black solid line) and in Case 2 (\(G_2(m;\Lambda _2)\), dashed-dotted line). The left panel corresponds to the weak Allee effect (\(M=-0.5\)), while the right panel to a strong Allee effect (\(M=0.02\)). (For interpretation of the color references, the reader is referred to the online version of this article)

Full size image
Fig. 2
figure 2

Regions in the \((M,\epsilon )\)-plane corresponding to the subcritical case (\(L<0\), grey region) and to the supercritical case (\(L>0\), black region). The value of \(\Lambda \) is set as in Case 1 (left) and as in Case 2 (right). The other parameter values are \(\tau =\beta =r=\delta =1\). The thick line on the left corresponds to the logistic case (independent of M). (For interpretation of the color references, the reader is referred to the online version of this article)

Full size image

Both cases show a significant reduction of the range of parameter \(\epsilon \) corresponding to the supercritical case (\(L>0\)), and hence to the existence of a stable equilibrium for the amplitude, especially when we consider a strong Allee effect. Moreover, by confining on realistic values of \(\epsilon \), namely \(0.5<\epsilon <1.5\) found in [7], we can observe only a subcritical region in presence of Allee effect; the third order equation cannot capture the amplitude of the pattern and the analysis should be pushed to higher order (\(O(\eta ^5)\)).

5 Bifurcation analysis and numerical simulations

In this Section, we investigate the influence of the Allee effect on stationary steady-state solutions, by combining numerical simulations and bifurcation diagrams. The numerical simulations are performed discretising system (1) via finite differences in space and a first-order explicit method in time, starting from a small random perturbation of the homogeneous steady state. For the bifurcation diagram, we exploited the continuation software pde2path. Thanks to the results of the weakly nonlinear analysis, we can perform the numerical continuation on the reduced system with only two equations for m and c, obtained from system (1), and considering the homogeneous solution \(d=d_*\). In fact, the ODE part is stable and the remaining eigenvalues are solely determined by the PDEs part, i.e., the instability only arises from the PDEs subsystem. In the following bifurcation diagrams, thick lines denote stable solutions while thin lines are unstable ones. The homogeneous branch is shown in black, while branches bifurcating from the homogeneous one having the same number of peaks are indicated with the same colours when changing the parameters, to facilitate the comparison (for interpretation of the colour references, the reader is referred to the web version of this article). The continuation software pde2path runs on Matlab and the scripts needed for the numerical continuation of system (1) are freely available in the GitHub folder [33]. In the following, we want to explore the influence of the Allee effect on the outcomes of the model, thus we vary parameters M and \(\Lambda \). Parameter \(\chi \) is taken as the bifurcation parameter, while we select two values for parameter \(\epsilon \), namely \(\epsilon =0.08\) and \(\epsilon =0.8\). The former corresponds to a case in which the weakly nonlinear analysis up to the third order is able to predict the amplitude of the pattern close to the critical value \(\chi _c\), while the latter is in the range of realistic values. All the other parameter values are taken as in [7] and reported here for convenience

$$\begin{aligned} \tau =1,\, \beta =1,\, r=1,\, \delta =1, \, \bar{L}=12\pi . \end{aligned}$$
(15)

We first want to highlight the influence of the Allee effect with respect to the case of logistic growth investigated in [7]. We fix \(\epsilon =0.08\) and we vary the Allee parameter M. The parameter \(\Lambda \) is taken as in Case 1, since this choice keeps the bifurcation points on the homogeneous branch fixed, allowing a direct comparison for selected values of the bifurcation parameter \(\chi \). The results are shown in Fig. 3. In the left panel, the bifurcation diagrams with respect to the parameter \(\chi \) are shown in the case of logistic growth (no Allee effect, for comparison) and with weak (\(M=-0.5\) and \(M=-0.1\)) and strong (\(M=0.02\)) Allee effect. The homogeneous branch is shown in black, while branches bifurcating from the homogeneous corresponding to the same wavenumber (leading to a specific number of peaks in the solution profile) are indicated with the same colours when changing the parameters, to facilitate the comparison. Dotted vertical lines at \(\chi =3.5\) are meant for reference and indicate the value used for the numerical simulations (right panel). We notice that the bifurcation points on the homogeneous branch remain fixed, as expected. We observe major changes induced by an increase of parameter M:

  • Loss of stability of the branches. With an Allee effect and also increasing the parameter M, we observe a progressive loss of stability of the branches. For instance, the magenta branch is stable with a logistic growth but not with an Allee effect, and also the red and blue branches become unstable.

  • More prominent bend for subcritically bifurcating branches. The branches bifurcate supercritically both with logistic growth and a weak Allee effect with \(M=-0.5\) (qualitatively similar comparing the growth functions, see Fig. 1). In particular, the first (blue) branch is stable for \(\chi >\chi _{k^2}\), with \(k^2=23\) (and corresponding solution presenting 11.5 peaks/bumps) while some of the others stabilise eventually. Increasing M, several branches become subcritical and present a wider stable region, namely for lower values of \(\chi \). This is clearly visible for \(M=0.02\).

  • Stabilisation of spatial-periodic solutions with lower number of peaks. For instance, setting \(\chi =3.5\) and performing numerical simulations starting from a random perturbation of the homogeneous state, the numerical solution shows 11 peaks (blue branch) without Allee effect, 11.5 peaks (red branch) with \(M=-0.5\), 10 peaks (cyan branch) with \(M=-0.1\), and 9.5 peaks (orange branch) with \(M=0.02\).

Fig. 3
figure 3

Influence of the parameter M on the bifurcation diagram of steady-state solutions, obtained considering \(\epsilon =0.08\) and \(\Lambda \) as in Case 1. Dotted lines serve as references and denote the values for which numerical simulations (right panel) were performed. (For interpretation of the color references, the reader is referred to the online version of this article)

Full size image

We now focus on the parameter \(\Lambda \) and in particular on the difference between the two possible ways of defining it (Case 1 and Case 2 in Sect. 4). Remember that the choice of \(\Lambda \) as in Case 2 allows maintaining the maximum growth rate as in the logistic growth (while in Case 1 we maintain the bifurcation point fixed, drastically reducing the maximum growth rate, see Fig. 1). To investigate the influence of this choice on the bifurcation diagram and steady-state solutions, we consider \(\epsilon =0.8\), corresponding to a realistic value for this parameter, and \(M=-0.5\) (all the other parameter values are listed in (15), except the bifurcation parameter \(\chi \) which is specified in the text for the numerical simulations). Results are illustrated in Fig. 4. In the bifurcation diagrams (left panel), the homogeneous branch is denoted in black and all the blue branches are unstable. Among the explored bifurcation points (with \(\chi _{k^2}<12\)), only two lead to regions of stable solutions and they are denoted in magenta and red. The magenta branch corresponds to solutions with 5.5 peaks, while the red branch with 5 peaks. The first difference comparing the bifurcation diagrams in Case 1 and Case 2 is, as expected, the position of the bifurcation points on the homogeneous branch: in particular, bifurcation points occur for lower values of \(\chi \) in Case 1. The magenta and red branches are subcritical in both cases, but they are not qualitatively the same. At the two reference values \(\chi =7.5\) and \(\chi =9\), indicated by dotted lines, we have different stable solutions. In particular, at \(\chi =7.5\), only the magenta branch turns out to be stable in Case 2, while in Case 1 only the red one is stable. On the contrary, at \(\chi =9\), both the magenta and the red branch are stable in Case 2, while again in Case 1 only the red one is stable. Moreover, we notice that to have a similar steady state solution (right panel), we have to select two different values of \(\chi \). In fact, we show in the right panel the numerical simulations obtained for \(\chi =9\) in Case 2 and \(\chi =7.5\) in Case 1 and the profile of the steady state solutions. Thus, the choice of \(\Lambda \) seems to have a qualitative and quantitative impact on the possible outcome of the system, which is crucial in applications.

Fig. 4
figure 4

Influence of the choice of parameter \(\Lambda \) as in Case 1 (upper panel) and as in Case 2 (lower panel) on the bifurcation diagram of steady state solutions, obtained considering a realistic value \(\epsilon =0.8\) and a weak Allee effect \(M=-0.5\). In the bifurcation diagrams (left panel), the magenta branch corresponds to solutions with 5.5 peaks, while the red branch with 5 peaks. The homogeneous branch is denoted in black and all the other unstable branches are shown in blue. Dotted lines serve as references and denote the values for which numerical simulations (right panel) were performed. (For interpretation of the color references, the reader is referred to the online version of this article)

Full size image

6 Conclusion and outlook

In this paper, we proposed a modification of the mathematical model describing inflammation and demyelination patterns in the brain caused by Multiple Sclerosis proposed in [7]. In particular, we hypothesized that a minimal amount of macrophages is needed to be able to start and sustain the inflammatory response, while if the concentration of macrophages is too low, then there is no inflammatory response. To model this, we introduced an Allee effect in the model function for macrophage activation. We investigated the emergence and the type of Turing patterns on a 1D domain by combining linearised and weakly nonlinear analysis, bifurcation diagrams and numerical simulations, focusing on the comparison with the previous model [7], formulated using a logistic growth for macrophage activation. We selected two cases for this comparison, one in which the bifurcation points on the homogeneous branch are the same in both the logistic and Allee case, the other in which the maximum macrophage activation rates are the same (even if it is attained for different values, see Fig. 1).

When an Allee effect is considered, we observed clear trends as the Allee parameter M increases. In detail, branches which are stable in the logistic case become unstable, leading to very few stable branches in the bifurcation diagrams. We also observed that the Allee effect leads more often to subcritical bifurcations or to a more prominent bend. Finally, we also observed the stabilisation of branches bifurcating for greater values of the bifurcation parameter (very far from the first bifurcation point and the critical value), leading to stable solutions with a lower number of peaks.

We plan to further investigate this model from several points of view.

In cases when the Landau coefficient L is negative and thus the third order Stuart–Landau equation is not able to capture the amplitude of the pattern, one should push the weakly nonlinear analysis to a higher order, obtaining a quintic Stuart–Landau equation for the amplitude of the pattern as in [7].

Another phenomenon worth to be investigated is the wavefront invasion, namely the existence of travelling wavefronts connecting two different steady solutions of the equations; this would require a weakly nonlinear analysis similar to that performed in Sect. 4, but with an additional expansion of the derivative with respect to the space variable, in order to take into account the slow modulation in space of pattern amplitude [6].

The analysis and the numerical simulations in this paper are performed on a one-dimensional domain. We plan to extend the results to two-dimensional domains, in order to check the formation of different types of spatial patterns in presence of Allee effect varying the parameters and to compare the results with the ones in [7] for the model without Allee effect. On the other hand, the presence of time-periodic spatial patterns was already noted in [7]. The effect of the Allee parameters on these patterns will also be a matter of future work.

Moreover, since we observed that the Allee effect sharpens the subcritical bifurcations, the bifurcation structure might show snaking branches leading to localised patterns. In particular, snaking branches might arise from secondary bifurcation points on branches bifurcating highly subcritically at \(\chi _{k^2}\gg \chi _c\). To identify such branches, the same weakly nonlinear analysis can be extended at the general bifurcation value \(\chi _{bif}\) (with similar but slightly different computations) and looking at the sign and value of the Landau coefficient L.