Keywords

1 Introduction

Peripheral Arterial Occlusive Diseases (PAODs) hold a high incidence worldwide with more than 200 million people affected annually [1]. Their etiology is mainly attributable to atherosclerosis, a chronic inflammation of the arterial wall that causes the narrowing of the lumen through the build-up of a fatty plaque and that preferentially develops in sites affected by a low or oscillatory Wall Shear Stress (WSS) profile [2].

The consequences of severe PAODs range from limb gangrene insurgence, which ultimately requires the amputation of the foot or leg, to the spreading of atherosclerosis to other parts of the body up to exposing the patient to high risk of heart attack and stroke [3].

Percutaneous Transluminal Angioplasty (PTA), with or without stent, and Vein Graft Bypass (VGB) are the preferred interventions aimed to restore the physiological circulation. These procedures suffer of a high rate of long-term failure, with a 3-years patency of 60% and a 5-years patency of 70% for PTA and VGB respectively [4].

Mathematical models, combined with computational simulations, are powerful tools that can offer a virtual environment to test clinical hypotheses and to drive the improvement of surgical interventions.

In this optical, many works simulated the post-surgical follow-up of both PTA and VGB procedures by using a heterogeneous thread of computational techniques [5,6,7,8,9,10,11,12,13]. Specifically, Agent-Based Models (ABMs) effectively describe pathophysiological processes in which spatial interactions play a major role [14], even though most of them perform their simulations by starting from healthy arteries [5, 7, 8, 10,11,12]. In some, the atherosclerotic plaque was introduced to provide a more realistic geometrical configuration for the pre-intervention condition, but without being generated as consequence of the pathological agents’ dynamics [13]. To the best of our knowledge, only few works addressed the modeling of atherosclerotic plaque formation and progression with an ABM approach, most of them by focusing on the inflammatory and immune events occurring in early atherosclerosis [15,16,17].

Our hypothesis is that an effective predictive model of surgical outcome must replicate the multiscale dynamics that lead to the formation of the atherosclerotic plaque itself, so that the simulation of an intervention is performed on a system whose current dynamics are already altered by the pathological status, which thus influences the outcome. The multiscale character of our approach is of crucial importance as triggering mechanisms and pathology development operate on different time and space scales, respectively seconds vs. weeks and tissue vs. cellular.

Accordingly, this work presents a computational model in which, starting from an artery in homeostatic conditions, the perturbation of specific cellular and extracellular dynamics led to the development of the pathology, ending up with a diseased artery.

The presented model was developed on an idealized portion of a Superficial Femoral Artery (SFA), where an ABM, locally replicating the plaque development, was coupled to Computational Fluid Dynamics (CFD) simulations defining the WSS profile at the wall.

2 Methods

Figure 1 shows the structure of the computational framework. Starting from an idealized 3D model of SFA, CFD simulations compute the WSS distribution along the entire lumen interface. Then, a discrete number (\( M = 10 \)) of 2D cross-sections is selected and an ABM is implemented for each of them to locally replicate the cellular and extra-cellular dynamics driving the wall remodeling. The resulting geometrical alterations, described by the ensemble of the M ABMs outputs, imply a variation of the flow/pressure pattern, and consequently of the WSS profile, which needs to be re-evaluated by coupling back the ABMs to the CFD simulation. For this purpose, a new 3D geometry is reconstructed starting from the M newly obtained 2D cross-sections. The CFD simulation is finally re-run with the same modalities and the loop re-starts.

Fig. 1.
figure 1

Computational Workflow. Starting from an idealized 3D model of SFA, the geometry is prepared for the CFD simulation that gives in input to the ABM the WSS profile at the lumen interface. The ABM performs the cellular and extracellular dynamics that modify the geometrical profile. A new 3D geometry is retrieved in output from the ABM module and the loop re-starts.

Full size image

The CFD simulation and the ABM, embedded in the dashed orange box in Fig. 1, constitute the core of the multiscale model. The first simulates at tissue level the average hemodynamics in a heartbeat, which is on the scale of seconds, while the second, replicating plaque growth, operates at cellular level on the scale of weeks.

