Vascular Adaptation: Pattern Formation and Cross Validation between an Agent Based Model and a Dynamical System.

Myocardial infarction is the global leading cause of mortality (Go et al., 2014). Coronary artery occlusion is its main etiology and it is commonly treated by Coronary Artery Bypass Graft (CABG) surgery (Wilson et al, 2007). The long-term outcome remains unsatisfactory (Benedetto, 2016) as the graft faces the phenomenon of restenosis during the post-surgery, which consists of re-occlusion of the lumen and usually requires secondary intervention even within one year after the initial surgery (Harskamp, 2013). In this work, we propose an extensive study of the restenosis phenomenon by implementing two mathematical models previously developed by our group: a heuristic Dynamical System (DS) (Garbey and Berceli, 2013), and a stochastic Agent Based Model (ABM) (Garbey et al., 2015). With an extensive use of the ABM, we retrieved the pattern formations of the cellular events that mainly lead the restenosis, especially focusing on mitosis in intima, caused by alteration in shear stress, and mitosis in media, fostered by alteration in wall tension. A deep understanding of the elements at the base of the restenosis is indeed crucial in order to improve the final outcome of vein graft bypass. We also turned the ABM closer to the physiological reality by abating its original assumption of circumferential symmetry. This allowed us to finely replicate the trigger event of the restenosis, i.e. the loss of the endothelium in the early stage of the post-surgical follow up (Roubos et al., 1995) and to simulate the encroachment of the lumen in a fashion aligned with histological evidences (Owens et al., 2015). Finally, we cross-validated the two models by creating an accurate matching procedure. In this way we added the degree of accuracy given by the ABM to a simplified model (DS) that can serve as powerful predictive tool for the clinic.


Introduction
Cardiovascular disease is the leading global cause of mortality and morbidity (Roger et al. 2012), accounting for more than 17.3 million deaths per year that represents the 31% of all global deaths (American Heart Association 2016). In 2013, nearly 801,000 people died in US from heart disease, stroke, or other heart related diseases (American Heart Association 2016). The direct and indirect costs of cardiovascular diseases and stroke total more than $316.6 billion (American Heart Association 2016). Coronary occlusion is the most common type of heart disease (Centers for Disease Control and Prevention 2015), and it prevents the blood from bringing nourishment to a portion of the myocardium, causing the heart attack (Reddy K 2015).
Coronary Artery Bypass Graft (CABG) surgery using an autologous vein graft (typically saphenous) is the most performed procedure in order to bypass the occlusion and to restore the physiological circulation (Alexander JH 2005).
Despite improvements in surgical techniques, the medium and long term efficacy of the procedure is far from being satisfying, also considering that between 3% and 12% of saphenous grafts occlude within 1 month of surgery (Joseph G 1998), a percentage that increases to 10-15% on a follow-up time of 1 year (Harskamp RE 2013).
The occlusion of the vein graft is due to a series of adaptation (or arterialization) events that takes place in the post-surgical period. The variation of hemodynamic conditions, from venous to arterial flow, plays a key role in the arterialization of the graft (de Vries MR 2016) as the exposure to the high pressure/high flow arterial circulation initiates two simultaneous processes: intimal hyperplasia and outward remodeling (Owens CD et al. 2015;Owens CD et al. 2010;Lytle BW et al. 1985).
Intimal hyperplasia is the universal response of a vessel to injury, for which a reduction of shear stress stimulates specific grow factors to switch their status from quiescent to active. Their activation causes a highlighted division of Smooth Muscle Cells (SMC) in the intimal layer, with subsequent synthesis of ExtraCellular Matrix (ECM). These two events combined lead the tunica intima to thicken and to narrow the lumen, and even though a moderate intimal hyperplasia formation is necessary for proper arterialization and long-term graft patency (Zwolak RM et al. 1987;Wallner K et al. 1999;Zhang L et al. 2004), the patency rates of vein grafts diminish immediately after surgery from 98% to 88% within the first month post-surgery owing to acute thrombosis (De Vries MR 2016).
In contrast, outward remodeling is characterized by preservation or loss of the lumen area through reorganization of the cellular and extracellular components in the media ).
Finally, the balance between intimal hyperplasia and outward remodeling determines the success or the failure of the bypass procedure Garbey et al. 2015;Szilagyi DE et al. 1973;Roddy SP et al. 2003).
In the present work we extensively used two mathematical models that replicate the phenomenon of restenosis and predict the graft adaptation outcome: a Dynamical System (DS) ) and an Agent Based Model (ABM) (Garbey et al. 2015). Furthermore, a primary goal was to prove how we were able to replicate experimental evidences that are commonly accepted by the scientific community, such as the fact that a drastic decrease in shear stress favors an uncontrolled proliferation of cells within the tunica intima and that an increase in wall tension is responsible for the augmented proliferation of cells within tunica media Garbey et al. 2015;Szilagyi DE et al. 1973;Roddy SP et al. 2003).
There are several advantages of having a tool able to replicate the main events that impact the final outcome of the procedure: i) to better understand the key events of restenosis; ii) to test in advance clinical hypotheses; iii) to anticipate the outcome of targeted therapies aimed to improve the durability of the graft.
The DS is a heuristic model and it was derived fundamentally from a conceptual diagram based on experimental observations, while the ABM is a stochastic model that starts from a bottom-up approach and implements biological knowledge at the level of the cells by using a cellular automaton principle. Cellular automata are used as mathematical models in order to investigate the selforganization of mechanical systems (Wolfram S 1983). A cellular automaton consists in a regular grid of cells, each of them in a specific state and each of them evolving according to defined rules (Deutsch A et al. 2005). The behavior of each cell is influenced by the particular states of its surrounding neighbors.
In this work, our goal was to:

1.
Replicate the clinical evidences with the ABM (Fig. 1), and retrieve the pattern formations related to intimal hyperplasia and outward remodeling.

2.
Define a matching procedure between the two models. This has a strong clinical interest. The DS can serve as a fast and user friendly predictive tool for the graft outcome, while the ABM is closer to the physiological realty, but it demands a high computational time and its setup is not simple. It is also easier to obtain patient specific data to set the DS than the ABM, which requires extensive biologic laboratory data. Fig. 2 shows the anatomy of a vein graft at the time of implantation.

Anatomy of a vein graft
As described by Garbey and Berceli (Garbey M et al. 2013), a graft suitable for implantation is composed of a thin intimal layer, and a thick medial layer. The first serves as blood-tissue interface, while the second provides structural support for the wall. The intimal layer is separated from the medial layer by a sheet of connective tissue, called Internal Elastic Lamina (IEL). A similar sheet separates the media from the external surface, the External Elastic Lamina (EEL). Fig. 3 shows the conceptual diagram the DS lays on. The aim of the model is to study the temporal dynamic of radius and thickness of the two innermost layers of the graft: the tunica intima and tunica media.

Dynamical System (DS)
Graft plasticity is described as a function of local biological mechanism and of the dynamic of SMCs and the driven mechanical stimuli are wall shear stress and intramural wall tension. Each of them promotes or prevents specific cellular activities and more in detail shear forces impact events occurring in the intimal layer, while tensile forces impact events within the medial layer.
The geometrical model used to describe the vein is a straight, thick, and circumferential symmetric cylinder with internal radius R 1 and external radius R 2 , and internal pressure P 1 and external pressure P 2 .
The mechanical parameters, being τ wall the shear stress at the wall and σ wall the wall tension, are defined according to the geometrical model chosen.
The wall shear stress is given by the formula: (1) U is the maximum velocity of the blood (recorded at the centerline), and μ is the dynamic viscosity of blood.
The wall tension is given by the formula: (2) ( 3) and σ θ (r) is the circumferential tension (or hoop stress) (Kleinstreurer C 2006), where For simplicity, we will abandon the subscript notation τ wall and σ wall in favor of τ and σ.
The dynamic of the cellular events, fully described in our previous publication , directly impacts the radius and thickness of each compartment and is provided by the following expressions (5) In Equation (5), the general variable represents the area of a wall compartment (I = intima, M = media) due to a specific cellular element density (SMC or ECM).
Δτ and Δσ are respectively the deviation of shear stress and tension from their baseline values. The model is formulated in such a way that the perturbations of shear stress and tension are the driving forces for the adaptation. If we let τ 0 be the shear stress and σ 0 be the tension at time t = 0, the growth rates can be expressed as a function of the difference between the mechanical condition at time t and a baseline setting: R IEL is the radius of the Internal Elastic Lamina, the layer of tissue that separates the intima from the media, and R ext is the radius of the graft. Finally are the constant parameters that drive the cellular events: • γ 5 : Outward remodeling driven by shear forces • γ 6 : Outward remodeling driven by tensile forces The dynamical system described by Equation (5) is non-linear because of the dependence of Δτ and Δσ on the areas of the lumen and of the two wall compartments, intima and media.
In order to accurately simulate the clinical evidences, the model implements systematically the following steps:

1.
Initial setup: geometrical properties of the vein graft are assigned.

2.
Basic solution generation: wall shear stress/wall tension baseline values are assigned such as the system is stable at an initial equilibrium point. With , the basic solution represents the healthy vein graft at the time of implantation.

3.
Perturbation of the system: shear/tensile forces and are modified in order to simulate the arterialization of the vein during the post-surgical follow up.

4.
Evolution of the system: the system evolves according to the cellular events activated by the correspondent parameters. The result of the evolution corresponds to the vein graft harvested at the end of the post-surgical follow up.

