Error analysis of finite difference/collocation method for the nonlinear coupled parabolic free boundary problem modeling plaque growth in the artery

TL;DR

This paper introduces a new finite difference and collocation method for solving a complex nonlinear free boundary problem modeling arterial plaque growth, demonstrating stability, convergence, and efficiency through numerical results.

Contribution

It presents a novel combination of front fixing, finite difference, and collocation methods for a nonlinear free boundary problem in atherosclerosis modeling, with proven stability and convergence.

Findings

Method is stable and convergent.

Numerical results confirm efficiency.

Effective in modeling plaque growth.

Abstract

The main target of this paper is to present a new and efficient method to solve a nonlinear free boundary mathematical model of atherosclerosis. This model consists of three parabolics, one elliptic and one ordinary differential equations that are coupled together and describe the growth of a plaque in the artery. We start our discussion by using the front fixing method to fix the free domain and simplify the model by changing the mix boundary condition to a Neumann one by applying suitable changes of variables. Then, after employing a nonclassical finite difference and the collocation method on this model, we prove the stability and convergence of methods. Finally, some numerical results are considered to show the efficiency of the method.

Tables

Table 1. Table 1: Maximum time-error with N=10 and various M by considering TFBM.

Error of		M=100		M=200		M=300		M=400		M=500		M=600
L		$1.832329e-05$		$7.33126e-06$		$4.11065e-06$		$2.58829e-06$		$1.70362e-06$		$1.12596e-06$
H		$2.57909e-07$		$1.03322e-07$		$5.79623e-08$		$3.65061e-08$		$2.40324e-08$		$1.58853e-08$
F		$6.60404e-09$		$4.21233e-09$		$3.16373e-09$		$2.26606e-09$		$1.60468e-09$		$1.11186e-09$

Table 2. Table 3: Maximum space-error with M=100 and various N by considering TFBM.

Error		N=2		N=4		N=6		N=8
L		0.399186e-05		0.286885e-05		0.190745e-05		0.119180e-05
H		0.504070e-06		0.359887-06		0.239578e-06		0.149847e-06
F		0.293787e-07		0.217007e-07		0.144762e-07		0.905871e-08

Table 3. Table 4: Condition number of the coefficient matrices (CN).

CN		N=4		N=6		N=8		N=10		N=12		N=14		N=16		N=18		N=20
Eq. (1)		272e+04		115e+06		529e+07		258e+09		131e+11		689e+12		369e+12		202e+15		singular
Eq. (2)		272e+04		115e+06		529e+07		258e+09		131e+11		689e+12		369e+12		202e+15		singular
Eq. (3)		367e+05		303e+07		229e+09		166e+11		117e+13		814e+14		558e+15		380e+17		singular

Table 4. Table 5: Maximum time-error with N=50 and various M by considering TFBL.

Error of		M=100		M=200		M=300		M=400		M=500		M=600
L		$2.21985e-06$		$9.74874e-07$		$5.66407e-07$		$3.63391e-07$		$2.41969e-07$		$1.61184e-07$
H		$8.08703e-06$		$3.67228e-06$		$2.13760e-06$		$1.37270e-06$		$9.14550e-07$		$6.09439e-07$
F		$3.89421e-06$		$1.68696e-06$		$9.75669e-07$		$6.24533e-07$		$4.15285e-07$		$2.76382e-07$

Table 5. Table 7: Maximum space-error with M=200 and various N by considering TFBL.

Error of		N=10		N=20		N=40		N=80		N=100
L		$1.53065e-10$		$7.04107e-12$		$6.98306e-13$		$1.21789e-13$		$7.34745e-14$
H		$2.96583e-11$		$8.76799e-13$		$2.18054e-13$		$8.28391e-14$		$5.54617e-14$
F		$1.44016e-11$		$6.17966e-12$		$3.09801e-12$		$1.69595e-12$		$1.34437e-12$

Table 6. Table 9: Variation of the radius of the plaque toward the various level of L 0 subscript 𝐿 0 L_{0} and H 0 subscript 𝐻 0 H_{0} in blood during the days.

$(L_0,H_0)$		$0$ day		$10$ day		$20$ day		$30$ day		$40$ day		$50$ day		$60$ day
$(150,45)$		$0.900000$		$0.897461$		$0.894874$		$0.892264$		$0.889634$		$0.886984$		$0.884316$
$(120,60)$		$0.900000$		$0.907537$		$0.914481$		$0.920822$		$0.926622$		$0.931935$		$0.936810$
$(46.5,139.5)$		$0.900000$		$0.900008$		$0.900018$		$0.900027$		$0.900037$		$0.900046$		$0.900055$

Equations

$$\frac{\partial \hat{L}}{\partial t}-\Delta \hat{L}=-k_1\frac{(M0-\hat{F})\hat{L}}{K_1+\hat{L}}-r_1\hat{L},R(t)<r<1,t>0,$$

$$\frac{\partial \hat{L}}{\partial n}+\alpha (\hat{L}-L_0)=0,\text{at}r=R(t),t>0,\frac{\partial \hat{L}}{\partial n}=0,\text{at}r=1,t>0,\hat{L}(r,0)=L_0,$$

$$\frac{\partial \hat{H}}{\partial t}-\Delta \hat{H}=-k_2\frac{\hat{H}\hat{F}}{K_2+\hat{F}}-r_2\hat{H},R(t)<r<1,t>0,$$

$$\frac{\partial \hat{H}}{\partial n}+\alpha (\hat{H}-H_0)=0,\text{at}r=R(t),t>0,\frac{\partial \hat{H}}{\partial n}=0,\text{at}r=1,t>0,\hat{H}(r,0)=H_0,$$

