A mathematical model of aortic aneurysm formation.

Abdominal aortic aneurysm (AAA) is a localized enlargement of the abdominal aorta, such that the diameter exceeds 3 cm. The natural history of AAA is progressive growth leading to rupture, an event that carries up to 90% risk of mortality. Hence there is a need to predict the growth of the diameter of the aorta based on the diameter of a patient's aneurysm at initial screening and aided by non-invasive biomarkers. IL-6 is overexpressed in AAA and was suggested as a prognostic marker for the risk in AAA. The present paper develops a mathematical model which relates the growth of the abdominal aorta to the serum concentration of IL-6. Given the initial diameter of the aorta and the serum concentration of IL-6, the model predicts the growth of the diameter at subsequent times. Such a prediction can provide guidance to how closely the patient's abdominal aorta should be monitored. The mathematical model is represented by a system of partial differential equations taking place in the aortic wall, where the media is assumed to have the constituency of an hyperelastic material.

Introduction

AAA is an abnormal dilatation most commonly of the infrarenal aorta. One definition of AAA is a diameter greater than 3 cm [1]. The clinical significance of AAA stems from the high mortality associated with rupture. Approximately 60% of patients with ruptured AAA die before reaching the hospital [2] and mortality rates of emergency surgical repair are as high as 35–70% [3]. The risk of rupture increases with AAA diameter. The pathogenesis of AAA is largely unknown and likely multifactorial. Diameter remains the only clinically useful and available marker for risk of rupture. Surgical repair is recommended for aneurysms measuring greater than 5.5 cm [4, 5] and there is insufficient evidence to recommend surgery for all patients with smaller AAA [2]. However, even elective surgical repair of AAA remains a major operation, incurring significant morbidity and mortality particularly in older, sicker patients where there is a higher prevalence of AAA. Furthermore, patients die every year from rupture of aneurysms smaller than 5.5 cm, and up to 60% of AAA larger than 5 cm remain stable [6]. Thus, some patients with smaller AAA are denied lifesaving surgery and others with larger AAA undergo unnecessary major surgery. The reasons why some smaller aneurysms go on to rupture while some larger ones remain stable are not understood. Techniques that provide early identification of small AAA with increased risk for rupture and large AAA with low risk for rupture will improve overall mortality by prompting personalized treatment plans for AAA.

There are a number of mathematical papers that describe the dynamics of the weakening and dilation of the arterial wall and the risk of rupture. The more comprehensive models include the nonlocal, nonlinear elastic nature of the arterial wall in response to wall shear stress [7-10]. However these models do not explain the clinical fact that some abdominal aortae rupture when the diameter of the cross section of the dilated aorta is less than 5.5 cm while others do not rupture even when much larger than 5.5 cm. There is a need to discover non-invasive biomarkers that will provide the following prognosis for a patient undergoing initial screening for AAA:

Given initial diameter of the cross section of the aorta, say R0, how will the diameter, say R(t), evolve with time?

How will the risk of rupture increase with time?

The present paper develops a mathematical model that addresses the first question, with IL-6 as the serum biomarker. IL-6 has already been suggested as prognostic biomarker for AAA in [11, 12]; additional potential biomarkers are reviewed in [13, 14].

The mathematical model consists of two parts. The first part focuses on the biology. It introduces a network of cells and cytokines which, in disease state, lead to the destruction of smooth muscle cells (SMCs) and elastin depletion in the media, and disruption of the extracellular matrix (ECM) in the adventitia. The biological predictions of this part of the model, in terms of the expression levels of the proteins involved, are in agreement with patients data [15]. The second part of the model develops the mechanics of diameter growth of the artery as a result of deficiency in elastin and disruptions in the ECM. Here, instead of the approach taken in [7-10], we use the fact that the arterial wall behaves like a hyperelastic material [16]. Hence its elastic strength and displacement under stress depend on the smooth muscle cells in the aortic wall. We can therefore relate the wall dilation to several cytokines whose expression is upregulated in disease state. Such cytokines could be used as potential biomarkers for predicting the arterial wall dilation. We shall focus on IL-6, which is upregulated in AAA patients [11, 12].

