Steady solution and its stability of a mathematical model of diabetic atherosclerosis

Abstract

Atherosclerosis is a leading cause of death worldwide. Making matters worse, nearly 463 million people have diabetes, which increases atherosclerosis-related inflammation. Diabetic patients are twice as likely to have a heart attack or stroke. In this paper, we consider a simplified mathematical model for diabetic atherosclerosis involving LDL, HDL, glucose, insulin, free radicals (ROS), β cells, macrophages and foam cells, which satisfy a system of partial differential equations with a free boundary, the interface between the blood flow and the plaque. We establish the existence of small radially symmetric stationary solutions to the model and study their stability. Our analysis shows that the plague will persist due to hyperglycemia even when LDL and HDL are in normal range, hence confirms that diabetes increase the risk of atherosclerosis.

Keywords:

• Atherosclerosis
• diabetes
• free boundary problem
• steady solution
• stability

Mathematic Subject classifications:

• 35Q92
• 92C50

1. Introduction

Atherosclerosis is a leading cause of death in the United States and worldwide. Atherosclerosis emerges as a result of multiple dynamical cell processes [8,24,38]. The atherosclerosis inflammatory processes have been modelled mathematically [1,27]. The governing equations used in LDL oxidation were either ordinary differential equation (ODE) [7,32] or partial differential equation (PDE) [6,10–12,17].The dynamics of cells were modelled by advection- diffusion equations that describe monocyte recruitment and chemoattractants, monocyte to macrophage differentiation, foam cell formation, T-cell recruitment and proliferation of SMCs [4,6,9,11,17,30,35]. The formation of a plaque was then modelled by a free boundary problem [11,17]. The results obtained from previous models are consistent with the guidelines issued by the American Heart Association periodically regarding to the risk of a heart attack associated with high level of cholesterol [11,17].

Nearly 463 million people worldwide have diabetes, a disease where patients' cells cannot efficiently take in dietary sugar, causing it to build up in the blood. When insulin secretion by the pancreas is insufficient or absent, due to (autoimmune) destruction of β-cells, the clinical picture of Type 1 Diabetes Mellitus (T1DM) results; when insulin is secreted in normal, or supernormal amounts, but it is ineffective in lowering glycemia to normal levels, Type 2 Diabetes Mellitus (T2DM) is said to be present. Mathematical modelling of the glucose-insulin feedback system has long been established [3,13,15,16,25,26]. Diabetic patients have an increased risk of cardiovascular disease (CVD), and have a more than twofold increase in the risk of dying from CVD [2].

Diabetes and atherosclerosis are chronic inflammatory conditions. Myeloid cells (neutrophils, monocytes, and macrophages) are involved in both atherosclerosis and diabetes [2,22,28,34,39]. The migration of circulating monocytes into the vessel wall is critical for the development of diabetic atherosclerosis. Moreover, macrophage migration inhibitory factor (MIF), which regulate the adhesion of monocytes, are dysregulated in hyperglycemia-induced atherosclerosis. Increased foam cells derived from macrophages promote the acceleration of atherosclerotic lesions in diabetic mice [21,33]. A more inflammatory monocyte/macrophage phenotype with secretion of higher levels of proinflammatory cytokines was detected in both animal models and patients with diabetes mellitus [23]. Diabetes cause endothelial cell dysfunction: Endothelial dysfunction due to inflammation and oxidative stress is a crucial characteristic in diabetes mellitus-linked atherosclerosis. Endothelial dysfunction is associated with decreased nitric oxide (NO) availability, either through loss of NO production or NO biological activity [28,37]. Excess glucose metabolites inhibits production of NO by blocking eNOS synthase activation and increase the production of ROS, and insulin resistance decrease endothelium-derived NO and increase the production of ROS. The excess generation of free oxygen radicals leads to apoptosis in endothelial cells [37]. In hyperglycemia, chronic inflammation increases vascular permeability, promotes the generation of adhesion molecules and chemokines, and stimulates accumulation of monocytes in the artery wall. Although the pathophysiology of diabetic vascular disease is generally understood, there is no mathematical model that includes the effect of diabetes on plaque growth until recent works [40,41]. In [40,41], governing equations were established to model the mechanism of how diabetes impairs endothelium-dependent (nitric oxide-mediated) vasodilation before the formation of atheroma and mechanism of how arteries affected by diabetes have altered vascular smooth muscle cell function. The local existence and uniqueness of solution to the diabetic atherosclerosis models were also proved.

In this paper, we consider a simplified mathematical model for diabetic atherosclerosis involving LDL, HDL, glucose, insulin, free radicals (ROS),β cells, macrophages and foam cells, which satisfy a system of partial differential equations with a free boundary, the interface between the blood flow and the plaque. We establish local existence of small radially symmetric stationary solutions to the model and study their stability. Our analysis shows that the plague will persist due to hyperglycemia even when LDL and HDL are in normal range, hence confirms that diabetes increase the risk of atherosclerosis.

2. Mathematical model

We will use a system of partial differential equations to develop a mathematical model of plaque formation in diabetic atherosclerosis. The mathematical model will be based on the network shown in [40,41]. The variables included in the model are listed in Table 1.

Table 1. The variables of the model: concentrations and densities are in units of ${\mathrm{g}/\mathrm{cm}}^3$.

L:	concentration of LDL	H:	concentration of HDL
G:	concentration of Glucose	I:	concentration of Insulin
R:	concentration of free radicals	M:	density of macrophages
β:	density of β-cells	F:	density of foam cells
$\mathbf{u}$:	fluid velocity (in $\mathrm{cm}/\mathrm{day}$)

Table 2. Steady concentrations of proteins and densities of cells are in units of ${\mathrm{g}/\mathrm{cm}}^3$.

Proteins and cells	Concentrations (${\mathrm{gcm}}^3$)	Explaination
L: LDL	$7\times {10}^{-4}-1.9\times {10}^{-3}$	Range is 70−190 mg/dl [17]
H: HDL	$4\times {10}^{-4}-1.6\times {10}^{-4}$	Range is 40−60 mg/dl [17]
G: Glucose	${10}^{-3}$	100 mg/dl [36]
I: Insulin	${10}^{-4}$	10 mg/dl [36]
β: β-cells	${10}^{-3}$	100 mg/dl [36]

We assume that all cells are moving with a common velocity $\mathbf{u}$ [11,17]; the velocity is the result of movement of macrophages and β cells into the intima. We also assume that all species are diffusing with appropriate diffusion coefficients. In this work, we consider only two dimensional radially symmetric plaque $\mathrm{\Omega }(t)$ which is given by $B(t)<r<1,r=\sqrt{x^2+y^2}$, where $r=B(t)$ is the free boundary between the intima and the lumen.

Following the charts for the variable interactions in [40,41], we assume that the variables satisfy the following PDEs in $\mathrm{\Omega }(t)$: (1) $\begin{array}{rl}& \frac{\mathrm{\partial }G}{\mathrm{\partial }t}-D_G\mathrm{\Delta }G={\tilde{G}}_0-(E_{G0}+S_{I}I)G-{\lambda }_{GH}(G-G_s)H,\end{array}$(1) (2) $\begin{array}{rl}& \frac{\mathrm{\partial }\beta }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}\beta )-D_{\beta }\mathrm{\Delta }\beta =(-d_0+r_{1}G-r_2{G}^2)\beta ,\end{array}$(2) (3) $\begin{array}{rl}& \frac{\mathrm{\partial }I}{\mathrm{\partial }t}-D_I\mathrm{\Delta }I=\frac{\sigma G^2\beta }{\alpha +G^2}-k_{I}I-{\lambda }_{RI}RI;\end{array}$(3) (4) $\begin{array}{rl}& \frac{\mathrm{\partial }L}{\mathrm{\partial }t}-D_L\mathrm{\Delta }L=-k_{L}RL+{\lambda }_{GL}(G-G_s)L-{\lambda }_{ML}\frac{ML}{K_{ML}+L},\end{array}$(4) (5) $\begin{array}{rl}& \frac{\mathrm{\partial }H}{\mathrm{\partial }t}-D_H\mathrm{\Delta }H=-k_{H}RH-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)},\end{array}$(5) (6) $\begin{array}{rl}& \frac{\mathrm{\partial }R}{\mathrm{\partial }t}-D_R\mathrm{\Delta }R=R_0-R(k_{L}L+k_{H}H)+{\lambda }_{GR}(G-G_s)R+{\lambda }_{IR}IR,\end{array}$(6) (7) $\begin{array}{rl}& \frac{\mathrm{\partial }M}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}M)-D_M\mathrm{\Delta }M=-{\lambda }_{ML}\frac{ML}{K_{ML}+L}+{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}+\lambda \frac{ML}{K_{MH}+H}-d_{M}M,\end{array}$(7) (8) $\begin{array}{rl}& \frac{\mathrm{\partial }F}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}F)-D_F\mathrm{\Delta }F={\lambda }_{ML}\frac{ML}{K_{ML}+L}-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}-d_{F}F,\end{array}$(8) Equation (1(1) $\begin{array}{rl}& \frac{\mathrm{\partial }G}{\mathrm{\partial }t}-D_G\mathrm{\Delta }G={\tilde{G}}_0-(E_{G0}+S_{I}I)G-{\lambda }_{GH}(G-G_s)H,\end{array}$(1) ) describes the dynamics of glucose and insulin, the third term on the right side of (1(1) $\begin{array}{rl}& \frac{\mathrm{\partial }G}{\mathrm{\partial }t}-D_G\mathrm{\Delta }G={\tilde{G}}_0-(E_{G0}+S_{I}I)G-{\lambda }_{GH}(G-G_s)H,\end{array}$(1) ) represents that HDL helps lower glucose and equation (2(2) $\begin{array}{rl}& \frac{\mathrm{\partial }\beta }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}\beta )-D_{\beta }\mathrm{\Delta }\beta =(-d_0+r_{1}G-r_2{G}^2)\beta ,\end{array}$(2) ) describes the formation and loss of β-cells [36]. The first term on the right-hand side of (3(3) $\begin{array}{rl}& \frac{\mathrm{\partial }I}{\mathrm{\partial }t}-D_I\mathrm{\Delta }I=\frac{\sigma G^2\beta }{\alpha +G^2}-k_{I}I-{\lambda }_{RI}RI;\end{array}$(3) ) is secretion of insulin [36], the second term is clearance of insulin. The third term on the right side of (3(3) $\begin{array}{rl}& \frac{\mathrm{\partial }I}{\mathrm{\partial }t}-D_I\mathrm{\Delta }I=\frac{\sigma G^2\beta }{\alpha +G^2}-k_{I}I-{\lambda }_{RI}RI;\end{array}$(3) ) models that ROS, in excess and over time, cause chronic oxidation stress, which in turn result in reduction of insulin secretion as well as increased apoptosis [29,31], where ${\lambda }_{RI}$ is the reduction rate of insulin due to ROS.

Equation (4(4) $\begin{array}{rl}& \frac{\mathrm{\partial }L}{\mathrm{\partial }t}-D_L\mathrm{\Delta }L=-k_{L}RL+{\lambda }_{GL}(G-G_s)L-{\lambda }_{ML}\frac{ML}{K_{ML}+L},\end{array}$(4) ) and (5(5) $\begin{array}{rl}& \frac{\mathrm{\partial }H}{\mathrm{\partial }t}-D_H\mathrm{\Delta }H=-k_{H}RH-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)},\end{array}$(5) ) describe the distribution of LDL, HDL. We merge LDL and HDL with oxidized LDL and oxidized HDL in these equations. LDL and HDL are lost by reaction of oxidation with free radicals, $k_L$ and $k_H$ are reaction rates of oxidization. LDL is ingested by macrophages and its production is enhanced by AGEs which is assumed to be proportional to the glucose. A reduction of ox-HDL through ingestion by foam cells is represented by the second term on right-hand side of (5(5) $\begin{array}{rl}& \frac{\mathrm{\partial }H}{\mathrm{\partial }t}-D_H\mathrm{\Delta }H=-k_{H}RH-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)},\end{array}$(5) ). ${\lambda }_{FH}$ is the reduction rate of ox-HDL due to ingestion by the foam cells.

Equation (6(6) $\begin{array}{rl}& \frac{\mathrm{\partial }R}{\mathrm{\partial }t}-D_R\mathrm{\Delta }R=R_0-R(k_{L}L+k_{H}H)+{\lambda }_{GR}(G-G_s)R+{\lambda }_{IR}IR,\end{array}$(6) ) models the concentration of radicals, where $R_0$ is the base line growth, the second represent the reduction of radicals due to the oxidization of LDL and HDL; the third term on the right hand of the equation models the mechanism whereby excess glucose metabolites inhibits production of NO by blocking eNOS synthase activation and increase the production of ROS [22,39]; where ${\lambda }_{GR}$ is the growth rate of ROS due to excess glucose. The fourth term on the right hand of the equation models the mechanism that insulin resistance decrease endothelium-derived NO and increase the production of ROS [18,19]; where ${\lambda }_{IR}$ is the growth rate of ROS due to insulin resistance.

