Pattern Formation in a Model of Inflammation

1. Introduction

1.1. Biological Background

Inflammation is an important mechanism participating in the immune response to harmful stimuli. Many cellular and molecular pathways are activated during the inflammatory process [1]. Inflammation has some disease-specific features, but the mechanism of the immune response is generic [2]. It starts with the identification of primary stimuli by cell surface patterns recognition receptors (PRRs). PRRs are proteins whose role is to sense a specific type of dangerous stimuli such as DAMPs (Damage-associated molecular pattern molecules) or PAMPs (Pathogen-associated molecular pattern molecules). The recognition process leads to the activation of intra-cellular signalling pathways such as NF-$\kappa$B, MAPK and JAK-STAT pathways which causes gene transcription activation on the cellular level. Therefore, numerous inflammatory agents such as inflammatory cytokines and chemokines are released [3]. Then, immune cells migrate from the bloodstream to the region of inflammation [4]. After the set-up of the immune response, inflammatory agents are mobilized in order to restore homeostasis and to resolve inflammation. In fact, an effective anti-inflammatory process is needed to promote the reduction of inflammation [5]. The sensors of blood-circulating monocytes do not detect stimuli, monocytes are not recruited to the site of inflammation. Conversely, in cases of uncontrolled inflammation, acute inflammation fails to stop tissue damage, leading to chronic inflammatory conditions. This imbalance contributes to the progression of various inflammatory diseases such as cancer [6], atherosclerosis [7] or skin diseases [8].

In the process of immune response, immune cells perform significant and crucial functions. A major step in the immune response is the recruitment of immune cells to the area of inflammation such as T cells, lymphocytes and monocytes. T cells are a type of leukocytes that differentiate into several subtypes such T helper Th1 and Th2 and regulate the maturation of monocytes into macrophages [9]. Th1 cells are considered pro-inflammatory because they secrete INF$\gamma$[10] to activate M1 type macrophages as shown in Figure 1. Classically activated macrophages (M1) produce pro-inflammatory cytokines, promote elimination of affected tissue cells and inhibits cell proliferation. In contrast, alternatively activated macrophages (M2 type) are induced by Il-4 and IL-13 produced by Th2 cells [11]. Unlike M1 macrophages, M2 macrophages promote cell repair and cells proliferation [12]. Similarly to macrophages, inflammatory cytokines can be categorized into two types: pro-inflammatory and anti-inflammatory. Pro-inflammatory cytokines, such as TNF-$\alpha$, IL-12, and IL-6, are produced to enhance inflammation. In contrast, IL-10, IL-4 and TGF-$\beta$ are anti-inflammatory cytokines produced to down-regulate inflammation [13] by inhibiting the production of pro-inflammatory cytokines by macrophages. For example, anti-inflammatory cytokines has an effect on the production of inflammatory cytokines such as IL-1$\beta$, TNF-$\alpha$, and IL-6 [14,15]. As macrophages produce more inflammatory cytokines, additional immune cells become activated causing the recruitment of even more immune cells to the inflammation site.

Inflammation plays a major role in skin diseases since it is essential for defending against infection and injury but can also lead to chronic skin conditions when uncontrolled. For example, acne results from inflammation of hair follicles and sebaceous glands [16]. Another skin disease example is eczema (atopic dermatitis) which causes the skin to become dry and itchy. Eczema is a type of dermatitis that can be caused by an overactive immune system to allergens or irritants. This disease is considered as an itchy inflammatory skin disease that touches one out of four school-aged children [17]. The emergence of such diseases starts by the disruption of the human skin [18]. The defense mechanism of the skin is composed of three layers that are the barrier of the skin, innate immunity and acquired immunity [19]. Human skin functioning is controlled by a complex network of cellular interactions that regulate essential physiological processes, including immune responses and inflammation. Several factors can cause a defect in the epidermal barrier such as genetics, obesity, excessive smoking and alcohol consumption [20,21]. The effect of such risk factors leads to genetic mutations and impairment of the epidermal barrier functions [22]. Epigenetic modification impacts the activity of the transcription of genes resulting to the translation of new mRNA and metylation of DNA. Therefore, they affect the functioning, polarization of cells and the secretion of cytokines [20]. As a result, skin inflammation can arise depending on the presence of immune cells, inflammatory cytokines profiles and the intensity of the immune response. Abnormal pro-inflammatory cytokines production by macrophages and dendritic cells leads to dysfunctional immune responses and emergence of several inflammatory skin diseases [23]. Some inflammatory skin diseases such as Psoriasis are depicted by an increase in the concentration of pro-inflammatory cytokines correlated with a severe chronic inflammation [24]. Atopic Dermatitis is highlighted by an abnormal immune response including enhanced type 2 inflammation correlated by overexpression of inflammatory cytokines (IL-4 and IL-13). Type 2 inflammatory cytokines inhibit the expression of epidermal proteins (FLG) and promote the itch-scratch cycle by activating neurons [25]. Hence, a dysfunction in the inflammatory processes induces loss of epidermal integrity, chronic inflammation of the skin, and increased sensitivity to infections. As a result, many symptoms appear on the skin such as red papules and plaques. Several approaches are considered in order to suppress the inflammatory process particularly the production of pro-inflammatory cytokines. For example, SOCS proteins regulates the signaling of pro-inflammatory cytokines through the modulation of the JAK/STAT pathway [24].

