A fully coupled computational fluid dynamics – agent-based model of atherosclerotic plaque development: Multiscale modeling framework and parameter sensitivity analysis

A fully coupled computational fluid dynamics – agent-based model of atherosclerotic plaque development: Multiscale modeling framework and parameter sensitivity analysis / Corti, Anna; Chiastra, Claudio; Colombo, Monika; Garbey, Marc; Migliavacca, Francesco; Casarin, Stefano. In: COMPUTERS IN BIOLOGY AND MEDICINE. ISSN 0010-4825. 118(2020), p. 103623. [10.1016/j.compbiomed.2020.103623] Original A fully coupled computational fluid dynamics – agent-based model of atherosclerotic plaque development: Multiscale modeling framework and parameter sensitivity analysis

here discussed in detail. Moreover, additional work was done to improve the model reliability, with peculiar attention directed to the sensitivity analysis of the ABM, performed to evaluate the model output in response to variations at the level of the input parameters. Indeed, due to the lack of a direct calibration of the driving coefficients of the model on experimental data, those coefficients were selected heuristically, leading to a certain level of uncertainty that needed to be quantified. The sensitivity analysis allowed us to understand the impact of the uncertain inputs on the model response, identifying the most influencing parameters, whose future calibration on experimental data will improve the model accuracy [29]. Finally, the results of this analysis provided further understandings of the model mechanisms, namely a verification of the ABM response with respect to the model laws, and the identification of an unexpected or not considered behavior. Figure 1 shows the structure of the complete multiscale framework that consists of four cyclically repeated steps [28]. First, an idealized 3D model of the lumen of a tortuous portion of healthy Superficial Femoral Artery (SFA) is built and a 3D mesh of the fluid dynamic domain is generated using ICEM CFD (v. 18 Specifically, depending on the WSS profile computed by the CFD simulation, the ABM replicates the physiologic or pathologic arterial wall remodeling occurring in a predefined time-frame (e.g. one week ABM simulated time). At the end of said period, here referred as coupling period, the ABM simulations are stopped to perform a new CFD simulation in the modified (i.e. updated) 3D geometry.

Multiscale framework
Indeed, the geometrical changes computed by the ABM in each 2D plane affect the fluid dynamic domain, implying the need to update the hemodynamics and the WSS distribution by coupling back the ABMs to the CFD module. To do that, a new 3D geometry of the lumen vessel is reconstructed by lofting the luminal curves of the M ABM outputs and the four main steps of the multiscale framework are then re-performed. The entire process stops at the end of a predefined follow-up period (e.g. at two months ABM simulated time). Each of the aforementioned steps, as well as the coupling process, require the user intervention. Fig. 1. Multiscale computational framework [28]. Starting from a 3D model of healthy artery, the physiologic/pathologic wall remodeling is simulated through a four-blocks scheme: (i) geometry preparation and meshing, (ii) CFD simulation, (iii) ABM simulation and (iv) new 3D geometry. The CFD and ABM modules constitute the multiscale core of the framework, acting on the second/tissue and weeks/cell scales, respectively.
The multiscale core is based on the CFD and the ABM modules, embedded in the dashed red box in A step-by-step extended description is provided below with the following order: i) geometry preparation and meshing, ii) CFD simulation, iii) ABM simulation, and iv) retrieval of the new 3D geometry.