The continuous interaction between wall remodeling and hemodynamics is captured by an iterative four-step cycle consisting in: (i) geometry preparation and meshing, (ii) CFD simulation, (iii) ABM simulation, and (iv) coupling and retrieval of the new 3D geometry. Cited steps are separately described below.

2.1 Geometry Preparation and Meshing

An idealized 3D geometry resembling a tortuous portion of SFA was initially built by using Rhinoceros® software (v. 6.0, Robert McNeel & Associates, Seattle, WA, USA).

The geometry presents a centerline with length \( L_{cent} = 84.023\, {\text{mm}} \), and circular inlet and outlet sections with diameters \( D_{inlet} = 4.1\,{\text{mm}} \) and \( D_{outlet} = 3.708\,{\text{mm}} \) respectively.

The model was imported in ICEM CFD® (v. 18.0, Ansys, Inc., Canonsburg, PA, USA) and discretized with tetrahedral elements and 5 boundary layers of prism elements close to the walls. The hybrid tetrahedral and prism mesh was created with the Octree method and globally smoothed by imposing 5 smoothing iterations and a quality criterion up to 0.4.

The 3D model was finally exported to Fluent® (v. 18.0, Ansys, Inc. Canonsburg, PA, USA) to perform the CFD simulations (Fig. 2).

Fig. 2.
figure 2

3D geometry preparation and meshing. A 3D portion of SFA is built in Rhinoceros® 6.0 (a) and M = 10 2D cross-sections selected for the future ABM analysis. The geometry is then discretized and meshed in preparation for the CFD simulation (b).

Full size image

2.2 CFD Simulation

Steady-state CFD simulations were performed on the 3D geometry presented in Sect. 2.1. The choice of a steady-state simulation is since at this developmental stage, an inlet pulsatile waveform-based analysis would have been too much time consuming. In addition, the replicated cellular events occur in the time scale of weeks, while the cardiac output waveform is in the order of the seconds. We can then accept that the current ABM implementation is insensitive to the transient of WSS within a single heartbeat.

At the inlet cross-section, a parabolic velocity profile, whose mean velocity was derived from the analysis of patient’s Doppler ultrasound image at the femoral artery level [18], was imposed. A reference zero pressure was applied at the outlet cross-section, while no-slip condition was chosen for the arterial walls, here supposed as a rigid structure. Blood was modeled as a non-Newtonian fluid with a density of 1060 kg/m3 and a viscosity defined by Caputo et al., as in [19].

A pressure-based couple method was adopted, with a least square cell-based scheme for the spatial discretization of the gradient, a second-order scheme for the pressure and a second-order upwind scheme for the momentum spatial discretization.

2.3 ABM Simulation

The following description refers to a single implementation of the ABM, whose basic principles are the same for each loop iteration and for each cross-section selected.

The simulation flowchart is showed in Fig. 3.

Fig. 3.
figure 3

ABM simulation flowchart.

Full size image

The ABM, implemented in Matlab® (v. 2016b, MathWorks, Natick, MA, USA), locally simulates the alterations due to atherosclerotic insurgence at cell-tissue level starting from a homeostatic condition. The model was developed assuming that the risk factors promoting the disease were already present. Accordingly, plaque growth in a specific region was exclusively associated with the alteration of WSS profile.

The main structure of the ABM follows the one developed by Garbey et al. [20] in the study of VGB post-surgical adaptation, differing however in agent types and simulated events.

The model is initialized to reflect the initial 2D geometry. Then, the cellular, extracellular, and lipid dynamics are performed. The latter are strongly influenced by the hemodynamic profile, which can promote either the maintaining of a homeostatic condition or the progression to a pathological one. Cited dynamics imply a certain level of perturbation within the arterial wall, which adjusts to the new geometrical/compositional conformation.

The simulation arrests either if a level of stenosis greater than 50% is recorded, condition after which a revascularization treatment is commonly performed, or at the end of the chosen follow-up time.

The dynamics of each agent is described with a probability law and the simulation is performed with a Monte Carlo method to introduce a level of randomness that simulates the noise characterizing a true experimental setup. The resulting stochastic nature of the ABM makes necessary to run N independent simulations starting from the same initial condition and taking as output the average of them. A reasonable tradeoff between computational time and minimization of the standard deviation was reached by setting N = 10.

Initialization.