$$\frac{\partial \hat{F}}{\partial t}-D\Delta \hat{F}+{\hat{F}}_r.v=k_1\frac{(M_0-\hat{F})\hat{L}}{K1+\hat{L}}-k_2\frac{\hat{H}\hat{F}}{K_2+\hat{F}}-$$

$$\lambda \frac{\hat{F}(M_0-\hat{F})\hat{L}}{M_0(\delta +\hat{H})}+\frac{{\mu }_1}{M_0}(M_0-\hat{F})\hat{F}-\frac{{\mu }_2}{M_0}\hat{F}(M_0-\hat{F}),R(t)<r<1,t>0,$$

$$\frac{\partial \hat{F}}{\partial n}+\beta \hat{F}=0,\text{at}r=R(t),t>0,\frac{\partial \hat{F}}{\partial n}=0,\text{at}r=1,t>0,F(r,0)=0,$$

$$M_0.v_r=\lambda \frac{(M_0-\hat{F})\hat{L}}{\delta +\hat{H}}-{\mu }_1(M_0-\hat{F})-{\mu }_2\hat{F},R(t)<r<1,t>0,$$

$$v(r,t)=0,\text{at}r=1,t>0,$$

$$\frac{dR(t)}{dt}=v(R(t),t),t>0,$$

$$R(0)=ϵ.$$

$$\frac{\partial \hat{L}}{\partial t}-\frac{1}{r^2}\frac{\partial }{\partial r}(r^2\frac{\partial \hat{L}}{\partial r})=-k_1\frac{(M_0-\hat{F})\hat{L}}{K_1+\hat{L}}-r_1\hat{L},R(t)<r<1,t>0,$$

$$\frac{\partial \hat{L}}{\partial r}+\alpha (\hat{L}-L_0)=0\text{at}r=R(t),t>0,\frac{\partial \hat{L}}{\partial r}=0\text{at}r=1,t>0,\hat{L}(r,0)=L_0,$$

$$\frac{\partial \hat{H}}{\partial t}-\frac{1}{r^2}\frac{\partial }{\partial r}(r^2\frac{\partial \hat{H}}{\partial r})=-k_2\frac{\hat{H}\hat{F}}{K_2+\hat{F}}-r_2\hat{H},R(t)<r<1,t>0,$$

$$\frac{\partial \hat{H}}{\partial r}+\alpha (\hat{H}-H_0)=0\text{at}r=R(t),t>0,\frac{\partial \hat{H}}{\partial r}=0\text{at}r=1,t>0,\hat{H}(r,0)=H_0,$$

$$\frac{\partial \hat{F}}{\partial t}-\frac{D}{r^2}\frac{\partial }{\partial r}(r^2\frac{\partial \hat{F}}{\partial r})+{\hat{F}}_r.v=k_1\frac{(M_0-\hat{F})\hat{L}}{K1+\hat{L}}-k_2\frac{\hat{H}\hat{F}}{K_2+\hat{F}}-$$

$$\lambda \frac{\hat{F}(M_0-\hat{F})\hat{L}}{M_0(\delta +\hat{H})}+\frac{{\mu }_1}{M_0}(M_0-\hat{F})\hat{F}-\frac{{\mu }_2}{M_0}\hat{F}(M_0-\hat{F}),R(t)<r<1,t>0,$$

$$\frac{\partial \hat{F}}{\partial r}+\beta \hat{F}=0\text{at}R(t)<r<1,t>0,\frac{\partial \hat{F}}{\partial r}=0\text{at}r=1,t>0,\hat{F}(r,0)=0,$$

$$M_0.v_r=\lambda \frac{(M_0-\hat{F})\hat{L}}{\delta +\hat{H}}-{\mu }_1(M_0-\hat{F})-{\mu }_2\hat{F},R(t)<r<1,t>0,$$

$$v(r,t)=0\text{at}r=1,t>0,$$

$$\frac{dR(t)}{dt}=v(R(t),t),t>0,$$

$$R(0)=ϵ.$$

$$\rho =\frac{r-R(t)}{1-R(t)},\hat{L}(\rho ,t)=\hat{L}(r,t),\hat{H}(\rho ,t)=\hat{H}(r,t),\hat{F}(\rho ,t)=\hat{F}(r,t),$$

$$\{(\rho ,t)\mid 0<\rho <1,t\ge 0\},$$

$$\frac{\partial \hat{L}}{\partial t}-\frac{1}{(1-R(t){)}^2}\frac{{\partial }^2\hat{L}}{\partial {\rho }^2}+\left(\frac{-2}{\rho (1-R(t){)}^2+R(t)(1-R(t))}+\frac{v(0,t)(\rho -1)}{(1-R(t){)}^2}\right)\frac{\partial \hat{L}}{\partial \rho }=f^{\hat{L}}(\hat{L},\hat{F}),$$

$$\frac{\partial \hat{L}}{\partial \rho }+\alpha (1-R(t))(\hat{L}-L_0)=0\text{at}\rho =0,t>0,\frac{\partial \hat{L}}{\partial \rho }=0\text{at}\rho =1,t>0,\hat{L}(\rho ,0)=L_0,$$

$$\frac{\partial \hat{H}}{\partial t}-\frac{1}{(1-R(t){)}^2}\frac{{\partial }^2\hat{H}}{\partial {\rho }^2}+\left(\frac{-2}{\rho (1-R(t){)}^2+R(t)(1-R(t))}+\frac{v(0,t)(\rho -1)}{(1-R(t){)}^2}\right)\frac{\partial \hat{H}}{\partial \rho }=f^{\hat{H}}(\hat{H},\hat{F}),$$