The evolution of macrophage density is modelled by (7(7) $\begin{array}{rl}& \frac{\mathrm{\partial }M}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}M)-D_M\mathrm{\Delta }M=-{\lambda }_{ML}\frac{ML}{K_{ML}+L}+{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}+\lambda \frac{ML}{K_{MH}+H}-d_{M}M,\end{array}$(7) ). Here the first and the second term on right-hand side accounts for the transition between M and L: when L molecules are ingested by a macrophage, the macrophage becomes a foam cell; when H combines with membrane protein on a foam cell in a process that clears it from the bad cholesterol, the foam cell turns into a macrophage. The third term is phenomenological [12]: the factor ML accounts for a sequence of reactions, stemming from oxidized LDL, which results in the activation of inflammatory macrophages, while the inhibition factor $\frac{1}{K_{MH}+H}$ is introduced to account for the fact that by oxidizing with radicals, HDL removes some of the radicals available to oxidize LDL. The fourth term accounts for the death of macrophages.

(8(8) $\begin{array}{rl}& \frac{\mathrm{\partial }F}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}F)-D_F\mathrm{\Delta }F={\lambda }_{ML}\frac{ML}{K_{ML}+L}-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}-d_{F}F,\end{array}$(8) ) is the equation for foam cells, where on the right side of the equation, a gain of foam cells from macrophages M ingesting L, and a loss of foam cells triggered by H. The death rate of foam cells is $d_F$. We assume that (9) $\begin{array}{r}M+F+\beta =M_0\end{array}$(9) By adding Eqs. (2(2) $\begin{array}{rl}& \frac{\mathrm{\partial }\beta }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}\beta )-D_{\beta }\mathrm{\Delta }\beta =(-d_0+r_{1}G-r_2{G}^2)\beta ,\end{array}$(2) ), (7(7) $\begin{array}{rl}& \frac{\mathrm{\partial }M}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}M)-D_M\mathrm{\Delta }M=-{\lambda }_{ML}\frac{ML}{K_{ML}+L}+{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}+\lambda \frac{ML}{K_{MH}+H}-d_{M}M,\end{array}$(7) ), and (8(8) $\begin{array}{rl}& \frac{\mathrm{\partial }F}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}F)-D_F\mathrm{\Delta }F={\lambda }_{ML}\frac{ML}{K_{ML}+L}-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}-d_{F}F,\end{array}$(8) ) and using (9(9) $\begin{array}{r}M+F+\beta =M_0\end{array}$(9) ), we get an equation for $\mathbf{u}=u(r){\overrightarrow{e}}_r$: (10) $\begin{array}{r}{M}_0\frac{1}{r}{\mathrm{\partial }}_r(ru)=(-d_0+r_{1}G-r_2{G}^2)\beta +\lambda \frac{(M_0-F-\beta )L}{K_{MH}+H}-d_{F}F-d_M(M_0-F-\beta ).\end{array}$(10) For simplicity, we consider only 2-dimensional plaques. Then the boundary of the plaque consists of i) ${\mathrm{\Gamma }}_M=\{r=1\}$, in contact with media; ii) a free boundary ${\mathrm{\Gamma }}_I=\{r=B(t)\}$, inside the lumen.

We assume flux boundary conditions of the form on the free boundary $r=B(t)$: (11) $\begin{array}{rl}& -\frac{\mathrm{\partial }L}{\mathrm{\partial }n}+{\alpha }_L(L-L_0)=0,-\frac{\mathrm{\partial }H}{\mathrm{\partial }n}+{\alpha }_H(H-H_0)=0,\end{array}$(11) (12) $\begin{array}{rl}& -\frac{\mathrm{\partial }G}{\mathrm{\partial }n}+{\alpha }_G(G-G_0)=0,-\frac{\mathrm{\partial }I}{\mathrm{\partial }n}+{\alpha }_I(I-I_0)=0,-\frac{\mathrm{\partial }R}{\mathrm{\partial }n}+{\alpha }_R(R-R_0)=0,\end{array}$(12) (13) $\begin{array}{rl}& -\frac{\mathrm{\partial }M}{\mathrm{\partial }n}+{\alpha }_M(M-{\overline{M}}_0)=0,-\frac{\mathrm{\partial }\beta }{\mathrm{\partial }n}+{\alpha }_M(\beta -{\beta }_0)=0,-\frac{\mathrm{\partial }F}{\mathrm{\partial }n}+{\alpha }_{M}F=0,\end{array}$(13) where we assume ${\overline{M}}_0+{\beta }_0=M_0$. Note that $G_0,I_0,L_0$ and $H_0$ are the glucose, insulin, LDL and HDL concentrations in the blood, so we shall be interested to see how these concentrations determine whether a small plaque will grow or shrink.

The equation for the free boundary is (14) $\begin{array}{r}\frac{dB(t)}{dt}=u(B(t)).\end{array}$(14) We assume no flux condition on the fixed boundary r = 1: (15) $\begin{array}{r}\frac{\mathrm{\partial }L}{\mathrm{\partial }n}=\frac{\mathrm{\partial }H}{\mathrm{\partial }n}=\frac{\mathrm{\partial }G}{\mathrm{\partial }n}=\frac{\mathrm{\partial }I}{\mathrm{\partial }n}=\frac{\mathrm{\partial }R}{\mathrm{\partial }n}=\frac{\mathrm{\partial }M}{\mathrm{\partial }n}=\frac{\mathrm{\partial }\beta }{\mathrm{\partial }n}=\frac{\mathrm{\partial }F}{\mathrm{\partial }n}=0.\end{array}$(15) and (16) $\begin{array}{r}u(1)=0.\end{array}$(16)

Parameter estimation: Most parameters on the right-hand side of (1(1) $\begin{array}{rl}& \frac{\mathrm{\partial }G}{\mathrm{\partial }t}-D_G\mathrm{\Delta }G={\tilde{G}}_0-(E_{G0}+S_{I}I)G-{\lambda }_{GH}(G-G_s)H,\end{array}$(1) )–(3(3) $\begin{array}{rl}& \frac{\mathrm{\partial }I}{\mathrm{\partial }t}-D_I\mathrm{\Delta }I=\frac{\sigma G^2\beta }{\alpha +G^2}-k_{I}I-{\lambda }_{RI}RI;\end{array}$(3) ) are given in [36] and most estimates of parameters on the right-hand side of (4(4) $\begin{array}{rl}& \frac{\mathrm{\partial }L}{\mathrm{\partial }t}-D_L\mathrm{\Delta }L=-k_{L}RL+{\lambda }_{GL}(G-G_s)L-{\lambda }_{ML}\frac{ML}{K_{ML}+L},\end{array}$(4) )–(13(13) $\begin{array}{rl}& -\frac{\mathrm{\partial }M}{\mathrm{\partial }n}+{\alpha }_M(M-{\overline{M}}_0)=0,-\frac{\mathrm{\partial }\beta }{\mathrm{\partial }n}+{\alpha }_M(\beta -{\beta }_0)=0,-\frac{\mathrm{\partial }F}{\mathrm{\partial }n}+{\alpha }_{M}F=0,\end{array}$(13) ) can be found given in [11,17]. See Table 3 and Table 4.

Table 3. Parameters' description and value.

Parameters	Description	Value
${\tilde{G}}_0$	net rate of production	$864{\mathrm{mg}\mathrm{dl}}^{-1}{\mathrm{d}}^{-1}$ [36]
	at zero glucose
$E_{G0}$	total glucose effectiveness	$1.44{\mathrm{d}}^{-1}$ [36]
	at zero insulin
$S_I$	total insulin sensitivity	$0.72\mathrm{ml}\mu U^{-1}{\mathrm{d}}^{-1}$ [36]
$d_0$	death rate at zero glucose	$0.06{\mathrm{d}}^{-1}$ [36]
	for β cell
$r_1$	rate constant	$0.84\times {10}^{-3}{\mathrm{mg}}^{-1}\mathrm{dl}{\mathrm{d}}^{-1}$ [36]
$r_2$	rate constant	$0.24\times {10}^{-3}{\mathrm{mg}}^{-1}{\mathrm{dl}}^2{\mathrm{d}}^{-1}$ [36]
σ	max secretion rate for insulin	$43.2\mu U{\mathrm{ml}}^{-1}{\mathrm{d}}^{-1}$[36]
α	constant in hill function	$20000{\mathrm{mg}}^2{\mathrm{dl}}^{-2}$ [36]
	with coefficient 2
$k_I$	clearance rate for insulin	$432{\mathrm{d}}^{-1}$ [36]
$D_G$	diffusion coefficient	$1.04\times {10}^{-1}{\mathrm{cm}}^2{\mathrm{d}}^{-1}$ [17]
	for glucose
$D_I$	diffusion coefficient	$1.042{\mathrm{cm}}^2{\mathrm{d}}^{-1}$ estimated
	for Insulin
$D_{\beta }$	diffusion coefficient	$8.64\times {10}^{-7}{\mathrm{cm}}^2{\mathrm{d}}^{-1}$ [17]
	for β cell
$D_R$	diffusion coefficient	$2.05\times {10}^{-1}{\mathrm{cm}}^2{\mathrm{d}}^{-1}$ [17]
	for radicals
$D_L$	diffusion coefficient	$29.89{\mathrm{cm}}^2{\mathrm{d}}^{-1}$ [17]
	for LDL
$D_H$	diffusion coefficient	$3.93{\mathrm{cm}}^2{\mathrm{d}}^{-1}$ [17]
	for HDL

Table 4. Parameters' description and value.

Parameters	Description	Value
$D_M$	diffusion coefficient	$8.64\times {10}^{-7}{\mathrm{cm}}^2{\mathrm{d}}^{-1}$ [17]
	for macrophage
$D_F$	diffusion coefficient	$8.64\times {10}^{-7}{\mathrm{cm}}^2{\mathrm{d}}^{-1}$ [17]
	for foam cells
$k_L$	reaction rate of	$2.35\times {10}^{-4}{\mathrm{g}}^{-1}{\mathrm{cm}}^3{\mathrm{d}}^{-1}$ [17]
	for LDL+ Radical
$k_H$	reaction rate of	$5.29\times {10}^{-6}{\mathrm{g}}^{-1}{\mathrm{cm}}^3{\mathrm{d}}^{-1}$ [17]
	for HDL+ Radical
${\lambda }_{GR}$	growth rate of radical	$124.69{\mathrm{d}}^{-1}$ estimated
	due to excess glucose
${\lambda }_{IR}$	growth rate of radical	$1253.055{\mathrm{d}}^{-1}$ estimated
	due to excess insulin
${\lambda }_{GH}$	rate of reduction	$860{\mathrm{cm}}^3{\mathrm{g}}^{-1}{\mathrm{d}}^{-1}$ estimated
	due to HDL
${\lambda }_{GL}$	rate of production	$860{\mathrm{cm}}^3{\mathrm{g}}^{-1}{\mathrm{d}}^{-1}$ estimated
	due to excess glucose
λ	parameter in Equation (7(7) $\begin{array}{rl}& \frac{\mathrm{\partial }M}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}M)-D_M\mathrm{\Delta }M=-{\lambda }_{ML}\frac{ML}{K_{ML}+L}+{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}+\lambda \frac{ML}{K_{MH}+H}-d_{M}M,\end{array}$(7) ) for M	$2.573\times {10}^{-3}{\mathrm{d}}^{-1}$ [12]
${\lambda }_{ML}$	rate of LDL ingestion	$144.5{\mathrm{d}}^{-1}$ [17]
	by macrophages
${\lambda }_{HF}$	rate of HDL ingestion	$10{\mathrm{d}}^{-1}$ [17]
	by foam cells
$d_M$	death rate of macrophages	$0.015{\mathrm{d}}^{-1}$ [17]
$d_F$	death rate of foam cells	$0.03{\mathrm{d}}^{-1}$ [17]
$K_{ML}$	LDL saturation for	${10}^{-2}{\mathrm{g}\mathrm{cm}}^{-3}$ [17]
	production of macrophages
$K_{FH}$	HDL saturation for	$0.5{\mathrm{g}\mathrm{cm}}^{-3}$ [17]
	production of foam cells
$K_{MH}$	parameter in Equation (7(7) $\begin{array}{rl}& \frac{\mathrm{\partial }M}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}M)-D_M\mathrm{\Delta }M=-{\lambda }_{ML}\frac{ML}{K_{ML}+L}+{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}+\lambda \frac{ML}{K_{MH}+H}-d_{M}M,\end{array}$(7) )	$-2.541\times {10}^{-3}{\mathrm{g}\mathrm{cm}}^{-3}$ [12]
$R_0$	source/influx of	$0.25\mathrm{g}{\mathrm{cm}}^{-3}{\mathrm{d}}^{-1}$ [17]
	radical into intima