The biological network introduced in the present paper includes macrophages, T cells, SMCs and fibroblasts, as well as cytokines such as IL-6, IL-10, IL-12, TNF-α and IFN-γ; MMP and TIMP are also included. The densities of the cells and the concentrations of the proteins satisfy a coupled system of partial differential equations (PDEs) in the arterial wall.

In the early stage of aneurysm, endothelial cells, under hemodynamic shear stress, secrete major attractant protein, MCP-1, and IL-6 [17]. IL-6 is also produced by SMCs (in conjunction with IL-1) [18]. This triggers recruitment of monocytes from the blood into the media layer of the arterial wall [19]. The monocytes mature into macrophages, while, at the same time, additional macrophages from the adventitia are also chemoattracted by MCP-1, IL-6 and IL-8 [20-22]. More precisely, IL6/MCP-1 combined, and MCP-1 and IL-8 separately chemoattract monocytes; for simplicity, we shall represent the chemoattractants by MCP-1 and IL-6. T cells are activated by contact with macrophages in the presence of IL-12, and macrophages are activated by IFN-γ produced by the T cells [23]. Macrophages produce TNF-α [24], MMP and TIMP [24, 25], and IL-10, IL-12 and IL-6 [26], but MMP is also produced by SMCs [27]. Fibroblasts produce collagen [28], and the collection of MMP, TIMP and collagen weaken the ability of the adventitia layer to withstand stress. Macrophages are known to cause apoptosis in SMCs [19], and this leads to reduction in elastin [29], thus weakening the elastic strength of the media; The apoptotic SMCs are known to produce MCP-1 [19].

Fig 1 shows SMCs residing in the media, fibroblasts residing in the adventitia, while macrophages and T cells are present in both layers of the aortic wall.

Fig 1

A schematic of the aortic wall including the aortic lumen, the intima, media, and adventitia.

Fig 2 is a schematics of the network within the arterial wall during aneurysm; all the cytokines and the ECM (elastin and collagen) are present in both the media and the adventitia. A schematics of the balloon-like bulge geometry in the aorta is shown in Fig 3.

Fig 2

Interaction network among cells and cytokines in their respective layers.

SMCs remain in the media, and fibroblasts remain in the adventitia.

Fig 3

Schematic representation of a 2D section of the computational domain.

Ω is blood vessel; Ω represents the media and Ω the adventitia; Γ, Γ and Γ are three free boundaries. Γ is the inner surface of the artery, Γ is the outer surface, and Γ is the surface between the media and adventitia.

1 Mathematical model of the biology

In this section we develop a mathematical model of aneurysm based on the diagram shown in Fig 2. The model, represented by a system of PDEs, includes the variables listed in Table 1. We assume that all cells are moving with a common velocity v; the velocity is the result of movement of macrophages, T cells and SMCs into the media and adventitia. We also assume that all species are diffusing with appropriate diffusion coefficients. The equation for each species of cells, X, has a form where the expression on the left-hand side includes advection and diffusion, and F accounts for various biochemical reactions, and chemotaxis. The equations for the chemical species are the same as for cells, but without the advection term, which is relatively very small compared to their large diffusion coefficients. Fig 3 shows a 2D projection of a blood vessel Ω with an aortic aneurysm Ω ∪ Ω.

Table 1

The variables of the model; concentration and densities are in units of g/cm3.

S:	density of SMCs	T:	density of T cells (Th1)
M:	density of macrophages	f:	density of fibroblasts
ρ:	density of ECM	P:	concentration of MCP-1
Tα:	concentration of TNF-α	Iγ:	concentration of INF-γ
I6:	concentration of IL-6	I10:	concentration of IL-10
I12:	concentration of IL-12	Q:	concentration of MMPs
Qr:	concentration of TIMP	p:	pressure in media (in g/cm/day2)
u:	displacement (in cm)	v:	cell velocity in media and adventitia (in cm/day)

Equation for MCP-1 (P)

