Abstract
Vascular calcification results in stiffening of the aorta and is associated with hypertension and atherosclerosis. Atherogenesis is a complex, multifactorial, and systemic process; the result of a number of factors, each operating simultaneously at several spatial and temporal scales. The ability to predict sites of atherogenesis would be of great use to clinicians in order to improve diagnostic and treatment planning. In this paper, we present a mathematical model as a tool to understand why atherosclerotic plaque and calcifications occur in specific locations. This model is then used to analyze vascular calcification and atherosclerotic areas in an aortic dissection patient using a mechanistic, multi-scale modeling approach, coupling patient-specific, fluid-structure interaction simulations with a model of endothelial mechanotransduction. A number of hemodynamic factors based on state-of-the-art literature are used as inputs to the endothelial permeability model, in order to investigate plaque and calcification distributions, which are compared with clinical imaging data. A significantly improved correlation between elevated hydraulic conductivity or volume flux and the presence of calcification and plaques was achieved by using a shear index comprising both mean and oscillatory shear components (HOLMES) and a non-Newtonian viscosity model as inputs, as compared to widely used hemodynamic indicators. The proposed approach shows promise as a predictive tool. The improvements obtained using the combined biomechanical/biochemical modeling approach highlight the benefits of mechanistic modeling as a powerful tool to understand complex phenomena and provides insight into the relative importance of key hemodynamic parameters.
Introduction
Atherogenesis is a complex, multifactorial and systemic process; the result of a number of factors, each operating simultaneously at several spatial and temporal scales. The bewildering molecular and cellular complexity is well-described in Lusis' classical review more than a decade ago (Lusis, ), which highlights a plethora of biological mechanisms and gene associations, revealing an incredible etiological complexity. In addition to the biological components of the disease, atherosclerosis is also known to be related to mechanical stimuli on the vessel wall and hemodynamic parameters (Suo et al., ). Experimental evidence indicates that hemodynamic stimuli influence mechanotransduction and affect permeability (Davies, ). Increased permeability can lead to penetration and accumulation of lipoproteins (e.g., Low Density Lipoproteins–LDL) in the arterial wall and thus initiation of atherosclerosis. In recent years, much research has been conducted in order to draw correlations between hemodynamics and the atherogenic process (Peiffer et al., ; Alimohammadi et al., ). Certain hemodynamic parameters have been identified as key; these include flow distribution, pressure and wall shear stress (WSS) indices. Nevertheless, given the incredible complexity of the atherogenesis process, these hemodynamic analyses, on their own, have been inconclusive (Peiffer et al., ) and a clear metric for plaque location remains elusive. Strong correlations between atherosclerotic disease and vascular calcification have been well-documented in the literature, including large cohort studies (Sangiorgi et al., ). Although the underlying molecular cause of calcification is unknown (Lanzer et al., ), the severity and extent of mineralization in calcification reflect atherosclerotic plaque burden (Demer and Tintut, ). Given the central role of inflammation in atherogenesis, an interesting possibility is that vascular mineral itself may initiate, promote, or perpetuate atherosclerosis by inducing inflammatory cytokines in monocytes that encounter and ingest hydroxyapatite crystals (Nadra, ).
Modeling and simulation have been used in a large number of studies in order to improve understanding of the role of hemodynamic variables in plaque formation. Although this work is mathematically elegant and can provide detailed insight into hemodynamics, it is disconnected from molecular research. Similarly, biomedical researchers often reduce the complexity of investigations of cardiovascular disease into manageable parts, for example, working on cell-lines or employing large-scale genome wide association studies (GWAS) to identify SNPs related to CVD (Tegner et al., ). However, such statistical genetic models have no mechanistic basis, and it is significant that Lusis highlighted the revival of functional studies in a relatively recent review (Lusis, ).
In this paper, a clear application of mathematics for healthcare will be made by unifying multi-mechanistic factors in the prediction of atherosclerosis location. The location of atherosclerotic plaque and vascular calcification will be investigated using a patient-specific biomechanical model of an aortic dissection (AD). From a physiological point of view, AD is a life threatening condition in which a tear forms in the wall of the aortic wall and blood splits the media layer, forming two lumina: the true lumen (TL) and false lumen (FL; Braverman, ). Blood flows from the TL into the FL via a primary tear and, in communicating dissections, returns to the TL via one or multiple tears downstream. The section of intima and media that separates the two lumina, called the intimal flap (IF), often stiffens over time due to fibrosis (Criado, ). Common comorbidities in patients with AD are atherosclerosis (Coady et al., ; Tsai et al., ) and inner wall calcification (de Jong et al., ). Additionally, patients suffering from AD usually have an elevated pulse pressure, which would likely be further increased in the presence of calcification, due to the reduction of vessel elastance (Demer and Tintut, ). The patient data used for this research showed an AD with significant atherosclerosis and calcification regions, consistent with the condition. Detailed hemodynamic characterization was achieved by using patient-specific dynamic boundary conditions representing the downstream vasculature, based on in vivo measurements collected for the same patient, treated in University College Hospital (UCH; Alimohammadi et al., ). Although complex flow simulations for this patient have been published previously, the work presented here, including the quantification and analysis of plaque and calcification areas as well as the multi-scale framework used in this context, is completely new.
The postulate of this research is that multiscale modeling and simulation can pave the way to study multi-mechanistic factors to explain disease in a cohesive modeling framework, which can integrate key markers at different biological scales and can provide insight into endothelial mechanotransduction, as well as potential predictive power in patient-specific analyses, compared to purely hemodynamic, biomechanical or biochemical approaches. We will use a virtual “follow-up” approach, combining a fluid-structure interaction (FSI) simulation model of a dissected aorta with a model of plaque formation, using a number of somewhat disparate indicators available in the literature, which will be described below. The patient-specific simulation results are then coupled to an endothelial permeability model following the three-pore approach (Olgac et al., ), also recently used by Kim and Giddens ().
This paper is organized as follows: Section Methods presents the methods and simulation details. Results are presented in Section Results, including a clear, interpretable metric for atherogenic potential, which shows a better performance when compared to others based on “established” descriptors in the current literature. Simulation results of the model will be presented in and compared to atherosclerotic plaques and calcifications indicated in the original CT scans. The discussion, limitations and conclusions of this work will be presented in Section Discussion.
Methods
Details of the simulation
The specifics of the FSI simulation used for the present study are presented briefly in this section and are described in detail in Alimohammadi et al. (). The patient gave informed oral consent and ethical approval was given by the National Research Ethics Service, UK, REC reference: 13/EM/0143. The fluid domain was reconstructed from CT scans of a female patient, starting at the ascending arch and extending down to the thoracic aorta, upstream of the iliac bifurcation. The supra-aortic branches were included in the model, but the visceral branches were not clearly resolvable from the CT data and so were omitted. The vessel wall was modeled by extruding the outer wall of the fluid geometry uniformly by 2.5 mm, based on reports of thickened aortic walls in hypertension (Malayeri et al., ), which occurs in ~70% of AD patients (Hagan et al., ; Khan and Nair, ). The intimal flap (IF), separating the TL and FL, was created by filling the gap between the two lumina and was 2.45 ± 0.34 mm (median ± median absolute deviation). Fluid and solid meshes were generated in ANSYS ICEM-CFD, and contained ~230,000 (with 7 prismatic layers at the wall) and 50,000 elements respectively.
Simulations were carried out using ANSYS mechanical and CFX. The vessel wall was modeled using the isotropic hyperelastic model of Raghavan and Vorp (). This model is comparable to a linear elastic model with a Young's modulus of 1 MPa, but displays a small amount of strain-stiffening. An external pressure of 52.5 mmHg (diastolic pressure in the descending aorta) was applied and at each of the outlets of the 3D domain, the solid geometry was restricted to planar motion about the center-point of the lumen, enabling expansion of the vessels.
At the inlet, a flow wave from another study of a patient with type-B aortic dissection was used (Karmonik et al., ), as such data was not available for the present patient. At the fluid outlets, three element Windkessel models were used, with parameters tuned to invasive pressure measurements on the patient, using an iterative technique described in a previous study (Alimohammadi et al., ). Blood was modeled as an incompressible fluid of density 1056 kg/m3. Turbulence was modeled using the hybrid k-ε, k-ω shear stress transport (SST) turbulence model, with a 1% turbulence level at the inlet. In order to account for the shear-thinning properties of blood, the Carreau-Yasuda (CY) viscosity model was utilized with the parameters reported by Gijsen et al. (). The CY viscosity model includes infinite and zero shear viscosities, a characteristic time constant and two exponents without direct physical meaning to describe the viscosity-shear response of blood. The parameters of Gijsen et al. () were fitted to a blood analog, which showed similar characteristics to blood samples. The viscosity predicted by this model is ~3.5 mPa s (a commonly used Newtonian viscosity) at 100/s, and decreases for higher shear rates to around 2.5 mPa s at 1000/s. At shear rates below 100/s, the viscosity increases, exceeding 10 mPa s for shear rates below 2/s.
Simulations were run with a time step of 5 ms, a periodic solution was achieved after two cycles and the third cycle was extracted for further analysis. In a previous study (Alimohammadi et al., ), it was shown that the low shear rates in the slow flow regions of the distal and proximal FL led to significantly increased viscosity.
Shear stress indices (SSI)
In the permeability calculations described in Section Endothelial Permeability Model, one of the inputs is an index indicative of the shear stress condition that leads to increased permeability. The time-average wall shear stress (TAWSS) and oscillatory shear index (OSI; Ku et al., ) are two commonly used indices that are considered important for plaque formation. TAWSS describes the average magnitude of the shear stress and the OSI gives an indication of the directionality of the shear stress, yielding 0 for uniaxial flows and 0.5 when there is no preferential direction. Typically, in permeability models for atherogenesis, TAWSS is used as the shear index related to permeability.
A number of studies (Malek et al., ; Xiang et al., ; Chiu and Chien, ; Meng et al., ) have shown that regions with low average shear stress combined with highly oscillatory shear stress have increased endothelial permeability along with other pathological responses. More recently, Sáez et al. () presented the 3D remodeling of endothelial cells as the combined effect of OSI and TAWSS in a computational framework, which fitted experimental works presented before in in vitro studies.
In light of the indications for increased permeability in low, oscillatory regions, we propose an index, HOLMES (Highly Oscillatory, Low MagnitudE Shear), given by:
This parameter is equivalent to half the reciprocal of relative residence time (RRT) which was previously identified as a potential index for combining these two characteristics (Himburg, ). The HOLMES indicator can be understood as a modified TAWSS, with the (0.5 − OSI) term further reducing the index in regions where the wall shear stress is both low in magnitude and oscillatory in nature. Additionally, HOLMES provides a conceptually alternative explanation, offering a linear (rather than reciprocal) index, proportional to WSS that intuitively corresponds to the observed effects of shear characteristics on endothelial permeability.
In the present study, we compare the efficacy of using the most-widely used hemodynamic/shear stress index (SSI) for atheroprone regions, i.e., TAWSS, with HOLMES, separate and in a mechanistic model, as shown in the sections below.
Endothelial permeability model
The endothelial permeability model proposed here is based on previous work describing the early stages of atherosclerosis, using a transport model of low density lipoprotein (LDL) from the artery lumen into the arterial wall, taking into account the effects of mechanical stimuli exerted by the blood flow on the endothelial cell layer and its pathways of volume and solute flux; see Díaz-Zuccarini et al. () for more details. An excellent and recent analysis and use of this model along the same lines using time-average wall shear stress (TAWSS) has been recently published by Kim and Giddens ().
The endothelial layer is described with a three-pore modeling approach considering the contributions of the vesicular pathway, normal junctions, and leaky junctions. The fraction of leaky junctions is calculated as a function of the mechanical stimuli and is used in conjunction with the pore theory to determine the transport properties of this pathway.
The LDL transport equations are decomposed using three main penetration pathways: leaky junctions, normal junctions and vesicular pathways; so the bulk of volume flux (Jv) through the endothelial membrane is given by:
where Jv, lj is the flux through leaky junctions and Jv, nj is the flux through normal junctions.
The volumetric flux through leaky junctions (Jv, lj) is calculated using a modified version of the Kedem-Ketchalsky equations for membrane transport:
where Lp, lj is the hydraulic conductivity, Δpend is the pressure difference through the endothelium, σd is the osmotic reflection coefficient and ΔΠ is the osmotic pressure. The value of Δpend is estimated by subtracting the externally applied pressure from either the average pressure throughout the domain (uniform pressure gradient) or the spatially varying time average pressure at each wall location (time-averaged pressure gradient).
According to the three pores theory, solute flux (LDL flux in this case) only occurs through endothelial leaky cell junctions and vesicles:
Assuming that the solute flux through the vesicular pathway (Js, v) is 10% of the solute flux through the leaky junction pathway (Js, lj) (Olgac et al., ).
Finally, in the interest of brevity, Equation (5) shows a general form of the function used to calculate Js, lj which is proportional to the magnitude of Jv, lj (more details about the calculation of Js, lj are shown elsewhere (Díaz-Zuccarini et al., ; Alimohammadi et al., ):
where Pi is the diffusive permeability, Pe the modified Peclet number, clum and cw, end the LDL concentration in the lumen and the sub-endothelial layer respectively, cave is the mean endothelial concentration and σ the solvent drag coefficient.
Endothelial cell shape will affect the amount of leaky junctions. Experimental findings have shown that in areas of altered hemodynamics, endothelial cells do not have a typical cobblestone shape, but rather exhibit a more circular shape as well as increased permeability. Previous models used a relationship between endothelial permeability and local WSS based on the Endothelial Cell Shape Index (ECSI; Levesque et al., ). ECSI is related to the cellular shape and takes values from zero to one, i.e., a circle has an ECSI of one whilst a straight line has an ECSI of zero.
Some approaches based on ECSI calculated this variable as a function of WSS (for instance Olgac et al., used WSS values in steady state simulations). However, as previously mentioned, recent work shows that ECSI can be affected by other indices as well, such as OSI. Based on the seminal paper by Levesque et al. (), excellent work from Sáez et al. (), presents how different OSI and TAWSS modify the endothelial cell shape. This combined relationship will be key for the modeling work presented here.
In the present study, we consider the role of the two different shear stress indices (SSI), i.e., TAWSS and HOLMES, described in Section Shear Stress Indices (SSI).
For each SSI, ECSI is defined according to
Figure 1 shows a comparison of normalized ECSI using Saez et al. results (Sáez et al., ), which are dependent on TAWSS and OSI, with the normalized ECSI using Equation (6) and HOLMES. Very similar behavior can be seen in both surfaces, showing how areas of low TAWSS and high OSI will increase the values of ECSI. Following this, we will use HOLMES as a combined index able to capture key haemodynamic features in the regions of interest for atherosclerotic plaque/calcification locations.
Figure 1
Leaky cells have high permeability to LDL, which can be linked to the magnitude of the SSI. Areas with high ECSI will be related to a higher number of mitotic cells (MC) which are calculated as follows:
Assuming that within the endothelium the quantity of leaky mitotic cells is ~80.5%, which represents ~45.3% of the total number of leaky cells (LC) in that area (Tedgui and Lever,
The ratio of endothelium (ϕ) covered by LCs is calculated using
where Rcell is the radius of a single cell. Finally, the total hydraulic conductivity of the endothelial leaky junctions (Lp, lj) is defined as:
where Lp, slj is the hydraulic conductivity of a single leaky junction calculated as follows:
with w and llj being the half-width (20 nm) and the length (2 μm) of the leaky junctions and μ the viscosity term used for the estimation of the LDL penetration.
Estimating plaque location and metrics evaluation
Previous work (Alimohammadi et al.,
Figure 2

Definition of the different regions in the aorta for comparison and analysis.
Comparison of the calcification/plaque area observed in the CT scans and those calculated by the chosen metric following the methodology shown here, was used to evaluate the efficacy of the different metrics and the role of some of the selected hemodynamic (viscosity and SSI) variables.
Results
Figure 3A shows an image of the 3D CT image reconstruction in Aquarius (TeraRecon, USA), in which the highest intensity regions indicate plaques (white) or vessel wall calcification. Figure 3B shows an image of the CT scans reconstructed in the LUT volume viewer in Fiji (Schindelin et al.,
Figure 3

The right anterior and left posterior views of the aorta. (A) 3D CT data visualized in Aquarius and (B) 3D LUT volume viewer in FIJI (Image J).
As seen in both panels, a number of completely formed plaques can be observed in the lower edge of the aortic arch, between the two tears in FL and around the iliac bifurcation. Calcifications can be seen along the aortic arch proximal to the supraaortic branches, and in the proximal and distal FL. The existence and locations of the plaques and calcifications were confirmed and approved by the vascular surgeon who managed this patient.
Wall shear stress indices
Figure 4 shows the distributions of the three SSI used in the present study. The TAWSS (Figure 4A) is moderate in the ascending aorta and in the FL between the two tears. Elevated regions of TAWSS can be observed in the visceral branches, at the coarctation and in the distal TL. The OSI (Figure 4B) shows scattered elevated regions throughout the domain, with particularly high values in the distal and proximal FL (which would not be captured using a rigid wall model; Alimohammadi et al.,
Figure 4

Wall shear indices from the FSI simulations. (A) TAWSS, (B) OSI, (C) HOLMES.
In the following sub-sections, we systematically compare the indices and parameters used in the literature to indicate or predict locations of calcification and/or atherogenesis with our integrated model, the use of HOLMES and spatially varying viscosity.
Analysis of the total hydraulic conductivity
As the volume flux is dependent on both hydraulic conductivity and pressure gradient across the endothelium, we first consider Lp alone.
Figures 4A,B show the total hydraulic conductivity of leaky junctions (Lp, lj) values calculated assuming a constant viscosity at the wall, equal to the plasma viscosity and TAWSS and HOLMES, respectively, as wall shear stress indices. Figures 4C,5D show the same wall shear stress indices, but with the viscosity at the wall predicted based on the continuum assumption using the Carreau-Yasuda viscosity model. Note that all panels use the same color scale. Tedgui and Lever (
Figure 5

Total hydraulic conductivity of leaky junctions (Lp, lj) values calculated using various parameter combinations from the FSI results shown in the right anterior view. (A) Constant plasma viscosity and TAWSS, (B) Constant plasma viscosity and HOLMES, (C) Non-Newtonian blood viscosity and TAWSS, (D) Non-Newtonian blood viscosity and HOLMES.
Under the assumption that plaque formation and calcification occur as a result of elevated Lp, lj, the correlation between Figures 6, 3B implies that this combination of variables is able to identify regions of calcification and plaque formation, and thus could be potentially used to predict further development of such pathologies.
Figure 6

Right anterior and left posterior views of Lp, lj, calculated using non-Newtonian viscosity and HOLMES.
The region of calcification with scattered small plaques observed in the FL between the two tears is captured in the Lp, lj prediction in Figure 6. Similarly, the calcification of the proximal FL (both views) and the distal FL (left posterior view), are predicted by elevated Lp, lj in these regions. The low Lp, lj in the distal TL and ascending arch in the left posterior view also correlate with minimal calcification in Figure 3B. The plaques at the lower edge of the aortic arch, best observed in Figure 3A in the left posterior view, correspond to regions of elevated Lp, lj in Figure 6.
Some regions of the image do not correlate, such as the low Lp, lj in the right anterior view of the ascending arch, where there is calcification in the CT image.
Analysis of the volume flux
Figure 6 provides images of the predicted volume flux across the vessel wall, Jv, lj, using various combinations of wall shear indices and viscosity assumptions, along with either a uniform pressure gradient or a time-averaged pressure gradient. These figures can be compared with Figure 3B to analyze the hypothesis that calcification and plaque formation are better correlated with volume flux, Jv, lj, rather than Lp, lj alone. Meyer et al. (
None of the constant plasma models (Figures 7A,B,E,F produced Jv, lj distributions that compared well with Figure 3B. Figures 7C,D using the uniform pressure gradient are qualitatively identical to Figures 4C,D as the uniform pressure distribution does not affect the relative values of Jv, lj as compared to Lp, lj. When using the time averaged pressure gradient (Figures 7G,H), the distributions do not significantly change in qualitative terms, indicating that the prediction is more sensitive to Lp, lj, and thus appropriate models for wall shear stress and blood viscosity at the wall, than to the pressure gradient. Figure 8, showing the distribution of Jv, lj using HOLMES, the non-Newtonian model and time averaged pressure distribution, is therefore qualitatively similar to Figure 6, although there are some differences in the magnitudes.
Figure 7

Volume flux (Jv, lj) calculated using various parameters. Wall shear indices (A) TAWSS, (B) HOLMES with constant plasma viscosity and uniform pressure gradient. Wall shear indices (C) TAWSS, (D) HOLMES with non-Newtonian blood viscosity and uniform pressure gradient. (E) TAWSS, (F) HOLMES with constant plasma viscosity and time averaged pressure gradient. (G) TAWSS, (H) HOLMES with non-Newtonian blood viscosity, and time-averaged pressure gradient.
Figure 8

Right anterior and left posterior views of Jv, lj, calculated using non-Newtonian viscosity, HOLMES and time-averaged pressure gradient.
Discussion
The results shown in Section Results are promising and thus warrant further investigation. They also open up a number of questions. Is the volume flux (which determines the magnitude of the solute flux) a metric for plaque location? If so, then what is required in order to estimate it correctly? As mentioned in the introduction, the question of plaque location remains elusive. There are presently no available tools to predict where this might occur, or indeed to sufficiently explain why. This may be because, if we take volume flux (Jv, lj) (or even hydraulic conductivity) as the metric for plaque location, its generation is multifactorial and complex and thus simplified models lack some of the key features that seem to impact plaque location (see Section Results). In the present study, strong correlations were only produced through the use of a multiscale model combining multiple hypotheses at the cellular and physiological level, as evidenced by the comparison of Figures 5, 7 with Figure 3. This complex approach involves the detailed formulation, extraction and combination of more than 11 different variables including the calculation of variable pressure differences along the arterial wall across the endothelium, the use of a non-linear model for blood viscosity and use of the HOLMES index to isolate regions of low, oscillatory shear. This formulation, making use of mechanistic models, also provides these metrics (Jv, lj and Lp, lj) with strong interpretability and physical meaning. As previously explained, the origin and interpretation of HOLMES as a combined index is strongly rooted in the most recent, relevant literature (Sáez et al.,
In addition to the clear qualitative improvement in the correlation between the CT scans and Lp, lj and Jv, lj distributions, we utilized a simple measure based on plaque areas, as shown in Figure 9 for the descending aorta, a key region in the case of AD (as explained in Section Estimating Plaque Location and Metrics Evaluation).
Figure 9

Example of atherosclerotic plaque identification from simulations using hydraulic conductivity and comparison with 3D CT scans. (A,B) Section of the artery showing plaques based on Lp, lj; (D,E) Plaque extracted from (A,B); (C) Section of the artery from CT scan; (F) Plaque extracted from (C).
When comparing Lp, lj values in specific regions of the aorta, using the results shown in Figure 5, the calcified regions predicted by the model overlapped 80% of the calcified regions observed in the CT scans for the HOLMES model. No calcification was predicted using TAWSS. The next step was to analyze the results of Jv, lj (given its stronger relationship with the LDL flux).
Figures 10A,C show Jv, lj in this region calculated with TAWSS and HOLMES, respectively and Figures 10B,D show the calcified region only. Figures 10E,F show the equivalent figures for the CT scans. Comparing the proportional area of the calcified region between the two cases, 26 and 93% of the regions overlap for TAWSS and HOLMES, respectively. When calculated using HOLMES in the selected region (Region 3 in Figure 2), 45% of the visible arterial wall is predicted to be prone to calcification, compared to 48% of the visible arterial wall observed in the CT scan image. Applying the same analysis using HOLMES to the whole aorta predicts between 80 and 95% of overlapping. This analysis supports the potential predictive power of the proposed model.
Figure 10

Example of atherosclerotic plaque identification from simulations using solute flux and comparison with 3D CT scans. (A,B) Section of the artery showing plaques based on Jv, lj; (D,E) Plaque extracted from (A,B); (C) Section of the artery from CT scan; (F) Plaque extracted from (C).
From a hemodynamics point of view, the heterogeneous distribution of calcification/plaque formation will be influenced by disturbed flow. Wall shear indices such as TAWSS and OSI capture only partial aspects of the complex flow fields. The wall shear stress vector has both direction and magnitude at each moment throughout the cardiac cycle. TAWSS captures the average magnitude, but contains no information about directionality. OSI characterizes the variability in the directionality. However, a number of indices (“emerging” multidirectional predictors, as defined by Gallo et al.,
In this paper, we have exclusively focused on the “established” predictors, for which the links between endothelial cell behavior and haemodynamics have been shown in the relevant literature as described above.
In the work presented here, the use of patient-specific data (including invasive haemodynamic measurements) and BCs, inclusion of wall motion and use of a non-Newtonian blood viscosity model are likely to improve the estimates of these hemodynamic variables and thereby increase the likelihood of an improved correlation. Importantly, as evidenced by Figures 4, 7, the interpretation of the role of these variables individually was unable to predict distributions of atherogenesis/calcification. It was only through a combination of HOLMES, a compound shear index, a non-Newtonian fluid viscosity and a sophisticated simulation model that this was achieved. In general, when extracting information related to values calculated at the wall, motion is extremely important. In this particular case and as shown previously (Alimohammadi et al.,
Limitations
One limitation of the fluid dynamics simulation is the absence of the visceral arteries, which would alter the flow in the descending aorta. As previously stated, it was not possible to resolve these vessels from the CT scans, so they were omitted. Nonetheless, even in their absence, we believe the usefulness of the model was demonstrated.
A patient-specific inflow waveform was not available for the present study, and as such, a waveform from the literature was selected from a patient suffering from a similar type-B AD. The absence of a patient-specific inflow reduces the specificity to the patient of the present results however, given that the boundary conditions were tuned to patient-specific values using the same inflow, the hemodynamic environment predicted by the simulations is expected to be an appropriate representation of a patient with type-B AD. The flow waveform at the inlet was applied as a uniform velocity, rather than mapping to parabolic or Womersley profiles. Whilst accurate axial velocity profile, as extracted using pcMRI can provide improved predictions of characteristics such as helical flow (Morbiducci et al.,
The isotropic hyperelastic model of Raghavan and Vorp (
Another limitation of this approach, which is fairly consistent throughout the related literature is the use of “visual maps,” in order to establish qualitative correlations between the in vivo data and the simulation results. Recent studies address the issue of quantification of these type of results by using statistical methods to establish quantitative correlations and statistical significance of the variables analyzed (simulated) with respect to plaque location, obtained from in vivo data (Morbiducci et al.,
The results presented in this paper when using the multi-criteria, patient-specific, multi-scale complex framework described here, compare well, qualitatively, to in vivo data and although there is no perfect match, this is to be expected, when considering the number of assumptions, simplifications and limitations as described (please see above). Moreover, as for patient-specific data, only anatomical information was fed into the model and hemodynamic variables were calculated based on invasive pressure measurements in order to properly characterize the flow, as described in Alimohammadi et al. (
Conclusions
In the present study, we have presented a mechanistic, mathematical model of endothelial mechanotransduction to understand plaque location/calcification. The model is tested on a patient-specific case for which in vivo measurements were obtained at University College Hospital and a patient-specific biomechanical model produced (Alimohammadi et al.,
As an input to the endothelial permeability model relating to shear stress, we hypothesized that regions of oscillatory, low magnitude shear stress would be susceptible to calcification, due to the known connection to increased permeability. We therefore, proposed the compound HOLMES shear index, which includes both magnitude and oscillatory characteristics and will thus emphasize oscillatory, low magnitude shear stress and found that it considerably improved the predictive power of the model over TAWSS-based analysis.
The role of the complex characteristics of the hemodynamics near the wall was investigated by hypothesizing that elevated blood viscosity near the wall would limit the convection of plasma into the vessel wall, and vice versa. It was found that when the spatially varying blood viscosity at the wall as estimated using an empirical non-Newtonian viscosity model, was used in the endothelial permeability model, the predicted Lp, lj and Jv, lj distributions bore increased resemblance to the observed regions of calcification.
This work provides a good example of the use of multiscale mathematical modeling to understand physiology. The promising results obtained from this approach warrant further investigation. Next steps will include studies of a larger number of patients to enable comparisons amongst patients and potentially statistical analyses, in order to investigate in detail the predictive power of the model.
Funding
EPSRC grant “Personalised Medicine Through Learning in the Model Space” (grant number EP/L000296/1). Leverhulme Trust Senior Research Fellowship “Exploring the Unknowable Using Simulation: Structural Uncertainty in Multiscale Models” (Fellowship number RF-446 2015-482).
Conflict of interest statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Statements
Author contributions
MA and CP conceived the study and carried out simulations and wrote the manuscript. OA acquired the clinical data and VD conceived the study and wrote the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
References
1
AlimohammadiM.AguO.BalabaniS.Díaz-ZuccariniV. (2014). Development of a patient-specific simulation tool to analyse aortic dissections: assessment of mixed patient-specific flow and pressure boundary conditions. Med. Eng. Phys.36, 275–284. 10.1016/j.medengphy.2013.11.003
2
AlimohammadiM.Pichardo-AlmarzaC.Giulia DiT.StavroulaB.AguO.Diaz-ZuccariniV. (2015a). Predicting atherosclerotic plaque locationin an iliac bifurcation using a hybrid CFD/biomechanical approach, in Part II, Bioinformatics and Biomedical Engineering: Volume 9044 of the Series Lecture Notes in Computer Science, eds OrtuñoF.RojasI. (Springer Link), 594–606.
3
AlimohammadiM.SherwoodJ. M.KarimpourM.AguO.BalabaniS.Díaz-ZuccariniV. (2015b). Aortic dissection simulation models for clinical support: fluid-structure interaction vs. rigid wall models. Biomed. Eng. Online14, 1661–1616. 10.1186/s12938-015-0032-6
4
ArzaniA.ShaddenS. C. (2016). Characterizations and correlations of wall shear stress in aneurysmal flow. J. Biomech. Eng.138, 014503–014510. 10.1115/1.4032056
5
BravermanA. C. (2010). Acute aortic dissection: clinician update. Circulation122, 184–188. 10.1161/CIRCULATIONAHA.110.958975
6
BrownA. G.ShiY.MarzoA.StaicuC.ValverdeI.BeerbaumP.et al. (2012). Accuracy vs. computational time: translating aortic simulations to the clinic. J. Biomech.45, 516–523. 10.1016/j.jbiomech.2011.11.041
7
CampbellI. C.RiesJ.DhawanS. S.QuyyumiA. A.TaylorW. R.OshinskiJ. N. (2012). Effect of inlet velocity profiles on patient-specific computational fluid dynamics simulations of the carotid bifurcation. J. Biomech. Eng.134, 051001–0510018. 10.1115/1.4006681
8
ChiuJ. J.ChienS. (2011). Effects of disturbed flow on vascular endothelium: pathophysiological basis and clinical perspectives. Physiol. Rev.91, 327–387. 10.1152/physrev.00047.2009
9
CoadyM. A.RizzoJ. A.ElefteriadesJ. A. (1999). Pathologic variants of thoracic aortic dissections. Penetrating atherosclerotic ulcers and intramural hematomas. Cardiol. Clin.17, 637–657. 10.1016/S0733-8651(05)70106-5
10
CriadoF. J. (2011). Aortic dissection: a 250-year perspective. Tex. Heart Inst. J.38, 694–700.
11
DaviesP. F. (1995). Flow-mediated endothelial mechanotransduction. Physiol. Rev. 75, 519–560.
12
de JongP. A.HellingsW. E.TakxR. A. P.IšgumI.van HerwaardenJ. A.MaliW. P. (2014). Computed tomography of aortic wall calcifications in aortic dissection patients. edited by Jeroen Hendrikse. PLoS ONE9:e102036. 10.1371/journal.pone.0102036
13
DemerL. L.TintutY.. (2008). Vascular calcification: pathobiology of a multifaceted disease. Circulation117, 2938–2948. 10.1161/CIRCULATIONAHA.107.743161
14
Díaz-ZuccariniV.Di TomasoG.AguO.Pichardo-AlmarzaC. (2014). Towards personalised management of atherosclerosis via computational models in vascular clinics: technology based on patient-specific simulation approach. Healthc. Technol. Lett.1, 13–18. 10.1049/htl.2013.0040
15
ErrenT. C. (2007). The case for a posteriori hypotheses to fuel scientific progress. Med. Hypotheses69, 448–453. 10.1016/j.mehy.2006.12.026
16
GalloD.SteinmanD. A.MorbiducciU. (2016). Insights into the co-localization of magnitude-based versus direction-based indicators of disturbed shear at the carotid bifurcation. J. Biomech.. [Epub ahead of print]. 10.1016/j.jbiomech.2016.02.010
17
GasserT. C.OgdenR. W.HolzapfelG. A. (2006). Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J. R. Soc. Interface3, 15–35. 10.1098/rsif.2005.0073
18
GijsenF. J.van de VosseF. N.JanssenJ. D. (1999). The influence of the non-newtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model. J. Biomech.32, 601–608. 10.1016/S0021-9290(99)00015-9
19
HaganP. G.NienaberC. A.IsselbacherE. M.BruckmanD.KaraviteD. J.RussmanP. L.et al. (2000). The International Registry of Acute Aortic Dissection (IRAD). JAMA283, 897–903. 10.1001/jama.283.7.897
20
HimburgH. A. (2004). Spatial comparison between wall shear stress measures and porcine arterial endothelial permeability. Am. J. Physiol. Heart Circ. Physiol.286, H1916–H1922. 10.1152/ajpheart.00897.2003
21
IsgumI.van GinnekenB.OlreeM. (2004). Automatic detection of calcifications in the aorta from CT scans of the abdomen. 3D computer-aided diagnosis. Acad. Radiol.11, 247–257. 10.1016/S1076-6332(03)00673-1
22
KarmonikC.BismuthJ. X.DaviesM. G.LumsdenA. B. (2008). Computational hemodynamics in the human aorta: a computational fluid dynamics study of three cases with patient-specific geometries and inflow rates. Technol. Health Care16, 343–354.
23
KhanI. A.NairC. K. (2002). Clinical, diagnostic, and management perspectives of aortic dissection. Chest122, 311–28.
24
KimS.GiddensD. P. (2015). Mass transport of low density lipoprotein in reconstructed hemodynamic environments of human carotid arteries: the role of volume and solute flux through the endothelium. J. Biomech. Eng.137:041007. 10.1115/1.4028969
25
KuD. N.GiddensD. P.ZarinsC. K.GlagovS. (1985). Pulsatile flow and atherosclerosis in the human carotid bifurcation. Positive correlation between plaque location and low oscillating shear stress. Arteriosclerosis5, 293–302. 10.1161/01.ATV.5.3.293
26
LanzerP.BoehmM.SorribasV.ThirietM.JanzenJ.ZellerT.et al. (2014). Medial vascular calcification revisited: review and perspectives. Eur. Heart J.35, 1515–1525. 10.1093/eurheartj/ehu163
27
LevesqueM. J.LiepschD.MoravecS.NeremR. M. (1986). Correlation of endothelial cell shape and wall shear stress in a stenosed dog aorta. Arteriosclerosis6, 220–229. 10.1161/01.ATV.6.2.220
28
LusisA. J. (2000). Atherosclerosis. Nature407, 233–241. 10.1038/35025203
29
LusisA. J. (2012). Life After GWAS: functional genomics in vascular biology. Arterioscler. Thromb. Vasc. Biol.32, 169–169. 10.1161/ATVBAHA.111.243543
30
MalayeriA. A.NatoriS.BahramiH.BertoniA. G.KronmalR.LimaJ. A.et al. (2008). Relation of aortic wall thickness and distensibility to cardiovascular risk factors (From the Multi-Ethnic Study of Atherosclerosis [MESA]). Am. J. Cardiol.102, 491–496. 10.1016/j.amjcard.2008.04.010
31
MalekA. M.AlperS. L.IzumoS. (1999). Hemodynamic shear stress and its role in atherosclerosis. JAMA282, 2035–2042.
32
MarzoA.SinghP.ReymondP.StergiopulosN.PatelU.HoseR. (2009). Influence of inlet boundary conditions on the local haemodynamics of intracranial aneurysms. Comput. Methods Biomech. Biomed. Engin.12, 431–444. 10.1080/10255840802654335
33
MengH.TutinoV. M.XiangJ.SiddiquiA. (2014). High WSS or Low WSS? Complex interactions of hemodynamics with intracranial aneurysm initiation, growth, and rupture: toward a unifying hypothesis. AJNR. Am. J. Neuroradiol.35, 1254–1262. 10.3174/ajnr.A3558
34
MeyerG.MervalR.TedguiA. (1996). Effects of pressure-induced stretch and convection on low-density lipoprotein and albumin uptake in the rabbit aortic wall. Circ. Res.79, 532–540. 10.1161/01.RES.79.3.532
35
MorbiducciU.GalloD.CristofanelliS.PonziniR.DeriuM. A.RizzoG.et al. (2015). A rational approach to defining principal axes of multidirectional wall shear stress in realistic vascular geometries, with application to the study of the influence of helical flow on wall shear stress directionality in aorta. J. Biomech.48, 899–906. 10.1016/j.jbiomech.2015.02.027
36
MorbiducciU.PonziniR.GalloD.BignardiC.RizzoG. (2013). Inflow boundary conditions for image-based computational hemodynamics impact of idealized versus measured velocity profiles in the human aorta. J. Biomech.46, 102–109. 10.1016/j.jbiomech.2012.10.012
37
NadraI. (2005). Proinflammatory activation of macrophages by basic calcium phosphate crystals via protein kinase C and MAP kinase pathways: a vicious cycle of inflammation and arterial calcification?Circ. Res.96, 1248–1256. 10.1161/01.RES.0000171451.88616.c2
38
OlgacU.KurtcuogluV.PoulikakosD. (2008a). Computational modeling of coupled blood-wall mass transport of LDL: effects of local wall shear stress. Am. J. Physiol. Heart Circul. Physiol.294, H909–H919. 10.1152/ajpheart.01082.2007
39
OlgacU.KurtcuogluV.PoulikakosD. (2008b). Computational modeling of coupled blood-wall mass transport of LDL: effects of local wall shear stress. Am. J. Physiol. Heart Circul. Physiol.294, H909–H919. 10.1152/ajpheart.01082.2007
40
PeifferV.SherwinS. J.WeinbergP. D. (2013a). Does low and oscillatory wall shear stress correlate spatially with early atherosclerosis? A systematic review. Cardiovasc. Res.99, 242–250. 10.1093/cvr/cvt044
41
PeifferV.SherwinS. J.WeinbergP. D. (2013b). Computation in the rabbit aorta of a new metric – the transverse wall shear stress – to quantify the multidirectional character of disturbed blood flow. J. Biomech.46, 2651–2658. 10.1016/j.jbiomech.2013.08.003
42
RaghavanM. L.VorpD. A. (2000). Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysm: identification of a finite strain constitutive model and evaluation of its applicability. J. Biomech.33, 475–482. 10.1016/S0021-9290(99)00201-8
43
ReymondP.CrosettoP.DeparisS.QuarteroniA.StergiopulosN. (2013). Physiological simulation of blood flow in the aorta: comparison of hemodynamic indices as predicted by 3-D FSI, 3-D rigid wall and 1-D models. Med. Eng. Phys.35, 784–791. 10.1016/j.medengphy.2012.08.009
44
RoyD.HolzapfelG. A.KauffmannC.SoulezG. (2014). Finite element analysis of abdominal aortic aneurysms: geometrical and structural reconstruction with application of an anisotropic material model. IMA J. Appl. Mathematics79, 1011–1026. 10.1093/imamat/hxu037
45
SáezP.MalvèM.MartínezM. A. (2015). A theoretical model of the endothelial cell morphology due to different waveforms. J. Theor. Biol.379, 16–23. 10.1016/j.jtbi.2015.04.038
46
SangiorgiG.RumbergerJ. A.SeversonA.EdwardsW. D.GregoireJ.FitzpatrickL. A.et al. (1998). Arterial calcification and not lumen stenosis is highly correlated with atherosclerotic plaque burden in humans: a histologic study of 723 coronary artery segments using nondecalcifying methodology. J. Am. Coll. Cardiol.31, 126–133. 10.1016/S0735-1097(97)00443-9
47
SchindelinJ.Arganda-CarrerasI.FriseE.KaynigV.LongairM.PietzschT.et al. (2012). Fiji: an open-source platform for biological-image analysis. Nat. Methods9, 676–682. 10.1038/nmeth.2019
48
SuoJ.FerraraD. E.SorescuD.GuldbergR. E.TaylorW. R.GiddensD. P. (2006). Hemodynamic shear stresses in mouse aortas: implications for atherogenesis. Arterioscler. Thromb. Vasc. Biol.27, 346–351. 10.1161/01.ATV.0000253492.45717.46
49
TedguiA.LeverM. J. (1984a). Filtration through damaged and undamaged rabbit thoracic aorta. Am. J. Physiol.247 (5 Pt 2), H784–H791.
50
TedguiA.LeverM. J. (1984b). Filtration through damaged and undamaged rabbit thoracic aorta. Am. J. Physiol. Heart Circ. Physiol.247, H784–91.
51
TegnerJ.SkogsbergJ.BjorkegrenJ. (2006). Thematic review series: systems biology approaches to metabolic and cardiovascular disorders. Multi-organ whole-genome measurements and reverse engineering to uncover gene networks underlying complex traits. J. Lipid Res.48, 267–277. 10.1194/jlr.R600030-JLR200
52
TsaiT. T.FattoriR.TrimarchiS.IsselbacherE.MyrmelT.EvangelistaA.et al. (2006). Long-term survival in patients presenting with type B acute aortic dissection: insights from the international registry of acute aortic dissection. Circulation114, 2226–2231. 10.1161/CIRCULATIONAHA.106.622340
53
XiangJ.NatarajanS. K.TremmelM.MaD.MoccoJ.HopkinsL. N.et al. (2010). Hemodynamic-morphologic discriminants for intracranial aneurysm rupture. Stroke42, 144–152. 10.1161/STROKEAHA.110.592923
Summary
Keywords
mathematical modeling, multiscale, atherosclerosis, patient-specific, aortic dissection, in vivo data
Citation
Alimohammadi M, Pichardo-Almarza C, Agu O and Díaz-Zuccarini V (2016) Development of a Patient-Specific Multi-Scale Model to Understand Atherosclerosis and Calcification Locations: Comparison with In vivo Data in an Aortic Dissection. Front. Physiol. 7:238. doi: 10.3389/fphys.2016.00238
Received
01 April 2016
Accepted
03 June 2016
Published
21 June 2016
Volume
7 - 2016
Edited by
Zbigniew R. Struzik, The University of Tokyo, Japan
Reviewed by
Umberto Morbiducci, Politecnico di Torino, Italy; Miguel Angel Martinez, University of Zaragoza, Spain
Updates

Check for updates
Copyright
© 2016 Alimohammadi, Pichardo-Almarza, Agu and Díaz-Zuccarini.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Mona Alimohammadi mona.alimohammadi.10@ucl.ac.uk
This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology
†These authors have contributed equally to this work.
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.