Multiscale Modeling of Vascular Remodeling Induced by Wall Shear Stress

Objective: Hemodynamics-induced low wall shear stress (WSS) is one of the critical reasons leading to vascular remodeling. However, the coupling effects of WSS and cellular kinetics have not been clearly modeled. The aim of this study was to establish a multiscale modeling approach to reveal the vascular remodeling behavior under the interaction between the macroscale of WSS loading and the microscale of cell evolution.

Methods: Computational fluid dynamics (CFD) method and agent-based model (ABM), which have significantly different characteristics in temporal and spatial scales, were adopted to establish the multiscale model. The CFD method is for the second/organ scale, and the ABM is for the month/cell scale. The CFD method was used to simulate blood flow in a vessel and obtain the WSS in a vessel cross-section. The simulations of the smooth muscle cell (SMC) proliferation/apoptosis and extracellular matrix (ECM) generation/degradation in a vessel cross-section were performed by using ABM. During the simulation of the vascular remodeling procedure, the damage index of the SMC and ECM was defined as deviation from the obtained WSS. The damage index decreased gradually to mimic the recovery of WSS-induced vessel damage.

Results: (1) The significant wall thickening region was consistent with the low WSS region. (2) There was no evident change of wall thickness in the normal WSS region. (3) When the damage index approached to 0, the amount and distribution of SMCs and ECM achieved a stable state, and the vessel reached vascular homeostasis.

Conclusion: The established multiscale model can be used to simulate the vascular remodeling behavior over time under various WSS conditions.

Introduction

Vascular remodeling is a process that the vessel changes its structure and function to adapt to the environmental alterations. This complex process involves endothelial hyperplasia, vascular smooth muscle cell (VSMC) proliferation/apoptosis, and extracellular matrix (ECM) generation/degradation. Vascular remodeling is a common pathology of various vascular disorders, such as atherosclerosis, bypass graft failure, in-stent restenosis, and arteriovenous fistula failure (Browne et al., 2015). It is affected by a variety of internal and external factors, for example, biology, chemistry, and physics, among which hemodynamics plays a particularly important role in vascular remodeling. Hemodynamics-induced low wall shear stress (WSS) is one of the critical reasons leading to vascular remodeling (Chatzizisis et al., 2008; Yu et al., 2015; Luong et al., 2016). WSS is a tangential friction force that blood flow exerts on the vessel wall. The pulsatility of blood flow, the rheological properties of blood flow, and the geometry of blood vessel together determine WSS, which is characterized by direction and magnitude (Chatzizisis et al., 2007). In the geometrically irregular region where disturbed flow occurs, it will produce low or oscillating WSS. The cutoff point of low WSS, physiological WSS, and high WSS varies with species and blood vessel types. Low WSS usually refers to WSS less than 1 Pa (Malek et al., 1999; Wentzel et al., 2012; Li et al., 2014; Wang et al., 2016).

Endothelial cell (EC) has WSS receptors, which can sense the changes of WSS and transduce them into biochemical signals (Baeyens et al., 2015). When VSMC receives biochemical signals, VSMC will produce vasoactive substances to act on VSMC and ECM, which leads to changes in the morphological structure of the vessel wall, ultimately adapts to mechanical environments. Therefore, we call vascular remodeling is a mechanochemical-biological process. This process involves the effect of a variety of vasoactive substances, such as nitric oxide (NO), platelet-derived growth factor, and matrix metalloproteinases (MMP). NO synthesized by NO synthase in EC has strong anti-inflammatory, antiapoptotic, antimitotic, and antithrombotic effects and can protect blood vessels. The activity of NO synthase is regulated by normal WSS, and the protein level of endothelial NO synthase is lower in low WSS, which weakens the protective effect on blood vessels (Yamamoto and David, 2010; Liu and Chen, 2020). Low WSS promotes the expression of platelet-derived growth factor, vascular endothelial growth factor, and pro-inflammatory cytokines, inhibits the expression of transforming growth factor, which contributes to the phenomena of overexpression of growth-promoting factor and the underexpression of growth-inhibiting factor, and ultimately promotes the phenotype transformation, migration, differentiation, and proliferation of VSMC (Grainger, 2004). The gene expression and activity of MMP are upregulated by low WSS, especially MMP-2 and MMP-9 (Cheng et al., 2006), which are the main proteases related to ECM degradation. ECM degradation will promote VSMC migration. In addition to promoting ECM degradation, low WSS also reduces ECM synthesis. Interferon is a pro-inflammatory cytokine produced by activated T lymphocytes, which responds to low WSS and is an effective inhibitor of ECM synthesis (Tousoulis et al., 2006).