$$\frac{\partial \hat{H}}{\partial \rho }+\alpha (1-R(t))(\hat{H}-H_0)=0\text{at}\rho =0,t>0,\frac{\partial \hat{H}}{\partial \rho }=0\text{at}\rho =1,t>0,\hat{H}(\rho ,0)=H_0,$$

$$\frac{\partial \hat{F}}{\partial t}-\frac{D}{(1-R(t){)}^2}\frac{{\partial }^2\hat{F}}{\partial {\rho }^2}+$$

$$\left(\frac{-2D}{\rho (1-R(t){)}^2+R(t)(1-R(t))}+\frac{v}{1-R(t)}+\frac{v(0,t)(\rho -1)}{(1-R(t){)}^2}\right)\frac{\partial \hat{F}}{\partial \rho }=f^{\hat{F}}(\hat{L},\hat{H},\hat{F}),$$

$$\frac{\partial \hat{F}}{\partial \rho }+\beta (1-R(t))\hat{F}=0\text{at}\rho =0,t>0,\frac{\partial \hat{F}}{\partial \rho }=0\text{at}\rho =1,t>0,\hat{F}(\rho ,0)=0,$$

$$\frac{1}{1-R(t)}\frac{\partial v}{\partial \rho }=f^v(\hat{L},\hat{H},\hat{F}),$$

$$v(\rho ,t)=0\text{at}\rho =1,t>0,$$

$$\frac{dR}{dt}=v(0,t),t>0,$$

$$R(0)=ϵ.$$

$$L(\rho ,t):=\exp (-\alpha (1-R(t))\frac{(1-\rho {)}^2}{2})(\hat{L}(\rho ,t)-L_0),$$

$$H(\rho ,t):=\exp (-\alpha (1-R(t))\frac{(1-\rho {)}^2}{2})(\hat{H}(\rho ,t)-H_0),$$

$$F(\rho ,t):=\exp (-\beta (1-R(t))\frac{(1-\rho {)}^2}{2})\hat{F}(\rho ,t).$$

$$\frac{\partial L}{\partial t}-\frac{1}{(1-R(t){)}^2}\frac{{\partial }^{2}L}{\partial {\rho }^2}+$$

$$\left(\frac{-2}{\rho (1-R(t){)}^2+R(t)(1-R(t))}+\frac{v(0,t)(\rho -1)}{1-R(t)}+\frac{2(1-\rho )\alpha }{1-R(t)}\right)\frac{\partial L}{\partial \rho }=f^L(L,F),$$

$$\frac{\partial L}{\partial \rho }=0\text{at}\rho =0,t>0,\frac{\partial L}{\partial \rho }=0\text{at}\rho =1,t>0,L(\rho ,0)=0,$$

Full text

1 Introduction

There are many phenomena that we have questions about and want to describe them or their behaviours. To find an answer we collect data and multiple sources of data support a rapid knowledge. However, our ability to analyze and interpret this data lags far behind data generation and storage capacity. The question that arises here is how can mathematics help us to solve this problem? The answer to this question is that mathematics helps us in this issue by modeling as described in Figure 1meerschaert2013mathematical . The first recognizable models were numbers since about 30.000 BC schichl2004models . The development of mathematical models continued in various fields. By the invention of calculus by Newton and Leibniz in 1671, differential equations came into existence newton1774methodus , and then, partial differential equations stand out in 1719 by Nicolaus Bernoulli cajori1928early . Early 16th century was the beginning of the modern interaction between mathematics and biology duo to the work of William Harvey harvey1737exercitatio . Thereafter, until the 20th century, biological modeling continued to advance by the work of Hardy in 1908 in population genetics hardy1908mendelian , Yula in 1925 in birth and death process yule1925growth , Luria and Delbruck in 1943 in estimating bacterial mutation rates luria1943mutations . The question that matters here is what biological issues to biological scientists are more valuable and selected for modeling. In reply to the question, according to the authors’ consideration, topics of the greatest importance to them are those that affect the human community more, which is the factor that causes a tendency to model for a more detailed examination of these issues. Due to human society conditions, studying diseases is the most important area of study. According to the ICD 10111International Classification of Diseasesworld2004international , the most important diseases that cause of death in the world are HIVmarinho2012model ; culshaw2003mathematical , Tumor rockne2009mathematical ; owen1999mathematical ; holmes2000mathematical ; chaplain1995mathematical , Cancer stephanou2006mathematical ; deng2004mathematical ; anderson2000mathematical ; franks2003modelling , Cardiovascular diseases (especially Atherosclerosis) yang2016mathematical ; eberhard2006numerical ; islam2016mathematical ; calvez2009mathematical ; calvez2010mathematical ; hao2014ldl ; 2015mathematical and Wound healing flegg2012wound ; terry2012mathematical ; paul2013mathematical ; javierre2009mathematical respectively. As described above, the heart attack or stroke that happens because of atherosclerosis diseases is one of the third leading cause of death in the world world2004international . There are two various perspectives to research in this area, ”modeling” and ”numerical analysis” point of view. In the case of modeling, the important point to note here is that there are different perspectives in studying Atherosclerosis seo2007mice ; ougrinovskaia2010ode . A very important perspective in the area of studying this disease is to study mathematical modeling in the forms of ordinary differential equation (ODE) and partial differential equation (PDE), that are distinguished from each other by the various factors such as choosing the region and the boundary, elements involved in biologically issues, boundary conditions and etc. For instance, in 2009 a 3-D nonlinear parabolic system of PDEs is considered as a mathematical model of Atherosclerosis involving the local blood flow dynamics by Calvezcalvez2009mathematical . In 2010 Calvez extended his previous work calvez2010mathematical . In both articles, Calvez assumed the artery to be an irregular 3-D cylinder. Friedman in 2014, 2015 presented mathematical models of plaque growth in the artery, which parabolic nonlinear 2-D system of PDEs with mixed boundary conditions and free boundary are considered for modeling friedman2015free ; hao2014ldl . It is worth noting that the difference between these models is that in friedman2015free the artery is assumed to be a very long circular cylinder but in hao2014ldl the artery is considered as an irregular cylinder. In other work, a mathematical model with the approach used in hao2014ldl by including the effect of reverse cholesterol transport of plaque growth which includes the (LDL, HDL) concentrations is developed by Friedman 2015mathematical . There are many other models that are associated with this disease eberhard2006numerical ; little2009model ; yang2016mathematical . In the case of numerical analysis, there are various mathematical methods and different perspectives in convergence and stability analyzing. For instance, in esmaili2018application mathematical modeling of a tumor is considered and numerical results with convergence and stability of numerical methods have been presented. In esmaili2017optimal mathematical modeling of optimal control of tumor with drug application has been presented and in esmaili2018application it has been numerically solved and analyzed. Additionally, in zhao2017convergence , a two-dimensional multi-term time fractional diffusion equation is solved numerically using a fully-discrete schem and convergence and superconvergence of the method is illustrated. For the mathematical analysis point of view, in ramezani2008combined a hyperbolic equation with an integral condition is solved using finite difference/spectral method.