Geometry preparation and meshing
A simplified 3D geometry of the lumen of a healthy artery with SFA-like features was initially built using the CAD software Rhinoceros (v. 6 The 3D geometry was imported into ICEM CFD to generate the 3D mesh of the lumen for the CFD simulation. A hybrid tetrahedral mesh with five boundary layers of prism elements was created with the Octree method [30]. As global mesh parameters, an element maximum size of 0.39 mm was set and the curvature/proximity based refinement was enabled with a minimum element size of 0.156 mm and a refinement of 20 edges along a radius of curvature. The five layers of prism elements were generated with an exponential growth law, setting 1.05 as height ratio. Finally, the mesh was globally smoothed by imposing five smoothing iterations and a quality criterion up to 0.4. The resulting mesh, shown in

CFD simulation
Steady-state CFD simulations were performed using Fluent to compute the hemodynamics in the 3D artery model. Since the arterial wall remodeling, computed by the ABM, occurs in the time scale of weeks, while the cardiac output waveform is in the order of the seconds, cellular dynamics were assumed to depend on the average WSS. Accordingly, to avoid excessive time consumption, a steady flow was imposed to approximate the average hemodynamics.
At the inlet cross-section, a constant parabolic-shaped velocity profile with physiological mean velocity was applied. The mean velocity was derived from the analysis of patient's Doppler ultrasound image at the SFA level [31], following a proper scaling to make it consistent with the current inlet area.
At the outlet cross-section, a reference zero pressure was imposed. No-slip wall boundary condition was specified at the arterial lumen, assumed as rigid. A density of ρ=1060 kg/m 3 was set for blood, modeled as a non-Newtonian Carreau fluid, as in [32]. The simulation was run using the pressure-based solver with coupled method as pressure-velocity coupling method, least square cell-based scheme for the spatial discretization of the gradient, second-order scheme for the pressure and second-order upwind scheme for the momentum spatial discretization [33].
At the end of the CFD simulation, WSS profiles were extracted at pre-selected M = 10 internal circular planes perpendicular to the centerline ( Fig. 2A) and used as hemodynamic input to the corresponding ABM. The choice of M = 10 planes guaranteed a reliable reconstruction of the 3D vessel geometry.

ABM simulation
In general, starting from an initial homeostatic condition, the 2D ABM replicates a physiologic or pathologic arterial wall remodeling (depending on the WSS profile) on a vessel wall cross-section, by locally simulating cell mitosis/apoptosis, ECM production/degradation and lipid infiltration in the intima. The model was developed assuming that the risk factors promoting the disease were already present. Accordingly, the process of plaque formation in a specific region is purely driven by the hemodynamic conditions, and specifically by the WSS profile.
The implemented ABM was inspired to the one developed by Garbey et al. [34,35] for the simulation of VGB post-surgical adaptation. However, different vessel structure and composition, as well as agent types and dynamics were implemented in the present work, which also deals with new cellular events. Figure 3 shows the ABM flowchart. Following the geometrical and hemodynamic initialization, at each time step of one hour, the model computes cell/ECM and lipid dynamics that drive the remodeling of the wall. Then, in order to retrieve smooth profiles and guarantee structural integrity, the lumen and external walls are regularized at each iteration until the end of the simulated period. To replicate the cellular and extracellular events, probabilistic behavioral rules were assigned to each agent and the simulations were performed with Monte Carlo method, allowing capturing the intrinsic variability of biological processes. Due to the stochasticity of the present ABM, the output of a single simulation cannot be considered as a representative solution. Thus, N = 10 independent simulations were run starting from the same initial condition and the average trend was evaluated. The choice of N was dictated by the need of a reasonable trade-off between computational time and minimization of the standard deviation.
A basic solution was generated by opportunely calibrating the agent dynamics in order to stabilize the system around an equilibrium working point. This condition, representative of the homeostatic state of a healthy artery, was then perturbed to simulate the process of atherosclerotic plaque formation.
Prior to the building of the fully coupled framework, the ABM behavior was verified both under physiologic and atherogenic conditions. For this purpose, a single CFD-ABM coupling was performed to initialize the ABM with the hemodynamic input and the ABM simulations were run for two months.
Differently, within the fully coupled framework, the ABM running time corresponded to the chosen coupling time between ABM and CFD modules.
The ABM simulations were run on a 16.00 GB RAM CPU, Intel® Core™ i7-4790, with 4 Cores and 8 Logical Processors.
Details on the ABM initialization, agent dynamics and geometrical regularization are provided below.
Initialization. The ABM was implemented on a 2D <130 x 130> hexagonal grid, representing a good compromise between affinity to the isotropic reality and level of complexity and computational efforts. The initial geometry is a 2D circular cross-section composed by 3 concentric layers, i.e. tunica intima, media and adventitia ( Fig. 4A) with the internal and external elastic laminae (IEL and EEL) separating the intima and the media and the media and adventitia, respectively (Fig. 4B). While sites in the lumen and in the external portion to the wall are initially empty, each site within the wall is randomly seeded with a cell or an ECM, coherently with the cell/ECM ratio of each layer [36][37][38].
Intima and media are initialized with SMCs and elastin and collagen as ECM, with a SMC/ECM ratio of 0.72 [36] and a collagen/elastin ratio of 0.63 [38], while fibroblasts and collagen fill the adventitial sites with a ratio of 0.43 [37]. Figures 4C and 4D show cellular and extracellular composition of intima, media and adventitia layers. For simplicity, no differentiated behaviors between SMC and fibroblasts or between elastin and collagen were implemented. Accepted manuscript at https://doi.org/10.1016/j.compbiomed.2020.103623

13
As previously mentioned, the 2D ABM is informed with an initial WSS profile, which can potentially trigger a pathologic vascular remodeling by perturbing the baseline cellular activity and favoring lipid infiltration and accumulation within the arterial wall [3,4]. Indeed, a low WSS affects the endothelial function by down-regulating atheroprotective genes and up-regulating the atherogenic ones, eventually promoting atherosclerotic plaque formation [4,39].
Accordingly, each site i of the lumen wall is initialized with a WSS value obtained from the 3D CFD simulation, and a level of endothelial dysfunction is computed as follows: is the WSS at site i and 0 = 1 is the assumed pathologic/physiologic WSS threshold. 0 was set in accordance with the study of Samady et al. [40], in which areas exposed to WSS lower than 1 Pa developed greater lumen area reduction. Moreover, the chosen threshold agrees with the physiological range of WSS in the SFA, identified to be between 1.5-2 Pa [41].
In the ABM, each dysfunctional endothelial site i, with ≠ 0, triggers a state of alteration that diffuses within the intima through isotropic diffusion, from a peak of intensity with a diffusion constant : from different endothelial sites are summed up to define the global level of inflammation of the k-th site , as follows: NL is the initial number of sites of the lumen wall (i.e. endothelial sites) and the resulting affects the agent dynamics, as described below, promoting atherosclerotic plaque formation.
Since the purpose of the present model was not to accurately replicate the mechanisms of endothelial dysfunction and the early inflammatory processes occurring during atherogenesis, the endothelium and inflammatory cells or molecules were not explicitly modeled. However, Eqs (1) - (3) were implemented to capture the key role of the hemodynamic input in the pathogenesis of atherosclerosis, whose effect, thanks to the mediation of the endothelial layer, is transferred to the interior sites. Specifically, if all the WSS values at the i-th sites are greater than the threshold, = 0 ∀ and = 0 everywhere. Under this condition, defined atheroprotective, the homeostasis of a healthy artery is replicated. On the contrary, if there is at least one site of the lumen wall exposed to a When a site containing cell/ECM is accessed (i.e. when it is in its potentially active state), a Monte Carlo simulation determines whether the potential event is happening or not, as shown in the "event assessment" phase of Fig. S1. The CPU generates a random number ∈ [0; 1] that is compared with the probability of the event itself, , and the event occurs if < .
The cellular events of interest are, as baseline of activity, mitosis/apoptosis for cells and deposition/degradation for ECM, and, only under atherogenic condition, the ABM also simulates the process of lipids infiltration.
In Fig. 5, the workflow adopted for the implementation of the agent dynamics and the corresponding parameter setting is shown. A detailed analysis of the biological processes occurring during atherosclerosis initiation and progression was performed, with a focus on the cellular, extracellular and lipid dynamics, which led to the definition of their probabilistic behavioral rules. The final probability equations depend on a set of coefficients α i , adopted to weigh a specific influencing factor in the global agent behavior or to set the probability in the interval (0:1).
The numerical value assumed by those coefficients was derived through an iterative process in which the output of the ABM was verified in terms of integrity and qualitative resemblance to histological or literature evidences. In the present section, the values obtained from the aforementioned process are listed below each equation and tagged as default values. This procedure allowed us to obtain a reasonable range for each coefficient. However, since they were not experimentally derived, they are associated with uncertainty, which is reflected in the model output. For this purpose, a sensitivity analysis of the α i parameters was performed and detailed in section 2.2.1. In future works, the parameters α i that emerged as driving ABM coefficients will be calibrated against experimental data and the ABM will be finally validated. Before implementing the pathological cellular dynamics, the baseline densities of probability were set for cell mitosis/apoptosis and ECM deposition/degradation to replicate the physiological conditions.
They were defined with Eq. (4) and Eq. (5), respectively: 1 , 4 and β were imposed to guarantee the maintenance of the physiologic cell/ECM ratio defined at the initialization phase for each tissue layer. 1  was introduced to obtain a baseline probability of said event within a realistic unit of measure [34].
While cell agents are responsible for cell mitosis and apoptosis and ECM production, ECM agents are involved in ECM degradation, meaning that the code scans the grid looking for cells or ECM, respectively. It results that, due to the prevalence of ECM on cells, the model has the tendency to preferentially degrade ECM, instead of producing it. Accordingly, to replicate a baseline condition where ECM production and degradation are averagely balanced, an adjusting coefficient β was introduced and calibrated for each layer. Specifically, since the intima and the media layer have the same cell/ECM ratio = 0.72, a single βint/med was defined for these two layers, while a different value, βadv, for the adventitia, being the adventitia composed by a cell/ECM ratio of 0.43. To this aim, ten ABM simulations were run under physiologic conditions with several tentative βint/med and βadv values and, at a 2-months follow-up, the ratio between final and initial ECM, , was computed for each layer. Considering the work of Garbey et al. [34], a first set of ten simulations was run with initial guesses of βint/med = 2.13 and βadv =2.5 and provided > 1 in the intima and media layers and = 1 in the adventitia layer. Other five values corresponding to 1, 1.5, 1.6, 1.75 and 2 were investigated to calibrate βint/med and, by interpolating the vs. plot in correspondence of = Fig. 6. Therefore, β = {1.57,1.57,2.5} were set for intima, media and adventitia, respectively, in order to guarantee stable trends of ECM in each layer under baseline conditions. With said calibrated coefficients, Eq. (4) and Eq. (5) drive the physiological wall remodeling, leading to the replication of the homeostatic state of a healthy artery. Differently, under atherogenic conditions (i.e. when at least one site of the lumen wall is associated with a < 1 ) cellular mitosis and ECM production in the intima are perturbed to model the increased cellular activity involving the intima layer during atherosclerosis. This translates into a modification of the probability densities defined with Eq. (4) and Eq. (5). Specifically, the probability of cell mitosis and ECM production in the intima increases with the inflammation level, the number of neighboring lipids and the closeness to the lumen [42], leading to the following: plaque formation under atherogenic condition, so that the pathological processes in the planes exposed to atherogenic WSS profile arose within two simulated months. This choice, although not realistic, was dictated by the need to reduce the elevated computational time.
Under atherogenic conditions, the ABM also implements the process of lipid infiltration in the intima. In order to simulate an earlier adaptive intimal thickening [3], lipid dynamics is activated once the intima thickens over a given threshold, here set as IT=6 sites. Since circulating low density lipoproteins were not explicitly modeled, the probability of lipid infiltration is computed as the probability of a site k at the lumen wall to allow lipids to invade the intima, expressed by: where 5 = 0.05 sets the event probability in the interval (0;1). The terms 6 exp(− ) and (1 + 7 ) promote lipid clustering, by increasing the probability of a lipid to occupy a site k close to another lipid, whose distance is and whose neighboring lipids is, in turn, . 6 = 10 weighs the distance term between k and its closest lipid, and 7 = 6 is a normalization constant to maintain the ratio in the interval (0;1). Also in this case, the terms and coefficients of Eq. (8) were set following the framework in Fig. 5 to obtain a lipid core resembling histological features [43].
At each time step only one lipid can enter the intima. To determine the site of access for the lipid, the ten sites of the lumen wall with higher are explored. Starting from the most probable site Monte Carlo simulation is applied and if > , then k is the designated site of access, otherwise the following site of the list is investigated, up to ten. This translates in assuming that a lipid has a total number of chances, = 10, to migrate into the intima.
In the present work, it was assumed that lipids might continue entering the intima until the lipid core potentially occupies maximum 15% of the lumen area. Although not representative of real biological mechanisms, this condition allowed replicating a progressive growth followed by a stabilization of the lipid core, which agrees with the choice not to simulate the phenomenon of plaque rupture, more likely associated with a continuous growth of the lipid-rich necrotic core [44]. Moreover, in PAD, which is the context of the current study, atherosclerotic plaques are usually characterized by a smaller lipid core than in coronary artery disease, and are less subjected to rupture, which further corroborates the assumption on the lipid core size [43].
Tissue plasticity and geometrical regularization. To accommodate the production or removal of an element while performing agent dynamics ("event manifestation" phase of Fig. S1), the tissue reorganizes by following a minimum energy principle, according to which agents move along the shortest path to the target site [34]. In case the active site is at the luminal or external wall border, the production of a cell/ECM results in the addition of a new agent, positioned in a random empty space surrounding the active agent itself. In same condition, if the active agent undergoes death or degradation, it is simply removed from the computational domain leaving an empty space. Differently, in case of an agent inside the arterial wall, a pushing or pulling movement of the surrounding elements allows the production or removal of an element. For example, in case of element production, a site adjacent to the mitotic/synthetic cell is freed thanks to the movement of the surrounding elements either towards the lumen or the exterior, respectively if the site is in the intima or in the media/adventitia.
Similarly, when an element is removed, its site is occupied by an inverse movement of the neighboring agents. Agent movement must always comply with the minimum energy principle, with the only exception constituted by the presence of lipid agents along the shortest path. Once they enter the intima layer, lipid agents must maintain their position throughout the entire simulation, thus constituting an obstacle to the movement of the surrounding elements. Consequently, the agent movement is performed along the shortest path that does not involve lipid agents, allowing preserving the lipid core.

Retrieval of the new 3D geometry
Within the fully coupled framework, at the end of the ABM simulation period, = 10 output solutions, different in terms of morphology, composition and plaque features, are obtained for each of the = 10 cross-sections of the artery. Accordingly, an innovative method was developed to select for each cross-section the output configuration that mostly resembled the corresponding average solution in terms of i) lumen radius, ii) external radius, and iii) plaque size. The procedure, described below, allowed building, at the end of each cycle of the framework in Fig. 1 where each j-th quantity is weighed by .
The same criterion was applied for all the = 10 cross-sections and the 3D geometry was finally reconstructed in Rhinoceros by lofting the lumen profiles of the selected configurations.