In the past three decades, the computational simulation method has emerged as a powerful tool to investigate biological phenomena, for example, studying causes, mechanisms, and therapeutic approaches of cardiovascular disease. By using the computational fluid dynamics (CFD) method, we could simulate blood flow in a special location of the cardiovascular system and have new insights into hemodynamics. CFD also shows hemodynamic parameters from the flow field, such as velocity, pressure, and WSS. The methods for vascular remodeling simulations can be divided into two types, namely, continuum models and discrete models (Garbey et al., 2015; Escuer et al., 2019). Continuum models usually use partial differential equations to describe the state changes of a system. One of the advantages of continuum models is that partial differential equations provide a clear mathematical relationship between the components in the system. In addition, the continuum models can be used to model all levels of a system. However, the continuum models do not reflect the state of discrete individuals and the interaction between discrete individuals in the system, while the discrete models such as the agent-based model (ABM) simulate each individual behavior and interactions among them. Therefore, the difference between the two types is that the former focuses on the overall state of the system, while the latter focuses on the interaction between individuals in the system. ABM is a discrete time, discrete event, rule-based modeling method. ABM includes various types, for example, lattice-based, lattice-free, Cellular Potts, lattice-gas, and subcellular element modeling methods (Wang et al., 2015). In ABM, each element usually represents an agent. Agents perform corresponding behaviors in the simulation process according to defined rules, including their own independent behaviors and interactions with other agents. Based on the characteristics of ABM, ABM is widely used in the simulation of biological phenomena, such as tumor microenvironment (Norton et al., 2019), bacterial activity in the soil habitat (Borer et al., 2019), and the transmission of new coronavirus disease 2019 (COVID-19) (Rockett et al., 2020). In addition, ABM is also widely used to simulate the activities of EC and VSMC. These studies mainly focus on in-stent restenosis (Boyle et al., 2010; Thorne et al., 2011; Tahir et al., 2014; Zahedmanesh et al., 2014; Zun et al., 2017; Keshavarzian et al., 2018, 2019; Li et al., 2019). For example, Zahedmanesh et al. (2014) investigated the role of MMP and ECM changes during in-stent restenosis by using a multiscale mechanobiological model. This multiscale mechanobiological model is 2D. Based on the 2D model, Zun et al. (2017) developed a 3D multiscale model of in-stent restenosis to study the effects of stent deployment depth, stent width, and reendothelization speed on the process of restenosis. A fully coupled ABM-finite element analysis (ABM-FEA) framework of in-stent restenosis was developed by Li et al. (2019) to overcome the limitation that the traditional ABM-FEA framework was often partially coupled. The abovementioned studies are based on the combination of ABM and FEA, which is called the multiscale simulation method. The FEA is for the second/organ scale, and the ABM is for the month/cell scale.

As mentioned above, WSS is one of the reasons that lead to vascular remodeling, which is based on experiments and statistics. However, the coupling effects of WSS and cellular kinetics have not been clearly modeled. The aim of this study was to establish a multiscale modeling approach to reveal the vascular remodeling behavior under the interaction between the macroscale of WSS loading and the microscale of cell evolution. CFD method and ABM, which have significantly different characteristics in temporal and spatial scales, were adopted to establish the multiscale model. The CFD method was used to simulate blood flow in a vessel and obtain the WSS in a vessel cross-section. The simulations of VSMC proliferation/apoptosis and ECM generation/degradation in a vessel cross-section were performed by using ABM.

Materials and Methods

Model Overview

The multiscale model of vascular remodeling contained two parts: the CFD model, which simulated blood flow in a vessel to obtain the WSS in a vessel cross-section, and the ABM, which simulated VSMC proliferation/apoptosis and ECM generation/degradation in a vessel cross-section. These two parts were coupled by a parameter named damage index, which was defined as deviation from the obtained WSS. The damage index was used to determine the probability of VSMC proliferation/apoptosis and ECM generation/degradation.

Computational Fluid Dynamics Model

The vascular remodeling of the proximal left anterior descending branch was investigated in this study to validate the proposed method. The 3D geometry of the proximal left anterior descending branch was a simplified conceptual model. It was reconstructed based on the actual radius and centerline which we got from the CT data, as illustrated in Figure 1A. The proximal left anterior descending branch is the middle part, and we extended both inlet and outlet to obtain a fully developed flow. The proximal left anterior descending branch has a radius of 1.535 mm and a length of 19.253 mm. Meshing and CFD simulations were carried out by using commercial software. The calculation domain was described by a nonstructural grid as shown in Figure 1B. Blood was considered as an incompressible Newtonian fluid with a density of 1,060 kg/m³ and a viscosity of 0.0035 Pa⋅s. Laminar flow at a speed of 50 ml/min was applied (Theodorakakos et al., 2008). The blood was set to be laminar flow because the Reynolds number was < 2,000 when the flow rate was set to 50 ml/min. The WSS was obtained by using steady CFD analysis. The wall of the model was assumed to be rigid with a nonslip boundary condition, and zero pressure was set at the outlet. The WSS of each element in a vessel cross-section, which would be utilized as the input for ABM simulations, was recorded at the end of the CFD simulation.

FIGURE 1

Figure 1. The 3D geometry of the proximal left anterior descending branch (A) and meshing result (B).

Agent-Based Model

An ABM was developed to simulate VSMC proliferation/apoptosis and ECM generation/degradation in a vessel cross-section. The ABM of a vessel cross-section was developed in MATLAB (R2014a). Figure 2 shows the flow diagram of ABM. As initial setup, the initialization of ABM included three parts, namely, geometry initialization, time initialization, and hemodynamic initialization. The simulation area was a 2D domain, which was discretized uniformly into 100 × 100 grids. The proximal left anterior descending branch was a simplified conceptual model, so the cross-section of the vessel is an ideal circle. We set a circular ring where the center was (50, 50) and radius was between 20 and 26 as the boundary of the vessel wall. The center of the lumen was set to be (50, 50), and its radius was 20. The proximal left anterior descending branch was selected as the vascular model, where the actual radius of the proximal left anterior descending branch was 1.535 mm, and the corresponding thickness of the vessel wall was 460.5 μm, which was consistent with the literature (Waller, 1989). Therefore, the rationality of the radius and wall thickness in ABM is acceptable. The ratio of VSMC/ECM was 1. The initial cellular pattern in a vessel cross-section was random according to the published literature (Hwang et al., 2013; Garbey et al., 2015; Casarin et al., 2018). Figure 3 shows the vessel cross-section after geometry initialization. Vascular remodeling was a process that evolves over time, so we set a time clock for VSMC and ECM agents. The cycle of VSMC was 12 h, and that of ECM was 4 h (Garbey et al., 2015). Each VSMC and ECM agent had a random age between 0 and cycle. As simulation began, the age increased according to the time step of 1 h. A function was used to define the damage index of VSMC and ECM as deviation from the obtained WSS as followed:

$$\text{D}=\{\begin{array}{c}\hfill \frac{{\text{WSS}}_0-\mathrm{WSS}}{{\text{WSS}}_0}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{WSS}<{\text{WSS}}_0\hfill \\ \hfill \mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\text{WSS}>{\text{WSS}}_0\hfill \end{array}(1)$$

FIGURE 2

Figure 2. Flow diagram of an agent-based model (ABM) in the vessel cross-section.

FIGURE 3

Figure 3. Geometry initialization.

where WSS₀ is a threshold for damage. The threshold for damage was 1 Pa, since low WSS usually refers to WSS less than 1 Pa (Malek et al., 1999; Wentzel et al., 2012; Li et al., 2014; Wang et al., 2016). That is, when WSS is greater than or equal to 1 Pa, its damage index is 0. When WSS is close to 0 Pa, its damage index is the largest, which is 1. When WSS is between 0 and 1, its damage index is linearly distributed (Boland et al., 2019).

The agent in the innermost layer of the vessel wall got the damage index, and then, the damage index of agents in the other layer was defined by the inner damage index and attenuation coefficient according to distance. That is, for an agent in the outer layer having the closest agent in inner layer, the damage index of the former one was defined as the latter one multiplied by the attenuation coefficient. After hemodynamic initialization, each agent in the vessel wall got a damage index. It showed that the lower the WSS, the greater the damage index. The damage index decreased gradually from the inner wall to the outer wall. During the vascular remodeling process, the lumen area decreased, which increased the WSS. Therefore, in each time step, the damage index was multiplied by a constant e^−0.0075 to replicate the WSS recovery. The damage index was used to calculate the probability induced by WSS, the functions of VSMC division probability, VSMC apoptosis probability, ECM generation probability, and ECM degradation probability were as follows:

$${\text{p}}_{\mathrm{div}-\mathrm{wss}}={\mathrm{\alpha }}_{\text{div}}\text{×D}(2)$$

$${\text{p}}_{\mathrm{apop}-\mathrm{wss}}={\mathrm{\alpha }}_{\text{apop}}\text{×D}(3)$$

$${\text{p}}_{\mathrm{gen}-\mathrm{wss}}={\mathrm{\alpha }}_{\text{gen}}\text{×D}(4)$$

$${\text{p}}_{\mathrm{deg}-\mathrm{wss}}={\mathrm{\alpha }}_{\text{deg}}\text{×D}(5)$$

where α_div, α_apop, α_gen, and α_deg are the parameters used to modify the probability.

Vascular remodeling is a mechanochemical-biological process. We considered the interaction between WSS, biochemical factors, VSMC, and ECM. These biochemical factors included Endothelin (ET), NO, and MMP-9. The mathematical relationships between WSS, biochemical factor content, VSMC, and ECM were established according to the literature (Guo and Kassab, 2009; Thorne et al., 2011; Nirala and Gohil, 2015). ET can be synthesized by a variety of cells, such as EC, VSMC, glial cell, and cardiomyocyte, mainly EC. ET can promote the proliferation and migration of VSMC. The mathematical relationship of WSS, EC, and ET was as follows (Thorne et al., 2011):

$${\text{C}}_{\text{ET}}=M(\mathrm{\delta }+\mathrm{\alpha }(1-{\text{e}}^{-k{\mathrm{\tau }}^n}))(6)$$

where C_ET represents the ET synthesized by a single EC per hour, the unit is pg/cell/h, τ represents WSS, the unit is Pa, and M=8×10⁻⁴pg/cell/h, δ=0.6, α=0.4, k=3.63, n = 1.68.

The NO in the vessel is mainly synthesized by EC, which mainly acts on VSMC and inhibits VSMC proliferation and migration. The mathematical relationship of WSS, EC, and NO was as follows (Guo and Kassab, 2009):

$${\text{C}}_{\text{NO}}=a{\mathrm{\tau }}^4+b{\mathrm{\tau }}^3+c{\mathrm{\tau }}^2+d\mathrm{\tau }+e(7)$$

where C_NO represents the NO synthesized by a single EC per hour, the unit is pg/cell/h, τ represents WSS, the unit is Pa, and a=4.365×10⁻⁷,b=−9.399×10⁻⁷, c=6.348×10⁻⁷, d=−9.939×10⁻⁸, e=9.333×10⁻⁹.

The main physiological role of MMPs is to degrade the ECM. More than 20 kinds of MMPs have been discovered so far. MMP-9 is one of them, which belongs to gelatinase and has a wide range of substrates. The mathematical relationship of WSS, EC, and MMP-9 was as follows (Nirala and Gohil, 2015):

$${\text{C}}_{\mathrm{MMP}-9}=\mathrm{\beta }\mathrm{\tau }+\mathrm{\gamma }(8)$$

where C_MMP−9 represents the MMP-9 synthesized by a single EC per hour, the unit is pg/cell/h, τ represents WSS, the unit is Pa, and β=4.939×10⁻⁸,γ=2.218×10⁻⁷.

After obtaining the content of ET, NO, and MMP-9 at a certain time, the mathematical relationship between the content of biochemical factors and the interaction of VSMC and ECM should be established. Substituting the content of biochemical factors into the following formulas, the content of biochemical factors could be converted into probability:

$${\text{p}}_{\text{ET}}=\frac{16.67}{-1.45\times {10}^{-9}{}^{*}\text{C}_{\text{ET}}+8\times {10}^4}(9)$$

$${\text{p}}_{\text{NO}}=\frac{16.67}{-1.45\times {10}^{-9}{}^{*}\text{C}_{\text{NO}}+8\times {10}^4}(10)$$

$${\text{p}}_{\mathrm{MMP}-9}=\frac{16.67}{-1.45\times {10}^{-9}{}^{*}\text{C}_{\mathrm{MMP}-9}+8\times {10}^4}(11)$$

where C_ET, C_NO, C_MMP−9 represent the content of ET, NO, and MMP-9, respectively.

The p_div−wss, p_apop−wss, p_gen−wss, and p_deg−wss were calculated based on the damage index, and the probability of the biochemical factor was superimposed on the abovementioned probability. Thus, the VSMC division probability P_div, VSMC apoptosis probability P_apop, ECM generation probability P_gen, and ECM degradation probability P_deg could be obtained. Considering that ET promotes the proliferation of VSMC and NO inhibits the proliferation of VSMC, the probability of VSMC division was modified as follows: P_div=p_div−wss + p_ET−p_NO. The probability of VSMC apoptosis was not modified, P_apop=p_apop−wss, because the concepts of promoting proliferation and inhibiting proliferation have nothing to do with apoptosis. Because MMP-9 can promote ECM degradation, the ECM degradation probability was modified as follows: P_deg=p_deg−wss + p_MMP−9, and the probability of ECM generation was not modified, P_gen=p_gen−wss.

Figure 4 shows the flowchart of the VSMC division. Simulating VSMC division mainly includes the following steps: (1) Find VSMC that may divide: select the VSMC that may divide according to the age, when the age of VSMC is an integral multiple of 12, the VSMC is likely to divide, (2) Determine whether the VSMC that may divide will divide: a random probability P_{div_rand} between 0 and 1 is generated for the VSMC that may divide. The random probability P_{div_rand} is compared with the probability division for VSMC. If the probability is larger than the random probability, then division for VSMC happened, otherwise it does not happen. The age of VSMC that has not been divided is increased by 1 h, and the abovementioned steps are repeated in the next time step. VSMC apoptosis and ECM generation/degradation are the same as VSMC division. The age of the new VSMC and ECM agents was 1, and the damage value was equal to the original agent.

FIGURE 4

Figure 4. Flow diagram of the vascular smooth muscle cell (VSMC) division.

In the simulation, we mimicked the process of positive and negative remodeling. Positive remodeling usually occurs at the early stage of vascular remodeling, which delays the degree of vascular stenosis. Only when the plaque burden exceeds 40%, the blood flow will be affected (Glagov et al., 1987). In the late stage of the vascular remodeling, as the vascular stenosis progresses, along with calcium salt deposition and fibrosis in the plaque, it causes further narrowing of the lumen, that is, negative remodeling. During the simulation process, when the plaque burden was less than 40%, the growth direction was outward of the wall; when the plaque burden was greater than 40%, the growth direction was inward of the wall.

Results

The first part of the multiscale model of vascular remodeling was the CFD model, which simulated blood flow in a vessel to obtain the WSS in a vessel cross-section. The WSS of the proximal left anterior descending branch was obtained by the CFD method, as shown in Figure 5A. A low WSS region was observed, and we randomly selected a cross-section in this region. The location indicated by the red arrow was the selected cross-section. Figure 5B shows the WSS of the selected cross-section, the red curve in the figure is drawn based on WSS, and the blue curve is the projection of the WSS data in the X-Y plane. Furthermore, WSS was converted into WSS damage index according to the damage index equation (1), as shown in Figure 5C. The red curve in the figure is drawn based on the damage index, and the blue curve is the projection of the damage index in the X-Y plane. When WSS was greater than the threshold WSS₀, the damage index is 0; when the WSS was less than the threshold WSS₀, the damage index showed the trend that the smaller the WSS was, the greater the damage index was.

FIGURE 5

Figure 5. The results of the computational fluid dynamics (CFD) model. (A) The wall shear stress (WSS) of the proximal left anterior descending branch. The location indicated by the red arrow is the selected cross-section. (B) The WSS of the selected cross-section. The red curve in the figure is drawn based on WSS, and the blue curve is the projection of WSS data in the X-Y plane. (C) The damage index of the selected cross-section. The red curve in the figure is drawn based on the damage index, and the blue curve is the projection of the damage index in the X-Y plane.

The second part of the multiscale model of vascular remodeling was the ABM, which simulated VSMC proliferation/apoptosis and ECM generation/degradation in a vessel cross-section. Since VSMC proliferation/apoptosis and ECM generation/degradation were probabilistic events, the simulation results of each run of the ABM were random. Figure 6 compares the vessel cross-section geometry over time obtained from a single run of the ABM. A total of 800 ticks were simulated, and the results were extracted every 100 ticks. The results from t = 100 to t = 800 are shown in Figures 6A–H, respectively. Lumen is represented in white, VSMC in red, and ECM in blue. In terms of geometry, the results of the cross-sections we obtained were consistent with the clinical histologic images (Varnava et al., 2002; Kroner et al., 2011). As shown in the figure, the vessel wall underwent positive remodeling from t = 1 to t = 400, the lumen area remained unchanged, and the vessel wall thickened outward. After t = 400, the vessel wall thickened inward, and lumen began to narrow, which was called negative remodeling. Comparing the geometry of a vessel cross-section with the WSS results, we can observe that the significant wall thickening region was consistent with the low WSS region. There was no evident change of wall thickness in the normal WSS region. The figure shows not only the change of geometry but also the change of plaque burden. The plaque burden can be calculated from the cross-section of vessels, and it increased gradually and ultimately reaches 51.19%. It can be observed that the change rate of plaque burden changed from fast to slow. This was because the damage index decreases continuously during the simulation, that is, the probability of VSMC proliferation/apoptosis and ECM degradation/generation decreased. At t = 800, the damage index approached to 0, the quantity and distribution of VSMC and ECM achieved a stable state, and the vessel reached vascular homeostasis.

FIGURE 6

Figure 6. The changes of geometry in the vessel cross-section over time were obtained from a single run of the ABM. A total of 800 ticks are simulated, and the results are extracted every 100 ticks. The results from t = 100 to t = 800 are shown in (A–H), respectively. Lumen is represented in white, and the vessel wall consists of VSMC and extracellular matrix (ECM): VSMC in red and ECM in blue. The figure not only shows the change of geometry but also the change of plaque burden.

As mentioned earlier, the simulation results of each run of the ABM were random. This was not only because VSMC proliferation/apoptosis and ECM degradation/generation were the probabilistic events, but also because the spatial distribution and age of VSMC and ECM were random during geometry initialization and time initialization. To minimize the influence of the probabilistic events and the random setting of initial conditions, more stochastic simulations need to be performed. Figure 7 shows 10 stochastic simulations of the ABM and mean trend. Figures 7A–D represents the quantity of VSMC, the quantity of ECM, the ratio of VSMC/ECM, and plaque burden, respectively. The quantity of VSMC increased according to the trend of the logic curve and gradually tended to be stable when t = 700. Under this WSS condition, the quantity of VSMC fluctuated around 1,282 (Figure 7A). The quantity of ECM also increased according to the trend of the logic curve and gradually tended to be stable when t = 700. Under the same WSS condition, the quantity of ECM fluctuated around 1,313 (Figure 7B). As shown in Figure 7C, the results of the 10 stochastic simulations showed that the ratio of VSMC/ECM was symmetrically distributed, and after reaching stability, the ratio of VSMC/ECM was slightly lower than the initial value and was 0.98. This was because the quantity of VSMC was slightly less than that of ECM. The plaque burden had the same trend as the quantity of VSMC and the quantity of ECM. The mean trend of plaque burden under this WSS condition was finally stable at 49%. Thus, although the results of each stochastic simulation were different, the trend of each simulation result remained the same. It also showed that the spatial distribution and age of VSMC and ECM were random during geometry initialization, and time initialization had less effect on the simulation results.

FIGURE 7

Figure 7. Temporal evolution of the quantity of VSMC (A), the quantity of ECM (B), the ratio of VSMC/ECM (C), and plaque burden (D). Notably, 10 stochastic simulations of the ABM are shown in color, and the mean trend is shown in black color.

Discussion

There is growing evidence that the hemodynamics-induced low WSS is one of the critical reasons leading to vascular remodeling (Roux et al., 2020; Mahmoudi et al., 2021). However, the coupling effects of WSS and cellular kinetics have not been clearly modeled. We established a multiscale modeling approach to reveal the vascular remodeling behavior under the interaction between the macroscale of WSS loading and the microscale of cell evolution. In the model of vascular remodeling induced by WSS, the simulations of VSMC proliferation/apoptosis and ECM generation/degradation in a vessel cross-section were performed.

The simulation results demonstrated that the multiscale model could simulate the mechanical-biological behavior between WSS, VSMC, and ECM. The significant wall thickening region was consistent with the low WSS region. There was no evident change of wall thickness in the normal WSS region. This phenomenon is consistent with the conclusion in the literature that low WSS induces vascular remodeling (Yu et al., 2015; Luong et al., 2016; Roux et al., 2020; Mahmoudi et al., 2021). Positive remodeling and negative remodeling were considered in the model of this study. It was observed that when the plaque burden was less than 40%, the vessel wall grew outward, and when the plaque burden was greater than 40%, the vessel wall grew inward, causing lumen stenosis. This phenomenon has been consistent with the previous study (Bhui and Hayenga, 2017). By comparing the results of multiple ABM simulations, it was found that although the results of each simulation were different, the trend was the same, which showed that some randomly distributed parameters in the simulation process have less influence on the results, and the main factor affecting the results was WSS. The WSS was utilized as the input for ABM simulations. Most studies have shown that low WSS is one of the reasons for vascular remodeling. Some studies suggest that vascular remodeling is related to WSS gradient (WSSG), oscillatory shear index (OSI), and relative residence time (RRT) (Gori and Boghi, 2011; Ameenuddin and Anand, 2020). WSSG, OSI, and RRT are the hemodynamic parameters derived from WSS. These parameters can provide further insight into the relationship between WSS and vascular remodeling. In this study, we only considered the effect of low WSS and did not consider the effects of high WSSG, high OSI, and high RRT. The next step of the study can be based on the other three factors to explore the differences in vascular remodeling under the action of these three factors.

Some studies have used the multiscale model to study in-stent restenosis, using Von Mises stress as input to study the effect of Von Mises stress on VSMC. Our model mainly takes WSS as input from the perspective of hemodynamics. Some models in the previous research did not consider ECM, which resulted in the simulation results being slightly smaller than the experimental data (Zun et al., 2017; Li et al., 2019). Aware of this, we considered the influence of ECM in our model.

There are several model parameters that can be adjusted to affect the dynamics of vascular remodeling. First, the degree of low WSS may be a key determinant of the individual natural history of vascular remodeling. As mentioned above, the cutoff point of low WSS, physiological WSS, and high WSS varies with species and blood vessel types. The threshold of low WSS determines the damage index and further affects the probability. Second, the periods of VSMC agent and ECM agent adopt different values in different studies, and we chose a value suitable for our model (Garbey et al., 2015; Zun et al., 2017; Corti et al., 2020). The best parameters in the model are expected to be verified by the experimental data.

Limitations

Although the multiscale model can be used to reveal the vascular remodeling behavior under the interaction between the macroscale of WSS loading and the microscale of cell evolution, this model still has some limitations. The CFD model and boundary conditions were simplified in this study. The patient-specific CFD model should be reconstructed, and boundary conditions should be obtained to ensure the accuracy of the CFD simulation results. In addition, the vessel wall was assumed to be a rigid wall, and the interaction between blood flow and the wall was not considered, which may have an impact on the calculated results. In the future study, the fluid-structure interaction methods should be considered. During geometry initialization, we simplified the vessel cross-section by using the ideal circle. Real vessel cross-section should be used in the future study to make the model more accurate. WSS updating plays an important role in coupling the CFD and ABM. In this study, we updated WSS by using some a priori functions rather than using the CFD model. This approach only described a change trend and did not realize the real WSS update, which leads to the present model being partially coupled. In the future study, we could select multiple cross-sections of the blood vessel for ABM. Then, a new 3D blood model can be reconstructed by the multiple ABM results. New WSS could be obtained as input after CFD simulation. We look forward to make our model fully coupled by using this method. In our multiscale model, only the simulations of VSMC proliferation/apoptosis and ECM generation/degradation were considered. There is a phenotypic transition of VSMC during vascular remodeling, and our model does not distinguish the phenotype of VSMC. Other factors such as leaky endothelial barrier, lipid deposition, and inflammatory response were excluded in this model, and these factors should be evolved during the next study. The role of three biochemical factors was added to the model, and the content of biochemical factors was used as a function of the number of EC and WSS. This method obviously simplifies the change process of the content of biochemical factor. It is possible to further improve the description of the changes and effects of biochemical factors by establishing convection-diffusion-reaction equations. In the future study, we hope to improve the model, carry out animal experiments, and verify the model by using the results of animal experiments.

Conclusion

With the usage of CFD and ABM, we were able to establish a multiscale modeling approach to reveal the vascular remodeling behavior under the interaction between the macroscale of WSS loading and the microscale of cell evolution. As expected, the thickening area of the vascular wall corresponded to the magnitude of WSS, showing that the lower the WSS, the easier the thickening of the vascular wall. The established multiscale model could be used to simulate the vascular remodeling behavior over time under various WSS conditions.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author Contributions

SC was responsible for modeling, simulation, data analysis, and manuscript preparation. HZ was responsible for language modification. QH assisted in the design of numerical simulation. YZ was responsible for language modification. AQ was responsible for supervision. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the National Natural Science Foundation of China (12172018) and the Joint Program of Beijing Municipal - Beijing Natural Science Foundation (KZ202110005004).

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.

Publisher’s Note

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.

References

Ameenuddin, M., and Anand, M. (2020). A mixture theory model for blood combined with low-density lipoprotein transport to predict early atherosclerosis regions in idealized and patient-derived abdominal aorta. J. Biomech. Eng. Trans. ASME 142:101008. doi: 10.1115/1.4047426

CrossRef Full Text | Google Scholar

Baeyens, N., Nicoli, S., Coon, B. G., Ross, T. D., Van den Dries, K., Han, J., et al. (2015). Vascular remodeling is governed by a VEGFR3-dependent fluid shear stress set point. eLife 4:e04645. doi: 10.7554/eLife.04645

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhui, R., and Hayenga, H. N. (2017). An agent-based model of leukocyte transendothelial migration during atherogenesis. PLoS Comput. Biol. 13:e1005523. doi: 10.1371/journal.pcbi.1005523

PubMed Abstract | CrossRef Full Text | Google Scholar

Boland, E. L., Grogan, J. A., and McHugh, P. E. (2019). Computational modelling of magnesium stent mechanical performance in a remodelling artery: effects of multiple remodelling stimuli. Int. J. Numer. Meth. Biomed. 35:e3247. doi: 10.1002/cnm.3247

PubMed Abstract | CrossRef Full Text | Google Scholar

Borer, B., Ataman, M., Hatzimanikatis, V., and Or, D. (2019). Modeling metabolic networks of individual bacterial agents in heterogeneous and dynamic soil habitats (IndiMeSH). PLoS Comput. Biol. 15:e1007127. doi: 10.1371/journal.pcbi.1007127

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyle, C. J., Lennon, A. B., Early, M., Kelly, D. J., Lally, C., and Prendergast, P. J. (2010). Computational simulation methodologies for mechanobiological modelling: a cell-centred approach to neointima development in stents. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 368, 2919–2935. doi: 10.1098/rsta.2010.0071

PubMed Abstract | CrossRef Full Text | Google Scholar

Browne, L. D., Bashar, K., Griffin, P., Kavanagh, E. G., Walsh, S. R., and Walsh, M. T. (2015). The role of shear stress in arteriovenous fistula maturation and failure: a systematic review. PLoS One 10:e145795. doi: 10.1371/journal.pone.0145795

PubMed Abstract | CrossRef Full Text | Google Scholar

Casarin, S., Berceli, S. A., and Garbey, M. (2018). A twofold usage of an agent-based model of vascular adaptation to design clinical experiments. J. Comput. Sci. 29, 59–69. doi: 10.1016/j.jocs.2018.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Chatzizisis, Y. S., Coskun, A. U., Jonas, M., Edelman, E. R., Stone, P. H., and Feldman, C. L. (2007). Risk stratification of individual coronary lesions using local endothelial shear stress: a new paradigm for managing coronary artery disease. Curr. Opin. Cardiol. 22, 552–564. doi: 10.1097/HCO.0b013e3282f07548

PubMed Abstract | CrossRef Full Text | Google Scholar

Chatzizisis, Y. S., Jonas, M., Coskun, A. U., Beigel, R., Stone, B. V., Maynard, C., et al. (2008). Prediction of the localization of high-risk coronary atherosclerotic plaques on the basis of low endothelial shear stress: an intravascular ultrasound and histopathology natural history study. Circulation 117, 993–1002. doi: 10.1161/CIRCULATIONAHA.107.695254

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, C., Tempel, D., van Haperen, R., van der Baan, A., Grosveld, F., Daemen, M. J., et al. (2006). Atherosclerotic lesion size and vulnerability are determined by patterns of fluid shear stress. Circulation 113, 2744–2753. doi: 10.1161/CIRCULATIONAHA.105.590018

PubMed Abstract | CrossRef Full Text | Google Scholar

Corti, A., Chiastra, C., Colombo, M., Garbey, M., Migliavacca, F., and Casarin, S. (2020). A fully coupled computational fluid dynamics - agent-based model of atherosclerotic plaque development: multiscale modeling framework and parameter sensitivity analysis. Comput. Biol. Med. 118:103623. doi: 10.1016/j.compbiomed.2020.103623

PubMed Abstract | CrossRef Full Text | Google Scholar

Escuer, J., Martinez, M. A., McGinty, S., and Pena, E. (2019). Mathematical modelling of the restenosis process after stent implantation. J. R. Soc. Interface 16:20190313. doi: 10.1098/rsif.2019.0313

PubMed Abstract | CrossRef Full Text | Google Scholar

Garbey, M., Rahman, M., and Berceli, S. A. (2015). A multiscale computational framework to understand vascular adaptation. J. Comput. Sci. 8, 32–47. doi: 10.1016/j.jocs.2015.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Glagov, S., Weisenberg, E., Zarins, C. K., Stankunavicius, R., and Kolettis, G. J. (1987). Compensatory enlargement of human atherosclerotic coronary arteries. N. Engl. J. Med. 316, 1371–1375. doi: 10.1056/NEJM198705283162204

PubMed Abstract | CrossRef Full Text | Google Scholar

Gori, F., and Boghi, A. (2011). Three-dimensional numerical simulation of blood flow in two coronary stents. Numer. Heat Tranf. A Appl. 59, 231–246. doi: 10.1080/10407782.2011.541147

CrossRef Full Text | Google Scholar

Grainger, D. J. (2004). Transforming growth factor beta and atherosclerosis: so far, so good for the protective cytokine hypothesis. Arterioscler. Thromb. Vasc. Biol. 24, 399–404. doi: 10.1161/01.ATV.0000114567.76772.33

CrossRef Full Text | Google Scholar

Guo, X., and Kassab, G. S. (2009). Role of shear stress on nitrite and NOS protein content in different size conduit arteries of swine. Acta Physiol. 197, 99–106. doi: 10.1111/j.1748-1716.2009.01999.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hwang, M., Garbey, M., Berceli, S. A., Wu, R., Jiang, Z., and Tran-Son-Tay, R. (2013). Rule-based model of vein graft remodeling. PLoS One 8:e57822. doi: 10.1371/journal.pone.0057822

PubMed Abstract | CrossRef Full Text | Google Scholar

Keshavarzian, M., Meyer, C. A., and Hayenga, H. N. (2018). Mechanobiological model of arterial growth and remodeling. Biomech. Model. Mechanobiol. 17, 87–101. doi: 10.1007/s10237-017-0946-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Keshavarzian, M., Meyer, C. A., and Hayenga, H. N. (2019). In silico tissue engineering: a coupled agent-based finite element approach. Tissue Eng. Part C Methods 25, 641–654. doi: 10.1089/ten.tec.2019.0103

PubMed Abstract | CrossRef Full Text | Google Scholar

Kroner, E. S., van Velzen, J. E., Boogers, M. J., Siebelink, H. M., Schalij, M. J., Kroft, L. J., et al. (2011). Positive remodeling on coronary computed tomography as a marker for plaque vulnerability on virtual histology intravascular ultrasound. Am. J. Cardiol. 107, 1725–1729. doi: 10.1016/j.amjcard.2011.02.337

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., Lei, L., Hu, Y., Zhang, Y., Zhao, S., and Zhang, J. (2019). A fully coupled framework for in silico investigation of in-stent restenosis. Comput. Methods Biomech. Biomed. Eng. 22, 217–228. doi: 10.1080/10255842.2018.1545017

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Yang, Q., Wang, Z., and Wei, D. (2014). Shear stress in atherosclerotic plaque determination. DNA Cell Biol. 33, 830–838. doi: 10.1089/dna.2014.2480

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C. D., and Chen, F. (2020). Increase of wall shear stress caused by arteriovenous fistula reduces neointimal hyperplasia after stent implantation in healthy arteries. Vascular 28, 396–404. doi: 10.1177/1708538120913748

PubMed Abstract | CrossRef Full Text | Google Scholar

Luong, L., Duckles, H., Schenkel, T., Mahmoud, M., Tremoleda, J. L., Wylezinska-Arridge, M., et al. (2016). Heart rate reduction with ivabradine promotes shear stress-dependent anti-inflammatory mechanisms in arteries. Thromb. Haemost. 116, 181–190. doi: 10.1160/TH16-03-0214

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahmoudi, M., Farghadan, A., McConnell, D. R., Barker, A. J., Wentzel, J. J., Budoff, M. J., et al. (2021). The story of wall shear stress in coronary artery atherosclerosis: biochemical transport and mechanotransduction. J. Biomech. Eng. Trans. ASME 143:041002. doi: 10.1115/1.4049026

CrossRef Full Text | Google Scholar

Malek, A. M., Alper, S. L., and Izumo, S. (1999). Hemodynamic shear stress and its role in atherosclerosis. JAMA J. Am. Med. Assoc. 282, 2035–2042. doi: 10.1001/jama.282.21.2035

PubMed Abstract | CrossRef Full Text | Google Scholar

Nirala, B. K., and Gohil, N. K. (2015). Effect of glycated serum albumin on functional markers in human umbilical vein endothelial cells in the presence of shear stress. J. Mech. Med. Biol. 15:1550026. doi: 10.1142/S0219519415500268

CrossRef Full Text | Google Scholar

Norton, K. A., Gong, C., Jamalian, S., and Popel, A. S. (2019). Multiscale agent-based and hybrid modeling of the tumor immune microenvironment. Processes 7:37. doi: 10.3390/pr7010037

PubMed Abstract | CrossRef Full Text | Google Scholar

Rockett, R. J., Arnott, A., Lam, C., Sadsad, R., Timms, V., Gray, K. A., et al. (2020). Revealing COVID-19 transmission in Australia by SARS-CoV-2 genome sequencing and agent-based modeling. Nat. Med. 26, 1398–1404. doi: 10.1038/s41591-020-1000-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Roux, E., Bougaran, P., Dufourcq, P., and Couffinhal, T. (2020). Fluid shear stress sensing by the endothelial layer. Front. Physiol. 11:861. doi: 10.3389/fphys.2020.00861

PubMed Abstract | CrossRef Full Text | Google Scholar

Tahir, H., Bona-Casas, C., Narracott, A. J., Iqbal, J., Gunn, J., Lawford, P., et al. (2014). Endothelial repair process and its relevance to longitudinal neointimal tissue patterns: comparing histology with in silico modelling. J. R. Soc. Interface 11:20140022. doi: 10.1098/rsif.2014.0022

PubMed Abstract | CrossRef Full Text | Google Scholar

Theodorakakos, A., Gavaises, M., Andriotis, A., Zifan, A., Liatsis, P., Pantos, I., et al. (2008). Simulation of cardiac motion on non-Newtonian, pulsating flow development in the human left anterior descending coronary artery. Phys. Med. Biol. 53, 4875–4892. doi: 10.1088/0031-9155/53/18/002

PubMed Abstract | CrossRef Full Text | Google Scholar

Thorne, B. C., Hayenga, H. N., Humphrey, J. D., and Peirce, S. M. (2011). Toward a multi-scale computational model of arterial adaptation in hypertension: verification of a multi-cell agent based model. Front. Physiol. 2:20. doi: 10.3389/fphys.2011.00020

PubMed Abstract | CrossRef Full Text | Google Scholar

Tousoulis, D., Antoniades, C., Koumallos, N., and Stefanadis, C. (2006). Pro-inflammatory cytokines in acute coronary syndromes: from bench to bedside. Cytokine Growth Factor Rev. 17, 225–233. doi: 10.1016/j.cytogfr.2006.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Varnava, A. M., Mills, P. G., and Davies, M. J. (2002). Relationship between coronary artery remodeling and plaque vulnerability. Circulation 105, 939–943. doi: 10.1161/hc0802.104327

PubMed Abstract | CrossRef Full Text | Google Scholar

Waller, B. F. (1989). Anatomy, histology, and pathology of the major epicardial coronary arteries relevant to echocardiographic imaging techniques. J. Am. Soc. Echocardiogr. 2, 232–252. doi: 10.1016/s0894-7317(89)80084-7

CrossRef Full Text | Google Scholar

Wang, Y., Qiu, J., Luo, S., Xie, X., Zheng, Y., Zhang, K., et al. (2016). High shear stress induces atherosclerotic vulnerable plaque formation through angiogenesis. Regen. Biomater. 3, 257–267. doi: 10.1093/rb/rbw021

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Butner, J. D., Kerketta, R., Cristini, V., and Deisboeck, T. S. (2015). Simulating cancer growth with multiscale agent-based modeling. Semin. Cancer Biol. 30, 70–78. doi: 10.1016/j.semcancer.2014.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Wentzel, J. J., Chatzizisis, Y. S., Gijsen, F. J. H., Giannoglou, G. D., Feldman, C. L., and Stone, P. H. (2012). Endothelial shear stress in the evolution of coronary atherosclerotic plaque and vascular remodelling: current understanding and remaining questions. Cardiovasc. Res. 96, 234–243. doi: 10.1093/cvr/cvs217

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamamoto, M., and David, T. (2010). “Endothelial nitric oxide concentration and its implications in carotid artery atherosclerosis - an integrated cell/haemodynamics approach,”in Proceedings of the IFMBE 6th World Congress of Biomechanics (WCB 2010). August 1-6, 2010 Singapore, PTS 1-3, Vol. 31, eds C. T. Lim and J. C. H. Goh (Berlin: Springer), 414–417. doi: 10.1007/978-3-642-14515-5_106

CrossRef Full Text | Google Scholar

Yu, Y., Cai, Z., Cui, M., Nie, P., Sun, Z., Sun, S., et al. (2015). The orphan nuclear receptor Nur77 inhibits low shear stress-induced carotid artery remodeling in mice. Int. J. Mol. Med. 36, 1547–1555. doi: 10.3892/ijmm.2015.2375

PubMed Abstract | CrossRef Full Text | Google Scholar

Zahedmanesh, H., Van Oosterwyck, H., and Lally, C. (2014). A multi-scale mechanobiological model of in-stent restenosis: deciphering the role of matrix metalloproteinase and extracellular matrix changes. Comput. Methods Biomech. Biomed. Eng. 17, 813–828. doi: 10.1080/10255842.2012.716830

PubMed Abstract | CrossRef Full Text | Google Scholar

Zun, P. S., Anikina, T., Svitenkov, A., and Hoekstra, A. G. (2017). A comparison of fully-coupled 3D in-stent restenosis simulations to in-vivo data. Front. Physiol. 8:284. doi: 10.3389/fphys.2017.00284

PubMed Abstract | CrossRef Full Text | Google Scholar