The equation for MCP-1 is given in the whole domain (media Ω, and adventitia Ω) by The first term on the right-hand side is the production of P by apoptotic SMCs, which is induced by macrophages [19]. The second term represents lost of MCP-1 as it binds to, and internalized by, macrophages which are chemoattracted to MCP-1; the internalization of MCP-1 may be limited due to the limited rate of receptor recycling. The last term in Eq (1) accounts for degradation of MCP-1.

Equation for macrophages (M)

According to [26] the two phenotype of macrophages, M1 and M2, are both present in aneurysm; M1 secretes inflammatory cytokines such as IL-12 and TNF-α, while M2 secretes immunoregulatory cytokines such as IL-10. For simplicity, we do not differentiate between the two phenotypes, but take into account that the productions of IL-12 and TNF-α are resisted by IL-10 [30].

The evolution of the macrophage density in Ω ∪ Ω is modeled by The first two terms on right-hand side account for the recruitment of macrophages by MCP-1 [19] and IL-6 [20-22]. The third term accounts for the activation of macrophages by IFN-γ [23], which is enhanced by TNF-α [31].

Equation for T cells (T)

The density of T cells in Ω ∪ Ω satisfies the equation We assume that T cells are activated by contact with macrophages in IL-12 environment, while IL-12 is inhibited by IL-10 [30].

Equation for SMCs (S)

The equation of the SMCs density in Ω is given by The second term of the right-hand side accounts for apoptosis of SMCs, caused by macrophages through a FasL/Fas-Caspase8-RIP1 mediated mechanism [19].

Equations for cytokines (in Ω ∪ Ω)

IL-6 is produced by macrophages [26], and by SMCs (in conjunction with IL-1) [18]; hence

IL-10 is produced by macrophages [26], so that

IL-12 is produced by macrophages [26], a process inhibited by IL-10. The equation for I12 is given by

Similarly, TNF-α is produced by macrophages, a process inhibited by IL-10 [32], so that

IFN-γ is produced by T cells [33], hence

MMP is produced by macrophages [25] and by SMCs [34], a process enhanced by TNF-α [24], and MMP is lost by binding with TIMP [25]. Accordingly, Q satisfies the following equation: TIMP is produced by macrophages [25], and is lost by binding to MMP, so that

Equation for ECM (ρ)

SMCs produce elastin [29] and fibroblasts produce collagen [28]. We assume that the density of fibroblasts is approximately constant, and that the ECM concentration is proportional to the combined concentrations of collagen and elastin, which are produced by fibroblasts and SMCs, respectively. Hence, where the term d Qρ represents ECM degradation by MMP. Here we used the notations: X+ = X if X > 0, X+ = 0 if X ≤ 0, and χ = 1 in Ω, χ = 0 in Ω.

Boundary conditions

Fig 3 shows a 2-D projection of the media and adventitia layers of the arterial wall. We shall take as our computational domain the right region lying between y = y1 and y = y2, and impose periodic boundary condition at y = y1 and y = y2.

T cells, which are primarily CD4+ Th1 cells [35], migrate to the adventitia and intima [36, 37], and so do macrophages [22, 37].

The concentration ρ is continuous across the media/adventitia membrane, while S satisfies no-flux boundary condition, , at this membrane. All the other variables satisfy the following flux conditions at the interface Γ: where X and X are the concentrations of cytokines or cell densities in media and adventitia, respectively. We take γ = 50 for cells, and γ = 500 for the cytokines [38].

On Γ and Γ, we have where .

Both MCP-1 and IL-6 are produced by endothelial cells [17], which lie near the inner boundary of the media, Γ. Hence where P0 and I60 are the concentrations of MCP-1 and I6 in the blood. We assume no-flux boundary conditions on Γ for all remaining variables. We also assume no-flux conditions on Γ for all variables (except M and T, as stated above)

Initial conditions

We take initially , where I60 = 6 × 10−9 gcm−3 in the source of influx of I6 from the endothelial cells [39], and as in [40], All other cell densities and cytokine concentrations are taken to be zero initially. We prescribe the initial geometry of aneurysm by

Mechanical hyperelasticity model