${\alpha }_G$	influx rate of	$1.0{\mathrm{cm}}^{-1}$ estimated
	glucose into intima
${\alpha }_I$	influx rate of	$1.0{\mathrm{cm}}^{-1}$ estimated
	insulin into intima
${\alpha }_L$	influx rate of	$1.0{\mathrm{cm}}^{-1}$ [17]
	LDL into intima
${\alpha }_H$	influx rate of	$1.0{\mathrm{cm}}^{-1}$ [17]
	HDL into intima
${\alpha }_{\beta }$	influx rate of	$0.2{\mathrm{cm}}^{-1}$ [17]
	β cell into intima
${\alpha }_M$	influx rate of	$0.2{\mathrm{cm}}^{-1}$ [17]
	Macrophages into intima

To estimate some of the parameters in the equations for proteins, we shall use the concept of accessible surface area [5,20] of a protein p, or briefly $A_p$, which is roughly the minimum surface area of the smooth shapes containing the protein. (17) $\begin{array}{r}{A}_p=11.1M_p^{2/3},\end{array}$(17) where $M_p$ is the molecular weight of the protein. The molecular weights of insulin, glucose and radical are: $M_I=5734\mathrm{g}/\mathrm{mol}(\mathrm{daltons}),M_G=180\mathrm{daltons},M_R=500\mathrm{daltons}.$From (18(18) $\begin{array}{rl}{A}_I& =3555.9974\times {10}^{-16}{\mathrm{cm}}^2,A_G=353.8653\times {10}^{-16}{\mathrm{cm}}^2,\\ A_R& =699.2562\times {10}^{-16}{\mathrm{cm}}^2.\end{array}$(18) ), we obtain (18) $\begin{array}{rl}{A}_I& =3555.9974\times {10}^{-16}{\mathrm{cm}}^2,A_G=353.8653\times {10}^{-16}{\mathrm{cm}}^2,\\ A_R& =699.2562\times {10}^{-16}{\mathrm{cm}}^2.\end{array}$(18) Diffusion coefficients. We assume that the diffusion coefficients of all the cells are the same, and take them to be $8.64\times {10}^{-7}{\mathrm{cm}}^2{\mathrm{d}}^{-1}$ [9,10]. In order to estimate the diffusion coefficients of the various proteins, we assume that the diffusion coefficient of protein p, $D_p$, is proportional to its area $A_p$, i.e. $D_p\sim K{A}_p$, where we take $K=2.93\times {10}^{12}{\mathrm{d}}^{-1}$ to be the same for all small molecules. Using (19(19) $\begin{array}{r}{\lambda }_{IR}=\frac{A_I}{A_G}{\lambda }_{GR}=10.049{\lambda }_{GR}.\end{array}$(19) ), we have $D_I=1.04{\mathrm{cm}}^2{\mathrm{d}}^{-1},D_G=0.104{\mathrm{cm}}^2{\mathrm{d}}^{-1}$.

Production rates: production rates in (7(7) $\begin{array}{rl}& \frac{\mathrm{\partial }M}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}M)-D_M\mathrm{\Delta }M=-{\lambda }_{ML}\frac{ML}{K_{ML}+L}+{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}+\lambda \frac{ML}{K_{MH}+H}-d_{M}M,\end{array}$(7) ) will be estimated from the steady state and the data in the existing literature. First from (17(17) $\begin{array}{r}{A}_p=11.1M_p^{2/3},\end{array}$(17) ), we have (19) $\begin{array}{r}{\lambda }_{IR}=\frac{A_I}{A_G}{\lambda }_{GR}=10.049{\lambda }_{GR}.\end{array}$(19) From (7(7) $\begin{array}{rl}& \frac{\mathrm{\partial }M}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}M)-D_M\mathrm{\Delta }M=-{\lambda }_{ML}\frac{ML}{K_{ML}+L}+{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}+\lambda \frac{ML}{K_{MH}+H}-d_{M}M,\end{array}$(7) ) we have (20) $\begin{array}{r}{R}_0={\lambda }_{IR}{I}_{\mathrm{steady}}+{\lambda }_{GR}{G}_{\mathrm{steady}},\end{array}$(20) where $I_{\mathrm{steady}}={10}^{-4}\mathrm{g}{\mathrm{cm}}^{-3},G_{\mathrm{steady}}={10}^{-3}{\mathrm{g}\mathrm{cm}}^{-3}$ from Table 2 and $R_0=0.25$ from Table 4.

Solving (19(19) $\begin{array}{r}{\lambda }_{IR}=\frac{A_I}{A_G}{\lambda }_{GR}=10.049{\lambda }_{GR}.\end{array}$(19) ) and (20(20) $\begin{array}{r}{R}_0={\lambda }_{IR}{I}_{\mathrm{steady}}+{\lambda }_{GR}{G}_{\mathrm{steady}},\end{array}$(20) ), we obtain (21) $\begin{array}{r}{\lambda }_{GR}=124.69{\mathrm{d}}^{-1},{\lambda }_{IR}=1253.055{\mathrm{d}}^{-1}.\end{array}$(21) Since the death rate of SMC is $0.86{\mathrm{d}}^{-1}$ [17], we estimate that ${\lambda }_{GL}=0.86\times {10}^3{\mathrm{d}}^{-1}{\mathrm{cm}}^3{g}^{-1}$ in Equation  (8(8) $\begin{array}{rl}& \frac{\mathrm{\partial }F}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}F)-D_F\mathrm{\Delta }F={\lambda }_{ML}\frac{ML}{K_{ML}+L}-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}-d_{F}F,\end{array}$(8) ).

We take the flux rates ${\alpha }_G={\alpha }_I={\alpha }_L={\alpha }_H=1.0{\mathrm{cm}}^{-1}$ and ${\alpha }_{\beta }={\alpha }_M=0.2{\mathrm{cm}}^{-1}$ [17].

3. Stationary solutions

In this section, we are going to prove that there exist stationary small plaque solutions $B(t)=1-ϵ$ for any sufficiently small ϵ. Dropping the time derivatives, the stationary solution satisfies (22) $\begin{array}{rl}& -D_G\mathrm{\Delta }G={\tilde{G}}_0-(E_{G0}+S_{I}I)G-{\lambda }_{GH}(G-G_s)H,\end{array}$(22) (23) $\begin{array}{rl}& \mathrm{\nabla }\cdot (\mathbf{u}\beta )-D_{\beta }\mathrm{\Delta }\beta =(-d_0+r_{1}G-r_2{G}^2)\beta ,\end{array}$(23) (24) $\begin{array}{rl}& -D_I\mathrm{\Delta }I=\frac{\sigma G^2\beta }{\alpha +G^2}-k_{I}I-{\lambda }_{RI}RI;\end{array}$(24) (25) $\begin{array}{rl}& -D_L\mathrm{\Delta }L=-k_{L}RL+{\lambda }_{GL}(G-G_s)L-{\lambda }_{ML}\frac{ML}{K_{ML}+L},\end{array}$(25) (26) $\begin{array}{rl}& -D_H\mathrm{\Delta }H=-k_{H}RH-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)},\end{array}$(26) (27) $\begin{array}{rl}& -D_R\mathrm{\Delta }R=R_0-R(k_{L}L+k_{H}H)+{\lambda }_{GR}(G-G_s)R+{\lambda }_{IR}IR,\end{array}$(27) (28) $\begin{array}{rl}& \mathrm{\nabla }\cdot (\mathbf{u}F)-D_F\mathrm{\Delta }F={\lambda }_{ML}\frac{(M_0-F-\beta )L}{K_{ML}+L}-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}-d_{F}F,\end{array}$(28) (29) $\begin{array}{rl}& M_0\frac{1}{r}{\mathrm{\partial }}_r(ru)=(-d_0+r_{1}G-r_2{G}^2)\beta +\lambda \frac{(M_0-F-\beta )L}{K_{MH}+H}-d_{F}F-d_M(M_0-F-\beta ).\end{array}$(29) No flux condition on the boundary r = 1: (30) $\begin{array}{r}\frac{\mathrm{\partial }L}{\mathrm{\partial }n}=\frac{\mathrm{\partial }H}{\mathrm{\partial }n}=\frac{\mathrm{\partial }G}{\mathrm{\partial }n}=\frac{\mathrm{\partial }I}{\mathrm{\partial }n}=\frac{\mathrm{\partial }R}{\mathrm{\partial }n}=\frac{\mathrm{\partial }\beta }{\mathrm{\partial }n}=\frac{\mathrm{\partial }F}{\mathrm{\partial }n}=u(1)=0.\end{array}$(30) The boundary conditions on $r=1-ϵ$ are (31) $\begin{array}{rl}& -\frac{\mathrm{\partial }L}{\mathrm{\partial }n}+\alpha (L-L_0)=0,-\frac{\mathrm{\partial }H}{\mathrm{\partial }n}+\alpha (H-H_0)=0,\end{array}$(31) (32) $\begin{array}{rl}& -\frac{\mathrm{\partial }G}{\mathrm{\partial }n}+{\alpha }_G(G-G_0)=0,-\frac{\mathrm{\partial }I}{\mathrm{\partial }n}+{\alpha }_I(I-I_0)=0,-\frac{\mathrm{\partial }R}{\mathrm{\partial }n}+{\alpha }_R(R-R_0)=0,\end{array}$(32) (33) $\begin{array}{rl}& -\frac{\mathrm{\partial }\beta }{\mathrm{\partial }n}+{\alpha }_M(\beta -{\beta }_0)=0,-\frac{\mathrm{\partial }F}{\mathrm{\partial }n}+{\alpha }_{M}F=0,\end{array}$(33) (34) $\begin{array}{rl}& u(1-ϵ)=0.\end{array}$(34)

Lemma 3.1

For sufficiently small $ϵ>0$, there exists a unique solution

$(G,I,\beta ,R,L,H,F,u)$ to the problem (22(22) $\begin{array}{rl}& -D_G\mathrm{\Delta }G={\tilde{G}}_0-(E_{G0}+S_{I}I)G-{\lambda }_{GH}(G-G_s)H,\end{array}$(22) )–(33(33) $\begin{array}{rl}& -\frac{\mathrm{\partial }\beta }{\mathrm{\partial }n}+{\alpha }_M(\beta -{\beta }_0)=0,-\frac{\mathrm{\partial }F}{\mathrm{\partial }n}+{\alpha }_{M}F=0,\end{array}$(33) ).

Proof.

We define a Banach space $\begin{array}{rl}\mathcal{B}& =\{(G,I,\beta ,R,L,H,F,u)\in {(C[1-ϵ,1])}^8,G_s\le G\le G_0,0\le I\le I_0,\\ 0& \le R\le R_0,0\le L\le L_0,0\le H\le H_0,0\le \beta \le {\beta }_0,0\le F\le M_0\}\end{array}$where $C[1-ϵ,1]$ is the space of continuous functions on $[1-ϵ,1]$.