ABM sensitivity analysis
Since the model output is largely affected by the parameter setting and none of the parameters was derived from experiments, a sensitivity analysis of the ABM parameters was performed. The goals were (i) to evaluate the oscillation of the model solution due to the uncertainty of the parameters or to the possible inter-subjects variability and (ii) to identify the parameters that mainly drive the ABM output.
For this purpose, first a mono-parametric and then a multi-parametric sensitivity analyses were carried out, as detailed below. In both analyses we defined v = {α2, α3, α5, α6, α7, IT, trylip} that consists in the parameter set under investigation. α1 and α4 were not included in the analysis because they were already calibrated in [34]. For each parameter, a triangular probability density function was defined, based on the parameter range and its most probable value, as shown in Tab Aligned with the purpose of the current analysis, all the ABM simulations were initialized with the same WSS profile, corresponding to the one computed at the 9 th plane of the SFA-like geometry ( Fig.   2A). Indeed, being the most critical hemodynamic scenario, it activates a prompt and intense atherogenic response, allowing appreciating the effects of parameter variation in the ABM response within just one month of follow-up. This is convenient, considering the high computational costs required by the sensitivity analysis. Moreover, since the focus of the analysis was the ABM, the hemodynamic update was not considered, but only 2D ABM simulations were run.
Mono-parametric sensitivity analysis. The probability density function of each parameter was divided in five equal probability intervals and the medium value for each interval was considered.
Moreover, two additional values were included in the analysis to explore the ABM behavior at the extremes of the parameter range, thus investigating seven values for each of the seven parameters, shown in Tab. 2. As mentioned above, in this analysis only one parameter at a time was varied, while keeping all the others at their default values, resulting in 49 cases, each with ten replicates.
The results were analyzed in terms of lumen area and intimal content of SMCs, ECM and lipids.
Indeed, the intima is the layer that is mostly affected by the pathologic wall remodeling occurring in atherosclerosis.
The statistical analysis of the results was performed in Matlab. Multi-parametric sensitivity analysis. Latin hypercube sampling (LHS) was adopted to randomly sample the triangular probability density function of each parameter and define the parameter set for the ABM simulations [45]. This method allows exploring the entire range of each parameter and achieves good accuracy with a limited number of simulations compared to simple random sampling [45]. In this study, the probability density functions of the j=7 {α2, α3, α5, α6, α7, IT, trylip} parameters were divided into k=10 equal probability intervals and the LHS matrix (k x j) was generated, identifying the k=10 ABM parameter combinations: For each k-th parameter set, ten simulations were run to account for the inherent stochasticity.
Partial Rank Correlation Coefficients (PRCC) were computed to quantify the correlation of the target outputs (i.e. lumen area and intimal content of SMCs, ECM and lipids) with each parameter, while removing the effect of the remaining parameters. To compute the PRCC, the average target outputs of the ten replicates for each k-th simulation was considered [45]. PRCC can span from -1 to +1, corresponding to a perfect negative/positive correlation, respectively, and a p-value is associated to each correlation to assess the statistical significance. Correlations were considered statistically significant if the corresponding p-value was lower than 0.05.