Agent Based Model (ABM)
An ABM is defined in a way similar to cellular automata (Hwang M 2011). An ABM is a computational model that simulates the actions and the interactions of single and autonomous agents aiming to investigate their action on the whole system. Each element, also called agent, individually assesses its own status and makes decisions on the basis of a set of pre-defined rules (Bonabeau E 2002;Gilbert N 2008).
The ABM implemented in this work lies on a cellular automata principle and studies the temporal evolution of the same variables tracked with the DS, such as thickness and radius of the wall compartments. Fig. 4 shows the conceptual flow chart that shows how the ABM evolves similarly to the DS, as described in section 2.1.
As initial setup, a hexagonal grid of 121×121 dimension is created. The fundamental unit of the grid is represented in Fig. 5. Every site of the ABM (in Fig. 5 labeled with letter i) is surrounded by six neighbors (labeled with letter j 1 , j 2 , …, j 6 ), which influence the state of the said site.
Given the initial lumen radius, wall thickness, intimal thickness, and SMC/ECM ratio, a random initial cellular pattern is generated. The model is built such that every site belonging to the graft wall is either part of the intimal layer or part of the medial layer and it is either occupied by an SCM or by an ECM discretized volume's element. The geometrical model chosen to describe the vein graft is a straight cylindrical pipe with properties similar to the model chosen for the DS. Among them it is important to underline how the assumption of cylindrical symmetry is still valid for the ABM.
A basic solution is generated within 2 months of simulated follow-up that is a time that ensures a stabilization of the initial conditions. The system is then perturbed and it evolves simulating the adaptation of the graft for the next 6 months of follow up.
The time stepping loop explores every single site of ABM. The dynamic of SMCs and ECMs is regulated by probability laws defined in Table 1. The simulation is based on a Monte Carlo method for which, at each cycle, the simulator generates a random number between 0 and 1 for each SMC or ECM that undergoes mitosis/apoptosis and synthesis/ degeneration. This random number is compared with the specific probability of proliferation for a SMC or synthesis/degeneration for an ECM. If the probability is higher than the random number, then the correspondent cellular event occurs, otherwise the cellular event does not happen. For every density of probability defined in Table 1, specific normalization constants need to be define to ensure that the probability lays between 0 and 1 to be comparable to the random number generated. The deviation of shear stress and wall tension from their baseline is again the key for the vascular adaptation exactly as for the DS.
In our implementation, we only considered the mitosis/apoptosis rate for SMCs and the deposition/degradation rate for ECM elements, ignoring any other additional cellular event that may affects the graft morphology. The same approach was followed in the implementation of the dynamical system. Our goal was to be able to replicate the two main phenomena responsible of restenosis by using a minimum set of cellular events and for each cellular event, by using a minimum set of parameters.

Basic Solution
Both intimal hyperplasia and outward remodeling share the same setup for the generation of the basic solution.
A fundamental step is the initialization of the system. We chose an initial lumen radius R l = 30 cells (~ 0.3 cm), an initial wall thickness th w = 20 cells (~ 0.2 cm), an initial intimal thickness th I = 3 cells (~ 0.03 cm), and finally an initial ratio . The parameters to generate the basic solution were chosen in order to compensate the different time scales of the cellular events. An SMC undergoes a potential mitotic/apoptotic event every 12 hours (T cell = 12 hours), while it potentially synthetizes/degrades a piece of ECM every 2 hours (T matrix = 2 hours).
Accordingly, we set and as the driving parameters respectively for the baseline rate of SMC mitosis/apoptosis and ECM synthesis/degradation.
All the parameters driving other cellular events are set to a null value in order to generate the basic solution.
Since the simulation is stochastic, a single result cannot be considered as a robust nor a reliable one. For this reason, the output of the model is systematically evaluated on the mean of 15 independent simulations run with the same conditions. This is valid both for the basic solution generation and for the post-surgical evaluation.

Perturbation and evolution of the system
Both intimal hyperplasia and outward remodeling are based on a feedback mechanism. Intimal thickness is caused by an initial reduction in shear stress. This latter promotes the division of SMC in intimal layer, which gives intimal thickening. However, the thickening of the intima increases blood flow and shear stress, saturating the driving force of the hyperplasia and driving the system to a stable plateau.
On the other hand, outward remodeling is caused by an increase in wall tension, which promotes SMC division in the medial layer, resulting in a thickening of the media. An increase in wall thickness lowers the tension saturating the driving force generating the outward remodeling and stabilizing the system to an equilibrium.