We define a map $\mathcal{M}$ from $\mathcal{B}\to \mathcal{B}:(G,I,\beta ,R,L,H,F,u)\to (\overline{G},\overline{I},\overline{\beta },\overline{R},\overline{L},\overline{H},\overline{F},\overline{u})$, where $(\overline{G},\overline{I},\overline{\beta },\overline{R},\overline{L},\overline{H},\overline{F},\overline{u})$ is given by (35) $\begin{array}{rl}& -D_G\mathrm{\Delta }\overline{G}={\tilde{G}}_0-(E_{G0}+S_{I}I)\overline{G}-{\lambda }_{GH}(\overline{G}-G_s)H,\end{array}$(35) (36) $\begin{array}{rl}& \mathrm{\nabla }\cdot (\mathbf{u}\overline{\beta })-D_{\beta }\mathrm{\Delta }\overline{\beta }=(-d_0+r_{1}G-r_2{G}^2)\overline{\beta },\end{array}$(36) (37) $\begin{array}{rl}& -D_I\mathrm{\Delta }\overline{I}=\frac{\sigma G^2\beta }{\alpha +G^2}-k_I\overline{I}-{\lambda }_{RI}R\overline{I};\end{array}$(37) (38) $\begin{array}{rl}& -D_L\mathrm{\Delta }\overline{L}=-k_{L}R\overline{L}+{\lambda }_{GL}(G-G_s)\overline{L}-{\lambda }_{ML}\frac{(M_0-F-\beta )\overline{L}}{K_{ML}+L},\end{array}$(38) (39) $\begin{array}{rl}& -D_H\mathrm{\Delta }\overline{H}=-k_{H}R\overline{H}-{\lambda }_{HF}\frac{\overline{H}F}{(K_{HF}+H)},\end{array}$(39) (40) $\begin{array}{rl}& -D_R\mathrm{\Delta }\overline{R}={\tilde{R}}_0-\overline{R}(k_{L}L+k_{H}H)+{\lambda }_{GR}(G-G_s)\overline{R}+{\lambda }_{IR}I\overline{R},\end{array}$(40) (41) $\begin{array}{rl}& \mathrm{\nabla }\cdot (\mathbf{u}\overline{F})-D_F\mathrm{\Delta }\overline{F}={\lambda }_{ML}\frac{(M_0-\overline{F}-\beta )L}{K_{ML}+L}-{\lambda }_{HF}\frac{H\overline{F}}{(K_{HF}+H)}-d_F\overline{F},\end{array}$(41) (42) $\begin{array}{rl}& M_0\frac{1}{r}{\mathrm{\partial }}_r(r\overline{u})=(-d_0+r_{1}G-r_2{G}^2)\beta +\lambda \frac{(M_0-F-\beta )L}{K_{MH}+H}-d_{F}F-d_M(M_0-F-\beta ).\end{array}$(42) No flux condition on the boundary r = 1: (43) $\begin{array}{r}\frac{\mathrm{\partial }\overline{L}}{\mathrm{\partial }n}=\frac{\mathrm{\partial }\overline{H}}{\mathrm{\partial }n}=\frac{\mathrm{\partial }\overline{G}}{\mathrm{\partial }n}=\frac{\mathrm{\partial }\overline{I}}{\mathrm{\partial }n}=\frac{\mathrm{\partial }\overline{R}}{\mathrm{\partial }n}=\frac{\mathrm{\partial }\overline{\beta }}{\mathrm{\partial }n}=\frac{\mathrm{\partial }\overline{F}}{\mathrm{\partial }n}=\overline{u}(1)=0.\end{array}$(43) The boundary conditions on $r=1-ϵ$ are (44) $\begin{array}{rl}& -\frac{\mathrm{\partial }\overline{L}}{\mathrm{\partial }n}+\alpha (\overline{L}-L_0)=0,-\frac{\mathrm{\partial }\overline{H}}{\mathrm{\partial }n}+\alpha (\overline{H}-H_0)=0,\end{array}$(44) (45) $\begin{array}{rl}& -\frac{\mathrm{\partial }\overline{G}}{\mathrm{\partial }n}+{\alpha }_G(\overline{G}-G_0)=0,-\frac{\mathrm{\partial }\overline{I}}{\mathrm{\partial }n}+{\alpha }_I(\overline{I}-I_0)=0,-\frac{\mathrm{\partial }\overline{R}}{\mathrm{\partial }n}+{\alpha }_R(\overline{R}-R_0)=0,\end{array}$(45) (46) $\begin{array}{rl}& -\frac{\mathrm{\partial }\overline{\beta }}{\mathrm{\partial }n}+{\alpha }_M(\overline{\beta }-{\beta }_0)=0,-\frac{\mathrm{\partial }\overline{F}}{\mathrm{\partial }n}+{\alpha }_M\overline{F}=0,\end{array}$(46) From Comparison Theorem, we can see that $G_s\le \overline{G}\le G_0,0\le \overline{I}\le I_0,0\le \overline{R}\le R_0,0\le \overline{L}\le L_0,0\le \overline{H}\le H_0,0\le \overline{\beta }\le {\beta }_0,0\le \overline{F}\le M_0$, hence $(\overline{G},\overline{I},\overline{\beta },\overline{R},\overline{L},\overline{H},\overline{F},\overline{u})\in \mathcal{B}.$ To prove that $\mathcal{M}$ is a contraction, we introduce function: (47) $\begin{array}{r}\eta (r)=\frac{1-r^2}{4}+\frac{1}{2}\log r;\end{array}$(47) then (48) $\begin{array}{r}\eta (1)={\eta }_r(1)=0,\eta <0;\mathrm{\Delta }\eta =-1,\eta =O({ϵ}^2),{\eta }_r=O(ϵ)\text{ as }ϵ\to 0.\end{array}$(48) Let $C(\gamma ,ϵ)$ be a number satisfying (49) $\begin{array}{r}-\frac{\mathrm{\partial }(\eta +C(\gamma ,ϵ)}{\mathrm{\partial }n}+\gamma (\eta +C(\gamma ,ϵ))=0,\text{ on }r=1-ϵ;\end{array}$(49) then (50) $\begin{array}{r}\eta +C(\gamma ,ϵ)>0,C(\gamma ,ϵ)=\frac{ϵ}{\gamma }+O({ϵ}^2).\end{array}$(50) Let $\begin{array}{rl}& ({\overline{G}}_j,{\overline{I}}_j,{\overline{\beta }}_j,{\overline{R}}_j,{\overline{L}}_j,{\overline{H}}_j,{\overline{F}}_j,{\overline{u}}_j)=\mathcal{M}(G_j,I_j,{\beta }_j,R_j,L_j,H_j,F_j,u_j),\\ & A=|G_1-G_2|+|I_1-I_2|+|{\beta }_1-{\beta }_2|+|R_1-R_2|\\ & +|L_1-L_2|+|H_1-H_2|+|F_1-F_2|;\\ & B=|{\overline{G}}_1-{\overline{G}}_2|+|{\overline{I}}_1-{\overline{I}}_2|+|{\overline{\beta }}_1-{\overline{\beta }}_2|+|{\overline{R}}_1-{\overline{R}}_2|\\ & +|{\overline{L}}_1-{\overline{L}}_2|+|{\overline{H}}_1-{\overline{H}}_2|+|{\overline{F}}_1-{\overline{F}}_2|;\end{array}$then for some constant C, $\begin{array}{rl}& |\mathrm{\Delta }({\overline{G}}_1-{\overline{G}}_2)|\le C(A+B);\\ & |\mathrm{\Delta }({\overline{I}}_1-{\overline{I}}_2)|\le C(A+B);\\ & |\mathrm{\Delta }({\overline{\beta }}_1-{\overline{\beta }}_2)-\mathrm{\nabla }\cdot (\mathbf{u}({\overline{\beta }}_1-{\overline{\beta }}_2))|\le C(A+B);\\ & |\mathrm{\Delta }({\overline{R}}_1-{\overline{R}}_2)|\le C(A+B);\\ & |\mathrm{\Delta }({\overline{L}}_1-{\overline{L}}_2)|\le C(A+B);\\ & |\mathrm{\Delta }({\overline{H}}_1-{\overline{H}}_2)|\le C(A+B);\\ & |\mathrm{\Delta }({\overline{F}}_1-{\overline{F}}_2)-\mathrm{\nabla }\cdot (\mathbf{u}({\overline{F}}_1-{\overline{F}}_2))|\le C(A+B);\end{array}$and use the maximum principle to get $\begin{array}{rl}& |({\overline{G}}_1-{\overline{G}}_2)|\le C(A+B)(\eta +C({\alpha }_G,ϵ));\\ & |({\overline{I}}_1-{\overline{I}}_2)|\le C(A+B)(\eta +C({\alpha }_I,ϵ));\\ & |({\overline{\beta }}_1-{\overline{\beta }}_2)|\le C(A+B)(\eta +C({\alpha }_M,ϵ));\\ & |({\overline{R}}_1-{\overline{R}}_2)|\le C(A+B)(\eta +C({\alpha }_R,ϵ));\\ & |({\overline{L}}_1-{\overline{L}}_2)|\le C(A+B)(\eta +C({\alpha }_L,ϵ));\\ & |({\overline{H}}_1-{\overline{H}}_2)|\le C(A+B)(\eta +C({\alpha }_H,ϵ));\\ & |({\overline{F}}_1-{\overline{F}}_2)-\mathrm{\nabla }\cdot (\mathbf{u}({\overline{F}}_1-{\overline{F}}_2))|\le C(A+B)(\eta +C({\alpha }_M,ϵ).\end{array}$Using the above, we have $B\le C(A+B)(C({\alpha }_G,ϵ)+C({\alpha }_I,ϵ)+2C({\alpha }_M,ϵ)+C({\alpha }_R,ϵ)+C({\alpha }_L,ϵ)+C({\alpha }_H,ϵ))$From (50(50) $\begin{array}{r}\eta +C(\gamma ,ϵ)>0,C(\gamma ,ϵ)=\frac{ϵ}{\gamma }+O({ϵ}^2).\end{array}$(50) ) and above, we can see $B<C^{\ast }A$ for some constant $C^{\ast }<1$, hence the mapping $\mathcal{M}$ is a contraction and the lemma holds.

As in [12], we assume that (51) $\begin{array}{r}\frac{\lambda L_0}{K_{MH}+H_0}-d_M=O(ϵ).\end{array}$(51) We assume that ${\tilde{G}}_s$ is the root of $-d_0+r_1{G}^2-r_2{G}^2=0$ such that ${\tilde{G}}_s>G_s$.

Theorem 3.2

For ϵ small, given $M_0>{\beta }_0>0$, $L_0$ and $H_0$ satisfy (51(51) $\begin{array}{r}\frac{\lambda L_0}{K_{MH}+H_0}-d_M=O(ϵ).\end{array}$(51) ). Then there exist $G_0$ satisfying $G_0={\tilde{G}}_s+O(ϵ)$ such that the solution in Lemma 1 satisfies $u(1-ϵ)=0$.

Proof.

Since $\begin{array}{rl}& |\mathrm{\Delta }(G-G_0)|\le C;\\ & |\mathrm{\Delta }(I-I_0)|\le C;\\ & |\mathrm{\Delta }({\beta }_1-{\beta }_0)-\mathrm{\nabla }\cdot (\mathbf{u}(\beta -{\beta }_0))|\le C;\\ & |\mathrm{\Delta }(R-R_0)|\le C;\\ & |\mathrm{\Delta }(L-L_0)|\le C;\\ & |\mathrm{\Delta }(H-H_0)|\le C;\\ & |\mathrm{\Delta }F-\mathrm{\nabla }\cdot (\mathbf{u}F)|\le C;\end{array}$and use the maximum principle to get $\begin{array}{rl}& |(G-G_0)|\le C(\eta +C({\alpha }_G,ϵ));\\ & |(I-I_0)|\le C(\eta +C({\alpha }_I,ϵ));\\ & |({\beta }_1-{\beta }_0)|\le C(\eta +C({\alpha }_M,ϵ));\\ & |(R-R_0)|\le C(\eta +C({\alpha }_R,ϵ));\\ & |(L-L_0)|\le C;(\eta +C({\alpha }_L,ϵ));\\ & |(H-H_0)|\le C(\eta +C({\alpha }_H,ϵ));\\ & |F|\le C(\eta +C({\alpha }_M,ϵ)).\end{array}$We obtain from (42(42) $\begin{array}{rl}& M_0\frac{1}{r}{\mathrm{\partial }}_r(r\overline{u})=(-d_0+r_{1}G-r_2{G}^2)\beta +\lambda \frac{(M_0-F-\beta )L}{K_{MH}+H}-d_{F}F-d_M(M_0-F-\beta ).\end{array}$(42) ) (52) $\begin{array}{r}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)=(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)+{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)+O(ϵ).\end{array}$(52) Using (51(51) $\begin{array}{r}\frac{\lambda L_0}{K_{MH}+H_0}-d_M=O(ϵ).\end{array}$(51) ) we have (53) $\begin{array}{r}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)={\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)+O(ϵ)\end{array}$(53) which implies that (54) $\begin{array}{r}{M}_0(1-ϵ)u(1-ϵ)=-{\int }_{1-ϵ}^{1}r{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)\mathrm{d}r+O({ϵ}^2).\end{array}$(54) The above implies that $u(1-ϵ)>0$ for $G_0>{\tilde{G}}_s$ while $u(1-ϵ)<0$ for $G_0<{\tilde{G}}_s$. So there is $G_0={\tilde{G}}_s+O(ϵ)$ such that $u(1-ϵ)=0$.