Sensitivity analysis of the coupling time
Back to the fully coupled CFD-ABM framework, an important decisional step was about the definition of the coupling time, namely at which time step the ABM simulations need to be paused to update the hemodynamics according to the geometrical changes. A short coupling period allows a better control of the model, but implies high computational time and efforts. Accordingly, a compromise between accuracy of the results and computational effort must be reached. To this aim, we developed an innovative technique based on a sensitivity analysis that was performed to assess the influence of the coupling time on the output, by testing three different cases on a 14 days follow-up period. In the first two cases, the ABM was coupled back to the CFD with a frequency of 7 and 3.5 days, respectively, while, in the third case, the first coupling was performed after 7 days, and then every 3.5 days. For each cross-section, the temporal evolution of the lumen area predicted in the three cases was evaluated, as well as the ABM simulation mode at each coupling interval, which can be either physiologic or pathologic, depending on the WSS profile computed at the corresponding coupling step.

ABM replication of homeostasis and atherosclerotic plaque generation
The ABM accurately and robustly replicated both the homeostatic condition of a healthy artery and the formation of an atherosclerotic plaque, when subjected to the hemodynamic stimuli. Figure 8 shows a qualitative comparison between histology and the ABM output, with the latter selected among the N independent runs for plane 1 (Fig. 8A) and plane 5 (Fig. 8B), for visualization purposes. When initialized with a physiologic WSS profile, the ABM output on a 2-months simulation did not show any substantial deviation from the initial configuration, as depicted in Fig. 8A. The slight alteration of the wall profile was within the physiological range, guaranteeing the preservation of the lumen and tissue layers areas. On the contrary, under atherogenic conditions, the ABM developed an atherosclerotic lesion with features resembling histological evidences, as shown in Fig. 8B. Both the ABM atherosclerotic output and the corresponding histology presented an asymmetric geometry, due to a focally localized thickening of the intima layer and the formation of a lipid-rich core (Fig. 8B).
However, in the ABM output (Fig. 8B) the thick layer of fibrous intimal tissue covering the lipid core was not present because lipids were still migrating in the intima at the stage of the ABM configuration in Fig. 8B. By implementation, such layer would form once the process of lipid infiltration arrests and SMCs and ECM remain the only active agents, coherently with the fact that lipid core formation precedes the increase of fibrous tissue [46]. However, in the present work, the tissue layer separating the lipid core from the lumen was not fibrotic but normal intima, namely SMCs and ECM. Finally, in the ABM solution, maintenance of baseline thickness and composition of the media and adventitia layer was in good agreement with the histological image. refers to a femoral artery of a 75-years old male subject [47], while in (B) to a coronary fibrous cap atheroma of a 24-years old man [48].
For each plane, the analysis of the temporal evolution of the ABM simulations and outputs allowed a further evaluation of the model dynamics and robustness, both in physiologic and pathologic mode.
Under physiologic condition, stable trends of total cells, ECM and wall area were observed, and final healthy configurations were generated (Fig. S2). As shown in Figs. S2-3, although the variability attributable to the inherent stochasticity, a good agreement among the outputs was appreciable and indicative of a robust replication of homeostasis.
Focusing on the atherogenic condition, further results are provided for the vessel cross-sectional plane 5 in Fig. 9. However, similar considerations apply for all the other cross-sections showing plaque formation. In Fig. 9, the results of the ten ABM simulations of arterial wall remodeling of plane 5 on a 2-months follow-up are illustrated in terms of temporal dynamics (Fig. 9A) and final ABM configurations (Fig. 9B). The normalized temporal trends of intimal, lumen, medial and adventitial area are provided, pointing out the monotonic decrease in lumen area due to the intimal thickening and the stability of media and adventitia layers. A considerable variation among the simulations involved the most active dynamics, i.e. the luminal and intimal areas, while little to negligible deviation was observed in adventitial and medial areas, respectively (Fig. 9A). The ABM output solutions at 2 months are provided in Fig. 9B for each of the a,...,j simulations. All the configurations replicated a pathologic wall remodeling with plaque generation, the latter shown in yellow. Although the intrinsic variability among the outputs due to the stochastic nature of the model, all the solutions agreed in terms of degree of stenosis and plaque size, location and morphology, as well as unaltered media and adventitia. Finally, the severity of the replicated pathology was proportional to the degree of atherogenic character of the WSS profile. In Fig. 10, the temporal evolution of an ABM cross-section out of ten is shown for planes 4, 5 and 9, providing an example of the ABM sensitivity to WSS. The percentage of lumen wall exposed to WSS < 1 Pa was 0.8, 16.5 and 52.6, respectively, while the lowest recorded WSS was 0.98 Pa, 0.69 Pa and 0.10 Pa. These WSS profiles triggered wall responses with different degree of intensity. Specifically, in plane 4, after two months the process of lipid infiltration was only at the beginning, with few lipid agents in the intimal layer, differently from plane 5 and 9, where lipids started migrating into the intimal layer within the first month, leading to the generation of a well discernible lipid core. In plane 9, the pathologic wall remodeling was faster than plane 5; moreover, although after one month the lipid core did not change significantly, the intima continued to grow, also thickening the layer between the lipid core and the lumen. At day 60, the configuration of plane 9 presented a more critical scenario compared to plane 5, in which the size of the lipid core was not stabilized yet and the lumen area was still largely preserved. Coherently with the definition and classification of advanced atherosclerotic plaque proposed by Stary et. al [46], the ABM configuration of plane 5 at day 60 (same as Fig. 8B) qualitatively resembles a type IV lesion, characterized by a dense accumulation of lipids without substantial lumen area change, and might progress to a condition of type V lesion, in which ECM is the major plaque component (as Fig. 10, plane 9 at day 60).
The lumen stenosis at the end of the two months follow-up were 10%, 20% and 80% for plane 4, 5 and 9, respectively.