The ABM lays on a <130 × 130> hexagonal grid such as each site is surrounded by 6 neighbors that, determining its local environment, affect its state.

According to the vessel’s anatomy, a 2D circular cross-section of the artery is generated with the wall composed by 3 concentric layers, i.e. tunica intima, media and adventitia. Each site within the wall is occupied by an agent randomly assigned to be a cell or Extra-Cellular Matrix (ECM), according to the relative tissue density. Sites belonging to the lumen or to the portion external to the wall are instead empty.

The process of atherosclerotic formation is triggered in regions interested by a low WSS profile where the altered flow damages the endothelial layer, which consequently allows for lipid infiltration and subsequent accumulation within the arterial wall [2].

Accordingly, the ABM is informed with an initial WSS distribution obtained from CFD simulations and used to associate to each site i at the lumen wall a level of endothelial dysfunction, defined as it follows:

$$ D\left( {WSS} \right)^{i} = D^{i} = \left\{ {\begin{array}{*{20}l} {1 - \frac{{WSS^{i} }}{{WSS_{0} }},} \hfill & { if\, WSS^{i} \, < \,WSS_{0} } \hfill \\ {0,} \hfill & { otherwise} \hfill \\ \end{array} } \right.. $$
(1)

\( WSS^{i} \) is the shear stress recorded at site i and \( WSS_{0} = 1\,Pa \) is the threshold of WSS below which the atherosclerotic process is initiated. \( WSS_{0} \) was set in accordance with findings of plaque localization in areas exposed to WSS < 1 [21] and the identification of the physiological range of WSS in the SFA between 1.5–2 Pa [22].

The level of endothelial dysfunction drives the state of alteration, labeled as Ai,k(Di, d), responsible for the atherosclerosis initiation. Ai,k(Di, d) represents the level of alteration originated by the site i and perceived by site k. It is computed for each site belonging to the wall and it propagates instantaneously across it from the i-th site of origin with isotropic diffusion and with a peak of intensity identified by Di, in formula:

$$ A^{i,k} \left( {D^{i} ,d} \right) = A^{i,k} = D^{i} \,*\,e^{{ - \frac{1}{2}\left( {\frac{d}{4\phi t}} \right)^{2} }} , $$
(2)

where d is the distance between sites i and k, and \( \phi \) the diffusion constant.

The ensemble of the alteration states defines the level of inflammation I that interests each site k of the wall. The resulting state of inflammation for a given k writes

$$ I^{k} = \sum\nolimits_{i = 1}^{{N_{L} }} {A^{i,k} } , $$
(3)

where NL is the number of sites initially belonging to the lumen wall.

Although the endothelial layer was not explicitly modeled, its key role in the pathogenesis of atherosclerosis was considered through the implementation of Eqs. (1)–(3). In presence of an atherogenic WSS profile, the mechanisms of plaque formation are activated, i.e. lipid dynamics and increased cellular activity in the intima, while, with an atheroprotective profile, the homeostatic condition of a healthy artery is maintained.

Agents Dynamics.

For designing and tuning purposes, the ABM was first run on a single cross-section for a 2 months follow-up and with a 1-h time step, in correspondence of which the grid is randomly scanned to compute the agents’ dynamics. Once the ABM will be coupled to the CFD module, the running time will be the chosen coupling time and the ABM analysis obviously extended to all the M cross-sections.

Every site containing a cell is accessed with a frequency of \( T_{cell} = 12 \,{\text{h}} \), while the ones containing ECM every \( T_{matrix} = 2\,{\text{h}} \). To determine whether the potential event is happening, a Monte Carlo method is applied, for which a random number \( test\, \in \,[0;1] \) is generated by the CPU and compared with the probability of the event itself, here labeled as \( p_{event} \). If \( p_{event} > test \), the event happens, otherwise the algorithm explores the next active site.

The investigated events were cellular mitosis/apoptosis, ECM deposition/degradation and lipids infiltration. Their probability densities are reported below.

First, a baseline density of probability was defined for both mitosis/apoptosis and ECM deposition/degradation with Eqs. (4) and (5) respectively:

$$ p_{mit} = p_{apop} = \alpha_{1} $$
(4)
$$ p_{prod} = \beta \,*\,p_{deg} = \alpha_{4} $$
(5)

\( \alpha_{1} \), \( \alpha_{4} \) and β were chosen to maintain a homeostatic balance between mitosis/apoptosis and ECM production/degradation. Accordingly, \( \alpha_{1} \) and \( \alpha_{4} \) were set to compensate the different cellular and extracellular timeframe, i.e. \( \alpha_{1} = \frac{1}{{T_{matrix} }} = 0.05 \) and \( \alpha_{4} = \frac{1}{{T_{cell} }} = 0.008 \). Finally, \( \beta \) was calibrated for each layer to compensate the tendency of the model to preferentially degrade ECM due to the prevalence of ECM on cells. To this aim, the model was run for a 2 months follow-up with several tentative β values and the ratio between final and initial ECM concentration, \( \frac{{ECM_{f} }}{{ECM_{i} }} \), was computed. As the choice must ensure the maintaining of a homeostatic condition, \( \beta \) was chosen by interpolating the \( \frac{{ECM_{f} }}{{ECM_{i} }} \) vs. \( \beta \) plot in correspondence of \( \frac{{ECM_{f} }}{{ECM_{i} }} = 1 \), obtaining \( \beta = \left\{ {1.57,1.57,2.5} \right\} \) for intima, media and adventitia, respectively.

Equations (4) and (5) represent the minimum set of equations to reproduce a homeostatic condition. Within the process of plaque formation, i.e. in presence of at least one site at the lumen wall exposed to \( WSS < WSS_{0} \), while media and adventitia maintain baseline activities, cellular mitosis is perturbed within the intima as follows:

$$ p_{mit} = \left\{ {\begin{array}{*{20}c} {\alpha_{1} \cdot\left( {1 + \alpha_{2} I^{k} } \right) if n_{lip} = 0} \\ {\alpha_{1} \cdot\left( {1 + \alpha_{2} I^{k} } \right)\left( {1 + \alpha_{3} n_{lip} } \right)\left\{ {1 + \exp \left( { - d_{lumen}^{k} } \right)} \right\} if n_{lip} \ne 0} \\ \end{array} } \right. $$
(6)

Similarly, the production of ECM in intima writes

$$ p_{prod} = \left\{ {\begin{array}{*{20}c} {\alpha_{4} \cdot\left( {1 + \alpha_{2} I^{k} } \right) if n_{lip} = 0} \\ {\alpha_{4} \cdot\left( {1 + \alpha_{2} I^{k} } \right)\left( {1 + \alpha_{3} n_{lip} } \right)\left\{ {1 + \exp \left( { - dist_{lumen}^{k} } \right)} \right\} if n_{lip} \ne 0} \\ \end{array} } \right. $$
(7)

\( \alpha_{2} = 1.5 \) drives the perturbation due to the inflammation state \( I^{k} \) of the specific site and \( \alpha_{3} = 0.1 \) accounts for the influence of the surrounding lipids, \( n_{lip} \). Finally, \( d_{lumen}^{k} \) is the distance of the site k from the lumen. With this formulation, the probability of mitosis and ECM production increases with the number of surrounding lipids and the closeness to the lumen.

The process of lipids infiltration only allows one lipid per time step to extravasate with a probability defined as it follows:

$$ p_{lipid} = \alpha_{5} \left( {1 + I^{k} } \right)\left\{ {1 + \alpha_{6} \cdot\exp \left( { - dist_{lip}^{k} } \right)} \right\}\left( {1 + \frac{{n_{lip} }}{{\alpha_{7} }}} \right), $$
(8)

where \( \alpha_{5} = 0.05 \) drives the probability of infiltration, while \( \alpha_{6} = 10 \) accounts for the increased probability of a lipid to occupy a site k near to another lipid agent in the wall, whose distance is \( dist_{lip}^{k} \). Finally, the term \( \left( {1 + \frac{{n_{lip} }}{{\alpha_{7} }}} \right) \) allows for the formation of a lipid cluster and subsequent buildup of the fatty plaque, with \( \alpha_{7} = 6 \) being a normalization constant to maintain the ratio in the interval [0; 1].

Tissue Plasticity and Geometrical Regularization.

The perturbation carried by the agents at the level of the wall is recovered by following a principle of free energy minimization. The sites belonging to the wall are re-arranged accordingly.

For example, to free a site for a newly introduced element, the mitotic/synthetic cell moves the surrounding elements either towards the lumen or the external support, respectively if the cell belongs to the intima or to the media/adventitia. The same principle applies in case of vacation of a site, i.e. cell apoptosis or ECM degradation, but with a reverse movement. Finally, the structure of the model is regularized and adjusted at each time step to obtain a smooth profile at each interface.

2.4 Retrieval of the New 3D Geometry

In an optical of an ABM-CFD coupling, \( N = 10 \) independent simulations were conducted for each of the \( M = 10 \) pre-selected cross-sections of the artery. The new 3D geometry was generated from the ensemble of the ABM solutions of each cross-section that mostly resembled the corresponding average solution about (i) lumen radius, (ii) external radius, and (iii) plaque size.

For each ABM output labeled as i = {1, …, N} lumen and external radii and plaque thickness as function of \( \vartheta \) were computed once assumed that they can be represented in polar coordinates. Being \( \vartheta \) the angle in degrees, the quantities were indicated as \( R_{j}^{i} (\vartheta ) \), with j = 1, 2, 3, respectively. The corresponding deviation, \( \Delta^{i} \), from the average configuration, \( \overline{{R_{j} (\vartheta )}} \), was defined in Eq. (9) as:

$$ \Delta^{i} = \mathop \sum \limits_{j = 1}^{3} \mathop \smallint \limits_{0}^{2\pi } w_{j} \sqrt {(R_{j} (\vartheta ) - \overline{{R_{j} (\vartheta ))^{2} }} } d\vartheta , $$
(9)

where \( w_{j} \) is the weight of the j-th measurement.

The cross-section in output of the ABM associated with the lowest \( \Delta \) was selected and the criterion repeated for all the \( M = 10 \) cross-sections to obtain the skeleton of the new geometry. The 3D geometry was finally reconstructed in Rhinoceros® by lofting the lumen profiles of said sections.

Finally, a sensitivity analysis was performed to define the best coupling time between CFD simulations and ABM, i.e. at which time step the ABM simulation needs to be paused to update the hemodynamic conditions according to the geometrical changes. A reasonable coupling time should constitute a good compromise between accuracy of the results and computational effort required.

Three different cases were tested on a 14 days follow-up period. In the first two cases, coupling times of 7 and 3.5 days were adopted for the whole follow-up. In the third case, an adaptive method was used, in which the first coupling was performed after 7 days, and then every 3.5 days.

3 Results and Discussion

3.1 Atherosclerotic Plaque Generation

The ability of the ABM to accurately replicate the plaque formation was evaluated on a single cross-section and on a 2 months follow-up. To reach a fast evaluation of the present model, the dynamics of plaque development was on purpose accelerated.

By using a 16.00 GB RAM CPU, Intel® Core™ i7-4790, with 4 Cores and 8 Logical Processors, the mean time of computation was \( t\bar{c} \) = 25.87 h for each simulation.

Basing on the ABM results, (i) output morphology of the model, and (ii) its sensitivity to the WSS profile were analyzed.

Output Morphology and Robustness.

Under atherogenic conditions, the ABM is able to generate an asymmetric lipid-rich atherosclerotic plaque with features resembling histological evidences, as shown in Fig. 4.

Fig. 4.
figure 4

ABM-Histology Comparison. (a) An ABM cross-section after 2 months of follow-up is compared with (b) a histology of coronary fibrous cap atheroma in a 24-years-old man [23].

Full size image

In the histology shown in Fig. 4(b), a thick layer of fibrous intimal tissue covers the lipid core, while it is not present in the cross-section generated by the ABM in Fig. 4(a). This is since said layer is thought to form once the lipid infiltration process stops and the only active dynamics gets back to be cells and ECM related only. However, in the presented simulations, at the end of the follow-up, lipids were still actively migrating into the intima and building the core up.

Figure 4(a) also highlights that, differently from the ABM developed by Bhui et al. [17], the plaque volumetric growth was due to lipid accumulation, SMCs proliferation, as also reported in [24], and ECM production, being ECM the major extracellular component of fibroatheromas after lipids [25].

In Fig. 5, the temporal dynamics of intimal, Fig. 5(a), lumen, Fig. 5(b), medial, Fig. 5(c), and adventitial area, Fig. 5(d), normalized on their initial value are provided. As expected, a decrease in lumen area is appreciable. Its monotonic trend is the consequence of having replicated only the inward thickening of the intima, neglecting the initial enlargement of the vessel in response to early stages of plaque formation [26], which, instead, was considered in [17]. Also, coherently with the implemented baseline cellular and ECM dynamics, medial and adventitial areas remained stable.

Fig. 5.
figure 5

Temporal dynamics of intima area (a), lumen area (b), medial area (c) and adventitial area (d). Trends are normalized on their respective initial value at time 0. The N = 10 independent simulations are reported in color, while, in bold black, the average trend associated with the respective standard deviation is shown.

Full size image

Sensitivity to WSS Profile.

The response of the ABM to 3 different WSS profiles, corresponding to an increasing level of atherogenic stimulus from case 1 to case 3, was investigated. The classification depends on the percentage of lumen wall exposed to low WSS values and on the magnitude of the lowest peak of WSS (case 1: 0.8% of lumen exposed with lowest peak \( WSS_{min} = 0.98 Pa \); case 2: 16.5% of lumen exposed with lowest peak \( WSS_{min} = 0.69 Pa \); case 3: 52.6% of lumen exposed with lowest peak \( WSS_{min} = 0.10 Pa \)). For each case, 10 simulations were run on a 2 months follow-up and the related average stenosis output was evaluated.

As expected, the ABM responded differently to the 3 WSS profiles, predicting a greater rate of lumen area reduction with an increased level of atherogenicity. Specifically, a 10%, 20% and 80% stenosis for a low, medium and high severity of WSS distribution were obtained as reported in Fig. 6, where case 1 is the least prone to generate the plaque, while case 3 the most prone.

Fig. 6.
figure 6

Sensitivity to WSS Profile. The geometrical evolution of an ABM cross-section on a 2 months follow-up with 15 days-time step is provided for an increasing level of atherogenicity. For each case, one representative cross-section is chosen among 10 independent simulations.

Full size image

3.2 ABM-CFD Coupling

The most suitable coupling profile between ABM and CFD was studied by running the complete framework described in Fig. 1 with a 14 days follow-up and with the modalities outlined in Sect. 2.4.

The solution of the complete framework at day 14 was obtained in about 20 h, 21.5 h and 20.5 h respectively when coupling every 7 days, 3.5 days and when applying the adaptive coupling time. Considering the average computation time of the ABM to be about 30 min per 1 simulated day, the different time consumption among the tested cases lies in the time required at each coupling step to generate the updated 3D geometry, mesh the model and perform the CFD simulation, estimated to be about 50 min.

From a preliminary analysis based on the temporal dynamics of lumen area, the choice of the coupling period did not imply any difference in the output. However, some differences were locally recorded in planes originally exposed to a WSS profile close to 1 Pa. In these planes, being the boundary between non-pathological and pathological condition very thin, a longer coupling time is not suitable to catch potential switching behaviors between the two conditions. Consequently, the choice of the coupling should follow the trends of these planes, identified as uncertain situations.

The latter suggests the use of a hybrid coupling time where the early phase of plaque development is simulated with a 3.5 days-coupling time, which is then extended to 7 days when the system has stabilized. The latter coupling scheme will be tested in future developments.

4 Conclusions

A coupled ABM-CFD multiscale computational framework of atherosclerotic plaque progression was developed.

The novelty of our approach is that the plaque formation is the result of the alteration of cellular and extracellular agents’ behavior at cell/tissue scale triggered by hemodynamic stimuli. The 2D ABM qualitatively replicates the morphological/compositional changes associated with the pathophysiological condition.

The strengths of the model are its multiscale nature and its modularity, which favor further improvements and investigation of additional phenomena. This will allow the model to serve as platform to test in advance the outcome of interventions and pharmacological therapies aimed to restore the physiological circulation on a long-term perspective.

A very first future step will regard the automation of the ABM-CFD coupling process to reduce time consumption during the coupling itself. Moreover, the influence of some settings of the CFD simulation, such as rigid walls and the steady-state analysis will be assessed.

Finally, a calibration of the model will be performed on clinical data by running the described framework starting from a patient-specific geometry.