1.2. Modelling of Inflammation

Theoretical modelling of inflammation is tackled by different approaches particularly the study of inflammation initiation due to respiratory infections [26]. Authors present a system of nonlinear ordinary differential equations to study the interactions between macrophages, pro-inflammatory, bacteria and anti-inflammatory cytokines. The study shows that an inflammatory state can set-up depending on the model parameters. According to [26], anti-inflammatory cytokines helps to reduce the intensity of inflammation. Mathematical modelling is a valuable tool to understand the interplay between pro-inflammatory and anti-inflammatory signaling [27]. The study presents a mathematical model to describe the interplay between Tumor Necrosis Factor (TNF) and IL-10 inflammatory cytokines in monocytes. Based on experimental data collection, authors show the role of IL-10 early feedback in switching the inflammation. Several studies describes the propagation of inflammation in tissue as a reaction-diffusion wave [28,29], in particular, in atherosclerosis models. The interplay between macrophages, pro-inflammatory cytokines and oxidized low-density lipoproteins (ox-LDL) triggers the set-up of inflammation.

1.3. Pattern Formation

Spatiotemporal patterns formation is widely manifested in various systems far from the thermodynamic equilibrium and exchanging matter and energy with the environment [30]. A. Turing suggested that organ formation in living organisms occurs due to the instability of a spatially homogeneous state in a system of interacting chemicals (morphogens) distributed over the growing embryo [31]. The non-homogeneous distribution of the concentrations of these chemicals influences the behavior of cells leading to their differentiation and organ formation. There are some experimental confirmations of this mechanism, in particular, for the growth of birds’ feathers [32], for the branching pattern of lungs [33], the formation of fingers [34], and the left-right asymmetry of body organs [35]. Turing instability has been also identified in many non-biological systems [36,37,38,39].

The conditions of the emergence of dissipative structures in a two-component systems can be formulated in terms of short-range activation and long-range inhibition, since one of the variables, referred to as a self-activator, should upregulate its own production at the uniform stationary state, and another variable should be a self-inhibitor, acting in the opposite way. Moreover, the diffusion coefficient of the self-inhibitor should be sufficiently higher than that of the self-activator. It is worth noting, that although this is the general concept for emergence of Turing patterns in the systems with greater number of variables as well [40,41], another possibilities for Turing patterns formation exist in multicomponent systems, that can bypass the requirement of short-range activator and long-range inhibitor [42]. However, this requirement is strictly necessary for two-component systems.

Pattern formation in biomedical models can be determined by different mechanisms. As such, the interaction of inflammation and chemotaxis can result in formation of fatty streaks in atherosclerosis [43]. Turing instability was suggested as the mechanism of pattern formation in Crohn’s disease [44]. In this work we study pattern formation in a generic model of inflammation with short-range activation (inflammation) and long-range inhibition (anti-inflammation) mechanism.

2. Model of Immune Response and Inflammation

Inflammation as a part of the immune response includes several steps common for different diseases. They are based on the interaction of different immune cells (monocytes, macrophages, lymphocytes) and cytokines as described in Figure 1. We formulate in this section a generic model of inflammation and reduce it to a simpler two-equation model in order to perform its stability analysis.

2.1. Generic Model of Inflammation

We formulate the model of inflammation for the concentrations of pro- and anti-inflammatory monocytes, macrophages, lymphocytes and cytokines. Equations for pro-inflammatory and anti-inflammatory monocytes $A_1$ and $A_2$ are considered in the following form:

$$\frac{d{A}_1}{dt}={\lambda }_{P{A}_1}\frac{C_3}{k_{1P}+k_{2P}{C}_3}{A}_1^0-{\lambda }_{C_1{A}_1}\frac{C_1}{k_{C_1}+C_1+k_4{C}_4}{A}_1-d_{A_1}{A}_1,$$