Mono-parametric sensitivity analysis
Among the investigated parameters, the performed analysis pointed out the presence of one most influent parameter, α2, whose primary effect on the SMC and ECM dynamics propagated to the lipid dynamics and largely affected the predicted lumen area reduction. Figure 11 provides details on the model sensitivity to α2 in terms of intimal SMC (Fig. 11A), ECM (Fig. 11B), and lipids (Fig. 11C) and lumen area (Fig. 11D) showing, for each studied output, the temporal trends along the simulation and seven box plots at one month of follow-up.
In the graphs reporting the time evolution of the variables, for each considered value of the parameter, the colored bold line represents the median trend and the associated band is the interquartile range (IQR 25 th -75 th percentiles). On the right, the corresponding box plots describe the data distribution of the specific variable at the end of the simulation obtained with the specific value of α2.
Acting on SMC and ECM dynamics, α2 indirectly affected also lipid dynamics by anticipating or delaying the process of lipid infiltration in the intima, whose starting moment is clearly evident in Fig.   11C and corresponds to the point on the x-axis at which a number of lipids greater than 0 is first observed. Specifically, the greater α2, the more SMC proliferation and ECM production were promoted, with an observed increase in such event rates and, obviously a higher intimal content of SMC and ECM at the end of the simulation (Figs. 11A-B). In turn, an augmented SMC proliferation and ECM production led to a faster thickening of the intima and, as consequence, to an earlier invasion of lipids in the wall (Fig. 11C). As expected, the starting point for lipid infiltration influenced the number of lipids in the intima observed at one month of follow-up. However, a saturation of the lipid content was observed with α2={2.436; 2.854}, due to a control on the lipid core size introduced in the lipid dynamics algorithm. Moreover, as natural consequence of the large effect that α2 has on the agent dynamics, the lumen area was considerably influenced by such parameter (Fig. 11D). An increase in α2 enhanced lumen area reduction rate, leading to the most critical scenario (i.e. smallest final lumen area) associated with the highest α2=2.854. Significant differences among the data distributions associated with different values of α2 were detected for all the studied outputs (p<0.05), and details on the multiple comparisons are provided in the supplementary material (Tab. S1). Finally, the results shown by the box plots pointed out a clearly monotonic relationship between α2 and the studied outputs. The graphs of the temporal dynamics and box plots of each considered output at the variation of the parameters {α3, α5, α6, α7, IT, trylip} are provided in the supplementary materials with the same modalities used for α2, see Figs. S4-9 and Tabs. S2-5. Within the studied range, parameters α3 and α7 did not show any significant influence neither in the agent dynamics, nor, as consequence, in the lumen area (Figs. S4 and S7).
While the global model output, represented by the lumen area, was affected almost exclusively by α2, which emerged as the driving parameter, the subset of parameters {α5, α6, IT, trylip} was identified to significantly influence lipid dynamics with minor or no effects in the other model outputs. Figure 12 provides the lipid temporal dynamics along one month of simulation and box plots of the final intimal lipid content to the variation of the parameters α5 (Fig. 12A), α6 (Fig. 12B), IT (Fig. 12C) and trylip (Fig. 12D). As regards α5, while in the range between 0.039 and 0.105, no significant differences in the final lipid content were detected, α5 significantly affected that output when decreased below 0.039, with almost inhibition of the lipid intimal infiltration for α5=0.001. The effect of α5 on the lipid dynamics slightly propagated to the ECM dynamics, producing a reduced final ECM content in the intima for α5=0.016 and α5=0.001 (p<0.05) (Fig. S5 and Tab. S2). Similarly to α5, also α6 had an effect on the lipid dynamics only when decreased to 0.245, although in this case a slight reduction of the lipid infiltration rate was produced with a minor consequence to the final lipid content, compared to α5 ( Fig.   S6 and Tab. S3). Furthermore, as expected, IT largely affected the final number of lipids, by controlling the starting point of the infiltration process. However, no propagation to the other dynamics was observed, except for a single significant difference recorded in the final content of ECM between IT=3.728 and IT=16.85, as shown in Fig. S8 and Tab. S4. Additionally, the number of chances for a lipid to invade the intima, trylip, acted on the infiltration rate, by increasing the probability that at each time step a lipid successfully enters. As for IT, only a significant difference was recorded in the final SMC content, between trylip=3 and trylip=20 ( Fig. S9 and Tab. S5). As previously mentioned for α2, a saturation of the lipid content may occur due to a control on the maximum number of lipids.
Finally, none of these parameters had an influence on the lumen area, meaning that lumen area reduction is mostly due to an augmented SMC proliferation and ECM production, rather than lipid accumulation in the intima.  Figure 13 shows the PRCCs between the target model outputs and each input parameters {α2, α3, α5, α6, α7, IT, trylip}. Coherently with the mono-parametric sensitivity analysis, the final intimal content of SMC, ECM and lipids and the final lumen area were the investigated outputs. Although only values of PRCC associated with a p<0.05 were considered as statistically significant, also PRCC values with p≈0.06 were taken into account as weakly significant. In accordance with the previous analysis, α2 was identified as the most influencing parameter with significant highly positive correlations with the final amount of ECM and lipids in the intima and a significant highly negative correlation with the lumen area (p<0.05). α2 highly correlated also with the final content of SMC, but associated with a weak significance (p=0.068). Moreover, the remaining parameters were not found to significantly correlate neither with SMC and ECM intimal content, nor with the final lumen area, although a slight influence was recorded in some cases by the mono-parametric analysis. Finally, high correlations with weak and high significance were detected between α6, IT, trylip and the final content of lipids, with α6 and trylip exhibiting positive correlations (weakly and highly significant, respectively), while IT a negative weakly significant one.