We set Ω = Ω ∪ Γ ∪ Ω and refer to this region, briefly, as the aortic wall. Aortic wall is assumed to have nonlinear, hyperelastic material properties [27, 41, 42]. Specifically, it is an ‘almost’ incompressible, homogeneous, and isotropic material with energy density function W of the form where β1 and β2 are forces that represent elastic coefficients; I is the first invariant of the Left Cauchy-Green tensor B (namely I = tr(B)), B = F F, while F is the deformation gradient tensor. In terms of the displacement vector u = (u), we can write , where I is the identity tensor.

The Cauchy stress tensor is given by where p is the hydrostatic pressure. The momentum equation is then where f is the body force and d is the tissue density. In the sequel we assume no body force and neglect tissue inertia, so that

The coefficients β1 and β2 are functions of the SMCs which we take to represent the elasticity coefficients. We shall use the following specific relations in Eq (14): The pressure, p, is a function of ECM, which we take to have the form for some parameters β, p*.

For boundary conditions, we take where σ is the stress tensor from the blood and κ is the mean curvature.

Results

Fig 4 shows the profile of average of cell densities and cytokines concentrations in the first 500 days, and Fig 5 shows how the aneurysm deforms at day 500.

Fig 4

Simulation results over 500 days.

The parameter used are as in Tables 2 and 3 with I60 = 6 × 10−9 g/ml.

Fig 5

Aneurysm deformation at day 500 with initial deformation shown in (A). The parameters used are as in Tables 2 and 3 with (B) I60 = 6 × 10−9 g/ml, and (C) I60 = 6 × 10−8 g/ml. The diameter grows from 2 cm to 2.5 cm in case (B), and to 4.5 cm in case (C).

Table 2

Parameters’ description and value.

Parameter	Description	Value
DM	dispersion coefficient of macrophages	8.64 × 10−7 cm2 day−1 [40, 46, 47]
DT	diffusion coefficient of T-cell	8.64 × 10−7 cm2 day−1 [40, 46, 47]
DIγ	diffusion coefficient of IFN-γ	1.08 × 10−2 cm2 day−1 [40, 46, 47]
DTα	diffusion coefficient for TNF-α	1.29 × 10−2 cm2 day−1 [47]
DS	diffusion coefficient of SMCs	8.64 × 10−7 cm2 day−1 [40, 46]
DP	diffusion coefficient of MCP-1	1.728 × 10−1 cm2 day−1 [40, 46]
DQ	diffusion coefficient of MMP	4.32 × 10−2 cm2 day−1 [40, 46]
DQr	diffusion coefficient for TIMP	4.32 × 10−2 cm2 day−1 [40, 46]
DI6	diffusion coefficient of IL-10	1.08 × 10−2 cm2 day−1 [40, 46, 47]
DI10	diffusion coefficient of IL-10	1.08 × 10−2 cm2 day−1 [40, 46, 47]
DI12	diffusion coefficient of IL-12	1.08 × 10−2 cm2 day−1 [40, 46, 47]
λTαM	activation rate of TNF-α due to macrophage	2.86 × 10−3 day−1 [47]
λMIγ	activation rate of macrophages by IFN-γ	0.005 day−1 [40, 46]
λT	activation rate of Th1 cells by IL-12	6 × 10-4 day−1 [47]
λI6M	production rate of IL-6 by macrophages	1.73 × 10−6 day−1 [16] & estimated
λI6S	production rate of IL-6 by SMCs	1.73 × 10−5 day−1 [18] & estimated
λIrT	production rate of IFN-γ by T cells	2.34 × 10−6 day−1 [47]
λI12M	production rate of IL-12 by macrophages	3.78 × 10−3 day−1 [47]
λI10M	production rate of IL-10 by macrophages	2 × 10−3 day−1 [47]
λQM	production rate of MMP by macrophages	3 × 10−4 day−1 [48]
λQS	production rate of MMP by SMCs	2.16 × 10−5 day−1 [40, 46]
λQTα	production rate of MMP by TNF-α	2 fitted
λQrM	production rate of TIMP by macrophages	6 × 10−5 day−1 [40, 46]
λPS	activation rate of MCP-1 due to SMCs	3 × 10−3 day−1 fitted
λρf	activation rate of ECM due to fibroblasts	3 × 10−4 day−1 [40]
λρS	activation rate of ECM due to SMCs	1 × 10−1 day−1 fitted
dM	death rate of macrophages	0.015 day−1 [40, 46, 47]
dρ	degradation rate of ECM	0.37 day−1 [40, 46, 47]
dP	degradation rate of MCP-1	1.73 day−1 [40, 46, 47]
dPM	internalization rate of MCP-1 by macrophages	2.08 × 10−4 day−1 [40, 46, 47]
dQQr	binding rate of MMP to TIMP	4.98 × 108 cm3 g−1 day−1 [40, 46, 47]
dQrQ	binding rate of TIMP to MMP	1.04 × 109 cm3 g−1 day−1 [40, 46, 47]
dQ	degradation rate of MMP	4.32 day−1 [40, 46, 47]
dQr	degradation rate of TIMP	21.6 day−1 [40, 46, 47]
dρQ	degradation rate of ECM due to MMP	2.59 × 107 cm3 g−1 day−1 [40, 46, 47]
dIγ	degradation rate of IFN-γ	0.69 day−1 [46]
dS	death rate of SMC	0.86 day−1 [40, 46, 47]
dSM	apoptosis rate of SMC by macrophages	1.72 day−1 fitted
dI12	degradation rate of IL-12	1.188 day−1 [40, 46, 47]
dTα	degradation rate of TNF-α	55.45 day−1 [40, 46, 47]
dI6	degradation rate of IL-6	0.173 day−1 [43]
dI10	degradation rate of IL-10	16.64 day−1 [47]