In this article, we want to solve a free boundary nonlinear system of coupled PDEs that model Atherosclerosis and consists of three parabolics, one elliptic and one ordinary differential equation which is introduced in friedman2015free . For the readers’ convenience, we highlight the main goals of this study as follows

• $\bullet$

We have fixed the domain using the front fixing method and simplified the model by changing the mix boundary condition to a Neumann one by applying a suitable change of variables to achieve more comfortable results for numerical analysis.

• $\bullet$

Applying the finite difference method, we have constructed a sequence, which converges to the exact solution of coupled partial differential equations.

• $\bullet$

In each time step, using Taylor theorem, the problem has changed to linear one (see (19)-(23)) and using the collocation method, equations (19)-(23) are solved numerically.

• $\bullet$

We have proved constructed sequence converges to the exact solution of the problem (see Theorem 5.2) and also the stability of the method is proven (see Theorem 5.3).

• $\bullet$

We have simulated the model using finite difference and collocation method for some pair of values $(L_{0},H_{0})$ to show the validity and efficiency of the presented method. It is critical to note that, construction of a new second-order non-classical discretization formula helps us to prove the convergence and stability Theorem appropriately.

We organize our paper as follows: In Section 2, we introduce the model of Atherosclerosis presented by Friedman in friedman2015free . We apply some changes to construct a more appropriate model for numerical and proof purposes in Section 3. In Section 4, we use finite difference and collocation method with convenience basis for approximating the solution of the problem. We discuss the stability and convergence of the method used for solving the model in Section 5. Finally, by presenting some numerical results, the theoretical statements are justified in Section 6.

2 Mathematical model

In this paper, we consider the following parabolic free boundary problem modeling plaque growth introduced in friedman2015free as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

There are several mathematical models that describe the growth of a plaque in the artery. All these models recognize the critical role of the “bad” cholesterols and the “good” cholesterols in determining whether a plaque, once formed, will grow or shrink.

In mathematical models, choosing the type of geometry of the problem is very important. Since the plaques grown in the artery are approximately spherically symmetric, in most of the models, it is assumed that the plaque grows radially-symmetric. In this model, it is assumed that the artery is a very long circular cylinder and a circular cross-section $0\leq r\leq 1$ is considered. The plaque is given by $R(t)<r<1$, where $r$ is measured in unit of cm, and $t$ is measured in unit of days. Also, The variables $\widehat{L}$, $\widehat{H}$ and $\widehat{F}$ are taken to be functions of $(r,t)$ only in the region $\{(r,t);R(t)<r<1,t>0\}$. In this model, $\widehat{L}$, $\widehat{H}$ and $\widehat{F}$ are the variables that illustrate the concentration of LDL and HDL and the density of foam cell in the plaque respectively and $v$ is the radial velocity. $k_{1}$ is the rate of ox-LDL ingestion by macrophages. $k_{2}$ is the rate of reverse cholesterol transport. $r_{1}$ and $r_{2}$ represent the degradation of the LDL and HDL caused by radicals respectively. $\mu_{1}$ and $\mu_{2}$ are the death rate of macrophages and foam cells respectively. Also, D is the diffusion coefficients of foam cells. Initially, in order to better understand the model, it is better to know a little how the plaque is calcified in the artery.

A plaque contains low density lipoprotein (LDL) or bad cholesterol, high-density lipoprotein (HDL) or good cholesterol, macrophages and foam cells (Figure 2.a). The process of plaque development begins with a lesion in the endothelial layer, penetration of low density lipoproteins in the intima (Figure 2.b) and by Free radicals in the initma (Figure 2.c) becoming oxidized LDL (Figure 2.d) and it is presented in the second term of (1). Notice that in this model $\widehat{L}$ (and $\widehat{H}$) and its oxidized form $\widehat{L}_{ox}$ (and $\widehat{H}_{ox}$) are merged. Free radicals are oxidative agents continuously released by biochemical reactions within the body, including the intima singh2015role ; bankson1993role . Endothelial cells, sensing the presence of ox-LDL, are activated and trigger monocyte chemoattractant protein, which triggers recruitment of monocytes into the intima (Figure 2.e) bankson1993role ; gerrity1981role . After entering the intima, monocytes differentiate into macrophages (M) (Figure 2.f). The ingestion of large amounts of ox-LDL (Figure 2.g), that shown in the first term of (1), transforms the fatty macrophages into the foam cells (Figure 2.h) bentzon2014mechanisms . Newly formed foam cells secrete chemokines which attract more macrophages, and the plaque is gradually calcified (Figure 2.i). At the same time that LDL enters the intima, high-density lipoprotein also enters the intima and becomes oxidized by free radicals (FR) as presented by the second term of (2). However, oxidized HDL ($H_{ox}$) is not ingested by macrophages. HDL helps to prevent Atherosclerosis by removing cholesterol from foam cells, and by the limiting inflammatory processes that underline Atherosclerosis, as shown in the first term of (2). Furthermore, HDL takes up free radicals that are otherwise available to LDL hao2014ldl . As the plaque continues to grow (Figure 2.i), the increased shear force may cause rupture of the plaque, possibly resulting in the formation of a thrombus (blood clot) and heart attack (for more information about the model see friedman2015free ).

