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, a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 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][8][9][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 R 0 , 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][8][9][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][21][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].   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 where the expression on the left-hand side includes advection and diffusion, and F X 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.

Equation for MCP-1 (P)
The equation for MCP-1 is given in the whole domain (media O M , and adventitia O A ) 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 first two terms on right-hand side account for the recruitment of macrophages by MCP-1 [19] and IL-6 [20][21][22]. The third term accounts for the activation of macrophages by IFN-γ [23], which is enhanced by TNF-α [31].

The density of T cells in
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 O M 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 Ω M [ Ω A )
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 @I 10 @t À D I 10 r 2 I 10 ¼ l I 10 M M |fflffl ffl{zfflffl ffl} production À d I 10 I 10 |fflfflffl ffl{zfflfflffl ffl} IL-12 is produced by macrophages [26], a process inhibited by IL-10. The equation for I 12 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 Qρ represents ECM degradation by MMP. Here we used the notations:

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 = y 1 and y = y 2 , and impose periodic boundary condition at y = y 1 and y = y 2 . 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, @S @n ¼ 0, at this membrane. All the other variables satisfy the following flux conditions at the interface Γ M : where X M and X A are the concentrations of cytokines or cell densities in media and adventitia, respectively. We take γ = 50 for cells, and γ = 500 for the cytokines [38].
where a M ðPÞ ¼ã M P PþK P . Both MCP-1 and IL-6 are produced by endothelial cells [17], which lie near the inner boundary of the media, Γ B . Hence @P @n þ a P ðP À P 0 Þ ¼ 0; where P 0 and I 60 are the concentrations of MCP-1 and I 6 in the blood. We assume no-flux boundary conditions on Γ B for all remaining variables. We also assume no-flux conditions on Γ A for all variables (except M and T, as stated above)

Initial conditions
We take initially I 6 ¼ I 60 2 , where I 60 = 6 × 10 −9 gcm −3 in the source of influx of I 6 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
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 B is the first invariant of the Left Cauchy-Green tensor B (namely I B = tr(B)), B = F F T , while F is the deformation gradient tensor. In terms of the displacement vector u = (u i ), we can write 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 , p Ã . For boundary conditions, we take where σ B is the stress tensor from the blood and κ is the mean curvature. 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).

Results
If at the initial screening for AAA we measure both the aorta diameter, R 0 , and the concentration of IL-6 in the blood, I 60 , then by simulating the model with these values of R(0) = R 0 and I 60 (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 I 60 (on the horizontal axis) and R 0 (on the vertical axis) the color on the column shows the diameter corresponding to the point (I 60 , R 0 ) at day T = 300 (upper) and T = 500 (lower).

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 (R 0 ) of a patient's  Tables 2 and 3 with I 60 = 6 × 10 −9 g/ml.  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?
The present paper addresses the first question with a mathematical model: Given R 0 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 R 0 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.

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

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, I 6 , I 10 , I 12 , T α , I γ , Q, Q r , ρ) T . 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 t u, where D t denotes the material derivative, In the equation for Z, F(Z) is a source term and D Z 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 O M (t) or O A (t). Since there is a jump for Z at the media/adventitia interface, Γ M (t), we split Z into two parts The boundary conditions can then be written in the form for an appropriate coefficientg, 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 sðuÞ ¼ À pðrÞI þ ð2b 1 ðSÞ þ 4b 2 ðSÞðtrðBÞ À 3ÞÞB: with the boundary conditions sn ¼ s B n À gkn; on G B ðtÞ; Multiplying the coupled System (22) and (23) by arbitrary test functions (W M , W A ,w), and performing integration by parts using the boundary conditions, we get the weak formulation: ðs; rwÞ OðtÞ À hs B n À gkn; wi G B ðtÞ þ hgkn; wi G M ðtÞ ¼ 0: Here we used (Á, Á) D and hÁ, Ái @D to denote the L 2 inner product on the domain D and on the boundary @D, respectively.

Lagrangian description for hyperelastic equation
To characterize the motion of the media and adventitia, we introduce a flow map xðx; tÞ for the position of the particlex at time t.
By change of variables, we can write the weak form for the equilibrium equation for w in Lagrangian coordinates, ðJŝF À T ;rŵÞÔ À hJŝ B F À Tn ;ŵi G B ð0Þ ¼ hgkn; wi G M ðtÞ[G B ðtÞ ð25Þ whereŝ ¼ À JpF À T þ Jðb 1 ðŜÞ þ 2b 2 ðŜÞðtrðBÞ À 3ÞÞF: As for the surface tension term, we can also write it in the Lagrangian coordinates as follows: hgkn; wi GðtÞ ¼ ZĜ gJjF À Tn j trðrwF À 1 Þ Àn T F À 1r wF À 1 F À Tn jF À Tn j 2 ! dŜ: By substituting this expression into the hyperelasitcity equation in Lagrangian coordinates Eq (25), we can get a nonlinear equation: a L ðû;ŵ;Ŝ;rÞ ¼ 0 ð27Þ 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.