Table 3

Parameters’ description and value.

Parameter	Description	Value
χC	chemotactic sensitivity parameter by MCP-1	10 cm5 g−1 day−1 [40, 46, 47]
KM	macrophage saturation	5 × 10−5 gcm−3 [40]
KP	MCP-1 saturation for influx of macrophages	5 × 10−9 gcm−3 [40]
KIγ	IFN-γ saturation for activation of macrophages	1 × 10−11 gcm−3 [40]
KI10	IL-10 saturation	2 × 10−7 g/cm3 [47]
KTα	TNF-α saturation	5 × 10−7 g/cm3 [47]
ρ0	ECM saturation	10−3 gcm−3 [40]
M0	source/influx of macrophages from blood	5 × 10−5 gcm−3 [40]
T0	source/influx of T cells into intima	1 × 10−3 gcm−3 [40]
P0	MCP-1 concentration	3 × 10−10 [40]
α˜	influx rate of macrophages into interstitium	0.2 cm−1 [40]
αT	influx rate of T cell into interstitium	0.2 cm−1 [40]
αI6	influx rate of IL-6 into interstitium	0.2 cm−1 [40]
αP	influx rate of MCP-1 into interstitium	0.2 cm−1 [40]
γ	influx rate of media into adventitia	0.1 cm−1 [40]
S0	initial steady state of SMCs	6 × 10−3 gcm−3 [40]
I60	source/influx of IL-6	6 × 10−9 ∼ 6 × 10−8 gcm−3 [39]
β10	Coefficient of AAA tissue material properties	17.4 N/cm2 [42] and estimated
β20	Coefficient of AAA tissue material properties	188.1 N/cm2 [42] and estimated
k1	Coefficient of AAA tissue material properties	30.4 N/cm2 [42] and estimated
k2	Coefficient of AAA tissue material properties	84 N/cm2 [42] and estimated
βp	pressure parameter	18.0 N/cm2 [44] and estimated

From Fig 4, we see that, as AAA progresses, the average densities of SMCs and ECM decrease. The concentration of the inflammatory cytokines TNF-α, IFN-γ and MCP-1 are increasing. On the other hand, the average concentration of MMP is increasing only for the first 200 days, and it then decreases as the concentration of TIMP keeps growing. The average concentration of IL-6 in the tissue decreases and nearly stabilizes after 500 days. Fig 5 shows growth of the aorta bulge after 500 days: from diameter of 2 cm to diameter of nearly 2.5 cm when I60 = 6 × 10−9 g/ml and nearly 4.5 cm when I60 = 6 × 10−8 g/ml.