Intimal hyperplasia-
The cellular event we used in order to replicate the hyperplasia of the tunica intima was the SMC mitosis in intimal layer. Its correspondent density of probability, with reference to Table 1, writes (7) (8) Δτ ABM (t) represents the biologic effect of shear stress in intima and it is assumed to decay exponentially with a time constant α 7 =20. y is the distance of the current ABM site from the center of the grid. d SMC is the standard dimension of an ABM site, assumed to be 0.01 cm.
is the difference between the shear stress recorded at time step t and the baseline τ 0 . In order to provide a global reduction of shear stress, the initial value τ 0 has been perturbed by a constant ε = 0.5, representing a 50% increase respect to its initial value.
is the normalization constant and finally α 3 = 0.2 is the constant parameter that drives the proliferation of SMC in intima. Since the model is driven by stochastic laws, the output of a singular simulation cannot be considered significant nor robust enough to represent a trustable result. In order to have a robust output, we run the model with the same conditions for 15 times, recording the state variables every month. The final output was the mean trend of the 15 simulations.
is the difference between the wall tension recorded at time step t and the baseline . We used the superscript notation "ABM" in order to distinguish the parameters of the ABM from the DS.
In order to provide a global increase of wall tension, the initial value has been perturbed by a constant χ = −0.005, i.e. a decrease of 0.5%.
is a normalization constant. Finally α 4 is the constant parameter driving the mitosis in the medial layer and it was set to 0.3.
The output of the system has been evaluated on the mean of 15 simulations, the same as for the intimal hyperplasia case.

Cross validation
We based the cross validation between the DS and the ABM on 3 points:

1.
Calibration of DS on the output of ABM run in intimal hyperplasia regime, i.e. evaluation of γ 1 from ABM output.

2.
Calibration of DS on the output of ABM run in outward remodeling regime, i.e. evaluation of γ 4 from ABM output.

3.
Qualitative comparison of competitive cellular events.
As the first goal of this work was to study separately intimal hyperplasia and wall remodeling, and consequently to setup the ABM with one cellular event at the time (SMC proliferation in intima and in media), we decided to cross-validate in parallel DS and ABM by first following the same separate methodology and then by adding a level of complexity with an analysis that involved two cellular events.
Furthermore, the two cellular events chosen act in a competitive fashion, meaning that they have opposite effects on a specific biological measure that is chosen as sample for the analysis. The goal of this third point is to lay the foundations for the development of a predictive model that can reproduce animal experimental evidences that were retrieved by our group in a recent publication (Zhang Y et al. 2016). The innovative predictive model will be characterized by a level of complexity between multiple scales and adaptation mechanisms that will go beyond our present controlled and simplified "theoretical biology" experiments. Supplemental Fig. 1 shows the general conceptual diagram followed to calibrate the DS, that can be applied indiscriminately both to calibration of γ 1 and γ 4 . Generally, a model calibration is the task of adjusting an already existing model to a reference system, and specifically referring to our case, we want to adjust the DS to the ABM output, by retrieving the value of the unknown variable of the DS.
The reference is the output of the ABM, and it varies depending on the event we are calibrating (intimal thickness or outward remodeling). The model to be adjusted is the DS, and its output is function of the unknown γ, the parameter driving the cellular event. For our scope, it is mandatory to run both models with the same initial conditions.
The calibration is performed by minimizing the "distance" between the output obtained from the ABM and the same output obtained from the DS.
The "distance" is evaluated using the Root Mean Square (RMS) deviation, which is also function of γ.
N is the number of time points, typically one every month. is the output of the dynamical system at the i-th time point, while is the output of the ABM at the same time point.
We minimized the RMS using a genetic algorithm from Matlab® Optimization Toolbox. Once the value of γ was retrieved, the DS was calibrated with a certain degree of error, and it shows the same output of the ABM.
The conceptual scheme is valid both for the calibration of intimal hyperplasia and of outward remodeling. The distinguishing factors are the output chosen as reference, and the unknown parameter to retrieve from the dynamical system.

Calibration of intimal hyperplasia-As
said, the intimal hyperplasia is generated by SMC division within the intimal layer. Consequently, the parameter to be retrieved from the dynamical system is γ 1 .
The reference output is the time dependent dynamic of lumen area obtained from the mean of 15 simulations of the ABM.

Calibration of outward remodeling-
On the other hand, outward remodeling is the result of SMC division within the medial layer. Accordingly, the parameter of the DS to be calibrated is γ 4 .
In this case the reference is the time dependent dynamic of the graft area, again obtained from the mean of 15 simulations of the ABM.

Competitive cellular events-
We calibrated the DS on the ABM output both in intimal hyperplasia and outward remodeling regime by using one single cellular event at a time.
The last stage of the cross validation was to compare qualitatively the output of the two models letting two competitive cellular events act simultaneously.
We set both models to activate at the same time SMC mitosis and ECM degradation in the medial layer. The output of the models that we tracked was the medial area. This resulted to be the most suitable variable to appreciate a competitive trend; the SMC mitosis enhances the thickness of the tunica media, while the ECM degradation prevents it.
Regarding the ABM, the initial conditions are the same already used to generate the intimal hyperplasia and the outward remodeling. The SMC mitosis in medial layer has already been studied by simulating the outward remodeling. The setup used to generate the competitive trend didn't change (Equation (5) and (6)).
With reference to Table 1, the density of probability of ECM degradation in the medial layer writes (13) and α 6 = 0.25 is the constant parameter driving the ECM degradation.
On the other hand, the DS was set up with following parameters: • γ 3 = −0.07; constant parameter γ 3 drives the ECM deposition in medial layer.
Obviously a negative value indicates a degradation of ECM, instead of a deposition.