For simplicity, we consider the model from cylindrical coordinates (1)-(5) to spherical coordinates as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the next section, we want to present a new reformulation of the model presented in (6)-(10).

3 A new reformulation of the model

Most of the mathematical models that originate from biological phenomena have some special features. Free or moving boundary condition, mix (robin) boundary condition and many other features of a mathematical model can cause some difficulties in applying classical numerical methods on them and researchers overcome these difficulties by applying suitable techniques. Since the model of Atherosclerosis (6)-10) is a free boundary model with mix boundary condition, to solve this model numerically, we need to remove the mentioned difficulties. To reach this aim, we must consider a new reformulation of this model for which the free boundary change to a fixe one and the robin (mix) boundary condition change to a Neumann one.

3.1 Fixing the domain

Due to the fact that the plaque grows radially symmetric with free boundary, using the front fixing method by the following variable changes

[TABLE]

the free boundary problem (6)-(10) is transformed into a problem with a fixed domain

[TABLE]

and the model becomes as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

3.2 changing the boundary condition

To obtain more comfortable results for numerical analysis, without loss of generality, we can change the mixed boundary condition of the model to a Neumann one by applying suitable variable changes as follows

[TABLE]

So, from (11)-(15) we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

4 Approximating the solution of the problem

In this section, we approximate the solution of the problem (19)-(23) for $0<\rho<1$ and $0<t<T$. Let $t_{i}:=ih\,\,(i=0,1,\ldots,M)$ be mesh points, where $h:=\dfrac{T}{M}$ is the time step and $M$ is a positive integer. The problem is solved using the finite difference-collocation method. In order to prove the stability and convergence of the method, we need to construct a non-classical discretization of second-order to approximate the time derivative. So, let us consider the following discretization formula for approximating the time derivative for a given function $u(\rho,t)$

[TABLE]

and the following approximation formula for linearizing the equations

[TABLE]

where

[TABLE]

and $C$ is a positive constant. In the following, we have assumed that

[TABLE]

Using finite difference method based on the approximation formula given by (24) and equations (11)-(13) we get

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $E_{t}^{u}$ is obtained by merging the errors of $E_{t}$ and $E^{u}$. From (26), there exists positive $C^{*}$ such that

[TABLE]

and using (15) we have

[TABLE]

Then, we conclude that there exists positive constant $C^{**}$ such that

[TABLE]

In the rest of this paper, the solution of the problem (6)-(7) is denoted by $(L_{n+1}^{ap},H_{n+1}^{ap},F_{n+1}^{ap})$, which is the approximated solution of the following problem

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $h^{*}=\dfrac{2h}{3}$ and $(L_{n+1}^{ap},H_{n+1}^{ap},F_{n+1}^{ap})$ is obtained as an approximation of $(L_{n+1},H_{n+1},F_{n+1})$ by solving (32)-(36) employing the collocation method. To implement this method, we employ $\{p_{i}(\rho)\}_{i=0}^{\infty}$ as trial functions as follows

[TABLE]

Then, we approximate $F_{n+1}$ by $F_{n+1}^{N}$ defined as follows

[TABLE]

Now, we consider the following equation

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

where $\Pi_{N}^{0,0}$ and $I_{N}^{0,0}$ are the orthogonal projection and Jacobi-Gauss-Lobatto interpolation operator with respect to $\rho$ on $[0,1]$, respectively. The approximation of the solution of (32) and (33) can be obtained by $L_{n+1}^{N}$ and $H_{n+1}^{N}$ in a similar manner which is used for $F_{n+1}^{N}$.

5 Convergence and stability

As we are aware, the first and most important step in numerical analysis is to prove the stability and convergence of the proposed numerical method. To do so, we need to introduce some mathematical preliminaries.

Definition 1

canuto2006spectral Suppose $\Omega=(a,b)^{d},~{}d\in\mathbb{N}$ and $L^{2}_{w^{\alpha,\beta}}(\Omega)$ is the space of square integrable functions in $\Omega$. Now, we can define the following inner product on $L^{2}_{w^{\alpha,\beta}}(\Omega)$

[TABLE]

[TABLE]

Definition 2

canuto2006spectral Suppose ${\mathbb{P}_{N}}^{d}$ is the space of all $d$ dimensional algebraic polynomials of degree at most $N$ in each variable. $\Pi_{N,w^{\alpha,\beta}}:L^{2}_{w^{\alpha,\beta}}(\Omega)\longrightarrow{\mathbb{P}}_{N}^{d}$ is an orthogonal projection if and only if for any $u\in L^{2}_{w^{\alpha,\beta}(\Omega)}$, we have

[TABLE]

Theorem 5.1

shen2011spectral * Let $\alpha,\beta>-1$. For any $u\in B^{m}_{\alpha,\beta}(I)$,*