We now proceed to prove some sharper estimates for the steady solution, these estimates are useful in the future sections. From (22(22) $\begin{array}{rl}& -D_G\mathrm{\Delta }G={\tilde{G}}_0-(E_{G0}+S_{I}I)G-{\lambda }_{GH}(G-G_s)H,\end{array}$(22) )–(28(28) $\begin{array}{rl}& \mathrm{\nabla }\cdot (\mathbf{u}F)-D_F\mathrm{\Delta }F={\lambda }_{ML}\frac{(M_0-F-\beta )L}{K_{ML}+L}-{\lambda }_{HF}\frac{HF}{(K_{HF}+H)}-d_{F}F,\end{array}$(28) ), we have (55) $\begin{array}{rl}& |\mathrm{\Delta }(G-[G_0+({\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GH}(G_0-G_s)H_0)\eta (r)])|\le O(ϵ);\\ & \left|\mathrm{\Delta }(I-[I_0+(\frac{\sigma G_0^2{\beta }_0}{\alpha +G_0^2}-k_I{I}_0-{\lambda }_{RI}{R}_0{I}_0)\eta (r)])\right|\le O(ϵ);\\ & |\mathrm{\Delta }({\beta }_1-[{\beta }_0+(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0\eta (r)])|\le O(ϵ));\\ & \left|\mathrm{\Delta }(R-[R_0+({\tilde{R}}_0-R_0(k_L{L}_0+k_H{H}_0)+{\lambda }_{GR}(G_0-G_s)R_0+{\lambda }_{IR}{I}_0{R}_0)\eta (r)])\right|\le O(ϵ);\\ & \left|\mathrm{\Delta }(L-[L_0+(-k_L{R}_0{L}_0+{\lambda }_{GL}(G_0-G_s)L_0-{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0})\eta (r)])\right|\le O(ϵ);\\ & \left|\mathrm{\Delta }(H-[H_0+(-k_H{R}_0{H}_0-{\lambda }_{HF}\frac{H_0{F}_0}{(K_{HF}+H_0)})\eta (r)])\right|\le O(ϵ);\\ & \left|\mathrm{\Delta }(F-\left({\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right)\eta (r))\right|\le O(ϵ).\end{array}$(55) Integrating above to obtain (56) $\begin{array}{rl}& |{\mathrm{\partial }}_r(G-[G_0+({\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GH}(G_0-G_s)H_0)\eta (r)])|\le O({ϵ}^2);\\ & \left|\mathrm{\Delta }(I-[I_0+(\frac{\sigma G_0^2{\beta }_0}{\alpha +G_0^2}-k_I{I}_0-{\lambda }_{RI}{R}_0{I}_0)\eta (r)])\right|\le O({ϵ}^2);\\ & |{\mathrm{\partial }}_r({\beta }_1-[{\beta }_0+(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0\eta (r)])|\le O({ϵ}^2));\\ & \left|{\mathrm{\partial }}_r(R-[R_0+({\tilde{R}}_0-R_0(k_L{L}_0+k_H{H}_0)+{\lambda }_{GR}(G_0-G_s)R_0+{\lambda }_{IR}{I}_0{R}_0)\eta (r)])\right|\le O({ϵ}^2);\\ & \left|{\mathrm{\partial }}_r(L-[L_0+(-k_L{R}_0{L}_0+{\lambda }_{GL}(G_0-G_s)L_0-{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0})\eta (r)])\right|\le O({ϵ}^2);\\ & \left|{\mathrm{\partial }}_r(H-[H_0+(-k_H{R}_0{H}_0-{\lambda }_{HF}\frac{H_0{F}_0}{(K_{HF}+H_0)})\eta (r)])\right|\le O({ϵ}^2);\\ & \left|{\mathrm{\partial }}_r(F-\left({\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right)\eta (r))\right|\le O({ϵ}^2).\end{array}$(56) Integrating the above and using (49(49) $\begin{array}{r}-\frac{\mathrm{\partial }(\eta +C(\gamma ,ϵ)}{\mathrm{\partial }n}+\gamma (\eta +C(\gamma ,ϵ))=0,\text{ on }r=1-ϵ;\end{array}$(49) ) and (50(50) $\begin{array}{r}\eta +C(\gamma ,ϵ)>0,C(\gamma ,ϵ)=\frac{ϵ}{\gamma }+O({ϵ}^2).\end{array}$(50) ), we can write (57) $\begin{array}{rl}G& =[G_0+({\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GH}(G_0-G_s)H_0)(\eta (r)+\frac{ϵ}{{\alpha }_G})]+G^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:G_0+G_1ϵ+G_2{ϵ}^2+O({ϵ}^3),\\ I& =[I_0+(\frac{\sigma G_0^2{\beta }_0}{\alpha +G_0^2}-k_I{I}_0-{\lambda }_{RI}{R}_0{I}_0)(\eta (r)+\frac{ϵ}{{\alpha }_I})]+I^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:I_0+I_1ϵ+L_2{ϵ}^2+O({ϵ}^3),\\ \beta & =[{\beta }_0+(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0(\eta (r)+\frac{ϵ}{{\alpha }_M})]+{\beta }^{\ast }{ϵ}^2+O({ϵ}^3));\\ R& =[R_0+({\tilde{R}}_0-R_0(k_L{L}_0+k_H{H}_0)+{\lambda }_{GR}(G_0-G_s)R_0+{\lambda }_{IR}{I}_0{R}_0)(\eta (r)++\frac{ϵ}{{\alpha }_R})]\\ & +R^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:R_0+R_1ϵ+R_2{ϵ}^2+O({ϵ}^3)\\ L& =[L_0+(-k_L{R}_0{L}_0+{\lambda }_{GL}(G_0-G_s)L_0-{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0})(\eta (r)+\frac{ϵ}{{\alpha }_L})]\\ & +L^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:L_0+L_1ϵ+L_2{ϵ}^2+O({ϵ}^3)\\ H& =[H_0+(-k_H{R}_0{H}_0-{\lambda }_{HF}\frac{H_0{F}_0}{(K_{HF}+H_0)})(\eta (r)+\frac{ϵ}{{\alpha }_H})]+H^{\ast }{ϵ}^2+O({ϵ}^3);\\ & =:H_0+H_1ϵ+H_2{ϵ}^2+O({ϵ}^3),\\ F& =\left({\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right)(\eta (r)+\frac{ϵ}{{\alpha }_M})+F^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:F_1ϵ+F_2{ϵ}^2+O({ϵ}^3),\end{array}$(57) where the constants $G^{\ast },I^{\ast },{\beta }^{\ast },L^{\ast },H^{\ast }$ and $F^{\ast }$ are of order $O(1)$.

Plugging (57(57) $\begin{array}{rl}G& =[G_0+({\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GH}(G_0-G_s)H_0)(\eta (r)+\frac{ϵ}{{\alpha }_G})]+G^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:G_0+G_1ϵ+G_2{ϵ}^2+O({ϵ}^3),\\ I& =[I_0+(\frac{\sigma G_0^2{\beta }_0}{\alpha +G_0^2}-k_I{I}_0-{\lambda }_{RI}{R}_0{I}_0)(\eta (r)+\frac{ϵ}{{\alpha }_I})]+I^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:I_0+I_1ϵ+L_2{ϵ}^2+O({ϵ}^3),\\ \beta & =[{\beta }_0+(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0(\eta (r)+\frac{ϵ}{{\alpha }_M})]+{\beta }^{\ast }{ϵ}^2+O({ϵ}^3));\\ R& =[R_0+({\tilde{R}}_0-R_0(k_L{L}_0+k_H{H}_0)+{\lambda }_{GR}(G_0-G_s)R_0+{\lambda }_{IR}{I}_0{R}_0)(\eta (r)++\frac{ϵ}{{\alpha }_R})]\\ & +R^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:R_0+R_1ϵ+R_2{ϵ}^2+O({ϵ}^3)\\ L& =[L_0+(-k_L{R}_0{L}_0+{\lambda }_{GL}(G_0-G_s)L_0-{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0})(\eta (r)+\frac{ϵ}{{\alpha }_L})]\\ & +L^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:L_0+L_1ϵ+L_2{ϵ}^2+O({ϵ}^3)\\ H& =[H_0+(-k_H{R}_0{H}_0-{\lambda }_{HF}\frac{H_0{F}_0}{(K_{HF}+H_0)})(\eta (r)+\frac{ϵ}{{\alpha }_H})]+H^{\ast }{ϵ}^2+O({ϵ}^3);\\ & =:H_0+H_1ϵ+H_2{ϵ}^2+O({ϵ}^3),\\ F& =\left({\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right)(\eta (r)+\frac{ϵ}{{\alpha }_M})+F^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:F_1ϵ+F_2{ϵ}^2+O({ϵ}^3),\end{array}$(57) ) into (29(29) $\begin{array}{rl}& M_0\frac{1}{r}{\mathrm{\partial }}_r(ru)=(-d_0+r_{1}G-r_2{G}^2)\beta +\lambda \frac{(M_0-F-\beta )L}{K_{MH}+H}-d_{F}F-d_M(M_0-F-\beta ).\end{array}$(29) ), we have (58) $\begin{array}{rl}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)& =(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)+{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)\\ & +K_1ϵ+K_2{ϵ}^2+O({ϵ}^3),\end{array}$(58) where (59) $\begin{array}{r}{K}_1={\lambda }_{HF}\frac{(M_0-{\beta }_0)L_1}{(K_{HF}+H_0)}-\frac{L_0({\beta }_1+F_1)}{(K_{HF}+H_0)}-\frac{(M_0-{\beta }_0)L_0}{(K_{HF}+H_0)}\frac{H_1}{(K_{HF}+H_0)}\\ -d_M(F_1+{\beta }_1)-d_F{F}_1+{\beta }_0(r_1{G}_1-2r_2{G}_0{G}_1)+{\beta }_1(-d_0+r_1{G}_0^2-r_2{G}_0^2),\end{array}$(59) and $G_1.I_1.L_1.H_1,{\beta }_1$ anf $F_1$ are defined in (57(57) $\begin{array}{rl}G& =[G_0+({\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GH}(G_0-G_s)H_0)(\eta (r)+\frac{ϵ}{{\alpha }_G})]+G^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:G_0+G_1ϵ+G_2{ϵ}^2+O({ϵ}^3),\\ I& =[I_0+(\frac{\sigma G_0^2{\beta }_0}{\alpha +G_0^2}-k_I{I}_0-{\lambda }_{RI}{R}_0{I}_0)(\eta (r)+\frac{ϵ}{{\alpha }_I})]+I^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:I_0+I_1ϵ+L_2{ϵ}^2+O({ϵ}^3),\\ \beta & =[{\beta }_0+(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0(\eta (r)+\frac{ϵ}{{\alpha }_M})]+{\beta }^{\ast }{ϵ}^2+O({ϵ}^3));\\ R& =[R_0+({\tilde{R}}_0-R_0(k_L{L}_0+k_H{H}_0)+{\lambda }_{GR}(G_0-G_s)R_0+{\lambda }_{IR}{I}_0{R}_0)(\eta (r)++\frac{ϵ}{{\alpha }_R})]\\ & +R^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:R_0+R_1ϵ+R_2{ϵ}^2+O({ϵ}^3)\\ L& =[L_0+(-k_L{R}_0{L}_0+{\lambda }_{GL}(G_0-G_s)L_0-{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0})(\eta (r)+\frac{ϵ}{{\alpha }_L})]\\ & +L^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:L_0+L_1ϵ+L_2{ϵ}^2+O({ϵ}^3)\\ H& =[H_0+(-k_H{R}_0{H}_0-{\lambda }_{HF}\frac{H_0{F}_0}{(K_{HF}+H_0)})(\eta (r)+\frac{ϵ}{{\alpha }_H})]+H^{\ast }{ϵ}^2+O({ϵ}^3);\\ & =:H_0+H_1ϵ+H_2{ϵ}^2+O({ϵ}^3),\\ F& =\left({\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right)(\eta (r)+\frac{ϵ}{{\alpha }_M})+F^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:F_1ϵ+F_2{ϵ}^2+O({ϵ}^3),\end{array}$(57) ).Using the fact (60) $\begin{array}{r}{\int }_{1-ϵ}^1{u}_r(s)\mathrm{d}s=0\end{array}$(60) and (58(58) $\begin{array}{rl}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)& =(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)+{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)\\ & +K_1ϵ+K_2{ϵ}^2+O({ϵ}^3),\end{array}$(58) ) we conclude that (61) $\begin{array}{r}(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)+{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)+K_1ϵ=O({ϵ}^2).\end{array}$(61) Now using (29(29) $\begin{array}{rl}& M_0\frac{1}{r}{\mathrm{\partial }}_r(ru)=(-d_0+r_{1}G-r_2{G}^2)\beta +\lambda \frac{(M_0-F-\beta )L}{K_{MH}+H}-d_{F}F-d_M(M_0-F-\beta ).\end{array}$(29) ), (57(57) $\begin{array}{rl}G& =[G_0+({\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GH}(G_0-G_s)H_0)(\eta (r)+\frac{ϵ}{{\alpha }_G})]+G^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:G_0+G_1ϵ+G_2{ϵ}^2+O({ϵ}^3),\\ I& =[I_0+(\frac{\sigma G_0^2{\beta }_0}{\alpha +G_0^2}-k_I{I}_0-{\lambda }_{RI}{R}_0{I}_0)(\eta (r)+\frac{ϵ}{{\alpha }_I})]+I^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:I_0+I_1ϵ+L_2{ϵ}^2+O({ϵ}^3),\\ \beta & =[{\beta }_0+(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0(\eta (r)+\frac{ϵ}{{\alpha }_M})]+{\beta }^{\ast }{ϵ}^2+O({ϵ}^3));\\ R& =[R_0+({\tilde{R}}_0-R_0(k_L{L}_0+k_H{H}_0)+{\lambda }_{GR}(G_0-G_s)R_0+{\lambda }_{IR}{I}_0{R}_0)(\eta (r)++\frac{ϵ}{{\alpha }_R})]\\ & +R^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:R_0+R_1ϵ+R_2{ϵ}^2+O({ϵ}^3)\\ L& =[L_0+(-k_L{R}_0{L}_0+{\lambda }_{GL}(G_0-G_s)L_0-{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0})(\eta (r)+\frac{ϵ}{{\alpha }_L})]\\ & +L^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:L_0+L_1ϵ+L_2{ϵ}^2+O({ϵ}^3)\\ H& =[H_0+(-k_H{R}_0{H}_0-{\lambda }_{HF}\frac{H_0{F}_0}{(K_{HF}+H_0)})(\eta (r)+\frac{ϵ}{{\alpha }_H})]+H^{\ast }{ϵ}^2+O({ϵ}^3);\\ & =:H_0+H_1ϵ+H_2{ϵ}^2+O({ϵ}^3),\\ F& =\left({\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right)(\eta (r)+\frac{ϵ}{{\alpha }_M})+F^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:F_1ϵ+F_2{ϵ}^2+O({ϵ}^3),\end{array}$(57) ) and (61(61) $\begin{array}{r}(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)+{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)+K_1ϵ=O({ϵ}^2).\end{array}$(61) ), we can write (62) $\begin{array}{r}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)=S_0\eta (r)+K_2^{\ast }{ϵ}^2+O({ϵ}^3),\end{array}$(62) where $K_2^{\ast }$ is a constant and $S_0$ is given by (63) $\begin{array}{rl}{S}_0& =(M_0-{\beta }_0)\frac{\lambda }{K_{MH}+H_0}[(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GL}(G_0-G_s)L_0]\\ & +\frac{\lambda L_0}{K_{MH}+H_0}[(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0-{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}]\\ & +{\lambda }_{HF}\frac{(M_0-{\beta }_0)L_0}{(K_{HF}+H_0{)}^2}[{\lambda }_{HR}{H}_0{R}_0]+(d_M-d_F){\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\\ & +d_M[(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0]+{[(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0]}^2\\ & +(r_1-2r_2)[(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GL}(G_0-G_s)L_0].\end{array}$(63) Now using (60(60) $\begin{array}{r}{\int }_{1-ϵ}^1{u}_r(s)\mathrm{d}s=0\end{array}$(60) ) and (62(62) $\begin{array}{r}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)=S_0\eta (r)+K_2^{\ast }{ϵ}^2+O({ϵ}^3),\end{array}$(62) ), we can write (64) $\begin{array}{r}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)=S_0[\eta (r)-{\eta }^{\ast }]+O({ϵ}^3),\end{array}$(64) where ${\eta }^{\ast }=\frac{1}{ϵ}{\int }_{1-ϵ}^1\eta (s)\mathrm{d}s=-\frac{1}{6}{ϵ}^2+O({ϵ}^3).$Taking $r=1-ϵ$ in (64(64) $\begin{array}{r}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)=S_0[\eta (r)-{\eta }^{\ast }]+O({ϵ}^3),\end{array}$(64) ), we get (65) $\begin{array}{r}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)(1-ϵ)=-\frac{1}{3}{S}_0{ϵ}^2+O({ϵ}^3),\end{array}$(65)

4. Linear stability

Let $(\tilde{G},\tilde{I},\tilde{\beta },\tilde{R},\tilde{L},\tilde{H},\tilde{F},\tilde{u})$ be the steady solution as in Theorem 1, we are going to study its linear stability. Let (66) $\begin{array}{rl}& G(r,t)=\tilde{G}(r)+\tau G^1(r,t)+O({\tau }^2),\\ & I(r,t)=\tilde{I}(r)+\tau I^1(r,t)+O({\tau }^2),\\ & L(r,t)=\tilde{L}(r)+\tau L^1(r,t)+O({\tau }^2),\\ & H(r,t)=\tilde{H}(r)+\tau H^1(r,t)+O({\tau }^2),\\ & R(r,t)=\tilde{R}(r)+\tau R^1(r,t)+O({\tau }^2),\\ & \beta (r,t)=\tilde{\beta }(r)+\tau {\beta }^1(r,t)+O({\tau }^2),\\ & F(r,t)=\tilde{F}(r)+\tau F^1(r,t)+O({\tau }^2),\\ & B(t)=1-ϵ+\tau B^1(t)+O({\tau }^2),\\ & u(r,t)=\tilde{u}(r)+\tau u^1(r,t)+O({\tau }^2).\end{array}$(66) The linearized equations are (67) $\begin{array}{rl}& \frac{\mathrm{\partial }{G}^1}{\mathrm{\partial }t}-D_G\mathrm{\Delta }{G}^1=S_I\tilde{I}{G}^1-{\lambda }_{GH}[(\tilde{G}-G_s)H^1+G^1\tilde{L}],\end{array}$(67) (68) $\begin{array}{rl}& \frac{\mathrm{\partial }{\beta }^1}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}{\beta }^1)-D_{\beta }\mathrm{\Delta }{\beta }^1=(-d_0+r_1\tilde{G}-r_2{\tilde{G}}^2){\beta }^1+\tilde{\beta }(r_1-2r_2\tilde{G})G^1,\end{array}$(68) (69) $\begin{array}{rl}& \frac{\mathrm{\partial }{I}^1}{\mathrm{\partial }t}-D_I\mathrm{\Delta }{I}^1=\frac{\sigma {\tilde{G}}^2{\beta }^1}{\alpha +{\tilde{G}}^2}-k_I{I}^1-{\lambda }_{RI}(\tilde{R}{I}^1+\tilde{I}{R}^1)\end{array}$(69) (70) $\begin{array}{rl}& +\frac{2\sigma G^1\tilde{\beta }}{\alpha +{\tilde{G}}^2}-\frac{2\sigma {\tilde{G}}^3\tilde{\beta }{G}^1}{(\alpha +{\tilde{G}}^2{)}^2};\\ & \frac{\mathrm{\partial }{L}^1}{\mathrm{\partial }t}-D_L\mathrm{\Delta }{L}^1=-k_L(\tilde{R}{L}^1+\tilde{L}{R}^1)+{\lambda }_{GL}[(\tilde{G}-G_s)L^1+\tilde{L}{G}^1]\end{array}$(70) (71) $\begin{array}{rl}& -{\lambda }_{ML}\frac{(M_0-\tilde{\beta }-\tilde{F})L^1}{K_{ML}+\tilde{L}}+{\lambda }_{ML}\frac{({\beta }^1+F^1)\tilde{L}}{K_{ML}+\tilde{L}}+{\lambda }_{ML}\frac{(M_0-\tilde{\beta }-\tilde{F})\tilde{L}{L}^1}{(K_{ML}+\tilde{L}{)}^2},\end{array}$(71) (72) $\begin{array}{rl}& \frac{\mathrm{\partial }{H}^1}{\mathrm{\partial }t}-D_H\mathrm{\Delta }{H}^1=-k_H(\tilde{R}{H}^1+\tilde{H}{R}^1)-{\lambda }_{HF}\frac{\tilde{H}{F}^1+\tilde{F}{H}^1}{(K_{HF}+\tilde{H})}+{\lambda }_{HF}\frac{\tilde{H}{F}^1\tilde{F}}{(K_{HF}+\tilde{H}{)}^2},\end{array}$(72) (73) $\begin{array}{rl}& \frac{\mathrm{\partial }{R}^1}{\mathrm{\partial }t}-D_R\mathrm{\Delta }{R}^1\\ & =-\tilde{R}(k_L{L}^1+k_H{H}^1)-R^1(k_L\tilde{L}+k_H\tilde{H})\\ & +{\lambda }_{GR}[(\tilde{G}-G_s)R^1+\tilde{R}{G}^1]+{\lambda }_{IR}(\tilde{I}{R}^1+\tilde{R}{I}^1),\end{array}$(73) (74) $\begin{array}{rl}& \frac{\mathrm{\partial }{F}^1}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}{F}^1)-D_F\mathrm{\Delta }{F}^1\\ & ={\lambda }_{ML}\frac{(M_0-\tilde{\beta }-\tilde{F})L^1}{K_{ML}+\tilde{L}}-{\lambda }_{ML}\frac{({\beta }^1+F^1)\tilde{L}}{K_{ML}+\tilde{L}}-{\lambda }_{ML}\frac{(M_0-\tilde{\beta }-\tilde{F})\tilde{L}{L}^1}{(K_{ML}+\tilde{L}{)}^2}\\ & -{\lambda }_{HF}\frac{\tilde{H}{F}^1+\tilde{F}{H}^1}{(K_{HF}+\tilde{H})}+{\lambda }_{HF}\frac{\tilde{H}\tilde{F}{H}^1}{(K_{HF}+\tilde{H}{)}^2}-d_F{F}^1.\end{array}$(74) The equation for $u^1$ is as follows from (10(10) $\begin{array}{r}{M}_0\frac{1}{r}{\mathrm{\partial }}_r(ru)=(-d_0+r_{1}G-r_2{G}^2)\beta +\lambda \frac{(M_0-F-\beta )L}{K_{MH}+H}-d_{F}F-d_M(M_0-F-\beta ).\end{array}$(10) ): (75) $\begin{array}{r}{M}_0\frac{1}{r}{\mathrm{\partial }}_r(r{u}^1)=(-d_0+r_1\tilde{G}-r_2{\tilde{G}}^2){\beta }^1+((r_1-2r_2{G}^1)\tilde{\beta }-d_F{F}^1-d_M(-F^1-{\beta }^1)\\ +\lambda \frac{(M_0-\tilde{F}-\tilde{\beta })L^1}{K_{MH}+\tilde{H}}-\lambda \frac{(M_0-\tilde{F}-\tilde{\beta })\tilde{L}{H}^1}{(K_{MH}+H{)}^2}-\lambda \frac{(F^1+{\beta }^1)\tilde{L}}{K_{MH}+\tilde{H}}.\end{array}$(75) No flux condition on the fixed boundary r = 1: (76) $\begin{array}{r}\frac{\mathrm{\partial }{L}^1}{\mathrm{\partial }n}=\frac{\mathrm{\partial }{H}^1}{\mathrm{\partial }n}=\frac{\mathrm{\partial }{G}^1}{\mathrm{\partial }n}=\frac{\mathrm{\partial }{I}^1}{\mathrm{\partial }n}=\frac{\mathrm{\partial }{R}^1}{\mathrm{\partial }n}=\frac{\mathrm{\partial }{M}^1}{\mathrm{\partial }n}=\frac{\mathrm{\partial }{\beta }^1}{\mathrm{\partial }n}=\frac{\mathrm{\partial }{F}^1}{\mathrm{\partial }n}=u^1(1)=0.\end{array}$(76) Flux boundary conditions on the free boundary $r=1-ϵ$: (77) $\begin{array}{rl}& -\frac{\mathrm{\partial }{L}^1}{\mathrm{\partial }n}+{\alpha }_L{L}^1=\frac{\mathrm{\partial }}{\mathrm{\partial }r}[\frac{\mathrm{\partial }\tilde{L}}{\mathrm{\partial }n}-{\alpha }_L\tilde{L}]B^1,\end{array}$(77) (78) $\begin{array}{rl}& -\frac{\mathrm{\partial }{H}^1}{\mathrm{\partial }n}+{\alpha }_H{H}^1=\frac{\mathrm{\partial }}{\mathrm{\partial }r}[\frac{\mathrm{\partial }\tilde{H}}{\mathrm{\partial }n}-{\alpha }_H\tilde{H}]B^1,\end{array}$(78) (79) $\begin{array}{rl}& -\frac{\mathrm{\partial }{G}^1}{\mathrm{\partial }n}+{\alpha }_G{G}^1=\frac{\mathrm{\partial }}{\mathrm{\partial }r}[\frac{\mathrm{\partial }\tilde{G}}{\mathrm{\partial }n}-{\alpha }_G\tilde{G}]B^1,\end{array}$(79) (80) $\begin{array}{rl}& -\frac{\mathrm{\partial }{I}^1}{\mathrm{\partial }n}+{\alpha }_I{I}^1=\frac{\mathrm{\partial }}{\mathrm{\partial }r}[\frac{\mathrm{\partial }\tilde{I}}{\mathrm{\partial }n}-{\alpha }_I\tilde{I}]B^1,\end{array}$(80) (81) $\begin{array}{rl}& -\frac{\mathrm{\partial }{R}^1}{\mathrm{\partial }n}+{\alpha }_R{R}^1=\frac{\mathrm{\partial }}{\mathrm{\partial }r}[\frac{\mathrm{\partial }\tilde{R}}{\mathrm{\partial }n}-{\alpha }_R\tilde{R}]B^1,\end{array}$(81) (82) $\begin{array}{rl}& -\frac{\mathrm{\partial }{\beta }^1}{\mathrm{\partial }n}+{\alpha }_M{\beta }^1=\frac{\mathrm{\partial }}{\mathrm{\partial }r}[\frac{\mathrm{\partial }\tilde{\beta }}{\mathrm{\partial }n}-{\alpha }_M\tilde{\beta }]B^1,\end{array}$(82) (83) $\begin{array}{rl}& -\frac{\mathrm{\partial }{F}^1}{\mathrm{\partial }n}+{\alpha }_M{F}^1=\frac{\mathrm{\partial }}{\mathrm{\partial }r}[\frac{\mathrm{\partial }\tilde{F}}{\mathrm{\partial }n}-{\alpha }_M\tilde{F}]B^1,\end{array}$(83) The equation for the free boundary leads (84) $\begin{array}{r}\frac{d{B}^1(t)}{dt}=u^1(1-ϵ,t)+{\tilde{u}}_r(1-ϵ).\end{array}$(84) Using (57(57) $\begin{array}{rl}G& =[G_0+({\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GH}(G_0-G_s)H_0)(\eta (r)+\frac{ϵ}{{\alpha }_G})]+G^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:G_0+G_1ϵ+G_2{ϵ}^2+O({ϵ}^3),\\ I& =[I_0+(\frac{\sigma G_0^2{\beta }_0}{\alpha +G_0^2}-k_I{I}_0-{\lambda }_{RI}{R}_0{I}_0)(\eta (r)+\frac{ϵ}{{\alpha }_I})]+I^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:I_0+I_1ϵ+L_2{ϵ}^2+O({ϵ}^3),\\ \beta & =[{\beta }_0+(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0(\eta (r)+\frac{ϵ}{{\alpha }_M})]+{\beta }^{\ast }{ϵ}^2+O({ϵ}^3));\\ R& =[R_0+({\tilde{R}}_0-R_0(k_L{L}_0+k_H{H}_0)+{\lambda }_{GR}(G_0-G_s)R_0+{\lambda }_{IR}{I}_0{R}_0)(\eta (r)++\frac{ϵ}{{\alpha }_R})]\\ & +R^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:R_0+R_1ϵ+R_2{ϵ}^2+O({ϵ}^3)\\ L& =[L_0+(-k_L{R}_0{L}_0+{\lambda }_{GL}(G_0-G_s)L_0-{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0})(\eta (r)+\frac{ϵ}{{\alpha }_L})]\\ & +L^{\ast }{ϵ}^2+O({ϵ}^3),\\ & =:L_0+L_1ϵ+L_2{ϵ}^2+O({ϵ}^3)\\ H& =[H_0+(-k_H{R}_0{H}_0-{\lambda }_{HF}\frac{H_0{F}_0}{(K_{HF}+H_0)})(\eta (r)+\frac{ϵ}{{\alpha }_H})]+H^{\ast }{ϵ}^2+O({ϵ}^3);\\ & =:H_0+H_1ϵ+H_2{ϵ}^2+O({ϵ}^3),\\ F& =\left({\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right)(\eta (r)+\frac{ϵ}{{\alpha }_M})+F^{\ast }{ϵ}^2+O({ϵ}^3)\\ & =:F_1ϵ+F_2{ϵ}^2+O({ϵ}^3),\end{array}$(57) ) and (77(77) $\begin{array}{rl}& -\frac{\mathrm{\partial }{L}^1}{\mathrm{\partial }n}+{\alpha }_L{L}^1=\frac{\mathrm{\partial }}{\mathrm{\partial }r}[\frac{\mathrm{\partial }\tilde{L}}{\mathrm{\partial }n}-{\alpha }_L\tilde{L}]B^1,\end{array}$(77) )–(83(83) $\begin{array}{rl}& -\frac{\mathrm{\partial }{F}^1}{\mathrm{\partial }n}+{\alpha }_M{F}^1=\frac{\mathrm{\partial }}{\mathrm{\partial }r}[\frac{\mathrm{\partial }\tilde{F}}{\mathrm{\partial }n}-{\alpha }_M\tilde{F}]B^1,\end{array}$(83) ), we have (85) $\begin{array}{rl}& -\frac{\mathrm{\partial }{L}^1}{\mathrm{\partial }n}+{\alpha }_L{L}^1=[k_L{R}_0{L}_0-{\lambda }_{GL}(G_0-G_s)L_0+{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}+O(ϵ)]B^1=:{\eta }_L{B}^1,\end{array}$(85) (86) $\begin{array}{rl}& -\frac{\mathrm{\partial }{H}^1}{\mathrm{\partial }n}+{\alpha }_H{H}^1=[k_H{R}_0{H}_0+{\lambda }_{HF}\frac{H_0{F}_0}{(K_{HF}+H_0)}+O(ϵ)]B^1=:{\eta }_H{B}^1,\end{array}$(86) (87) $\begin{array}{rl}& -\frac{\mathrm{\partial }{G}^1}{\mathrm{\partial }n}+{\alpha }_G{G}^1=-[{\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GH}(G_0-G_s)H_0+O(ϵ)]B^1=:{\eta }_G{B}^1,\end{array}$(87) (88) $\begin{array}{rl}& -\frac{\mathrm{\partial }{I}^1}{\mathrm{\partial }n}+{\alpha }_I{I}^1=[\frac{\sigma G_0^2{\beta }_0}{\alpha +G_0^2}-k_I{I}_0-{\lambda }_{RI}{R}_0{I}_0+O(ϵ)]=:{\eta }_I{B}^1{B}^1,\end{array}$(88) (89) $\begin{array}{rl}& -\frac{\mathrm{\partial }{R}^1}{\mathrm{\partial }n}+{\alpha }_R{R}^1\\ & =-[{\tilde{R}}_0-R_0(k_L{L}_0+k_H{H}_0)+{\lambda }_{GR}(G_0-G_s)R_0+{\lambda }_{IR}{I}_0{R}_0+O(ϵ)]B^1=:{\eta }_R{B}^1,\end{array}$(89) (90) $\begin{array}{rl}& -\frac{\mathrm{\partial }{\beta }^1}{\mathrm{\partial }n}+{\alpha }_M{\beta }^1=-[(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0+O(ϵ)]B^1=:{\eta }_{\beta }{B}^1,\end{array}$(90) (91) $\begin{array}{rl}& -\frac{\mathrm{\partial }{F}^1}{\mathrm{\partial }n}+{\alpha }_M{F}^1=-[{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}+O(ϵ)]B^1=:{\eta }_F{B}^1.\end{array}$(91) Now using (75(75) $\begin{array}{r}{M}_0\frac{1}{r}{\mathrm{\partial }}_r(r{u}^1)=(-d_0+r_1\tilde{G}-r_2{\tilde{G}}^2){\beta }^1+((r_1-2r_2{G}^1)\tilde{\beta }-d_F{F}^1-d_M(-F^1-{\beta }^1)\\ +\lambda \frac{(M_0-\tilde{F}-\tilde{\beta })L^1}{K_{MH}+\tilde{H}}-\lambda \frac{(M_0-\tilde{F}-\tilde{\beta })\tilde{L}{H}^1}{(K_{MH}+H{)}^2}-\lambda \frac{(F^1+{\beta }^1)\tilde{L}}{K_{MH}+\tilde{H}}.\end{array}$(75) ) and (84(84) $\begin{array}{r}\frac{d{B}^1(t)}{dt}=u^1(1-ϵ,t)+{\tilde{u}}_r(1-ϵ).\end{array}$(84) )–(91(91) $\begin{array}{rl}& -\frac{\mathrm{\partial }{F}^1}{\mathrm{\partial }n}+{\alpha }_M{F}^1=-[{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}+O(ϵ)]B^1=:{\eta }_F{B}^1.\end{array}$(91) ), we can write (92) $\begin{array}{r}{M}_0{\mathrm{\partial }}_r{u}^1=J_1+J_2{B}^1,\end{array}$(92) where (93) $\begin{array}{rl}{J}_1& =(M_0-\tilde{\beta }-\tilde{F})\frac{\lambda }{K_{MH}+\tilde{H}}(L^1-\frac{{\eta }_L}{{\alpha }_L}{B}^1)\\ & -(M_0-\tilde{F}-\tilde{\beta })\frac{\lambda \tilde{L}}{K_{MH}+\tilde{H}}[{\beta }^1-(\frac{{\eta }_{\beta }}{{\alpha }_M}{B}^1+F^1-\frac{{\eta }_F}{{\alpha }_M}{B}^1)]\\ & -(M_0-\tilde{\beta }-\tilde{F})\frac{\lambda }{(K_{MH}+\tilde{H}{)}^2}(H^1-\frac{{\eta }_H}{{\alpha }_H}{B}^1)+d_M({\beta }^1-\frac{{\eta }_{\beta }}{{\alpha }_M}{B}^1+F^1-\frac{{\eta }_F}{{\alpha }_M}{B}^1)\\ & -d_F(F^1-\frac{{\eta }_F}{{\alpha }_M}{B}^1)+{\beta }_0(r_1(G^1-\frac{{\eta }_G}{{\alpha }_G}{B}^1)-2r_2\tilde{G}(G^1-\frac{{\eta }_G}{{\alpha }_G}))\\ & +(-d_0+r_1\tilde{G}-r_2{\tilde{G}}^2)({\beta }^1-\frac{{\eta }_{\beta }}{{\alpha }_M}{B}^1);\end{array}$(93) (94) $\begin{array}{rl}{J}_2& =(M_0-{\beta }_0)\frac{\lambda }{K_{MH}+H_0}\frac{{\eta }_L}{{\alpha }_L}-(M_0-{\beta }_0)\frac{\lambda L_0}{K_{MH}+H_0}(\frac{{\eta }_{\beta }}{{\alpha }_M}+\frac{{\eta }_F}{{\alpha }_M})\\ & -(M_0-{\beta }_0)\frac{\lambda }{(K_{MH}+H_0{)}^2}\frac{{\eta }_H}{{\alpha }_H}+d_M(\frac{{\eta }_{\beta }}{{\alpha }_M}+\frac{{\eta }_F}{{\alpha }_M})\\ & -d_F\frac{{\eta }_F}{{\alpha }_M}+{\beta }_0(r_1\frac{{\eta }_G}{{\alpha }_G}-2r_2{G}_0\frac{{\eta }_G}{{\alpha }_G})+(-d_0+r_1{G}_0-r_2{G}_0^2)\frac{{\eta }_{\beta }}{{\alpha }_M}.\end{array}$(94) Integrating (92(92) $\begin{array}{r}{M}_0{\mathrm{\partial }}_r{u}^1=J_1+J_2{B}^1,\end{array}$(92) ) we obtain: (95) $\begin{array}{r}{u}^1(1-ϵ,t)=-\frac{1}{M_0}{\int }_{1-ϵ}^1{J}_1\mathrm{d}r-\frac{B^1}{M_0}{\int }_{1-ϵ}^1{J}_2\mathrm{d}r\end{array}$(95) Set (96) $\begin{array}{rl}{S}_1& =-\frac{(M_0-{\beta }_0)}{M_0}\frac{\lambda }{K_{MH}+H_0}\frac{[k_L{R}_0{L}_0-{\lambda }_{GL}(G_0-G_s)L_0+{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}]}{{\alpha }_L}\\ & -\frac{(M_0-{\beta }_0)}{M_0}\frac{\lambda L_0}{K_{MH}+H_0}(\frac{[(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0]}{{\alpha }_M}+\frac{\left[{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right]}{{\alpha }_M})\\ & +\frac{(M_0-{\beta }_0)}{M_0}\frac{\lambda }{(K_{MH}+H_0{)}^2}\frac{[k_H{R}_0{H}_0+{\lambda }_{HF}\frac{H_0{F}_0}{(K_{HF}+H_0)}]}{{\alpha }_H}\\ & +\frac{d_M}{M_0}(\frac{[(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0]}{{\alpha }_M}+\frac{\left[{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right]}{{\alpha }_M})\\ & -\frac{d_F}{M_0}\frac{\left[{\lambda }_{ML}\frac{(M_0-{\beta }_0)L_0}{K_{ML}+L_0}\right]}{{\alpha }_M}+\frac{{\beta }_0}{M_0}\left(r_1\frac{[{\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GH}(G_0-G_s)H_0]}{{\alpha }_G}\right)\\ & -\frac{{\beta }_0}{M_0}\left(\frac{2r_2{G}_0}{M_0}\frac{[{\tilde{G}}_0-(E_{G0}+S_I{I}_0)G_0-{\lambda }_{GL}(G_0-G_s)L_0]}{{\alpha }_G}\right)\\ & +\frac{(-d_0+r_1{G}_0-r_2{G}_0^2)}{M_0}\frac{[(-d_0+r_1{G}_0-r_2{G}_0^2){\beta }_0]}{{\alpha }_M}.\end{array}$(96) Then from (84(84) $\begin{array}{r}\frac{d{B}^1(t)}{dt}=u^1(1-ϵ,t)+{\tilde{u}}_r(1-ϵ).\end{array}$(84) ), (95(95) $\begin{array}{r}{u}^1(1-ϵ,t)=-\frac{1}{M_0}{\int }_{1-ϵ}^1{J}_1\mathrm{d}r-\frac{B^1}{M_0}{\int }_{1-ϵ}^1{J}_2\mathrm{d}r\end{array}$(95) ) and (65(65) $\begin{array}{r}\frac{M_0}{r}{\mathrm{\partial }}_r(ru)(1-ϵ)=-\frac{1}{3}{S}_0{ϵ}^2+O({ϵ}^3),\end{array}$(65) ), we have (97) $\begin{array}{r}\frac{d{B}^1(t)}{dt}=(S_1ϵ+O({ϵ}^2))B^1+a_1.\end{array}$(97) where (98) $\begin{array}{r}{a}_1=-\frac{1}{M_0}{\int }_{1-ϵ}^1{J}_1\mathrm{d}r.\end{array}$(98) Using (93(93) $\begin{array}{rl}{J}_1& =(M_0-\tilde{\beta }-\tilde{F})\frac{\lambda }{K_{MH}+\tilde{H}}(L^1-\frac{{\eta }_L}{{\alpha }_L}{B}^1)\\ & -(M_0-\tilde{F}-\tilde{\beta })\frac{\lambda \tilde{L}}{K_{MH}+\tilde{H}}[{\beta }^1-(\frac{{\eta }_{\beta }}{{\alpha }_M}{B}^1+F^1-\frac{{\eta }_F}{{\alpha }_M}{B}^1)]\\ & -(M_0-\tilde{\beta }-\tilde{F})\frac{\lambda }{(K_{MH}+\tilde{H}{)}^2}(H^1-\frac{{\eta }_H}{{\alpha }_H}{B}^1)+d_M({\beta }^1-\frac{{\eta }_{\beta }}{{\alpha }_M}{B}^1+F^1-\frac{{\eta }_F}{{\alpha }_M}{B}^1)\\ & -d_F(F^1-\frac{{\eta }_F}{{\alpha }_M}{B}^1)+{\beta }_0(r_1(G^1-\frac{{\eta }_G}{{\alpha }_G}{B}^1)-2r_2\tilde{G}(G^1-\frac{{\eta }_G}{{\alpha }_G}))\\ & +(-d_0+r_1\tilde{G}-r_2{\tilde{G}}^2)({\beta }^1-\frac{{\eta }_{\beta }}{{\alpha }_M}{B}^1);\end{array}$(93) ) and (67(67) $\begin{array}{rl}& \frac{\mathrm{\partial }{G}^1}{\mathrm{\partial }t}-D_G\mathrm{\Delta }{G}^1=S_I\tilde{I}{G}^1-{\lambda }_{GH}[(\tilde{G}-G_s)H^1+G^1\tilde{L}],\end{array}$(67) )–(74(74) $\begin{array}{rl}& \frac{\mathrm{\partial }{F}^1}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\mathbf{u}{F}^1)-D_F\mathrm{\Delta }{F}^1\\ & ={\lambda }_{ML}\frac{(M_0-\tilde{\beta }-\tilde{F})L^1}{K_{ML}+\tilde{L}}-{\lambda }_{ML}\frac{({\beta }^1+F^1)\tilde{L}}{K_{ML}+\tilde{L}}-{\lambda }_{ML}\frac{(M_0-\tilde{\beta }-\tilde{F})\tilde{L}{L}^1}{(K_{ML}+\tilde{L}{)}^2}\\ & -{\lambda }_{HF}\frac{\tilde{H}{F}^1+\tilde{F}{H}^1}{(K_{HF}+\tilde{H})}+{\lambda }_{HF}\frac{\tilde{H}\tilde{F}{H}^1}{(K_{HF}+\tilde{H}{)}^2}-d_F{F}^1.\end{array}$(74) ), we can show that (99) $\begin{array}{r}|a_1|\le C{ϵ}^{3/2}{({\int }_{1-ϵ}^1[|G_{rr}^1{|}^2+|I_{rr}^1{|}^2+|R_{rr}^1{|}^2+|L_{rr}^1{|}^2+|H_{rr}^1{|}^2+|{\beta }_{rr}^1{|}^2+|F_{rr}^1{|}^2])}^{1/2}.\end{array}$(99) We conclude from (97(97) $\begin{array}{r}\frac{d{B}^1(t)}{dt}=(S_1ϵ+O({ϵ}^2))B^1+a_1.\end{array}$(97) ) (98(98) $\begin{array}{r}{a}_1=-\frac{1}{M_0}{\int }_{1-ϵ}^1{J}_1\mathrm{d}r.\end{array}$(98) ) and (99(99) $\begin{array}{r}|a_1|\le C{ϵ}^{3/2}{({\int }_{1-ϵ}^1[|G_{rr}^1{|}^2+|I_{rr}^1{|}^2+|R_{rr}^1{|}^2+|L_{rr}^1{|}^2+|H_{rr}^1{|}^2+|{\beta }_{rr}^1{|}^2+|F_{rr}^1{|}^2])}^{1/2}.\end{array}$(99) ) that

Theorem 4.1

For fixed ϵ, the steady solution in Theorem 1 is linearly asymptotically stable if $S_1<0$, and linearly unstable if $S_1>0$.

5. Plaque persistence due to hyperglycemia

For patients without diabetes, i.e. $G_0=G_s$, it is shown in [12] that an ϵ initial plague shrinks to zero if (100) $\begin{array}{r}\frac{\lambda L_0}{K_{MH}+H_0}-d_M<0\end{array}$(100) In this section, we are going to show that for patients with hyperglycemia, where the following inequality (101(101) $\begin{array}{r}(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)+{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)>0.\end{array}$(101) ) holds, an initial plaque will persist even if (100(100) $\begin{array}{r}\frac{\lambda L_0}{K_{MH}+H_0}-d_M<0\end{array}$(100) ) holds. (101) $\begin{array}{r}(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)+{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)>0.\end{array}$(101)

Theorem 5.1

If (101(101) $\begin{array}{r}(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)+{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)>0.\end{array}$(101) ) holds, assume that (102) $\begin{array}{rl}& 0<L(r,0)\le L_0+Aϵ,H(r,0)<H_0+Aϵ,F(r,0)\le Aϵ,\\ & G_s<G(r,0)<G_0+Aϵ,\beta <{\beta }_0+Aϵ,B(0)=1-ϵ/2\end{array}$(102) for some constant $A\ge 0$, then (103) $\begin{array}{r}\underset{t\to \mathrm{\infty }}{lim sup}B(t)<1.\end{array}$(103)

Proof.

Using strong maximum principle, we can show that (104) $\begin{array}{rl}& L_0-k^{\ast }ϵ\le L(r,t)\le L_0+A^{\ast }ϵ+max|L(r,0)|e^{-c_{1}t},\end{array}$(104) (105) $\begin{array}{rl}& H_0-k_1ϵ\le H(r,t)\le H_0+A_1ϵ+max|H(r,0)|e^{-c_{1}t},\end{array}$(105) (106) $\begin{array}{rl}& 0<F(r,t)<A^{\ast }ϵ,\end{array}$(106) (107) $\begin{array}{rl}& G_0-k_2ϵ\le G(r,t)\le G_0+A_2ϵ+max|H(r,0)|e^{-c_{2}t},\end{array}$(107) (108) $\begin{array}{rl}& {\beta }_0-k_3ϵ\le \beta (r,t)<{\beta }_0+A_3ϵ;\end{array}$(108) where $k^{\ast },k_1,k_2,k_3$ and $A_1,A_2,A_3,A^{\ast }$ are positive constants.

Using (10(10) $\begin{array}{r}{M}_0\frac{1}{r}{\mathrm{\partial }}_r(ru)=(-d_0+r_{1}G-r_2{G}^2)\beta +\lambda \frac{(M_0-F-\beta )L}{K_{MH}+H}-d_{F}F-d_M(M_0-F-\beta ).\end{array}$(10) ) and the above inequalities, we have for sufficiently large t (109) $\begin{array}{r}{M}_0{\mathrm{\partial }}_{r}u(r,t)=(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)+{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)+O(ϵ).\end{array}$(109) Integrating above to obtain (110) $\begin{array}{rl}u(r,t)& =-\frac{1}{M_0}{\int }_r^1[(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)\\ & +{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)+O(ϵ)\phantom{(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)}]<0.\end{array}$(110) Now (14(14) $\begin{array}{r}\frac{dB(t)}{dt}=u(B(t)).\end{array}$(14) ) and (110(110) $\begin{array}{rl}u(r,t)& =-\frac{1}{M_0}{\int }_r^1[(M_0-{\beta }_0)(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)\\ & +{\beta }_0(-d_0+r_1{G}_0^2-r_2{G}_0^2)+O(ϵ)\phantom{(\frac{\lambda L_0}{K_{MH}+H_0}-d_M)}]<0.\end{array}$(110) ) imply that $B(t)$ is decreasing, so (103(103) $\begin{array}{r}\underset{t\to \mathrm{\infty }}{lim sup}B(t)<1.\end{array}$(103) ) holds.

6. Conclusion and discussion

Atherosclerosis is a leading cause of death worldwide; it emerges as a result of multiple dynamical cell processes including hemodynamics, endothelial damage, innate immunity and sterol biochemistry. As atherosclerosis is a cardiovascular condition that affects critical circulatory systems, studying human atheroma poses logistical and ethical problems, as access to live atherosclerotic tissue is limited and disturbances risk triggering plaque rupture. As a result, the appropriate framework to consider emergent dynamical behaviour of this type of disease is mathematical and computational modelling. The existing mathematical models [10–12,17] that describe the growth of a plaque in the artery recognize the critical role of low density lipoprotein (LDL) and high-density lipoprotein (HDL), in determining whether a plaque, once formed, will grow or shrink. Making matters worse, nearly 463 million people worldwide have diabetes. In part because diabetes increases atherosclerosis-related inflammation, diabetic patients are twice as likely to have a heart attack or stroke. Past work has shown that hyperglycemia and insulin resistance alter function of multiple cell types, including endothelium, smooth muscle cells and platelets, indicating the extent of vascular disarray in this disease. Although the pathophysiology of diabetic vascular disease is generally understood, there is no mathematical model to date that includes the effect of diabetes on plaque growth until recent mathematical models for diabetic atherosclerosis [40,41]. In [40,41], governing equations were established to model the mechanism of how diabetes impairs endothelium-dependent (nitric oxide-mediated) vasodilation before the formation of atheroma and mechanism of how arteries affected by diabetes have altered vascular smooth muscle cell function. The local existence and uniqueness of solution to the diabetic atherosclerosis models were also proved. In the present work, we consider a simplified mathematical model for diabetic atherosclerosis involving LDL, HDL, glucose, insulin, free radicals (ROS),β cells, macrophages and foam cells, which satisfy a system of partial differential equations with a free boundary, the interface between the blood flow and the plaque. We establish local existence of small radially symmetric stationary solutions to the model and study their stability. Our analysis shows that the plague will persist due to hyperglycemia even when LDL and HDL are in normal range, hence confirms that diabetes increase the risk of atherosclerosis.

There are still questions to be answered regarding diabetic atherosclerosis. One of the questions is to study the risk to plague growth associated with hyperglycemia by performing numerical simulation of mathematical models. We resorted to estimating parameters for the model, based on existing literature or inference from other cell processes. Animal studies of diabetic atherosclerosis do exist for mouse, rabbit and pig [14,21]. As atherosclerosis is a cardiovascular condition that affects critical circulatory systems, studying human atheroma poses logistical and ethical problems, as access to live atherosclerotic tissue is limited and disturbances risk triggering plaque rupture. Consequently, data are limited. As a result, establishing biologically relevant kinetic parameters that can be used to simulate pathway dynamics is challenging, and comprehensive parameterizations and validation of the mathematical model are difficult.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Research reported in this paper was supported by National Institute of General Medical Sciences of the National Institutes of Health under Award Number UL1GM118973.