A mathematical model of the atherosclerosis development in thin blood vessels and its asymptotic approximation

TL;DR

This paper introduces a new mathematical model for atherosclerosis in thin blood vessels, incorporating recent experimental findings on macrophages, and demonstrates that a simplified 2D model can accurately approximate the complex 3D model.

Contribution

It proposes an improved mathematical model of atherosclerosis considering macrophage roles and provides an asymptotic approximation to simplify the complex 3D model.

Findings

Existence and uniqueness of the model's positive solution are proven.

Asymptotic approximation justifies replacing 3D with a 2D model.

Simplified 2D model maintains sufficient accuracy.

Abstract

Some existing models of the atherosclerosis development are discussed and a new improved mathematical model, which takes into account new experimental results about diverse roles of macrophages in atherosclerosis, is proposed. Using technic of upper and lower solutions, the existence and uniqueness of its positive solution are justified. After the nondimensionalisation, small parameters are found. Then asymptotic approximation for the solution is constructed and justified with the help of asymptotic methods for boundary-value problems in thin domains. The results argue for the possibility to replace the complex $3D$ (dimensional) mathematical model with the corresponding simpler $2D$ model with sufficient accuracy measured by these small parameters.

Equations

$$\left\{\begin{array}[]{rcl}\partial_{t}u_{m}&=&d_{m}\partial^{2}_{xx}u_{m}+f_{1}(u_{c})-\beta u_{m},\\[5.69054pt] \partial_{t}u_{c}&=&d_{c}\partial^{2}_{xx}u_{c}+f_{2}(u_{c})u_{m}-\gamma u_{c}+b,\end{array}\right.$$

$$f_1(u_c)=\frac{{\alpha }_1+{\beta }_1{u}_c}{1+u_c/{\tau }_1},f_2(u_m)=\frac{{\alpha }_2{u}_m}{1+u_m/{\tau }_2}.$$

$$\left\{\begin{array}[]{rcl}\partial_{t}u_{m}&=&d_{m}\Delta_{x_{1}x_{2}}u_{m}-\beta u_{m},\\[5.69054pt] \partial_{t}u_{c}&=&d_{c}\Delta_{x_{1}x_{2}}u_{c}+f_{2}(u_{c})u_{m}-\gamma u_{c}+b,\end{array}\right.$$

$${\Omega }_{\epsilon }=\{(x_1,x_2)\in {\mathbb{R}}^2:x_1\in (0,L),x_2\in (0,\epsilon )\}$$

$$d_m{\partial }_{x_2}{u}_m=\epsilon f_1(u_c)\text{at}{x}_2=\epsilon ,x_1\in (0,L),$$

$${\partial }_{x_2}{u}_m=0,{\partial }_{x_2}{u}_c=0\text{at}{x}_2=0,$$

$${\partial }_{x_2}{u}_m=f_1(u_c),{\partial }_{x_2}{u}_c=0\text{at}{x}_2=h,$$

$$f_2(u)>0\text{for}u>0,f_2(0)=0,f_2(u)\to f_2^{+}\text{as}u\to +\infty ;$$

$$f_1(u)>0\text{for}u>c_0,f_1(c_0)=0,f_1(u)\to f_1^{+}\text{as}u\to +\infty ,$$

$$S+E{\rightleftarrows }_{k_{-1}}^{k_1}SE\stackrel{k_2}{⟶}P+E$$

$$\frac{ds}{dt}=-k_{1}es+k_{-1}c,\frac{de}{dt}=-k_{1}es+(k_{-1}+k_2)c,$$

$$\frac{dc}{dt}=k_{1}es-(k_{-1}+k_2)c,\frac{dp}{dt}=k_{2}c,$$

$$\frac{ds}{dt}=-f(s),\frac{dp}{dt}=f(s),$$

$$f(s)=\frac{\lambda s^p}{\mu +s^p}\text{or}f(s)=\frac{\alpha +\lambda s^p}{\mu +s^p}(p>0,\mu >\frac{\alpha }{\lambda },\lambda >0);$$

$$f(s)=\frac{\lambda g(s)}{\mu +g(s)}(g^{\prime }(s)\ge 0,g(0)=0,\mu >0).$$

$$C_{R,\rho_{0}}:=\Big{\{}x=(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}:\ \ x_{1}\in(0,l),\quad\sqrt{x_{2}^{2}+x_{3}^{2}}\,=:\rho\in(R,R+\rho_{0})\Big{\}}$$

$$\Gamma_{0,l}:=\partial C_{R,\rho}\cap\big{(}\{x:\ x_{1}=0\}\cup\{x:\ x_{1}=l\}\big{)}$$

$${\rm Y}:=\partial C_{R,\rho }\setminus {\Gamma }_{0,l}.$$

$$d_L{\partial }_{\nu }L={\alpha }_0\frac{L}{k_0+L}+{\text{c◯}}_1\text{on}\omega ,$$

$$d_m{\partial }_{\nu }m=\frac{{\alpha }_1{L}_{ox}}{k_1+L_{ox}}+{\text{c◯}}_2\text{on}\Omega \supset \omega ;$$

$$\begin{array}[]{rcl}&&M_{1}\ \ (\text{classically activated;}\ \ 85\%\sim 95\%)\\ &\nearrow&\\ m&&\\ &\searrow&\\ &&M_{2}\ \ (\text{alternatively activated})\end{array}$$

$$d_{M_2}{\partial }_{\nu }{M}_2=-\frac{{\alpha }_2{M}_2^{p_2}}{k_2^{p_2}+M_2^{p_2}}\text{on}\Omega \supset \omega (p_2>1);$$

$${\partial }_{t}L-d_L{\Delta }_{x}L=-\Theta (x,t)L+\psi (x,t),$$

$${\partial }_t{L}_{ox}-d_{L_{ox}}{\Delta }_x{L}_{ox}=\Theta (x,t)L-\underset{intakeof{L}_{ox}by{M}_1}{\underbrace{{\lambda }_{L_{ox}{M}_1}\frac{L_{ox}}{K_{M_1}+L_{ox}}{M}_1}}-\underset{intakeof{L}_{ox}by{M}_2}{\underbrace{{\lambda }_{L_{ox}{M}_2}\frac{L_{ox}}{K_{M_2}+L_{ox}}{M}_2}}$$

$${\partial }_{t}m-d_m{\Delta }_{x}m=-\underset{transformationinto{M}_1}{\underbrace{{\lambda }_{m{M}_1}\frac{m}{k_{M_1}+m}}}-\underset{transformationinto{M}_2}{\underbrace{{\lambda }_{m{M}_2}\frac{m}{k_{M_2}+m}}}\underset{deathofm}{\underbrace{-{\beta }_{m}m}},$$

$${\partial }_t{M}_1-d_{M_1}{\Delta }_x{M}_1={\lambda }_{m{M}_1}\frac{m}{k_{M_1}+m}+{\lambda }_{L_{ox}{M}_1}\frac{L_{ox}}{K_{M_1}+L_{ox}}{M}_1-\underset{transformationintoF}{\underbrace{{\lambda }_{M_{1}F}\frac{M_1}{K_F+M_1}}}-{\beta }_{M_1}{M}_1,$$

$${\partial }_t{M}_2-d_{M_2}{\Delta }_x{M}_2={\lambda }_{m{M}_2}\frac{m}{k_{M_2}+m}+{\lambda }_{L_{ox}{M}_2}\frac{L_{ox}}{K_{M_2}+L_{ox}}{M}_2-{\beta }_{M_2}{M}_2,$$

$${\partial }_{t}F-d_F{\Delta }_{x}F={\lambda }_{M_{1}F}\frac{M_1}{K_F+M_1}-{\beta }_{F}F$$

$$d_m{\partial }_{\nu }m=\frac{{\alpha }_1{L}_{ox}}{k_1+L_{ox}}+\frac{{\alpha }_3{M}_1^{p_1}}{(k_3+M_2)(k_4^{p_1}+M_1^{p_1})}\text{on}\Omega (p_1>1).$$

$$d_L{\partial }_{\nu }L={\alpha }_0\frac{L}{k_0+L}+\frac{{\alpha }_5{F}^{p_3}}{k_5^{p_3}+F^{p_3}}\text{on}\omega (p_3>1).$$

$$\frac{d{u}_1}{dt}=g(u_n)-k_1{u}_1,$$

$$\frac{d{u}_i}{dt}=u_{i-1}-k_i{u}_i,i=2,\dots ,n(k_i>0),$$

$$g(u)=\frac{\nu +\alpha u^p}{K+u^p}(K>0,\nu ,\alpha \ge 0,p>0).$$

Full text

    MOS subject classification:   35B40, 92C20, 35K57, 35K50, 74K10

1. 1 Introduction

2. 2 Existing mathematical models of the atherosclerosis development

3. 3 Improved mathematical model

4. 3.1 Modelling the penetration of LDLs and monocytes

5. 3.2 Modelling the efflux of macrophages

6. 3.3 Processes in the intima

7. 3.4 Modelling of the stage 3 (chronic inflammatory reaction)

8. 3.5 Nondimensionalization

9. 3.6 Existence and uniqueness of the solution to problem (3.15)

10. 4 Asymptotic analysis of problem (3.15)

11. 4.1 Construction of the asymptotic approximation

12. 4.2 Justification

13. 5 Conclusions

1. Introduction

Cardiovascular diseases occupy a leading place in mortality in the world. The main cause of these diseases is atherosclerosis. Therefore the development of atherosclerosis is intensively investigated in last time. There are several theories about atherosclerosis (e.g. see Refs. [21, 13, 37, 42]), but none can explain the whole process because this disease is associated with many risk factors. The atherosclerotic process is not fully understood till now. However, many researchers agree that the damage or dysfunction of the arterial endothelium and high level of low density lipoproteins in blood vessels play the main role in the development of atherosclerosis.

In existing mathematical models, researchers are trying to take into account various factors and different types of molecules involved in the development of this illness. As a result we have models containing from two to more then thirty differential equations. However, even a such cumbersome system cannot take into account all features of atherosclerosis.

Starting from a review of current modeling approaches of atherosclerosis, this article aims at specifying the optimal prospects for research on the mathematical study of atherosclerosis involving special rigorous asymptotic method that enables the reasonable approximation of the original model with its features. The main idea is that atherosclerosis should be looked at as complex system of enzyme reactions that are able to choose their individual dynamics.

The paper is organized as follows. Section 2 describes the basic stages of mechanism of the atherosclerosis development. Moreover, it presents some preliminary models focusing on their shortcomings. Taking into account those shortcomings, Section 3 offers a new mathematical model of atherosclerosis. Using technic of upper-lower solution, the existence and uniqueness of its positive solution are justified. After the nondimensionalisation, we find small parameters of the mathematical model and make the asymptotic analysis as those small parameters tend to zero in Section  4. Namely, we find the corresponding limit problem, construct the asymptotic approximation, find its residuals, estimate them and prove asymptotic estimates for the difference between the solution and the approximating function. In Section 5, we discuss several generalizations and research perspectives.

2. Existing mathematical models of the atherosclerosis development

For convenience of readers we present the following short glossary:

• •

LDLs and HDLs – low and high density lipoproteins; they transport lipid and cholesterol to cells and tissues around the body;

• •

Free radicals – extremely reactive molecules which scavenge electrons causing oxidation of the target;

• •

Monocytes – the largest of the white blood cells (leukocytes); they are able to get in places of inflammation or tissue damage;

• •

Cytokines – small proteins that are important in cell signaling; they are secreted by certain cells of the immune system and have an effect on other cells including the stimulation or inhibition of their growth and other functional activities;

• •

Intima – the first thin layer of the blood vessel wall after the endothelium.

The mechanism of the atherosclerosis development can be shortly sketched as follows (see e.g. Ref. [13]):

Stage 1

1. 1.

After damage of the vessel endothelium, circulating lipids (mostly (LDLs)) begin to accumulate at the site of the injury. 2. 2.

The process of atherosclerosis begins when LDLs penetrate into the intima of the arterial wall where they are oxidized (ox-LDLs).

Stage 2

1. 3.

Ox-LDLs in the arterial intima is considered by the immune system as a dangerous substance, hence an immune response is launched: monocytes circulating in the blood adhere to the endothelium and then they penetrate to the intima. 2. 4.

Once in the intima, these monocytes are converted into macrophages.

Stage 3

1. 5.

The macrophages phagocytose (ingest) the ox-LDLs and become foam cells. Simultaneously chronic inflammatory reaction is started:

• •

macrophages secrete pro-inflammatory cytokines that promote the recruitment of new monocytes and thereby support the production of new pro-inflammatory cytokines;

• •

macrophages with a large intake of the ox-LDLs become foam cells that cause the growth of the intimal layer and thereby amplify the endothelial dysfunction;

• •

this auto-amplification phenomenon is compensated by an anti-inflammatory response mediated by anti-inflammatory cytokines.

This is the short sketch of the atherosclerosis development. In reality, many complex biochemical reactions are hidden behind of these steps. Some of them will be discussed later.

Now let us consider existing mathematical models. The first models of penetration of cholesterol in the arterial wall (the first step in the stage 1) were considered in Refs. [10, 36]. These models are one dimensional systems of differential equations with respect to the radial coordinate (the vessel wall was assumed to be uniform along its length).

In Ref. [26] the authors have developed simple models of the reactions arising in the arterial intima. The first model represents the following one-dimensional reaction-diffusion system:

[TABLE]

for $x\in(0,L).$ Here $\partial_{t}u=\frac{\partial u}{\partial t},$ $\partial^{2}_{xx}u=\frac{\partial^{2}u}{\partial x^{2}},$ the value $u_{m}$ is the concentration of monocytes, macrophages and foam cells together in the intima, $u_{c}$ is the concentration of cytokines. The function $f_{1}(u_{c})$ describes the recruitment of monocytes from the blood flow, $f_{2}(u_{c})u_{m}$ is the rate of production of the cytokines which depends on their concentration and on the concentration of the blood cells:

[TABLE]

The negative terms “$-\beta u_{m}$” and “$-\gamma u_{c}$” correspond to the natural death of the cells and chemical substances, and the last term $b$ in the right-hand side describes the ground level of the cytokines in the intima.

The second model deals with the system

[TABLE]

in the two-dimensional thin rectangle

[TABLE]

(here $\varepsilon$ is a small parameter,  $\Delta_{x_{1}x_{2}}u:=\partial^{2}_{x_{1}x_{1}}u+\partial^{2}_{x_{2}x_{2}}u$) with the nonlinear boundary condition that takes into account the recruitment of monocytes from the blood flow:

[TABLE]

and the homogeneous Neumann boundary condition on the rest part of the boundary $\partial\Omega_{\varepsilon};$ the homogeneous Neumann boundary conditions are imposed everywhere on $\partial\Omega_{\varepsilon}$ for the cytokine concentration $u_{c}$.

The authors have analyzed one-dimensional model depending on the parameters $\alpha_{1},\alpha_{2},\beta_{1},\tau_{1},\tau_{2}$, proved the existence of travelling wave solutions and with the help formal asymptotic analysis (see the appendix B in Ref. [26]) showed that the system (2.1) can obtained from (2.3)–(2.4) as $\varepsilon\to 0.$

In Ref. [27] the authors continued to study two-dimensional system (2.3) in the strip $\Omega_{h}=\{(x_{1},x_{2})\in\mathbb{R}^{2}:\ \ -\infty<x_{1}<+\infty,\ \ 0<x_{2}<h\}$ with the following boundary conditions:

[TABLE]

where the functions $f_{1}$ and $f_{2}$ are sufficiently smooth and satisfy the following conditions:

[TABLE]

and $f_{1}^{\prime}(u)>0.$ The positivity of the solution of the system (2.3)–(2.4) and the existence of travelling waves in the monostable case are proved. Results of numerical simulations for those systems were also obtained in Refs. [26, 27].

These models include only two variables $u_{m}$ and $u_{c}$ (it is essentially for the method proposed by the authors). Obviously, it is not enough to describe the full picture of atherosclerosis. In addition, these models omit some essential features of atherosclerosis development, namely, LDLs penetration and oxidation (see the stage 1);  transformation of monocytes into macrophages and then into foam cells, the chronic inflammatory reaction (see the stage 3);  the diffusion coefficients of macrophages and cytokines are quite different $(8.64\times 10^{-7}$ and $1.08\times 10^{2}$ $cm^{2}\,day^{-1},$ see e.g. Ref. [22]) and cytokines are active in very low concentrations and their secretion is short and strictly regulated (see e.g. Refs. [1, 43]). Also, the size of the damaged endothelium and cylindrical shape of vessels were not taken into account (boundary conditions (2.4) and (4.18) are set on the whole side both of the rectangle and strip).

At first time, all three stages and in addition the formation of a plaque were mathematically simulated in Ref. [22]. The model includes the following key variables: LDLs and HDLs, free radicals and ox-LDLs,  six types of cytokines (MCP-1, IFN-$\gamma$, IL-12, PDGF, MMP and TIMP), and four types of cells (macrophages, foam cells, T-cells and smooth muscle cells). This model consists of eighteen partial differential equations in a plane domain. Unfortunately, the existence and uniqueness of the solution of this system were not justified, the nondimensionalization was not made, monocytes were not involved in the model. In addition, there are many special assumptions, e.g.:  all cells are of the same volume and surface area, so that the diffusion coefficients of the all cells have the same coefficient; the cells are moving with a common velocity,  the domain has polygonal boundaries.

But, even those equations are not enough to describe a complete model of atherosclerosis. For instance, in Refs. [11, 44] the authors showed that more than twenty substances (different cells, various cytokines, hormones, chemical mediators and effectors) are involved in the pathogenesis of atherosclerosis.

Now researchers study different stages of atherosclerosis in more detail. In Ref. [12] the preliminary lipoprorein oxidation model includes nine ordinary equations, the extended model has twelve equations. The paper[20] discusses the central roles of macrophages in different states of atherosclerosis, focusing on the role of inflammatory biomarkers in predicting primary cardiovascular events. There (Section 4) it was showed how many different pro- and anti-inflammatory cytokines, chemokines, mediators, enzymes and biomarkers related only to macrophages are involved in the progression of atherosclerosis.

So, the following question arises. How should we mathematically simulate the atherosclerosis development? Should we collect all those equations together? Then we get a system with about forty equations.

3. Improved mathematical model

It is known that many disease are a set of biochemical reactions. Most of these reactions are enzymatic reactions that take into account the influence of smaller molecules (like cytokines and enzymes). Therefore, the answer on the question above is the balance between biological realism and simplicity which can be described by known scenarios of enzymatic reactions and their mathematical models.

The mechanism of the most basic enzymatic reaction (a substance S reacting with an enzyme E to produce a substance P), first proposed by Michaelis and Menten[31], is represented schematically by

[TABLE]

and described by the system of differential equations (see e.g. Ref. [34] [Sec. 6.1])

[TABLE]

with the initial conditions: $s(0)=s_{0},\ e(0)=e_{0},\ c(0)=0,\ p(0)=0.$ Here $s=[S],$ $e=[E],$ $c=[SE],$ $p=[P]$ are corresponding concentrations, $k$’s are the rate constants.

Taking into account the second and third initial conditions, we can reduce this system to the following two differential equations:

[TABLE]

where $\displaystyle{f(s)=\frac{k_{2}\,e_{0}\,s}{K_{m}+s}}$ and $\displaystyle{K_{m}=\frac{k_{-1}+k_{2}}{k_{1}}}$ is the Michaelis constant.

Thus, we don’t need any differential equation for the enzyme $E$ to describe the enzymatic reaction that converts the substance $S$ to $P.$

Typical choices of such functions in biological applications are

[TABLE]

and more general case

[TABLE]

In the present paper this approach will be used to describe biochemical reactions in the development of atherosclerosis.

From the review above, we see that the following molecules and cells

• •

LDLs $(L),$ ox-LDLs $(L_{ox}),$ monocytes $(m),$ macrophages $(M),$ foam cells $(F)$

qualitatively reproduce the atherosclerosis development. Therefore, we consider their concentrations $L,L_{ox},m,M,F$ as key variables in our model and do not include other molecules with very low concentration.

Now let us go to modelling. At first, we define the following pipe domain (see Fig. 1):

[TABLE]

that will be a prototype of an artery. Its boundary consists of two bases

[TABLE]

and two cylindrical parts

[TABLE]

In the inner cylindrical part \Upsilon_{R}:=\big{\{}x\in\mathbb{R}^{3}:\ \ x_{1}\in(0,l),\ \rho=R\big{\}}, we consider two smooth surfaces $\omega$ and $\Omega$ such that $\omega\subset\Omega\ \subseteq\ \Upsilon_{R}.$ The surface $\omega$ is a prototype of the damage or dysfunction of the arterial endothelium, and $\Omega$ is a surface of the penetration of monocytes into the intima.

3.1. Modelling the penetration of LDLs and monocytes

The penetration of LDLs was simulated by usual linear Robin boundary condition (see equastion (23) in Ref. [22]). However, as has been found in Ref. [13], the damaged endothelium is involved in the activation of endothelial adhesion molecules:  (ICAM-1), (VCAM-1),  P-selestin,  and E-selestin that attract LDLs to endothelial cells. Thus, the penetration of LDLs is not pure diffusion. We admit that the penetration is obeyed to the Michaelis-Menten rule (when enzyme reactions appear at the boundary surface of the diffusion medium it leads to nonlinear boundary conditions [42, 38]. So, we propose the following boundary condition:

[TABLE]

and $\partial_{\nu}L=0$ on $\Upsilon\setminus\omega$ (the non-flux boundary condition);  $\partial_{\nu}L$ is the outward normal derivative. Hereinafter all coefficients in differential equations and boundary conditions are positive (biological meaning of the coefficients $\alpha_{0}$ and $k_{0}$ before (4.1)). The second term “${\copyright_{1}}$” in (3.3) simulates another process and will be presented later.

The key to the early inflammatory response is the activation of the endothelial cells through special enzyme (apoB-LPs), which is intensified by ox-LDLs, that leads to the recruitment of monocytes from the lumen into the intima [33]. Therefore, we model this enzyme reaction by the following boundary condition:

[TABLE]

in addition the non-flux boundary condition on $\Upsilon\setminus\Omega.$ The term “${\copyright_{2}}$” will be presented later.

3.2. Modelling the efflux of macrophages

Recently, there has been a great deal of interest in the role of macrophages in atherosclerotic disease (see e.g. Refs. [24, 30, 39, 19, 20, 14]). It was discovered that monocytes are transformed into two different kind of macrophages $M_{1}$ and $M_{2}:$

[TABLE]

Classically activated macrophages $(M_{1})$ promote the inflammation and product cytokines $IL$-23, $IL$-6, $IL$-12, $IL$-1, $TNF$-$\alpha;$ while the second type macrophages $(M_{2})$ secrete anti-inflammatory cytokines $IL$-4, $IL$-13, $IL$-1, $IL$-10 (see Ref. [20]). In addition, it is turned out that $M_{2}$-macrophages with the ingested ox-LDLs inside can return (efflux) to the blood flow in the vessel [14]. We model this by

[TABLE]

and the non-flux boundary condition on $\Upsilon\setminus\Omega.$

Remark 3.1**.**

Hereinafter, parameters $p_{1},p_{2},p_{3}$ are greater than $1.$ This means some delay in time and less mobility of macrophages and foam cells because of their size.

3.3. Processes in the intima

Evolution and oxidation of LDLs are modelled by the following two reaction-diffusion equations:

[TABLE]

in the domain $C_{R,\rho},$ where $d_{L}$ and $d_{L_{ox}}$ are the diffusion coefficients of LDLs and ox-LDLs in the intima respectively. The function $\psi$ is smooth, nonnegative and its compact support is contained in $\omega$ (its biological meaning in Remark 3.5).

In papers Refs. [22, 12], the term “$k_{L}rL$” is used to simulate the production of ox-LDLs in the reaction of LDLs with free radicals $r$ whose baseline growth is a bifurcation parameter $r_{0}.$ However, as was noted in Refs. [12], it is impossible to measure $r_{0}$ since many social factors such as drinking, smoking and ecological factor highly and unpredictably increase the parameter $r_{0}.$ Therefore, we use a smooth function $\Theta(x,t),$ $(x,t)\in C_{R,\rho_{0}}\times[0,T],$ that can depend both on the concentration of free radicals $r$ and on other values (like HDLs, vitamins $C$ and $E$) and assume that $\Theta|_{\Gamma_{0,l}}=0,$ $\Theta|_{t=0}=0,$ and $0\leq\Theta\leq c_{2}.$

The last two terms in the right-hand side of equation (3.8) are reduction terms that represents the phagocitosis (intake) of ox-LDLs by macrophages $M_{1}$ and $M_{2}$ respectively. By virtue of the diagram (3.5) we regard that the reduction rate $\lambda_{L_{ox}M_{2}}$ is much lower than $\lambda_{L_{ox}M_{1}}$ and $K_{M_{2}}<K_{M_{1}}.$

In Ref. [22]   the authors did not take into account two different types of macrophages and have chosen the term “$\lambda_{L_{ox}M}L_{ox}M$” to simulate the ingestion of ox-LDLs. However, we prefer terms like $\lambda\frac{L_{ox}}{K_{M}+L_{ox}}\,M$ since they appear in many modern prey-predator models and other biochemical reactions [40].

The equation (3.8) is supplemented with the boundary condition $\partial_{\nu}L_{ox}=0$ on $\Upsilon$.

Equations, which model the evolution of “$m$”, “$M_{1}$”, and “$M_{2}$” in the intima, are as follows:

[TABLE]

[TABLE]

[TABLE]

where $\lambda_{mM_{2}}<\lambda_{mM_{1}}$ and $k_{M_{2}}<k_{M_{1}}$ because of (3.5). We would like to note once again that monocytes were not involved in the model of the paper Ref. [22].

The remained macrophages “$M_{1}$” with a large intake of “$L_{ox}$” become foam cells “$F$”. This reaction can be simulated by the following equation:

[TABLE]

in the intima. Obviously, $\partial_{\nu}F=\partial_{\nu}M_{1}=0$ on $\Upsilon.$

Remark 3.2**.**

Since the transformation of a substance $S$ into a substance $P$ is a set of enzyme reactions, we use terms like $\lambda\frac{S}{k_{P}+S}$ in equations (3.9) – (3.12). However, we can take more general terms (3.1) and (3.2) without any restrictions. In reality, when selecting the view of each term, experimental study should be taken into account.

3.4. Modelling of the stage 3 (chronic inflammatory reaction)

As was mentioned above $M_{1}$-macrophages secrete pro-inflammatory cytokines that promote the recruitment of new monocytes; on the other hand, $M_{2}$-macrophages secrete anti-inflammatory cytokines that inhibit the recruitment of new monocytes. We model this by the term “${\copyright_{2}}$” in the boundary condition (3.4), namely

[TABLE]

Remark 3.3**.**

The activator-inhibitor system was firstly considered by Gierer and Meinhardt [16] in a pattern formation model (the $``u``$-activator – $``v``$-inhibitor term in this model is ${\frac{u^{2}}{v(1+Ku^{2})}})$ and now it has been used in many applications (see e.g. Ref. [35]).

The foam cells cause the growth of the intimal layer, thereby they amplify the endothelial dysfunction and additional penetration of LDLs. This feedback can be simulated by the term “${\copyright_{1}}$” in the boundary condition (3.3), namely

[TABLE]

Remark 3.4**.**

Such kind of feedback loop is common in the biochemical control circuit of substances $u_{1},\ldots,u_{n},$ which is described by the system of equations

[TABLE]

where $g$ is the feedback function. Typical choice of $g$ in applications is

[TABLE]

The first model of this type was proposed by Goodwin [18] with the following feedback function: $\displaystyle{g_{G}(u)=\frac{\nu_{0}}{1+(\frac{u}{k})^{p}}}.$

3.5. Nondimensionalization

The equations (3.7), (3.8), (3.9)–(3.12) with the boundary conditions mentioned above, with the Dirichlet conditions

[TABLE]

and with the initial conditions $L=L_{ox}=m=M_{1}=0=M_{2}=F=0$ at $t=0$ form the improved mathematical model of atherosclerosis development.

Remark 3.5**.**

It is easy to verify that the compatibility condition at $t=0$ is satisfied. General case $L|_{t=0}=L_{0},$ which describes the LDL concentration at the beginning of the atherosclerosis development, can be reduced with the substitution $L-L_{0}$ to the homogeneous one. This substitution explains the appearance of the function $\psi$ in (3.7) that depends on the LDL concentration in the blood.

Now we nondimensionalise the equations by setting

[TABLE]

[TABLE]

Then the pipe domain $C_{R,\rho_{0}}$ is transformed into a new pipe domain

[TABLE]

Its boundary consists of the bases (we drop the asterisks for algebraic simplicity)

[TABLE]

and two cylindrical parts $\Upsilon_{1,\varepsilon}:=\Upsilon_{1}\cup\Upsilon_{1+\varepsilon}.$ The surfaces $\omega$ and $\Omega$ are transformed into two cylindrical smooth surfaces $\omega_{1}$ and $\Omega_{1}$ respectively that belong to $\Upsilon_{1}.$

Let us denote by ${\bf u}:=(u_{1},\ldots,u_{6}).$ Then the nondimensional system is as follows:

[TABLE]

where the matrix $\mathfrak{D}={\rm diag}\left(1,d_{2},\ldots,d_{6}\right)$ introduces the diffusion constants

[TABLE]

the reaction terms

[TABLE]

and

[TABLE]

Since the boundary conditions (3.6), (3.14) and (3.13) are localized on $\omega$ and $\Omega$ respectively, in (3.17) we introduce special smooth functions $\varphi_{1}$ and $\varphi_{2}$ such that $\text{supp}(\varphi_{1})\subset\omega_{1},$ $0<\varphi_{1}\leq 1$ on $\omega_{1},$ $\text{supp}(\varphi_{2})\subset\Omega_{1},$ $0<\varphi_{2}\leq 1$ on $\Omega_{1}.$ All constants $\{\lambda_{i}\},$ $\{\mu_{i}\},$ $\{\delta_{i}\},$ and $\{\gamma_{i}\}$ are positive; $\mu_{3}<\mu_{2},$ $\lambda_{2}<\lambda_{1},$ $\mu_{5}<\mu_{4},$ $\mu_{8}<\mu_{6},$ $\lambda_{3}<1,$ $\lambda_{5}<1.$

3.6. Existence and uniqueness of the solution to problem (3.15)

For this we will use the method of upper and lower solutions, which was developed in Ref. [38]. This method leads not only to the basic results of existence and uniqueness of solutions but also to some their qualitative properties. For the convenience of readers we present some definitions from this book in adapting to problem (3.15).

Definition 3.1**.**

A vector-function ${\bf f}=(f_{1},\ldots,f_{6})$ $(f_{i}\equiv f_{i}(t,x,u_{1},\ldots,u_{6}))$ is said to possess a quasi-monotone property if for any $i\in\{1,\ldots,6\}$ there exist $a_{i},b_{i}\in\mathbb{N}_{0},\ a_{i}+b_{i}=5,$ such that $f_{i}(\cdot,\cdot,u_{i},[{\bf u}_{a_{i}}],[{\bf u}_{b_{i}}])$ is monotone nondecreasing in $[{\bf u}_{a_{i}}]$ and is monotone nonincreasing in $[{\bf u}_{b_{i}}].$ Here $[{\bf u}_{a_{i}}]$, $[{\bf u}_{b_{i}}]$ denote the $a_{i}$-components and $b_{i}$-components of the vector ${\bf u}\in\mathbb{R}^{6}_{+},$ respectively.

For example, the function

[TABLE]

is monotone nondecreasing in $u_{1},$ $u_{3},$ $u_{6}$ and is monotone nonincreasing in $u_{4},$ $u_{5}$ (recall that the function $\Theta$ is nonnegative (see subsection 3.1)); thus $a_{2}=3,$ $b_{2}=2.$

It is easy to verify that the vector-functions ${\bf F}$ and ${\bf G}$ are quasi-monotone.

Definition 3.2**.**

A pair of functions ${\bf{\tilde{u}}}=(\tilde{u}_{1},\ldots,\tilde{u}_{6})$, ${\bf{\hat{u}}}=(\hat{u}_{1},\ldots,\hat{u}_{6})$ from the space $\mathcal{C}(\overline{C_{\varepsilon}}\times[0,T])\cap\mathcal{C}^{0,1}((C_{\varepsilon}\cup\Upsilon_{1,\varepsilon})\times[0,T])\cap\mathcal{C}^{1,2}(C_{\varepsilon}\times(0,T))$ are called coupled upper and lower solutions to the problem (3.15) if ${\bf\tilde{u}}\geq{\bf\hat{u}}$ (component-wise) and if they satisfy the differential inequalities

[TABLE]

the boundary inequalities

[TABLE]

and the initial inequalities $\tilde{u}_{i}(0,x)\geq 0\geq\hat{u}_{i}(0,x),$ and $\tilde{u}_{i}=\hat{u}_{i}=0$ on $\Gamma^{\varepsilon}_{0,\ell}\times(0,T)$ for every $i=1,\ldots,6.$

Since f_{i}\big{(}t,x,0,[{\bf 0}]_{a_{i}},[{\bf\tilde{u}}]_{b_{i}}\big{)}\geq 0  if $[{\bf\tilde{u}}]_{b_{i}}\geq[{\bf 0}]_{b_{i}}$  and g_{i}\big{(}x,0,[{\bf 0}]_{c_{i}},[{\bf\tilde{u}}]_{d_{i}}\big{)}\geq 0 if $[{\bf\tilde{u}}]_{d_{i}}\geq[{\bf 0}]_{d_{i}}$ for $i=1,\ldots,6,$  ${\bf\hat{u}}\equiv\bf{0}$ is the lower solution to the problem (3.15).

Then the requirement for an upper solution is reduced to inequalities ${\bf\tilde{u}}(x,0)\geq{\bf 0},$

[TABLE]

Obviously, that the solution to the linear system

[TABLE]

is an upper solution to problem (3.15). In (3.18)

[TABLE]

and it is easily seen that for each $i\in\{1,\ldots,6\}$ the corresponding components of the vector-functions ${\bf\tilde{F}}=(\widetilde{f}_{1},\ldots,\widetilde{f}_{6}),$ ${\bf\tilde{G}}=(\widetilde{g}_{1},\ldots,\widetilde{g}_{6})$ satisfy the inequalities

[TABLE]

In fact, the system (3.18) is split into a system with respect $\widetilde{u}_{1}$ and $\widetilde{u}_{2},$ and four differential equations. Using classical results for linear parabolic boundary-value problems, we can state that there exists a unique nonnegative solution ${\bf\tilde{u}}$ to the problem (3.18) (this means the existence of an upper solution of (3.15)) and

[TABLE]

With the help of (3.19) it is easy to prove that for every $i\in\{1,\ldots,6\}$ there exit positive constants $K_{i}$ such that for the corresponding components of the vector-functions ${\bf F}$ and ${\bf G}$ the following inequalities hold:

[TABLE]

for all ${\bf 0}\leq{\bf u}\leq{\bf\tilde{u}},$ ${\bf 0}\leq{\bf v}\leq{\bf\tilde{u}}$ and $(x,t)\in\overline{C}_{\varepsilon}\times[0,T].$

Hence all conditions of Theorem 7.2 from Ref. [38][Chapter 9] are satisfied and we can state that the problem (3.15) has a unique solution ${\bf u}$ and

[TABLE]

Remark 3.6**.**

The existence and uniqueness of the solution to problem (3.15) can be also obtained from Ref. [29][Charpter VII].

4. Asymptotic analysis of problem (3.15)

In the boundary condition for the component $u_{1}$ (see (3.17)) the parameter $\gamma_{1}$ is equal to

[TABLE]

where $d_{L}$ is the diffusion coefficient of LDLs in the intima (it is equal to $29.89\,cm^{2}\,day^{-1}$ (see Refs. [12, 22])),  $k_{0}$ is the concentration of LDLs in the blood (estimated k_{0}\in\big{(}7\times 10^{-4}\,,\,1.9\times 10^{-3}\big{)}\ g\,cm^{-3}) at which the penetration rate is half of $V_{max},$  $\alpha_{0}$ is the penetration rate of LDLs from the blood into the intima through the damaged endothelium $\omega$ (estimated range is $(10^{-7},\ 3\times 10^{-4})\,g\,cm^{-2}day^{-1}),$  $R$ is the radius of a blood vessel. The size of blood vessels varies from $4\times 10^{-4}\,cm$ (radius of capillaries) to $1.2\,cm$ (radius of the aorta). Atherosclerosis usually affects arteries of large and medium caliber. Therefore, we regard $R\in(10^{-2}\,,\,1)\,cm.$ As a result,

[TABLE]

is a small parameter that can be compared with $\varepsilon$ (nondimensional thickness of the intima). According to Ref. [2], the ratio between the wall width and the radius of blood vessels is relatively constant and it lays between $0.15$ and $0.20$. Then, $\varepsilon$ varies from $0.05$ to $0.07$. Hence, we can introduce new values

[TABLE]

similarly in the boundary conditions for $u_{3}$ and $u_{4}$:

[TABLE]

where $\varrho_{i}\geq 1,$ and parameters $\eta_{1},\ldots,\eta_{5}$ can be involved in the study of active biochemical phenomena on the vessel wall like the endothelium dysfunction, adhesion of LDLs and monocytes, and their penetration or efflux.

For further asymptotic analysis we assume that $\varrho_{1}=\ldots=\varrho_{5}=:\varrho\geq 1.$ This assumption is made only for the sake of simplicity and it will be clear how to proceed in general case from further calculations.

Also small parameters are appeared through the coefficients

[TABLE]

[TABLE]

Because of $\varepsilon\in(0.05,0.07),$ we regard that $3\leq\tau_{1}<\tau_{2}<2\tau_{1}.$ Thus, we get the following parabolic system perturbed by the small parameter $\varepsilon:$

[TABLE]

where $\mathfrak{D}_{\varepsilon}={\rm diag}\left(1,d_{2},d_{3},d_{4},{\varepsilon^{\tau_{1}}},{\varepsilon^{\tau_{2}}}\right),$ {\bf G}({\bf u})=\big{(}g_{1}({\bf u}),0,g_{3}({\bf u}),g_{4}({\bf u}),0,0\big{)} and

[TABLE]

Next our aim is to construct the asymptotic approximation for the solution to the problem (4.2) as the small parameter $\varepsilon\to 0,$  to derive the limit problem $(\varepsilon=0),$ and to prove the corresponding asymptotic estimates.

It should be noted that the limit process in the thin tube domain $C_{\varepsilon}$ is accompanied by the perturbed coefficients both in the differential equations and in the boundary conditions.

4.1. Construction of the asymptotic approximation

At first we rewrite the system (4.2) in the cylindrical coordinates $(x_{1},r,\vartheta):$

[TABLE]

Here

[TABLE]

is the Laplace operator in the cylindrical coordinates, the thin tube domain $C_{\varepsilon}$ is transformed into the thin plate

[TABLE]

with the following parts of the boundary: {\it\Gamma}^{\varepsilon}_{0,\ell}:=\partial\mathfrak{C}_{\varepsilon}\cap\big{(}\{x_{1}=0\}\cup\{x_{1}=\ell\}\big{)},

[TABLE]

Two cylindrical smooth surfaces $\omega_{1}$ and $\Omega_{1}$ are transfigured into smooth plane surfaces ${\it\omega}_{1}$ and ${\it\Omega}_{1}$ that belong to ${\it\Upsilon}_{1}.$

Using the asymptotic approach to construct approximations for solutions of boundary-value problems in thin domains (see e.g. Refs. [17, 32, 28]), we propose the following approximation for the solution to problem (4.6):

[TABLE]

where ${\bf v}=(v_{1},\ldots,v_{6})$ and ${\bf W}=(W_{1},\ldots,W_{4},0,0).$

Let us substitute ${\bf R}^{\varepsilon}$ into (4.6) instead of ${\bf u}^{\varepsilon}.$ Since the operator $\Delta_{x_{1},r,\vartheta}$ in the variables $(x_{1},\xi_{2},\vartheta)$ takes the form

[TABLE]

where $\xi_{2}=\frac{r-1}{\varepsilon},$ the collection of coefficients at the same power of $\varepsilon$ with regards (3.20) and (3.21) gives the following four $(i=1,\ldots,4)$ one-dimensional Neumann boundary-value problems with respect to $\xi_{2}:$

[TABLE]

and two ordinary differential equations

[TABLE]

In (4.9), $\delta_{1,\varrho}$ is Kronecker’s symbol (recall that $\varrho\geq 1$), $\Delta_{x_{1},\vartheta}v=\frac{\partial^{2}v}{\partial x_{1}^{2}}+\frac{\partial^{2}v}{\partial\vartheta^{2}},$ and the variables $x_{1},\vartheta,t$ are regarded as parameters.

The solvability condition for the problem (4.9) at a fixed index $i\in\{1,\ldots,4\}$ is given by the differential equation

[TABLE]

Let $\{v_{i}\}_{i=1}^{6}$ be solutions of the respective differential equations (boundary and initial conditions for them will be determined later). Then solutions of the problems (4.9) exist and the additional relations

[TABLE]

supply their uniqueness. Since $g_{2}\equiv 0,$ the solution ${W}_{2}\equiv 0.$ Thus, ${\bf W}=(W_{1},0,W_{3},W_{4},0,0).$

Equations (4.11) and (4.10) form the limit system that is supplied by corresponding boundary and initial conditions. As a result, we get the coupled parabolic-ordinary system

[TABLE]

with the reaction terms

[TABLE]

[TABLE]

Here we have written down the reaction terms of the limit problem without the assumption that $\varrho_{1}=\ldots=\varrho_{5}=:\varrho.$

The method of upper and lower solutions can be used to the system (4.13) as well (see Ref. [38][§ 9.7]). Similar as in subsection 3.6, we can verify that the reaction terms ${\bf F}^{(0)}$ are quasi-monotone, ${\bf\hat{v}}\equiv{\bf 0}$ is a lower solution of (4.13), and there exists an upper solution ${\bf\tilde{v}}.$ Then it follows from Theorem 7.3 (Ref. [38][Chapter 9]) that the limit problem (4.13) has a unique solution ${\bf v}$ and

[TABLE]

From positivity lemma (see Ref.[38][Ch. 2]) it follows that either ${\bf v}>{\bf 0}$ (component-wise) in $(0,\ell)\times[0,2\pi]\times(0,T)$ or ${\bf v}\equiv{\bf 0}$ that is impossible.

To be specific we will carry out further investigation for the more interesting case $\varrho=1$ (if $\varrho>1,$ then ${\bf W}\equiv 0).$

Thanks to the boundary conditions for the functions $v_{1},\ldots,v_{4},$ the initial value problem for $v_{5}$ at $x_{1}=0$ has the view

[TABLE]

Due to Cauchy-Lipschitz theorem on existence and uniqueness of solutions to first-order equations with given initial conditions, $v_{5}(0,\vartheta,t)=0$ for $\vartheta\in[0,2\pi]$ and $t\in[0,T].$ Similarly, we verify that $v_{6}$ and $v_{5}$ vanish at $x_{1}=0$ and $x_{1}=\ell$ as well. Thus, ${\bf R}^{\varepsilon}|_{x_{1}=0}={\bf R}^{\varepsilon}|_{x_{1}=\ell}.$

The approximation function ${\bf R}^{\varepsilon}$ takes into account the inhomogeneity of the right-hand side ${\bf F}$ and the boundary conditions on ${\it\omega}_{1}$ and ${\it\Omega}_{1}$ that belong with their closures in the horizontal side ${\it\Upsilon}_{1}$ of the thin domain $\mathfrak{C}_{\varepsilon}.$ Since these boundary conditions are localized inside of ${\it\Upsilon}_{1},$ the functions $W_{1},W_{3},W_{4}$ vanish outside of $\overline{{\it\omega}_{1}}$ and $\overline{{\it\Omega}_{1}}$ respectively. Therefore, \partial_{r}\big{(}v_{1}+\varepsilon^{2}W_{1}\big{)}|_{r=1}=0 on ${\it\Upsilon}_{1}\setminus{\it\omega}_{1},$ and $\partial_{r}{\bf R}^{\varepsilon}|_{r=1}=0$ on ${\it\Upsilon}_{1}\setminus{\it\Omega}_{1}.$ Clearly that $\partial_{r}{\bf R}^{\varepsilon}|_{r=1+\varepsilon}=0$ and $\partial_{\vartheta}^{k}{\bf R}^{\varepsilon}|_{\vartheta=0}=\partial_{\vartheta}^{k}{\bf R}^{\varepsilon}|_{\vartheta=2\pi}$ for $k=0,1.$

Now due to (4.11) we can rewrite problems for $W_{1},W_{3},W_{4}$ as follows

[TABLE]

Since $g_{i}({\bf v})|_{t=0}=0,$ we have $W_{i}|_{t=0}=0$ at $i\in\{1,3,4\}.$ Therefore, ${\bf R}^{\varepsilon}|_{t=0}={\bf 0}.$

4.2. Justification

In virtue that $W_{1},W_{3},W_{4}$ vanish in neighbourhoods of the sides $x_{1}=0$ and $x_{1}=\ell$ and ${\bf v}>{\bf 0}$ (see before), we can regard that ${\bf 0}\leq{\bf R}^{\varepsilon}\leq{\bf\tilde{u}}$ for $\varepsilon$ small enough. Therefore, due to (3.20) and (3.21) we have

[TABLE]

in $\overline{\mathfrak{C}}_{\varepsilon}\times[0,T].$

Remark 4.1**.**

Here and in what follows all constants $\{c_{i}\}$ in inequalities are independent of the parameter $\varepsilon.$

Substituting ${\bf R}^{\varepsilon}$ in the equations of the system (4.6) and the respective boundary conditions on ${\it\omega}_{1}$ and ${\it\Omega}_{1},$ and thanks to (4.9) and (4.13), we find

[TABLE]

in $\mathfrak{C}_{\varepsilon}\times(0,T),$

[TABLE]

respectively on ${\it\omega}_{1}$ and ${\it\Omega}_{1}.$ In (4.17) {\bf\Psi}^{\varepsilon}:=\big{(}{\bf F}({\bf v})-{\bf F}({\bf R}^{\varepsilon})\big{)}+{\bf\Theta}^{\varepsilon} and the components of ${\bf\Theta}^{\varepsilon}$ are as follow

[TABLE]

for $i=1,\ldots,4$ and

[TABLE]

Taking into account that $r\in(1,1+\varepsilon)$ and $\Psi^{\varepsilon}_{i}=\Theta^{\varepsilon}_{i}$ for $i\in\{2,5,6\},$ there exist positive constants $c_{3}$ and $\varepsilon_{0}$ such that for all $\varepsilon\in(0,\varepsilon_{0})$

[TABLE]

Thus, the difference between the approximation function ${\bf R}^{\varepsilon}$ and the solution ${\bf u}^{\varepsilon}$ to the problem (4.6) satisfies the following relations:

[TABLE]

or

[TABLE]

in the variables $x=(x_{1},x_{2},x_{2}).$ Here ${\bf\Phi}^{\varepsilon}:={\bf G}({\bf v})-{\bf G}({\bf R}^{\varepsilon}),$ ${\Phi}^{\varepsilon}_{i}=0,\ i=2,5,6,$ and

[TABLE]

To deduce an estimate for the difference ${\bf R}^{\varepsilon}-{\bf u}^{\varepsilon}$ we use Theorem 7.3 in Chapter 5 of Ref. [29], where an apriori estimate for the solution to the semilinear parabolic initial-boundary value problem

[TABLE]

was proved under the following conditions for any $u$:

[TABLE]

Here $\nu_{1},\mu_{1}>0,$ $\beta_{i}\geq 0,\ i\in\{0,1,\ldots,4\}.$ As a result, for any solution $u$ from $C^{2,1}(\Omega\times(0,T])\cap C(\overline{\Omega}\times[0,T])\cap C^{1,0}(\overline{\Omega}\times(0,T])$ the following estimate holds:

[TABLE]

where the constants $\lambda_{1}$ and $\lambda$ are defined with the constants $\nu_{1},\mu_{1},\beta_{0},\beta_{1},\beta_{3}.$

With small changes the estimate carries over to every component of the solution ${\bf R}^{\varepsilon}-{\bf u}^{\varepsilon}$ to the problem (4.22). For each $i\in\{1,\ldots,4\}$ the condition (4.23) is satisfied. Then using Young’s inequality $(ab\leq\varsigma a^{2}+\tfrac{1}{4\varsigma}b^{2}$ for all $a,b\in\mathbb{R}$ and $\varsigma>0)$ with an suitable $\varsigma,$ (4.20), (4.16), properties of the components of (3.16), (4.3)-(4.5) and nonnegativity of ${\bf u}^{\varepsilon}$ and ${\bf R}^{\varepsilon}$, we derive analogues of the inequality (4.24)

[TABLE]

and analogues of the second inequality in (4.25) for $i\in\{1,3,4\}$

[TABLE]

Thus, due to (4.26) we have

[TABLE]

For $i\in\{5,6\}$ $v_{i}\equiv R^{\varepsilon}_{i}$ and $\Phi^{\varepsilon}_{i}\equiv 0,$ and we can repeat the proof of Theorem 7.3 [29] under conditions (4.23) and (4.24) only. In this case the constant $\lambda_{1}$ in (4.26) is independent of $\varepsilon.$ As a result, similarly as before we prove that

[TABLE]

Sequentially substituting inequality (4.34)–(4.39) one in the other and taking

[TABLE]

we get

[TABLE]

Substituting (4.42) in (4.34), (4.41) and then obtained inequalities sequentially in (4.40), (4.35) and (4.36), we get estimates for the other components. Thus, we have proved the following theorem.

Theorem 4.1**.**

There exist positive constants $C_{0},\varepsilon_{0}$ such that for all values $\varepsilon\in(0,\varepsilon_{0})$ the difference between the solution ${\bf u}^{\varepsilon}$ to the problem (4.2) and the solution ${\bf v}$ to the the limit coupled parabolic-ordinary system (4.13) satisfies the following estimate

[TABLE]

Remark 4.2**.**

We think that the estimate (4.43) could be obtained with the help of special integral representations using Green’s functions for the corresponding linear initial-boundary value problems (see e.g. Ref. [38][§9.6]). But the main difficulty will be in establishing the dependence of the constant on the small parameter $\varepsilon$ in the proof of Lemma 9.6.1 [38].

Multiplying each parabolic equation in (4.22) by the respective difference $R_{i}^{\varepsilon}-u^{\varepsilon}_{i},$ integrating over $C_{\varepsilon}\times(0,T)$ by parts and using the estimates (4.27)–(4.33) and (4.43), we find that for any $t\in(0,T)$

[TABLE]

Hence, the following estimate for gradients holds.

Corollary 4.1**.**

[TABLE]

5. Conclusions

Overview mathematical approaches presented in Section 2 and the results obtained in this paper suggest that the development of new models of the atherosclerosis should be based on known scenarios of enzymatic reactions, namely on the corresponding differential equations describing these reactions. Furthermore, those mathematical models should be nondimensionalised to introduce small (or large) scales in the system dynamics with a multiscale perspective to derive the respective simplified limit problem. In our case this is the small parameter $\varepsilon$ (see Section 4).

It is very important for a proposed multi-scale method to justify its stability and accuracy. The proof of the error estimate between the constructed approximation and the exact solution is a general principle that has been applied to the analysis of the efficiency of a multi-scale method. In the present paper, the asymptotic approximation ${\bf R}^{\varepsilon}$ (see (4.7)) for the solution ${\bf u}^{\varepsilon}$ to problem (4.2) is constructed and justified. The results obtained in Theorem 4.1 and Corollary 4.1 argue that it is possible to replace the perturbed parabolic system (4.2) in the thin tube domain $C_{\varepsilon}$ with the corresponding limit problem (4.13) in the rectangle $(0,\ell)\times(0,T)$ with the sufficient accuracy measured by the parameter $\varepsilon$ characterizing the non-dimensional thickness of the intima layer.

From the viewpoint of mathematical modelling, this work is only an early start to the study of the diverse roles of macrophages in atherosclerosis, because there remain a few important problems unexplored. Among these open problems the following problems are most important and interesting in author’s opinion.

• •

From a practical point of view it is very important to stabilize the development of atherosclerosis. A basic mathematical question about this problem is whether the solution ${\bf u}^{\varepsilon}$ to the problem (4.2) converges to the steady state as $t\to+\infty.$ Obviously, first of all this will depend on the parameters $\eta_{1},\ldots,\eta_{5}$ and $p_{1},$ $p_{2},$ $p_{3}$ that are responsible for active biochemical phenomena on the vessel wall (the endothelium dysfunction, adhesion of LDLs and monocytes, penetration of monocytes, and efflux of $M_{2}$-macrophages), and also on the other parameters in the problem (4.2).

Due to the asymptotic results obtained in this paper, the problem is reduced to the study of the time evolution of the solution to the limit problem (4.13). This will be a simpler problem, since there are no nonlinear boundary conditions. It is known (see Ref. [38], Chapter 10) that in the case of a homogeneous Neumann boundary condition the asymptotic behavior (as $t\to+\infty)$ is compared with the behavior of the solution to the corresponding ordinary differential system.

• •

To study the influence of the parameters of the limit problem (4.13) on the velocity of the atherosclerosis development. In this connection, it will be interesting to examine the unboundedness of the solution to the problem (4.13) with emphasis on its blowing-up behavior in finite time.

• •

During the atherosclerosis development, the intimal thickness and the endothelial damage area grow and this should be taken into account. We did this in the boundary condition (3.14) through the special feedback term. However, it will be more appropriate to consider the corresponding free boundary problem modeling both the growth of the intimal tissue and damage area. For this the methods proposed in Refs. [4, 5, 6, 7, 8] will be helpful.

Of course, the introduced model is an incomplete representation of all the processes of the development of atherosclerosis. However, we believe that it reflects many of the main features of the atherosclerosis and can give new understanding into the complex interactions of the disease. In addition, we can complicate the model taking into account other features. For instance, when the effect of diffusion and convection are both taken into consideration then the new term ${\bf b}\cdot\nabla u_{i}$ appears in the corresponding differential equation (see Refs. [3, 15]). In this case the method of upper and lower solutions can also apply to reaction diffusion-convection system with quasimonotone functions (see Ref. [38], Chapter 12). Also it would be interesting to include some chemotaxis terms in the model (see Refs. [9, 23, 25, 35]).

Acknowledgment

This research was begun at the University of Stuttgart under support of the Alexander von Humboldt Foundation in the summer of 2015. The author is grateful to Prof. Christian Rohde for the hospitality and wonderful working conditions. Then this research was continued in the framework of the Marie Curie IRSES project ”EU-Ukrainian Mathematicians for Life Sciences” (FP7-People-2011-IRSES Project number 295164) and of the Marie Curie RISE project ”Approximation Methods for Modelling and Diagnosis Tools”. Some results were reported on the workshops of these projects.