[TABLE]

5.1 Convergence

In this section, we want to prove the convergence of the above method under the title of Convergence Theorem. For this purpose, using the principle of mathematical induction, we need to show that for $k=0,1,\ldots,M$, (for an arbitrary $M$) there exist positive constants $L^{*}$, $H^{*}$, $F^{*}$ and $R^{*}$ such that

[TABLE]

where $L_{k}$, $H_{k}$, $F_{k}$ and $R_{k}$ are the exact solutions of (32)-(36) in $t=t_{k}$, respectively. First, we assume that

[TABLE]

In the following, to prove the convergence theorem, we need to present the following lemma.

Lemma 1

Let $F$ be the exact solution of (21) on $[0,1]\times[0,T]$, $F_{n+1}^{ap}=F_{n+1}^{N},~{}\forall n\geq 1$ and $\dfrac{\partial^{2}F}{\partial\rho^{2}}$ be $C^{1}$-smooth function. Then, for each $0<\rho<1$, there exist positive constants $c_{1}$, $c_{2}$ and $c$ such that

[TABLE]

[TABLE]

where $K(N)$ is the error generated by the spectral method and

[TABLE]

Proof

For every $N\in\mathbb{N}$, there exists a polynomial $F^{N}_{1}$ such that

[TABLE]

and

[TABLE]

Now, by taking the inner product of both sides of (39), we conclude that

[TABLE]

[TABLE]

where $F_{n+1,1}^{N}(\rho):=F_{1}^{N}(\rho,t_{n+1})$ and for each $u\in C^{2}[0,1]$, $\mathtt{L}$ is defined as follows

[TABLE]

where $G(\rho)$ is defined in (41) and is a bounded continuous function. Therefore we have

[TABLE]

[TABLE]

[TABLE]

Then, by using Cauchy-Schwarz inequality there exists a constant $c_{5}$ such that

[TABLE]

By combining (40) with (48) we obtain that

[TABLE]

Therefore, from (48) and (49) and using Cauchy-Schwarz inequality and Young inequality, we get that

[TABLE]

[TABLE]

[TABLE]

So there exists a positive constant $c_{4}$ such that

[TABLE]

where

[TABLE]

and $\mathtt{L}$ is defined in (46). On the other hand, from (39) we have

[TABLE]

where

[TABLE]

Then, by applying the recurrence relation (52) repeatedly in (50) we can conclude that

[TABLE]

where $\|E_{F}^{*}\|_{\infty}$ and $K_{1}(N)$ are as in (44).

$\blacksquare$

Lemma 2

Let $F_{n+1}^{ap}=F_{n+1}^{N},~{}\forall n\geq 1$ and $\dfrac{\partial^{2}F}{\partial\rho^{2}}$ be $C^{1}$-smooth function. Then, for each $0<\rho<1$, there exist positive constants $M^{*}$ and $C^{*}$ such that the following inequality holds

[TABLE]

where $L^{*},H^{*}$ and $F^{*}$ are obtained from (42) and $\|E_{F}^{*}\|_{\infty}$ and $K_{1}(N)$ are as in (44).

Proof

By taking the inner product of both sides of (39), we can deduce that

[TABLE]

[TABLE]

Using (53), (42), Theorem 5.1, Definition 2, Cauchy-Schwarz inequality and Young inequality, we have

[TABLE]

where $c_{5}$ is obtained from (48). $\blacksquare$

Finally, using Lemmas 1 and 2 we have

[TABLE]

where $\|E_{F}^{*}\|_{\infty}$ and $K_{1}(N)$ are as in (44).

Similar to the proof of lemma 1, we can show that there exist positive constants $c_{6},\,c_{7},\,c_{8},\,c_{9}$ such that for $0\leq t\leq T,$

[TABLE]

[TABLE]

where

[TABLE]

and $c$ is a positive constant.

In the following theorem, the convergence of the proposed method is proved.

Theorem 5.2

(Convergence Theorem) Let $L_{n+1}^{ap}=L_{n+1}^{N}$, $H_{n+1}^{ap}=H_{n+1}^{N}$ and $F_{n+1}^{ap}=F_{n+1}^{N}$. Under the assumption of Lemma 1 there exist positive constants $M_{1},M_{2}$ such that

[TABLE]

where

[TABLE]

[TABLE]

Proof

From (56), (57) and (58) we can conclude that there exists a constant $M_{4},$ such that

[TABLE]

where

[TABLE]

[TABLE]

and $\|E_{F}^{*}\|$ is as in (44). Therefore, there exists a constant $M_{3}$ that

[TABLE]

Then, by applying the above recurrence relation we have

[TABLE]

Finally, it may be concluded that there exist constants $M_{1}$ and $M_{2}$ such that

[TABLE]

$\blacksquare$

Employing the General Sobolev inequalities, there exists a positive constant $M_{5}$ such that for $0\leq t\leq T$,

[TABLE]

Using (60) and (56), we can choose proper $N$ and $h$ such that

[TABLE]

Similar to the proof of the Theorem 5.2 and (61), we can show that

[TABLE]

and

[TABLE]

Thus, using (61), (62), (63) and Theorem 5.2, we can deduce that the sequence $\{L_{n},H_{n},F_{n},R_{n}\}_{n=0}^{\infty}$ converges to the exact solution of the problem (6)-(10) on $[0,1]\times[0,T]$. Now, using the principle of mathematical induction, Theorem (5.2), (42), (43), (61), (62) and (63), we have

[TABLE]

5.2 Stability

In this section, we want to prove the stability of the presented method. For this purpose, first, we consider the perturbed problem as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the following theorem, the stability of the proposed method is proved.