(1)

$$\frac{d{A}_2}{dt}={\lambda }_{P{A}_2}\frac{C_3}{k_{1P}+k_{2P}{C}_3}{A}_2^0-{\lambda }_{C_2{A}_2}\frac{C_2}{k_{C_2}+C_2}{A}_2-d_{A_2}{A}_2.$$

(2)

The constants ${\lambda }_{P{A}_1}$ and ${\lambda }_{P{A}_2}$ characterize the rate of monocytes penetration through the vessel wall, $A_1^0$ and $A_2^0$ are the densities of the corresponding monocytes in blood, $d_{A_1}$ and $d_{A_2}$ is death rate of $A_1$ and $A_2$ monocytes. Cell influx depends on the concentration of inflammatory cytokines $C_3$ with saturation. [45]. Thus, the first term in the right-hand side of this equation represents the influx of monocytes from the blood [46], the second term describes the differentiation of monocytes into macrophages [47,48,49], and the third term corresponds to the death of monocytes.

Next, in the equations for macrophages $M_1$ and $M_2$:

$$\frac{d{M}_1}{dt}={\lambda }_{C_1{A}_1}\frac{C_1}{k_{C_1}+C_1+k_4{C}_4}{A}_1-{\lambda }_{L_{ox}{M}_1}\frac{L_{ox}}{k_{L_{ox}}+L_{ox}}{M}_1-d_{M_1}{M}_1,$$

(3)

$$\frac{d{M}_2}{dt}={\lambda }_{C_2{A}_2}\frac{C_2}{k_{A_2}+C_2}{A}_2-{\lambda }_{L_{ox}{M}_2}\frac{L_{ox}}{k_{L_{ox}}+L_{ox}}{M}_2-d_{M_2}{M}_2,$$

(4)

the first term characterizes the differentiation of monocytes into macrophages [47,48,49], the second term describes transformation of macrophages into foam cells [50] in the case of atherosclerosis but it can be canceled for other diseases, and the third term corresponds to the death of macrophages.

Equations for T-helper cells $T_1$ and $T_2$:

$$\frac{d{T}_1}{dt}={\lambda }_{P{T}_1}\frac{C_3}{k_{1P}+k_{2P}{C}_3}{T}_1^0-d_{T_1}{T}_1,$$

(5)

$$\frac{d{T}_2}{dt}={\lambda }_{P{T}_2}\frac{C_3}{k_{1P}+k_{2P}{C}_3}{T}_2^0-d_{T_2}{T}_2.$$

(6)

model the influx of T-cells from the blood (similar to monocytes) [46] and death of T-cells. Here ${\lambda }_{P{T}_1}$ and ${\lambda }_{P{T}_2}$ are the rates of cell penetration, $T_1^0$ and $T_2^0$ are the densities of $T_1$ and $T_2$ cells in blood, $d_{T_1}$ and $d_{T_2}$ are the death rates of $T_1$ and $T_2$ cells.

Equations for T-helper cytokines $C_1$ (IFN-$\gamma$, IL-1$\beta$) and $C_2$ (IL-4, IL-13):

$$\frac{d{C}_1}{dt}={\lambda }_{C_1}{T}_1-k_{C_1{A}_1}{C}_1{A}_1-d_{C_1}{C}_1,$$

(7)

$$\frac{d{C}_2}{dt}={\lambda }_{C_2}{T}_2-k_{C_2{A}_2}{C}_2{A}_2-d_{C_2}{C}_2.$$

(8)

describe productions of cytokines $C_1$ and $C_2$ by T-lymphocytes with the rates ${\lambda }_{C_1}$ and ${\lambda }_{C_2}$. They promote the differentiation of $A_1$ and $A_2$ monocytes into $M_1$ and $M_2$ macrophages with rate constants $k_{C_1{A}_1}$ and $k_{C_2{A}_2}$ and consumed (internalized) due to this interaction [48]. The third term corresponds to the degradation of cytokines with the rate constants $d_{C_1}$ and $d_{C_2}$.

Equation for pro-inflammatory cytokines (IL-6, IL-12 and TNF-$\alpha$):

$$\frac{d{C}_3}{dt}={\lambda }_{C_4{M}_1}\frac{1}{k_{C_4}+C_4}{M}_1-d_{C_3}{C}_3.$$

(9)