CFD-ABM coupling period
As regards the temporal trend of the lumen area along the total 14 days of simulation, no differences resulted from the choice of a smaller or larger coupling period. Indeed, for each plane, the standard deviations of the curves obtained with the three case studies were connected or partially overlapped, meaning that the error committed by adopting the greatest coupling time (i.e. 7 days) was in the range of the noise of the stochastic ABM simulations. However, two different scenarios emerged in terms of influence of the coupling period on the ABM simulation mode, which can be physiologic or pathologic depending on the WSS profile computed at each coupling step. Indeed, while for planes where the minimum WSS values were far from 1 Pa, there was complete agreement in the ABM simulation mode along the 14 days (Fig. 14A), discordance was observed in planes where the lowest WSS values were close to 1 Pa (Fig. 14B). In the last scenario, the shortest coupling period (i.e. 3.5 days) allowed catching switches between the physiologic and pathologic mode, ignored by the other two cases (see  For each plane the output is analyzed in terms of normalized lumen area over time (curves) and simulation mode, i.e. physiologic/pathologic, along the 14 days of follow-up, for each case study a) coupling every 7 days b) coupling every 3.5 days and c) first coupling after 7 days and then every 3.5 days.

Discussion
Several computational studies have already used ABMs to simulate vascular adaptation processes in response to the alteration of the baseline working conditions [16,17,26,27,34,[18][19][20][21][22][23][24][25]. In particular, in [16][17][18][19][20][21][22][23][24] ABMs of the arterial wall remodeling following stenting procedure were implemented to investigate the mechanisms of in-stent restenosis. However, in those works, the intervention procedure was simulated on a healthy artery, thus neglecting the underlying pathology which, instead, we consider to have a role in the final outcome. Herein we presented a novel framework that simulates the process of atherosclerotic plaque development in relation to the hemodynamics by coupling 2D ABM to CFD simulations in a 3D vessel geometry. Furthermore, we performed a detailed sensitivity analysis to identify the driving coefficients of the coupled model, thus laying the foundations for a future quantitative calibration and validation based on experimental and/or clinical data, which were not addressed in the present work.
The strengths of the present framework are (i) the inclusion of important biological aspects related to the disease, i.e. cellular and extracellular dynamics, and (ii) the definition of a loop where molecular and tissue levels are strictly interconnected, and the effect of a perturbation applied to one node is directly reflected on the other ones. The framework is modular and versatile. After a proper experimental calibration and subsequent validation, it might allow the implementation of additional phenomena, such as drug therapies or intervention procedures, on a model of diseased artery. For instance, the effect of anti-proliferative or LDL-lowering drugs (e.g. statins) might be investigated, either by introducing advection-diffusion-reaction equations for those species, or by including them in the model as agents. In both cases, the current agent rules should be modified to take into account the mutual influence between their dynamics and the newly introduced factor/agent. In this context, it might be possible to study a potential plaque stabilization or regression [49]. Moreover, the model might serve to create a virtual population of patients with different patterns (e.g. size and location) of lipid core, degree of stenosis and, once implemented, fibrotic cap and calcifications. Indeed, although until now the model only includes SMCs, ECM and lipid agents, implementation of calcifications is currently under investigation. This, together with further improvements and an experimental validation, might allow investigating the effects of such patterns on the lesion progression and on the intervention outcomes. Finally, by using previously developed in house techniques [14,50], the framework might be also informed with monocyte related gene expression data following a specific treatment. All these aspects are thought to represent a turning point in terms of ability to predict the treatment outcome.
To the best of our knowledge, the only work that partially resembles the proposed framework is that by Bhui et al. [27] in which a 3D ABM was coupled to CFD simulations to simulate the WSS-driven leukocyte trans-endothelial migration and subsequent plaque formation. However, differently from their model, in our ABM plaque volumetric growth was mainly due to SMCs proliferation, ECM production and lipid accumulation (Figs. 8B, 9B and 10), rather than leukocytes infiltration. This is in accordance with the study of Doran et al. [51], in which a key role to SMCs was recognized, and, the research of Stary et al. [46], who identified ECM as the major extracellular component of fibroatheromas after lipids. Finally, while in the model of Bhui et al. [27] wall remodeling was simulated according to Glagov's phenomenon [52], in the present work, for simplicity, only the stenotic effect of plaque growth was considered. Consequently, a monotonic trend of lumen area reduction was observed in our model (Fig. 9A), while, according to [52], a compensatory enlargement of the vessel wall should preserve lumen area in the early stages of plaque formation.
Although being able to capture some major aspects of the pathologic wall remodeling associated with atherosclerosis, as shown in Figs. 8, 9 and 10, the present ABM does not replicate the formation of a fibrotic cap. However, instead of fibrous tissue, it simulates a thickening of the intima (e.g. increased SMCs and ECM) between the lipid core and the lumen in advanced stages of the plaque, as in Fig. 10.
In view of a future implementation of structural (i.e. mechanical) aspects, it might be important to distinguish the fibrotic cap, since it is a key factor in determining plaque stability/instability [3].
Finally, the ABM successfully captured wall response to different WSS stimuli (Fig. 10). Although endothelial cells were not included in the ABM, the contribution of the WSS on the endothelial dysfunction, subsequently triggering the process of plaque formation, was considered through the implementation of a diffusion equation for the inflammation. Similarly, in previous works aimed at investigating the phenomenon of in-stent restenosis [18,19], the endothelium was not explicitly modeled, but a method to estimate the nitric oxide was used, thus considering its influence in SMCs activity and the effect WSS exerts on it.
The sensitivity analysis performed on the ABM provided insights in the working mechanisms of the model, by identifying the most influencing parameters, the interactions among the agent dynamics and the contribution of SMC, ECM and lipid dynamics in the lumen area change. Both the monoparametric and multi-parametric sensitivity analyses recognized α2 (i.e. weight of the effect of inflammation in cell/ECM dynamics) as the driving parameter. From a designing perspective, this is the most important finding, suggesting that a future calibration of α2 will reduce most of the epistemic uncertainty associated with the model, thus improving the accuracy of the results. Moreover, in a general view, knowing in advance which is the most important parameter to be calibrated allows planning experiments optimized for the specific purpose, thus avoiding a waste of time and resources.
The reasons why α3 did not show any contribution probably lie in the considered range and in the fact that the term (1 + 3 ) only involves a minor portion of SMC/ECM agents, thus constituting a local factor that does not produce a net effect in the cell/ECM dynamics. As regards the lipid dynamics and its driving parameters, the performed sensitivity analysis pointed out some crucial aspects, previously unknown. First, the saturation of the lipid infiltration rate for α5>0.039 and α6>3.488 while keeping all the other parameter fixed, is due to the fact that, with trylip=10 chances, the probability is already enough to allow a lipid to enter at each time step. As a consequence, if only one parameter at time is varied, trylip is the one that mostly controls the rate of lipid accumulation in the intima. This parameter may represent the global endothelium permeability to lipids. Indeed, while the probability of a single endothelial site to allow a lipid to enter the wall is computed as expressed by Eq. 8, the number of sites at each time step potentially favorable for lipid entry constitutes a global measure that largely control the phenomenon. However, when the combined effects of the parameters were investigated, also the contribution of α6 was identified, confirming what previously stated and recognizing the potentialities of multi-parametric sensitivity analysis.
Due to the high computational cost of the ABM (mean computation time of 25.87 hours for a 2months simulation) the LHS/PRCC sensitivity analysis was performed on a small sample size and results were obtained by running 100 ABM simulations. Consequently, only few correlations were identified as statistically significant. However, we decided to take into account also high PRCC associated with p≈0.06 because we attribute the inability to get more significance for those high correlations to the small sample size. Indeed, such correlations corresponded to parameters that were found to significantly affect that specific output in the mono-parametric sensitivity analysis. On the contrary, high p-values were associated with low PRCC values and results of the mono-parametric analysis revealing no relationship between said input and output. A low efficacy of the PRCC is, indeed, associated with non-monotonic input-output relationships [45].
To obtain more reliable correlation estimates between each output and input, it is necessary to hugely increase the sample size, implying large computational effort. For this purpose, the Matlab code of the ABM might be converted in C, which is thought to extremely reduce the computation time for each ABM simulation. This will also allow extending the sensitivity analysis to more input parameters that were not considered in the present work, namely those related to the diffusion equation for the inflammation and the threshold on the WSS condition. In particular, the WSS threshold of 1 Pa is a strong assumption which is thought to largely affect the output of the model. In fact, in case the threshold was lowered to 0.5 Pa, for example, plane 5 in Fig. 10, whose minimum WSS is 0.69 Pa, would remain in physiologic condition as plane 1, and plane 9 would develop a less severe plaque.
Moreover, in a future perspective, the combination of sensitivity analysis and inverse problem solution proposed by Casarin et al. [53,54] can be a valuable tool to further narrow the range of optimal setting of the ABM model coefficients.
The sensitivity analysis on the coupling period revealed that a shorter coupling time should be preferred until the ABM simulation mode for each plane stabilizes, namely until the WSS profile become clearly pathologic/physiologic (i.e. far from 1 Pa) or switching behaviors have not been detected for enough time. However, in this work only three cases of coupling period were investigated.
The automation of the fully coupled CFD-ABM framework will reduce the user time consumption in the coupling processes, namely for the generation of the CFD model, the CFD simulation settings and the initialization of the subsequent ABM simulations. Once automated, a more extensive sensitivity analysis on the coupling period will be possible, as well as a sensitivity analysis on the number M of 2D cross-sectional planes, which was not addressed in the present work, although considered important to determine the effect of the spatial resolution on the results. Finally, the automation of the framework and the conversion of the Matlab code in C, will facilitate the calibration of the ABM parameters on patient-specific geometries and the future validation, which otherwise would require excessive computational efforts and time.

Conclusions
In this methodological work, we developed a multiscale CFD-ABM framework, able to capture the mutual interaction between hemodynamics and arterial wall remodeling in atherosclerosis. The framework successfully simulated plaque formation in areas affected by disturbed hemodynamics of an idealized SFA model and updated the fluid dynamics following plaque growth.
Qualitatively, the ABM replicated the main morphological and compositional changes involved in atherosclerosis, generating a pathologic arterial wall configuration coherent with histological images.
Quantitatively, the output of the model was associated with uncertainty, which was related to its stochasticity and the input parameters. Replicating the ABM simulations N=10 times keeping all the parameters fixed allowed having an estimation of the aleatory uncertainty while the combination of mono-parametric and multi-parametric sensitivity analyses provided an estimation of the output oscillation due to uncertainty in the input parameters.
The sensitivity analysis of the ABM parameters revealed that the lumen area reduction, which is the most clinically relevant effect of plaque formation, was exclusively governed by the weight of the WSS-induced inflammation, represented by α2, which acts on the SMC proliferation and ECM production in the intima layer. As a consequence, α2 was responsible for most of the uncertainty of the model output. This finding suggests that the identification of the exact value for α2 will be a turning point towards the definition of a simplified, but reliable model. Other parameters were found to influence the process of formation of the lipid core, without affecting neither SMC/ECM dynamics, nor the lumen area change. These parameters, having a local effect, are less important to be calibrated, but may express the inter-variability of the lipid core size among individuals.
In conclusion, the results of the sensitivity analysis lay the foundations for a future parameter calibration and model validation based on experimental and/or clinical data, which is required for a more systematic assessment of the reliability and usability of the present multiscale CFD-ABM framework of atherosclerosis.