Theorem 5.3

Let $\epsilon_{1}$ be a positive constant and $|p_{i}|<\epsilon_{1}\,\,(i=1,\ldots,4)$. Then, under the assumptions of Lemma 1, there exist positive constants $M_{1}^{*},\,M_{2}^{*}$such that

[TABLE]

where

[TABLE]

Proof

If we solve the perturbed problem (65)-(70) using the presented method, one can conclude that there exist a positive constant $M_{4}^{*}$ such that

[TABLE]

where $\|E\|_{\infty}\leq M^{*}{h^{*}}^{2}$ and $M^{*}$ is a positive constant and also,

[TABLE]

[TABLE]

So we get

[TABLE]

then, by applying the above recurrence relation we have

[TABLE]

Therefore, there exist costats $M_{1}^{*}$ and $M_{2}^{*}$ such that we have

[TABLE]

$\blacksquare$

6 Numerical experiments

The main target of this section is to investigate the numerical solution of the model of Atherosclerosis from two perspectives: 1-numerical solution point of view 2-biological simulation point of view. We first examine the numerical results from the first perspective. We solve the model of Atherosclerosis by applying the finite difference/collocation method. The most important question to ask is how to construct trial functions which satisfy the boundary condition to establish the collocation method. The first and simple way that comes to mind is to use a linear combination of monomial polynomials. So, we approximate the functions $L(\rho,t),H(\rho,t),F(\rho,t)$ in the form of (37) as follows

[TABLE]

where

[TABLE]

which will denote by ”TFBM” (Trial Functions Based on Monomials). Also, the Gauss quadrature points $\{x_{i}^{0,0}\}_{i=1}^{N}$ (i.e., the zeros of Legendre polynomial of degree $N+1$) are considered as collocation points. The typical parameter values and the initial conditions for our numerical simulations are: $k_{1}=10,\,k_{2}=10,\,K_{1}=10^{-2},\,K_{2}=0.5,\,D=8.64\times 10^{-7},\,\mu_{1}=0.015,\,\mu_{2}=0.03,\,r_{1}=2.42\times 10^{-5},\,r_{2}=5.45\times 10^{-7},\,\lambda=2.573\times 10^{-3},\,\delta=-2.541\times 10^{-3},\,M_{0}=5\times 10^{-5},\,\alpha=1,\,\beta=0.01,\,L_{0}=14\times 10^{-4},\,H_{0}=6\times 10^{-4},\,F_{0}=0$, which are taken from hao2014ldl ; friedman2015free . Notice that in our simulations, the radius of initial plaque $\epsilon$ in (5) is considered to be 0.9 cm.

Now, in order to verify our numerical results, we need to present the following definition.

Definition 3

A sequence $\{x_{n}\}_{n=1}^{\infty}$ is said to converge to $x$ with order $p$ if there exists a constant $C$ such that $|x_{n}-x|\leq Cn^{-p},~{}~{}\forall n$. This can be written as $|x_{n}-x|=\mathcal{O}(n^{-p})$. A practical method to calculate the rate of convergence for a discretization method is to use the following formula

[TABLE]

where $e_{n_{1}}$ and $e_{n_{2}}$ denote the errors with respect to the step sizes $\dfrac{1}{n_{1}}$ and $\dfrac{1}{n_{2}}$, respectively gautschi1997numerical .

It is valuable to point out that our numerical calculations are carried out using the MATLAB 2018a program in a computer with the Intel Core i7 processor (2.90 GHz, 4 physical cores).

Lies in the fact that it is hard or sometimes impossible to reach the exact solution of most of the coupled nonlinear models analytically, and because we have shown that the presented method for solving the model is stable and convergent (See Theorem 5.2), so, we consider numerical results for the large $M=300$ and $N=10$ as an exact solution and for other values of $M$, we report the time-error in Table 1. To better see the time-error of the presented approach numerically, we report the obtained results in Figure 3. As we see in Table 2 and Figure 4, the non-classical finite difference method presented in (24) has almost $\mathcal{O}(h^{2})$ error, which verifies our theoretical results presented in Theorem 5.2.

[TABLE]

In Table 3, we have presented the maximum space-errors with $M=100$ and various values of $N$. To better see the space-error of numerical results, Figure 5 is presented.

It is worth to note that in the numerical results, the discrepancy between an exact value and some approximation to it, is called ”maximum error” and is denoted by $E_{\infty}$. In Table 4, the condition numbers (CN) of the coefficient matrices of the collocation method are shown. In this table, high values of condition numbers are highlighted in bold. Because of the fact that the condition number of the coefficient matrices grows very fast when $N>10$, unfortunately, the Matlab software can not accurately extract the results. To overcome this difficulty, we need to propose proper trial functions which reduce the condition number significantly. In this case, we decide to consider Legendre polynomials (Jacobi polynomials with $\alpha=\beta=0$) to construct trial functions for the collocation method doha2018fully . So, let $L_{n}(x)$ be the Legendre polynomial of degree $n$ and set

[TABLE]

where the constants $c_{n}$ and $d_{n}$ are uniquely determined in such a way that $p_{n}(x)$ satisfies the boundary conditions; in other word, $\dfrac{\partial p_{n}}{\partial x}(\pm 1)=0,~{}~{}\forall n\geq 0$.

According to the features of Legendre polynomials we have

[TABLE]

therefore, by solving (73) one can easily conclude that

[TABLE]

Thus (72) becomes as follows

[TABLE]

So we can approximate the functions $L(\rho,t),H(\rho,t),F(\rho,t)$ in the form of (37) as follows

[TABLE]

where

[TABLE]

which we will denote by ”TFBL” (Trial Functions Based on Legendre polynomials).