Generalization of the ABM
As mentioned in section 2, the DS is based on the assumption of circumferential symmetry, and for calibration purposes, the ABM has been used with the same assumption.
With the thickening of the intima we assist to an invasion of the lumen by the intima itself. With a circumferential symmetry assumption, this invasion is symmetric as well. We used the circumferential symmetry condition in order to keep the structure of the graft intact. Indeed, we were more interested in the time dependent evolution of targeted events than in the specific SMC/ECM distribution within the various graft compartments.
Although, by abating this hypothesis, we were able to study the occlusion of the lumen in a way closer to the physiological realty. More important, in this way we were able to add a certain degree of accuracy to our analysis using the ABM, something that was not possible to achieve with the DS.

Results and Discussions
Using the ABM, we retrieved the known patterns of intimal hyperplasia and outward remodeling, and depending on the event, we tracked the temporal evolution of different biological variables. However, we saw in the previous section how we need to establish first the basic solution, in order to replicate the condition of the healthy vein at the implantation before being able to simulate the surgical follow up.

Basic Solution generation
The aim of the generation of a basic solution is to stabilize the system around an equilibrium point that does not differ too much from the initialization of the system itself. To verify it, we tracked the number of SMCs and ECMs in the wall for a time t = 2 months, starting from t = 0 and recording the status of the variables every month. Supplemental Fig. 2 shows how, with the parameters chosen to generate the basic solution (stated in section 3.1), we were able to obtain a stable trend for both SMCs (Supplemental Fig. 2a) and ECMs (Supplemental Fig. 2b).
Furthermore, Supplemental Fig. 2 proves that t = 2 months is indeed a suitable time window in order to generate a stable system.
Once stabilized around an equilibrium point, the flexibility of the ABM allows us to conduct several virtual experiments by simply tuning the parameters that regulates the cellular event to be replicated. Table 2 shows the complete list of all the simulated events with the correspondent parameters setup.

Intimal Hyperplasia pattern formation
As reported in Table 2, the intimal hyperplasia was simulated by setting the parameter regulating the SMC proliferation in intima (α 3 ) to 0.2 and by choosing a global reduction of shear stress equal to the 50% of its initial value (ε = 0.5). All the other parameters were held to the basic solution configuration.
The hyperplasia of the intima is well represented after 6 months. Fig. 6 compares the cellular pattern at the time of the implant (Fig. 6a) with the pattern recorded after a follow up of 6 months ( Fig. 6b) obtained from a single run of the ABM. As expected, the tunica intima has thickened and the lumen has narrowed.
A separate consideration has to be done for the SMC activity within the 2 layers of the wall: it is clear how an accentuated division of SMC has led the intima to grow, augmenting sensibly its size. While this event was to be expected, it looks like an accentuated cellular activity (SMC division again) has been recorded in the tunica media too, even without apparently modify the thickness of the media itself.
This can be due either to the stochastic nature of the simulation (Fig. 6 represents the output of one single ideal simulation), or to the presence of a secondary/indirect activity in the media.
To clarify this point, we studied the temporal evolution of lumen area, intimal area, and medial area on the average of 15 runs of the ABM with reference to Fig. 7.
All the trends have been normalized on their initial values and the post-surgical follow up time is 6 months as mentioned when we described the principles of the ABM. Fig. 7a and 7b confirm our observation on lumen and intimal area retrieved from the comparison of cellular pattern seen in Fig. 6. Fig. 7b shows a quasi-logistic growing trend for the intimal area, while Fig. 7a shows a quasi-logistic decrease for the lumen area. Indeed, both of them present an initial exponential growth (or a decrease, depending on the variable studied), followed by an inflection point, and finally by a plateau, that represents the new equilibrium of the system after having recovered the perturbation.
Both figures well represent the feedback mechanism governing the whole system. At the time of implant, Δτ ABM (t) has its maximum value and it drives the SMC division in intimal layer to make the intimal area grow exponentially. This is well appreciable in Fig. 7b between month 0 and month 4. In parallel, as the intima thickens, it starts to invade the lumen, diminishing the lumen area. This is clear from Fig. 7a between month 0 and month 4. However, the encroachment of the lumen causes a progressive shear stress reduce, which in turn causes a reduce of Δτ ABM (t). Consequently, the rate of SMC division starts to slow down and the thickening of the intima diminishes as well. The latter is clearly appreciable between months 4 and 5. The plateau is reached once the difference of shear stress has been completely saturated, SMC division ceases, and the solution stabilizes to an equilibrium point at month 6.
Finally, we can see from Fig. 7c how, on an average of 15 simulations, the medial thickness substantially didn't change. Basing on this, we can state that the increase of SMC activity recorded in medial thickness (Fig. 6b) was only due to the stochastic nature of the model, and it doesn't cause any significant variation to the thickness of the media.
The modularity of the ABM allowed us to solely focus on the effect that a reduction of shear stress across the wall employs at cellular level and especially at intimal level. As largely described in this section, it is clear how low shear stress promotes the neointimal generation fostering intimal hyperplasia, and favoring in this way the restenosis of the graft. This outcome is largely confirmed with experimental evidences presented in literature, where many works address that the uncontrolled proliferation of SMC within tunica intima is favored by low shear stress (Heise M et al. 2003;Ojha M 1984;Salam TA et al. 1996;Jackson M et al. 2009;Huynh TTT et al. 1999).