Macrophages $M_1$ produce pro-inflammatory cytokines $C_3$ with the rate constant ${\lambda }_{C_4{M}_1}$ activating endothelial cells [48]. The second term corresponds to the degradation of cytokines with the rate $d_{C_3}$. Anti-inflammatory cytocines $C_4$ downregulate $C_3$ production. Equation for anti-inflammatory cytokines (IL-10 and TGF-$\beta$):

$$\frac{d{C}_4}{dt}={\lambda }_{C_4}{M}_2-d_{C_4}{C}_4.$$

(10)

The $M_2$ macrophages produce anti-inflammatory cytokines $C_4$ with rate ${\lambda }_{C_4}$ [47,48,51]. They prevent cell penetration in the tissue. The second term corresponds to the degradation of these cytokines.

2.2. Reduced inflammation model

In order to study pattern formation in the model of inflammation, we will formally reduce system (1)-(10) to a system of two equations. This reduction is based on the approximations of detailed equilibrium for some variables (fast and slow variables) and on the assumption that some constants are small and the corresponding terms can be neglected. The simplified model provides a starting point for further investigations of the complete model without these assumptions.

Equating zero the right-hand side of equation (3) for ${\lambda }_{L_{ox}{M}_1}=0$, we get

$$M_1=\frac{{\lambda }_{C_1{A}_1}}{d_{M_1}}\frac{C_1{A}_1}{k_{C_1}+C_1+k_4{C}_4}.$$

(11)

Similar, from equation (4),

$$M_2=\frac{{\lambda }_{C_2{A}_2}}{d_{M_2}}\frac{C_2{A}_2}{k_{A_2}+C_2}.$$

(12)

From equations (7), (8) for $k_{C_1{A}_1}=k_{C_2{A}_2}=0$,

$$C_1=\frac{{\lambda }_{C_1}}{d_{C_1}}{T}_1,C_2=\frac{{\lambda }_{C_1}}{d_{C_1}}{T}_2.$$

(13)

Next, approximating P by $C_3$, we get from equations (5) and (6) for $k_{2P}=0$:

$$T_1=\frac{{\lambda }_{P{T}_1}{T}_1^0}{k_{1P}{d}_{T_1}}{C}_3,T_2=\frac{{\lambda }_{P{T}_2}{T}_2^0}{k_{1P}{d}_{T_2}}{C}_3.$$

(14)

As before, we find from equations (1), (2):

$$A_1=\frac{{\lambda }_{P{A}_1}{A}_1^0}{k_{1P}{d}_{A_1}}{C}_3,A_2=\frac{{\lambda }_{P{A}_2}{A}_2^0}{k_{1P}{d}_{A_2}}{C}_3.$$

(15)

Finally, substituting all these relations in equations (9), (10) we get for ${\lambda }_{C_{4}F}=0$:

$$\frac{d{C}_3}{dt}=\frac{a_1{C}_3^2}{1+a_2{C}_3+a_3{C}_4}-{\sigma }_1{C}_3,$$

(16)

$$\frac{d{C}_4}{dt}=\frac{b_1{C}_3^2}{1+b_2{C}_3}-{\sigma }_2{C}_4,$$

(17)

where parameters $a_i$, $b_i$, ${\sigma }_i$ can be expressed through the original parameters. This model will be analyzed in the next section.

3. Pattern Formation

3.1. Analysis of Non-Spatially-Distributed System

Denoting for convenience $A=C_3$, $B=C_4$, we obtain the non-spatially-distributed system, considered in this section:

$$\left\{{\displaystyle \begin{array}{l}\frac{\partial A}{\partial t}=\frac{a_1{A}^2}{1+a_{2}A+a_{3}B}-{\sigma }_{1}A,\\ \frac{\partial B}{\partial t}=\frac{b_1{A}^2}{1+b_{2}A}-{\sigma }_{2}B,\end{array}}\right.$$

(18)

all the parameters of which are positive. The nullclines of the system (18) are

$$\begin{gathered} {\displaystyle \begin{array}{l}\end{array}A=0,B=\frac{1}{a_3}\left(\left(\frac{a_1}{{\sigma }_1}-a_2\right)A-1\right), \\ B=\frac{b_1{A}^2}{{\sigma }_2(1+b_{2}A)}.} \end{gathered}$$

(19)

The first nullcline, referred to as A-nullcline, is a pair of lines, intersecting at $(0,-1/a_3)$. The second nullcline, B-nullcline, is a hyperbola with the asymptotes $A=-1/b_2$, always situated in the left half-plane, and $B=(b_1/b_2{\sigma }_2)A-b_1/b_2^2{\sigma }_2$, always has a positive slope and intersects the ordinate axis in the lower half-plane. B-nullcline goes through the point $(0,0)$, touching the abscissa axis at it, and it has a positive second derivative equal to $2b_1/{\sigma }_2$.

The system (18) has a stationary state $(0,0)$ for all values of parameters. As Figure 2 illustrates, depending on the values of parameters, the system can have one or two stationary states in the first quadrant, where the solutions have physical meaning. Note that the phase vectors are not defined at the lines $A=-1/b_2$ and $B=-(1+a_{2}A)/a_3$ where one of the denominators of system (18) vanish, but they do not belong to the first quadrant.

In order to identify the parameter region, where the system has one or two stationary states in the first quadrant, it is useful first of all to retrieve the conditions, under which the nullclines touch each other in the first quadrant. They can be obtained by consideration of the following system, which respresent the conditions, that the nullclines have the same first derivative at their common point:

$$\left\{{\displaystyle \begin{array}{l}\frac{1}{a_3}\left(\left(\frac{a_1}{{\sigma }_1}-a_2\right)A-1\right)=\frac{b_1{A}^2}{{\sigma }_2(1+b_{2}A)},\\ \frac{1}{a_3}\left(\frac{a_1}{{\sigma }_1}-a_2\right)=\frac{b_{1}A(2+b_{2}A)}{{\sigma }_2{(1+b_{2}A)}^2}.\end{array}}\right.$$

Its straightforward transformation gives the following conditions for the touching of the nullclines in the first quadrant (see example in Figure 2a):

$$\left\{{\displaystyle \begin{array}{l}\frac{a_1}{{\sigma }_1}>a_2+b_2,\\ \frac{a_1}{{\sigma }_1}=2\sqrt{\frac{a_3{b}_1}{{\sigma }_2}}+a_2-b_2.\end{array}}\right.$$

(20)

From the geometrical reasoning it is clear that there are two stationary states in the first quadrant if the slope of the second line of A-nullcline is greater than the slope in the case of touching the B-nullcline, but not greater than the slope of its asymptote (see Figure 2b):

$$\left\{{\displaystyle \begin{array}{l}\frac{a_1}{{\sigma }_1}>a_2+b_2,\\ \frac{a_3{b}_1}{b_2{\sigma }_2}+a_2>\frac{a_1}{{\sigma }_1}>2\sqrt{\frac{a_3{b}_1}{{\sigma }_2}}+a_2-b_2.\end{array}}\right.$$

(21)

The case of one stationary state in the first quadrant, apart from the case of touching of nullclines, described by Eqs. (20), is realised in the situation, when the slope of the second line of A-nullcline is greater than the slope of the asymptote of the B-nullcline (see Figure 2c):

$${\displaystyle \frac{a_1}{{\sigma }_1}>\frac{a_3{b}_1}{b_2{\sigma }_2}+a_2,}$$

(22)

or in a borderline situation, when their slopes are equal, but the second line of A-nullcline lies higher than the asymptote of the B-nullcline (see Figure 2d):

$$\left\{{\displaystyle \begin{array}{l}\frac{a_1}{{\sigma }_1}=\frac{a_3{b}_1}{b_2{\sigma }_2}+a_2,\\ \frac{a_3{b}_1}{b_2^2{\sigma }_2}>1.\end{array}}\right.$$

(23)

The Jacobian matrix at the stationary state $(0,0)$ is as follows:

$$\left(\begin{array}{ll}-{\sigma }_1& 0\\ 0& -{\sigma }_2\end{array}\right),$$

(24)

which indicates that this point is always a stable node. According to the second Poincaré Index Theorem [52], the neighbouring simple stationary state along the B-nullcline can be only a saddle. Therefore, in the cases of one stationary state in the first quadrant, described by Eqs. (22) and (23), this state in always unstable. In the case of two stationary states in the first quadrant, described by Eqs. (21), the state with smaller coordinates is also unstable. The stability of the stationary state with greater coordinate is a less trivial question, and it seems to be dependent on the specific values of parameters. Further we will consider the case, corresponding to Figure 2, fixing the values $a_2=a_3=b_2={\sigma }_1={\sigma }_2=1$, $b_1=4$, and varying the value of $a_1$. In this case, there exist a stable stationary state in the region $4<a_1<5$:

$$\begin{gathered} {\displaystyle A_S\left(a_1\right)=\frac{a_1-2+\sqrt{a_1^2-16}}{2(5-a_1)}, \\ {B}_S\left(a_1\right)=\frac{a_1^2-a_1-8+[a_1-1]\sqrt{a_1^2-16}}{2(5-a_1)},} \end{gathered}$$

(25)

that can be checked by straightforward calculations. This state corresponds to persistent inflammation under strong enough activation of corresponding cytokines.

3.2. Turing Structures in a Spatially-Distributed System

In this section, we consider a spatially-distributed system

$$\left\{{\displaystyle \begin{array}{l}\frac{\partial A}{\partial t}=\frac{a_1{A}^2}{1+A+B}-A+\Delta A,\\ \frac{\partial B}{\partial t}=\frac{4A^2}{1+A}-B+D_B\Delta B,\end{array}}\right.$$

(26)

focusing on the formation of stationary non-uniform structures in it, i.e., the Turing patterns. Their emergence is determined by the system linearized at a corresponding stationary state [31]. For simplicity of calculations, we set here $a_2=a_3=1,b_1=4,b_2=1,{\sigma }_1={\sigma }_2=1$, and consider $a_1$ as a parameter of the problem. In general case of a system

$$\left\{{\displaystyle \begin{array}{l}\frac{\partial \tilde{A}}{\partial t}=a_{11}\tilde{A}+a_{12}\tilde{B}+D_1\Delta \tilde{A},\\ \frac{\partial \tilde{B}}{\partial t}=a_{21}\tilde{A}+a_{22}\tilde{B}+D_2\Delta \tilde{B},\end{array}}\right.$$

where $\tilde{A}$ and $\tilde{B}$ denote the perturbations of the stationary state, the conditions for Turing instability can be formulated in the following way:

$$\left\{\begin{array}{l}{a}_{11}{a}_{22}-a_{12}{a}_{21}>0,\\ a_{11}+a_{22}<0,\\ D_2{a}_{11}+D_1{a}_{22}-2\sqrt{D_1{D}_2(a_{11}{a}_{22}-a_{12}{a}_{21})}>0.\end{array}\right.$$

(27)

The first two inequalities ensure the stability of a corresponding stationary state in the absence of diffusion, while the third condition indicates that in its presence the modes with some wavelengths become unstable, eventually producing a non-uniform stationary pattern. Of note, the formation of Turing structures in the considered system (26) is impossible around the uniform stationary state $(0,0)$ due to the fact, that both variables there act as self-inhibitors, downregulating their own production at the uniform stationary state, which follows from Eq. (24).

As it was mentioned in Section 3.1, in the region $4<a_1<5$ there exists a stationary state with coordinates $A_S>0$, $B_S>0$, which is stable in the absence of diffusion implying the fulfillment of the first two conditions for Turing instability around this state. The third condition, applied to the system (26), turns into the inequality

$$\frac{a_1-1}{a_1}{D}_B-2\sqrt{\frac{-a_1^2+16+[8-a_1]\sqrt{a_1^2-16}}{8a_1}{D}_B}-1>0.$$

Figure 3 depicts the part of the region in the parametric space where this condition is also fulfilled, and the formation of the Turing patterns can take place. Note that the diffusion coefficients should not differ very much, i.e., the value of $D_B=2.64$ already results in the Turing instability for every considered value of $a_1$. This is crucial from the biological point of view, since diffusivities of different cytokines can hardly differ an order of magnitude or more under identical conditions.

Turing patterns formation for the considered system was studied numerically in the two-dimensional square domain. The size of the domain $L=200$ was chosen to be sufficiently large in order for several of the unstable wavelengths to fit it under any of the used parameter sets. As initial conditions we consider a sufficiently strong local perturbation of the form

$$A(x,y)=B(x,y)=10{\left(1+10({(x-L/2)}^2+{(y-L/2)}^2)\right)}^{-1}$$

which corresponds to the local initiation of the plaque formation process. The simulations were performed with the C++ implemented code. Zero flux conditions were considered at the boundaries. The equations were solved using the method of splitting into physical processes, i.e., at each time step, the reaction and diffusion terms were calculated sequentially. The space step was chosen to be sufficiently small to correctly reproduce the emerging structures; the time step was chosen in order to minimize the computation time while maintaining the functionality of the implicit Crank-Nicholson method, used to solve the diffusion equations, and the fourth order Runge-Kutta method, used to solve the reaction equations. It was verified that the refinement of discretization does not qualitatively affect the result.

An important moment is that numerical simulations of the considered system (26) demonstrated unrestrained growth of solutions under sufficiently high values of $D_B$ and $a_1$. This is related to the fact that the model describes only the initial stage of inflammation progression. There are various other anti-inflammatory mechanisms initiated at more advanced stages of inflammation [53,54]. Moreover, there are some simplifying assumptions in the model, such as cell influx is not limited and cytokines production is a linear function of cell concentration without limitation mechanisms. In order to implicitly account for the fact of restrained inflammation growth, the following restrictions were introduced and implemented within the program code:

$$\begin{gathered} {\displaystyle \frac{\partial A}{\partial t}=\left\{\begin{array}{l}{\displaystyle \frac{a_1{A}^2}{1+A+B}-A+\Delta AifA<2A_S\left(a_1\right),}\hfill \\ 0ifA\ge k_{lim}{A}_S\left(a_1\right);\hfill \end{array}\right. \\ \frac{\partial B}{\partial t}=\left\{\begin{array}{c}{\displaystyle \frac{4A^2}{1+A}-B+D_B\Delta BifB<2B_S\left(a_1\right),}\hfill \\ 0ifB\ge k_{lim}{B}_S\left(a_1\right);\hfill \end{array}\right.} \end{gathered}$$

(28)

the functions of $A_S\left(a_1\right)$ and $B_S\left(a_1\right)$ being defined in Eqs. (25). The value of the constant $k_{lim}$ is considered as parameter. We will return to this question in the discussion, and we will prove that the solution is bounded under some conditions on parameters (Appendix).

Figure 4 shows the dynamics of Turing patterns formation under $k_{lim}=2$ for different values of parameters $D_B$ and $a_1$, the most right snapshots displaying visually stabilized patterns. Depending on the parameter values, the considered system can exhibit several qualitatively different ways of self-completion of initial local excitation into Turing patterns. Some aspects of the shown dynamics correspond to the ones, previously reported for other classical nonlinear systems [55,56,57].

Close to the boundary of the Turing instability region under low values of $D_B$, the spot-like structures eventually appear. In this case, as shown in Figure 4a, initial excitation results in the propagation of a circular wave, and the values of model variables in its interior region oscillate around the stationary state $A_S,B_S$. The spots close to the initial excitation arise in result of the rotational symmetry breaking of these oscillatory tails. The positions of the spots, situated further from the initial excitation, are strongly affected by their proximity to the boundaries of the computational domain. Under higher values of $a_1$, as Figure 4b demonstrates, the rotational symmetry breaking does not happen and the resulting patterns represent concentric circles, deformed due to the boundary effects.

Deeper in the Turing instability region, the nature of the patterns changes drastically. Under higher values of $D_B$ and sufficiently low values of $a_1$, this process begins with the elongating of the initial spot in one direction, along which it has a slightly increased ratio of activator to inhibitor. The evolution of the developed snakelike structures, that can tear, corrugate and fuse together, may result in stripes, aligned in one direction, as illustrated by Figure 4c, or in labyrinth patterns, as exemplified by Figure 4d, as well as in various transient patterns. Under sufficiently high values of $a_1$ the initial excitation transforms into thick spot, which develops into Turing pattern via consecutive symmetry breaking, as Figure 4e shows.

Interestingly, the choice of the value of $k_{lim}$ by itself affects the form of the emerging patterns. For example, Figure 5a shows the evolution of the considered system under the values of model parameters, corresponding to Figure 4c, but under greater value of $k_{lim}=3$. In this case, the spot-like structures begin to break apart after the initial elongation, eventually resulting in a domain covered in spots, rather than parallel stripes. Under further increase of $k_{lim}$ for the same values of $D_B$ and $a_1$, the initially formed spot-like structure becomes much more stable and retains its shape during the numerical simulation (at least until $t=35000$). However, such localized structure can be perturbed by external noise. Sufficiently strong noise can result in displacement of the structure, and even stronger noise can lead also to its self-completion. Figure 5b illustrates the corresponding numerical simulation, where the random number in the range $[-0.5;0.5]$ is added to the values of both variables at every grid point at every time step. Note that during the simulation these perturbations once add up to disturbance sufficiently strong to initiate a new localized structure near the left boundary.

4. Discussion

We suggest in this work a generic model of inflammation, which includes the main cell types involved in this process, and pro- and anti-inflammatory cytokines. After some simplifications, we reduce this model to a two-equation model since the analysis of the complete model would be excessively difficult, and it would not give an analytical result.

Inflammatory diseases. Inflammation is a complex physiological process including many cell types, cytokines and chemokines with multiple interactions between them and various positive and negative feedbacks. We can represent these multiple factors as two big groups, pro-inflammatory and anti-inflammatory factors. Pro-inflammatory factors stimulate their own production and the production of anti-inflammatory factors, while the latter downregulate the production of both of them. This kind of interaction together with diffusion (or random motion for cells) can result in a non-homogeneous spatial distribution of these factors leading to the formation of lesions (or plaques). Inflammatory skin diseases manifest spectacular examples of different plaques. Among them, linear, target, multiform and serpiginous plaques, which are also observed in modelling results (Figure 6 and Figure 7).

We hypothesize that the interaction of pro- and anti-inflammatory factors determine the lesion shape, size, and growth rate. Mathematical modelling shows that in the presence of only pro-inflammatory factors, the lesion grows as a single circular plaque, which is not the case for the majority of patients. The presence of different plaque forms (Figure 6 and Figure 7) confirms the hypothesis about the role of anti-inflammatory factors in their formation.

We consider as example skin lesions observed in eczema. Figure 6 shows linear and target lesions and the corresponding patterns obtained in numerical simulations. Figure 7 presents drawings of some other lesions with different structures. Most of them have some similarities with the simulated structures. Cases 1 and 7 correspond to stripe patterns c, case 3 to single dot and case 5 to multiple dot patterns a and b. Case 6 resembles the transient cross pattern e (right), and case 8 the labyrinth pattern d.

Short range activation - long range inhibition. Conditions of the emergence of dissipative structures imply that the diffusion coefficient of inhibitor is greater than the diffusion coefficient of activator. Activating and inhibiting pathways include cells (monocytes, macrophages, foam cells, smooth muscle cells, T-helper cells) and cytokines. Molecular weights of pro-inflammatory cytokines (IL-6 - 26 kDa, IL-12 - 70 kDa, TNF-$\alpha$ - 17.3 kDa) is an average greater that of anti-inflammatory cytokines (IL-10 - 18 kDa, TGF-$\beta$ (active form) - 12 kDa). Since the diffusion coefficient is in inverse proportion to the power $1/3$ of molecular weight, then anti-inflammatory cytokines diffuse faster.

The motion of pro- and anti-inflammatory monocytes and macrophages is similar to each other. However, foam cells, being among the main sources of pro-inflammatory cytokines, do not move. On the other hand, smooth muscle cells migrate to form fibrous cap and to protect the plaque from further influx of monocytes. Though we do not explicitly consider smooth muscle cells in the model, the cytokines which determine their migration, such as platelet-derived growth factor (PDGF) and heparin binding epidermal growth factor (HB-EGF) produced by macrophages, TGF-$\beta$ and so on [58,59], indirectly downregulate the influx of monocytes. In the framework of our model they can be considered as anti-inflammatory factors.

Thus, in spite of the complexity of atherogenesis and numerous factors involved in this processes, faster diffusion of anti-atherogenetic factors seems plausible.

Wave propagation. Reduced inflammatory model developed in this work has from one to three non-negative stationary points. Their number and stability depend on the parameters of the model. Besides stable zero stationary point, the model can have one more stable and unstable points. Therefore, it is a bistable case, and the transition from the disease-free stationary point to a disease stationary point occurs if the initial perturbation is sufficiently large.

If we consider a spatially distributed system with a localized in space initial perturbation of the disease-free equilibrium, then transition to the disease equilibrium occurs in the form of travelling waves. Existence and stability of such waves in the two-equation pro-inflammatory system was studied in [28]. Transition waves for more complete pro- and anti-inflammatory systems will be studied in the subsequent works. Let us only mention here that in the pro-inflammatory model the disease equilibrium is stable. Wave propagation describes growth of a single circular plaque. In the case of pro- and anti-inflammatory model, the disease equilibrium can become unstable. In this case, the solution behind the wave is non-homogeneous in space, and we obtain multiple plaques with different shapes.

Random perturbations of cell and cytokine concentrations specific for physiological systems can be essential at the initial stage of plaque growth. After the initiations, small randomly located lesions can merge and form a single plaque, as it is the case of the process is dominated by pro-inflammatory factors, or they can evolve to a multiple plaque with a complex lesion structure under the influence of anti-inflammatory factors.

Limitations and perspectives. Studying inflammatory diseases, we encounter the usual difficulties of the biomedical modelling: physiological processes are too complex, they contain many unknown or partially known factors, the values of kinetic parameters are basically unknown. In this case, we have to rely on qualitative models which are in agreement with established physiological data. We have to simplify the model grouping them into two classes, pro- and anti-inflammatory.

Some of the kinetic constants are estimated in the literature, some other are considered as free parameters and varied in numerical simulations. Generic inflammatory model, obtained by a formal reduction from a more complete model, contains a relatively small number of parameters. Their values are chosen from the conditions of the emergence of spatial structures. The application of these qualitative ideas and results to specific inflammatory diseases will require further investigations.