In Figs 4 and 5, the initial values of the pro-inflammatory cytokines were taken to be below their steady state. On the other hand the SMCs and ECM densities were taken to be above their steady state. Other choices of the initial conditions yield similar profiles (not shown here).

If at the initial screening for AAA we measure both the aorta diameter, R0, and the concentration of IL-6 in the blood, I60, then by simulating the model with these values of R(0) = R0 and I60 (in Eq (13)), we can predict the AAA diameter at any future time T. Figs 6 and 7 simulate the cases T = 300 and T = 500 days. Thus, for every pair of I60 (on the horizontal axis) and R0 (on the vertical axis) the color on the column shows the diameter corresponding to the point (I60, R0) at day T = 300 (upper) and T = 500 (lower).

Fig 6

AAA growth.

The x-axis scales I60 from 6 × 10−9 g/ml to 6 × 10−8 g/ml; the y-axis scales R0 from 2 cm to 3 cm; The color represents the diameter of aortic bulge at 300 day.

Fig 7

AAA growth.

The x-axis scales I60 from 6 × 10−9 g/ml to 6 × 10−8 g/ml; the y-axis scales R0 from 2 cm to 3 cm; The color represents the diameter of aortic bulge at 500 day.

Conclusions

Abdominal aortic aneurysm is a localized enlargement of the abdominal aorta, which may lead to rupture of the aorta. The disease is asymptomatic until rupture, which is nearly always fatal. The risk of rupture associated with the increased diameter of the aorta varies among people. Hence it is difficult to determine, on the basis of just the current diameter (R0) of a patient’s aortic wall, how closely this patient’s abdominal aorta should be monitored. In order to make such a determination, we need address the following questions:

how fast the diameter will grow over time?

how does the risk of rupture depend on the growing diameter?

Several non-invasive prognostic biomarkers for AAA have been suggested [13, 14]; specifically, IL-6 [11, 12].

The present paper addresses the first question with a mathematical model: Given R0 and the serum concentration of IL-6, the model predicts the dynamical growth of the diameter, R(t), for any future time t. The mathematical model includes the basic biology underlying AAA formation as a disease that degrades the elastic strength of the aortic wall. Simulations of the model can show how the initial diameter R0 and the serum concentration of IL-6 determine the diameter R after a given period of time T; This is done in Fig 6 in the cases of T = 300 days and T = 500 days. Thus, the level of IL-6 concentration in the blood can be used to suggest how soon to schedule the next screening of the patient’s abdominal aorta.

The model is represented by a system of partial differential equations, based on some simplified assumptions. In particular, the model assumes that the disease is associated with initial inflammation, and that the aortic wall is a ‘pure’ hyperelastic material. Furthermore, some of the parameters were estimated crudely, for the lack of available data. The authors plan to follow up this study with patients data, which are unavailable at this time. The patients will be followed every 6 months with CT scan and measurement of their aneurysm, and of course any patient who comes for surgery or rupture will be recorded. When patients data become available, the model could then be refined, in particular by adjusting some parameters for specific groups of patients. The model could then be used as a predictive tool for monitoring the enlargement of the abdominal aorta.

Parameter estimates

All the parameter values are listed in Tables 2 and 3. Most parameters are taken from the literature while some are estimated below. A few remaining parameters are fitted.

λ and λ

According to [16], 106 macrophages produced 10 ng/ml of IL-6. Assuming steady state, we have Since d = 0.173 day−1 [43], we get λ = 1.73 × 10−6 day−1.

According to [18], 4 × 104 cells/ml SMCs produced 4 ng/ml of IL-6. Assuming steady state, we have Since d = 0.173 day−1, we get λ = 1.73 × 10−5 day−1.

β10, β20, k1 and k2

By [42], β1 = 0.022 N/mm2 and β10 = 0.174 N/mm2. We assume that in aneurysm S is approximately half its normal density, i.e., . Then by Eq (17), which gives k1 = 30.4 N/cm2.