Outward remodeling pattern formation
The setup of the parameters to generate the outward remodeling is also listed in Table 2. Specifically, the parameter regulating the SMC division in medial layer (α 4 ) was set to 0.3, and χ = −0.005, that corresponds to an increase in wall tension of the 0.5% of its initial value.
As done for the intimal hyperplasia case, we studied the cellular pattern formation for the outward remodeling, represented in Fig. 8.
The principle is the same seen in Fig. 6 for intimal hyperplasia, i.e. Fig. 8a shows the cellular pattern at the time of implant, while Fig. 8b details the pattern after a follow up of 6 months. Also for the outward remodeling the analysis is based on a single ideal run of the ABM.
The expected outcome of the outward remodeling is well represented. Comparing the 2 different stages, it appears how an augmented SMC activity in the medial layer brought the media to thicken, while no significant SMC activity is appreciable in the intimal layer, which maintain the thickness recorded at the basic solution state.
It has to be noticed that the relative variation of medial area is much lower than the relative intimal area variation recorded in the intimal hyperplasia regime. This is only due to the limited grid size chosen for the simulation. Indeed, acting on a grid of 121×121 size, an Intel(R) Xeon(R) CPU E3-1270 V2 @ 3.50GHz machine takes 23 hours in order to simulate 8 months of, comprehensive of 2 months of basic solution generation and 6 months of postsurgical follow up. With a bigger grid, the computational time is forecasted to grow accordingly.
The observations deduced from Fig. 8 were confirmed by studying the trends of the average of 15 independent simulations of the ABM run in outward remodeling regime (Fig. 9) Similarly to what done for the intimal hyperplasia case, Fig. 9a shows the mean trend of intimal area, 9b the mean trend of medial area, and finally 9c the mean trend of graft area.
It is interesting to appreciate how the feedback mechanism is visible even in outward remodeling and almost within the same time frame recorded in intimal hyperplasia regime. Also, the three phases appreciated in intimal hyperplasia are still well appreciable even for medial hyperplasia: an initial exponential growth, followed by an inflection point, and finally by the plateau.
The driving force of the outward remodeling is the difference of wall tension, , which is maximum at the time of implantation. Tensile forces enhance division of SMC in the medial layer, which causes the medial area to increase exponentially. This is visible from month 0 and month 4 from Fig. 9b. Consequently, the same trend is recorded for the graft area (Fig. 9c), which increases linearly with the medial area. It is not surprising that the intrinsic growth of medial area is bigger than the one of the entire wall, considering that the intimal area remains substantially unvaried during outward remodeling, and this latter trend is finely appreciable in Fig 9a, which also confirmed our previous considerations on Fig. 8.
The increase of medial area decreases the current wall tension, reducing the driving force . As for the intimal hyperplasia, the system starts to stabilize around month 4, until the difference of wall tension is totally saturated (month 6). Here SMC division arrests, medial and graft area stops growing and the system reaches equilibrium.
Within this section the effect of an increase of wall tension was studied by solely turning on the parameter that drives the proliferation of SMC in the tunica media. Basing on the considerations retrieved from Figs. 8-9 one can conclude that the outward remodeling is favored by an increase in wall tension. In addition, this consideration is largely corroborated by experimental evidences in literature, both retrievable by previous works from our group (Jiang Z et al. 2004;Hwang M et al. 2013) and by works of other investigators (Glagov S 1997;Huynh TTT et al. 1999).

Cross validation and calibration
4.4.1 Calibration-The principle of the calibration of the DS on the output of the ABM was detailed in 3.3.1 and 3.3.2 and summarized in Table 2. The ABM was initially perturbed by an increment of wall tension of the 0.5% of its basic solution condition value in order to activate the cellular events that depend from wall tensile forces, and the parameters leading the SMC mitosis and the ECM degradation in medial layer were respectively triggered to the value of 0.3 and 0.25. On the other hand, we set the correspondent parameters of the DS to be −0.07 for the ECM degradation (the negative sign simply reverses the effect of ECM deposition) and 0.05 for the SMC mitosis.
In both the analysis we took as reference the average of 15 independent simulations run with the ABM. Specifically, the lumen area and the graft area are the most suitable variable to calibrate the DS on the ABM output respectively in intimal hyperplasia and in outward remodeling regime. For each of them, the outputs were normalized on their initial values. The goal of the calibration was to retrieve the value of the parameters γ 1 and γ 4 , respectively for SMC division in intima and media, able to drive the DS to produce an output reasonably close to the one given by the ABM. From our analysis γ 1 = 2.82 and γ 4 = 3.14.
We evaluated the goodness of the calibration by calculating the relative percentile error with the Percentile Root Mean Square (PRMS) deviation: All the variables of (14) have been already defined with (12). From a qualitative evaluation of the plots it is clear how we were able to match the two models with a high level of accuracy. This is confirmed by the percentile errors recorded. For the intimal hyperplasia, PRMS = 0.58%, and for the outward remodeling, PRMS = 0.47%.