Remark 1

The scaling factors in the trial functions (76) play the role of precondition factor for the collocation matrices and reduce the condition number. huang2018spectral ; khosravian2017generalized .

[TABLE]

It should be noted that in order to use Legendre polynomials, we map the domain of the problem (19)-(23) to $[-1,1]$.

We also point out that the Gauss quadrature points $\{x_{i}^{0,0}\}_{i=1}^{N}$ are considered as the collocation points. As mentioned in TFBM case, because of the stability and convergence of the presented method, we consider the solution of the problem with $N=100$ and $M=1000$ as an exact solution. In Table 5, we have presented the maximum time-error $N=51$ and various values of $M$. It is observed that the time-error of the presented approach numerically, we report the obtained results in Figure 6. As we can see from Table 6 and Figure 7, the non-classical finite difference method presented in (24) has almost $\mathcal{O}(h^{2})$ error, which verifies our theoretical results presented in convergence Theorem 5.2. In Table 7, we illustrate the maximum space-errors $M=200$ and various values of $N$ by considering the solution of the problem with $N=200$ and $M=200$ as an exact solution. To better observe the space-error of numerical results, Figure 8 is presented. Also, the condition number of the coefficient matrices constructed by collocation method for equations (1)-(4) are shown in Table 8 and Figure 9; which shows that, as we expect from TFBL, the condition number of matrices do not increase significantly by increasing $N$ even for $N=150$ . For a better comparison of condition numbers in both TFBM and TFBL cases, Figure 10 is presented. One can easily see that using Legendre polynomials, the condition number of the coefficient matrices in the collocation method is significantly less than TFBM.

Now, in this position, the examination of the numerical solutions from the perspective of biology and simulation is reported by presenting the rate of tumor growth with three various values of pairs $(L_{0},H_{0})$. The results are presented in Table 9 and Figure 11, and they are compared to the risk map illustrated in figure 12 and Retrieved from friedman2015free . As we expect, the level of $L_{0}$ and $H_{0}$ in blood directly affects the growth and shrink of the plaque, that means for the values of $(L_{0},H_{0})$ below the ”zero growth”, the plaque grows as shown in Figure 13, and for the values of $(L_{0},H_{0})$ above the ”zero growth” the plaque shrinks, as shown in Figure 14. To better see the radius changes, a small part of the plaque in the vessel is magnified and is presented in the figures mentioned above. It is noteworthy that arrows in these two figures indicate the direction of growth or shrink of the plaque.

N Eq. (1) Eq. (2) Eq. (3)

10 3.245 3.245 4.096

20 17.475 17.475 7.230

40 1.302e+02 1.302e+02 13.61

80 1.022e+02 1.022e+02 26.14

100 1.990e+02 1.990e+02 32.05

Table 8: Condition number of the coefficient matrices.

7 Conclusion

There are many mathematical methods for solving biological models. However, mathematical modeling often produces nonlinear differential equations. Therefore, we cannot always obtain the exact solution of these equations; so, developing numerical techniques to solve these equations is required. In this study, a mathematical model of Atherosclerosis is solved numerically and the convergence and stability analysis are presented. For the reader’s convenience, we give the main contributions of this study as follows

$\bullet$ In this article, we use the front fixing method to convert the moving boundary problem (6)-(10) to a fix one (11)-(15), because classical numerical methods are not effective to solve free and moving boundary problems and moreover, because of the suitability of the front fixing method to apply to problems with regular geometries along with the mesh-based methods.

$\bullet$ To achieve more comfortable results for numerical analysis, we have simplified the model by changing the mix boundary condition of the equations (11)-(13) to a Neumann one by applying a suitable change of variables (17)-(18).

$\bullet$ To solve nonlinear systems, one way is to linearize the equations of the system and then, solve the linear one by classical methods. Since the model studied in this article is nonlinear, we have used Taylor theorem simultaneously both to linearize the equations and for constructing new second-order non-classical discretization formula to approximate time discretization (Finite difference method).

$\bullet$ In this article, we use spectral collocation method in space. To construct trial functions which satisfy the boundary conditions, the first way that comes to mind is to use a linear combination of monomial polynomials. But, there are some problems to use this kind of polynomials. some of these problems are as follows

    1. The condition number of collocation matrices increases significantly by increasing the size of the matrices for $N\geq 10$ and because of that, we are limited to increase the collocation points and achieve arbitrary errors. (See Table 4 and Table 3)

    2. The high error of collocation method affects the error of the finite difference method (See Table 1).

To the above reasons, we use a linear combination of classical orthogonal polynomials or orthogonal functions to construct trial functions.

So, using orthogonal polynomials, the condition number of collocation matrices does not increase significantly by increasing the collocation points even for $N=150$ (See Table 8, Figure 9 and Figure 10), and because of that, the error of the collocation method decreases compared to the TFBM case (See Table 7 and Figure 8).

$\bullet$ Moreover, the convergence and stability of the presented method is proved (See Theorem 5.2 and Theorem 5.3) and the order of convergence is presented. Numerical results in both TFBM and TFBL cases show that the finite difference method displays an $\mathcal{O}(h^{2})$ order of convergence, as we expect from convergence Theorem 5.2 (See Figure 4 and Figure 7), and the space-error shows that using the collocation method, the results are converging to the exact solution.

$\bullet$ As illustrated in Figure 11 and Table 9, we present the effect of the level of $L_{0}$ and $H_{0}$ in blood on the growth and shrink of the plaque. It is easily can be seen that for $(L_{0},H_{0})=(150,45)$ (below the ”zero growth”), the plaque grows and for $(L_{0},H_{0})=(120,60)$ (above the ”zero growth”) the plaque shrinks.