Similarly, by [42], β2 = 1.461 N/mm2 and β20 = 1.881 N/mm2 so that hence k2 = 84 N/cm2.

β

By [44], the peak wall stress measured for normal people and AAA patients were 30-44 N/cm2 and 28-37 N/cm2, respectively. We assume that ECM in AAA patients is reduced to 75% of its value in healthy individuals. Hence which implies, by Eq (18), that β = 18 N/cm2.

1.1 σ

We assume σ = −p I, where p is average central aortic systolic pressure (100-120 mm Hg) [45]; we take p = 110 mm Hg.

Numerical methods

Weak formulation

In this subsection, we rewrite in compact form the coupled biological and mechanical systems. We denote by Z the vector of the densities or concentrations of cells and cytokines, i.e. Z = (P, M, T, S, I6, I10, I12, T, I, Q, Q, ρ). The entries of Z satisfy advection-diffusion equations with nonlinear source term, as seen from Eqs (1)–(12): where the velocity v is defined by v = D u, where D denotes the material derivative, In the equation for Z, F(Z) is a source term and D is a diagonal matrix with the diagonal terms for the diffusion coefficients of the biological species.

Some of the components of Z live only in either Ω(t) or Ω(t). Since there is a jump for Z at the media/adventitia interface, Γ(t), we split Z into two parts The boundary conditions can then be written in the form for an appropriate coefficient , some of its component being zero.

The governing equation for the displacement u can be written as where the stress σ is is defined in terms of the displacement by with the boundary conditions

Multiplying the coupled System (22) and (23) by arbitrary test functions (W, W,w), and performing integration by parts using the boundary conditions, we get the weak formulation: Here we used (⋅, ⋅) and 〈⋅, ⋅〉∂ to denote the L2 inner product on the domain D and on the boundary ∂D, respectively.

Lagrangian description for hyperelastic equation

Ω(t), Ω(t) and their boundaries are moving in time. We set Ω(t) = Ω(t) ∪ Ω(t) ∪ Γ(t) and , , . To characterize the motion of the media and adventitia, we introduce a flow map for the position of the particle at time t.

Then each function can be described in Lagrangian coordinates, namely for each function f(x, t), its Lagrangian description is We also introduce the following variables in Lagrangian coordinates: The displacement , the deformation tensor , and its determinant .

By change of variables, we can write the weak form for the equilibrium equation for w in Lagrangian coordinates, where

As for the surface tension term, we can also write it in the Lagrangian coordinates as follows:

By substituting this expression into the hyperelasitcity equation in Lagrangian coordinates Eq (25), we can get a nonlinear equation:

Finite element discretization on the moving mesh

In order to deal with the free boundary and the jump at the interface, we use the moving mesh method.

The initial domains and are discretized by triangulations, denoted by and , respectively. We consider the case that and are matching on the interface . We use piecewise linear and continuous function to approximate the displacement. Then the triangulation for the current domain can be obtained by moving the grid nodes of according to the flow map .

We next introduce the following time discretization for the material derivative and for where we use the superscript n to indicate the function at time nk, where k is the time step. Note that no interpolation is needed for evaluating the material derivative D at grid points, since the grid points are moving according to the flow map.

For the space discretization, we use piecewise linear and continuous function spaces and to approximate the solutions. Let , and be the bases of and , respectively, where N, N and N are the numbers of basis functions of the corresponding spaces. The approximation solutions at time nk can be written as summations with coefficients to be determined by the discrete system, corresponding to the System (24), for any and . Here we use the explicit discretization of the coefficients for the elasticity, and the source terms for the biological species. We denote the first two equations in Eulerian coordinates of the System (28) by

By this discretization, the elastic equation and the biological equation are decoupled. Given the solutions in last time step, we use the Newton’s nonlinear iteration method to solve the system a = 0 in the initial domain. Then we update the triangulations for the current domain by the approximate flow map . Following the standard technique of finite element method, we then solve the variational problem a = 0 for the biological species. We thus get the solution of current time step.