Qualitative cross validation-
Through an examination of competitive cellular events, we evaluated the flexibility and complex dynamics within the ABM that have previously been observed in the DS . The goal of this qualitative cross-validation is to show how the ABM is able to replicate peculiar dynamics that have been already achieved with the DS in published works .
We run the ABM plugging in both the SMC division and the ECM degradation at the same time. We chose to run our analysis in the medial layer and not in the intimal layer because an Garbey et al. Page 15 J Theor Biol. Author manuscript; available in PMC 2018 September 21. initial loss of mass was possible during the simulation. While the medial layer is originally thick enough to handle a negative oscillation of mass, the intimal layer is too thin at the time of implant (~ 0.03 cm, approximately the thickness of 3 sites of ABM), and in presence of an early loss of mass, there would have been the risk to totally make the intimal layer disappearing, losing in this way the integrity of the graft wall.
The results are reported in Fig. 11a and 11b, which show the dynamic of the two biologic processes separately. Again, having to deal with a stochastic simulation, a robust output is given only considering the mean trend of 15 independent simulations. A direct comparison between the two cellular events is shown in the normalized plot of Supplemental Fig. 3. Both the cellular events have their own influence on the medial area. As detailed in the method section, we chose these two events for their opposite effect on the area of the media. This is a conclusion easily retrievable from Fig. 11c and Fig. 11d that show the dynamic of the medial area as an output of both the ABM and DS.
Two separate considerations are retrievable from the Fig. 11c and Fig 11d: i) competitive biologic processes lead to an oscillatory response of the models; ii) we are in effect able to retrieve the same outputs from both models, proving the idea of the cross validation.  Fig. 9b. In it, the same biological measure of Fig. 11c (medial area) was retrieved with the ABM but only by simulating SMC mitosis within medial layer. Indeed, in Fig. 9b the rate of growth of medial area is not affected by any initial decay and the system is driven toward the plateau directly through an exponential rate of growth.
A wide consideration out of Fig. 11 allows us to deeply study the nature of the oscillations recorded. In the early stage of vascular adaptation, the ECM degradation prevails on the rates of SMC division, as it is clear from Supplemental Fig. 3, where the negative slope of ECM dynamic is higher than the positive slope of the SMC dynamic, and leads to a decrease in the medial area (Fig. 11c). After 1 month, the ECM degradation has saturated its effect. The sign of the slope inverts, turning positive, while the SMC dynamic curve maintains its positive slope. Now both of the cellular events enhance the medial area growth at the same time. After an additional month, the medial area has returned to its initial size. The remaining trend is rather similar to the one appreciated for the outward remodeling simulation. The increase in medial area causes a decrease of wall tension, which saturates the driving forces of the remodeling. A plateau is again reached at the 5 th month of follow up.
Finally, a strong proof of the cross validation between the two models is given by the qualitative trend reported in Fig. 11d compared to the mean trend reported in Fig. 11c. The latter clearly show how we can retrieve the same non-monotonic temporal trend of medial area by using the DS instead of the ABM.

Details on the endothelium dynamic
Finally, we focused on the encroachment of the lumen during hyperplasia of the intima, by adding a certain level of detail to the process and by getting closer to the experimental observations on the vascular endothelium structure at the harvesting of the vein graft.
The vascular endothelium is a thin layer that serves as functional barrel between the circulating blood and the vessel wall. Its integrity is essential for the maintenance of the homeostasis, and consequently it is directly implicated in the development of vascular pathologies.
Indeed, according with the work of Cox JL et al. (1991), in the early stages of adaptation the leading mechanism is graft thrombosis, which includes disruption of the endothelium. This results to be one of the key events for the enhancement of intimal hyperplasia. According with Nick Roubos et al. (Roubos et al. 1995), the endothelial loss causes denudation of the surface of intima, causing the deposition of fibrin and platelets and increasing the intimal hyperplasia rate.
This thesis was confirmed in the work of Davies MG et al. (Davies et al. 1993), experiments conducted on rabbit models showed how, in order not to experience an accentuated neointimal hyperplasia, a post-surgical therapy should be aimed to mitigate the endothelial injury and its consequences during the first 5 days following the implantation. Other works confirmed the importance to preserve the endothelial integrity not to incur in an uncontrolled cellular growth within the intimal layer, such as the one of LoGerfo et. al (LoGerfo FW et al. 1981), where it is shown that endothelial morphology is best preserved when the vein grafts were pretreated with papaverine before excision, thesis consolidated in a later work of the same team (LoGerfo et. al 1983). It is interesting to notice how they proved that veins, voluntarily injured to stimulate hyperplasia, showed an accentuated loss of endothelial integrity that can be certainly considered as one of the main trigger events of the restenosis phenomenon. To model this, we run the ABM in intimal hyperplasia regime neglecting the hypothesis of circumferential symmetry. Again 2 months of basic solution generation were followed by 6 months of hyperplasia. We compared experimental evidences of vein graft narrowing with the results of our computation. Our goal was to modify the original model in order to obtain an output as close as possible to the physiological realty. Our hypothesis is that some evidences are only retrievable with the usage of the ABM, while they do not result from the DS. Fig. 12 shows a comparison between the output of the ABM run in the same conditions, but with (Fig. 12a) and without (Fig. 12b) circumferential symmetry assumption.
Even though both the simulations show an evident invasion of the lumen by the intimal layer, the irregular injury of the endothelium appreciated from histologic analysis in the works previously cited (Davies et al. 1993;Lo Gerfo et al. 1981;Lo Gerfo et al. 1983), is retrieved only in Fig 12b, that so represents a closer similarity to the experimental evidences. Therefore, through a very simple modification of our ABM, we were able to finely simulate one of the trigger events of the restenosis, making our model closer to the physiological vascularization of the graft.
Another proof of closeness with the physiological reality is given by the study of Fig. 13, which shows histological evidences of several experimental vein grafts harvested 8 months after the original surgical intervention (Owens CD et al. 2015). Especially from panels F-H it results clear how the invasion of the lumen does not respect any criterion of symmetry. The loss of integrity of the endothelium is not regulated by any recurrent or specific pattern. These experimental evidences are comparable with the outcome obtained in Fig 12b with a high level of confidence. Indeed they both shows how the encroachment of the lumen is characterized by a non-uniform and non-symmetric trend of invasion by SMC and ECM.

Conclusions
With the extensive usage of an ABM we were able to replicate clinical and experimental evidences of the two main events leading to the coronary vein graft restenosis phenomenon.
With our cross validation we proved that a simplified model, the DS, can be easily calibrated in order to be closer to the physiological reality than what it was in its original definition. In this way it may better serve as a powerful predictive tool for the clinic.
We also proved to be able to get closer to the physiological reality with the ABM by abandoning the initial hypothesis of circumferential symmetry of the graft. As shown in the previous section, one of the main cause of the accentuated is the loss of regularity in the endothelium structure, a feature that we were able to finely replicate with a simple modification of our model.
As future developments, we will add an additional level of complexity to the ABM by better discerning the entity of the elements of the wall. So far, the wall has been made either by SMC or ECM. We are interested in how the presence of collagen and elastin modifies the dynamic of the system and how it influences the progression of the restenosis.
The impact of macrophage activity will be considered in future developments of the ABM. The intrinsic SMC and ECM activities will be weighted by a time factor in order to simulate the augmentation of biologic activity occurring right after the injury caused by surgical operation. ).
In real-life application of cellular automata, the traditional definition is considered too restrictive, and there have been many relaxations (Hwang M 2011;O' Sullivan D 2001). In our specific case a future development will consist in abating the hypothesis of hexagonal grid initially assumed and in choosing a more suitable unstructured representation.
The original frameworks described in previous works published by our group Garbey M et al. 2015) are both based on a 2D geometry simplification. The study of a 2D cross-section of the graft is indeed suitable enough to show the dynamics of interest. Once the current methodology will have been quantitatively validated thanks to a large data base of histology cross section images and ad hoc spatial statistics, a future perspective will be to move forward the model for 3D animal specific anatomy in order to address further biological questions on the impact of local wall curvature, anastomosis and bifurcations. Clinical observation of a vein graft success and failure 6 months after surgery(left) replicated with the ABM(right). The black portion represents the lumen, the pink one the tunica intima, and the beige one the tunica media.    Conceptual scheme of the Dynamical System, illustrating the primary interacting elements in vein graft adaptation ). On one hand, a variation of inlet flux provokes a variation of Wall Shear stress that impacts specific cellular events involving both SMC and ECM dynamics, which contribute to modify the thickness of the wall. On the other hand, a variation of transmural pressure dictates the wall tension to drive SMC to proliferate and to synthetize ECM within the medial layer. This process augments the thickness of the media, which, combined with the increase of intimal thickness, impact both the lumen radius and the wall thickness. The variation in geometry of wall and lumen is linked with a feedback system to the hemodynamic conditions creating an interacting network         Example of lumen encroachment simulated a with and b without circumferential symmetry hypothesis.  Histological evidences for negative remodeling and intimal hyperplasia as a cause of late lumen loss in a great saphenous vein bypass graft. Significant restenosis is recorded in graft A-B and in graft F-H, in which the stenotic areas demonstrate loss of total vessel area (Owens CD et al. 2015).  Table 2 Parameters setting and conditions for virtual experiments