A Computational Model on Cartesian Adaptive Grid for Thrombosis Simulation

,


I. INTRODUCTION
Thrombosis is a pathological intrinsic coagulation attributed to blood stasis, tissue injury, and hypercoagulability [1]. In view of the fact that thrombosis is a primary cause of major cardiovascular and cerebrovascular diseases such as cerebral infarction, and heart attack, the in-depth studies of thrombosis prevention and treatment depending on thrombogenesis mechanism have been carried out in a variety of ways, such as animal models [2], clinic [3], vitro models [4] and computational simulations [5]. As explored by several studies [6]- [9], the geometry shapes of vascular lumen in the vicinity of a developing clot determine the hemodynamic, which impacts thrombus formation.
The applications of computer technology have been widespread in every area. The mechanism of thrombus formation can also be studied through model based computer simulation experiments. Compared with traditional experimental techniques, computer simulations are superior in low cost, easier operation, and various experimental design [10]. On the ground that the computational model used in a simulation decides the value of the experiment, great progress has been made by numerous researchers in The associate editor coordinating the review of this manuscript and approving it for publication was Stavros Souravlas . developing computational models for thrombosis simulation in the past few decades [11]. Among the numerous studies presented in recent years, hybrid particle-continuum models [12]- [16] demonstrate the advantages of alterable precision and microstructural retention.
In the hybrid models, the fluid-solid interactions on the surface of thrombin will become intractable if thrombin and plasma are modelled in different forms. In existing studies, two approaches are implemented to simplify the fluid-solid interactions: particle-based approach and continuum-based approach. In the models proposed in Tosenberger's [12] and AL-Saad's [16] papers, the thrombin is modelled in particle clusters, and the plasma is modelled with smoothed particles hydrodynamic method. The interactions are reduced to particle-particle interaction, yet the computation of hemodynamics is laborious in the artery which contains a large amount of plasma particles. On the contrary, in the model proposed by Mukherjee [13], both plasma and thrombin are modelled in continuum. A scaffold grid is created based on the discrete particles in the thrombus and then embedded in the continuum field as fictitious domains. Similar strategies are used in the phase-field based model proposed by Hosseinzadegan [14] and the arbitrary Lagrangian-Eulerian formulation based model proposed by Zheng [15]. Even the continuum-based strategies capture macroscale dynamics better than particle-based strategies. There are still problems on: (i) retaining the complex morphological characteristics of thrombus when converting particles into continuum field; (ii) reducing the burden of the data and the calculations acting on the thrombus surface.
It is the main factor of the above problems that the surface of thrombus is a free surface. In the researches of fluid simulation, the adaptive meshes are used as the efficient schemes for the numerical solution of surface moving. In conclusion, the motivations of this study are the application value of hybrid computational models in thrombus researches and the above problems to be solved in the existing models. A computational model for thrombus growth is constructed with Cartesian Adaptive Grid in this paper. The baseline model is a particle-continuum hybrid model, in which plasma, thrombus, and vessel are described in continuum, and the platelets are described in particles. To solve the problems (i), multiphase Level Set Method (LSM) is used to define the interface of two different materials in the continuum field, and the process of thrombus growth is a series of interface expansions with certain speed functions depending on the shape of platelet. Then, an improved Adaptive Mesh Refinement (AMR) strategy is used in LSM to solve the problem (ii). At last, the proposed model is applied to simulating the primary thrombus formation in stenotic arteries. The main contributions of the paper are as following: First, other than covering the whole thrombus into continuum at each time step in simulations, the thrombus and its growth process are implemented in continuum directly, and the conversion is only operated on the platelets aggregated into the thrombus.
Second, to obtain an accuracy thrombus shape and keep an acceptable data amount, a Cartesian adaptive grid is built for the conversion implemented by Local Level Set Method. Besides, to improve the of computing and data management, the refined nodes are stored by red-black tree.
This paper is organized as follows. Section II provides biological background on thrombus formation. Section III presents a detailed description of the proposed hybrid model. Section IV contains the demonstration and the experiment results of the numerical simulations. Section V gives a brief discussion and conclusion of this paper.

II. BIOLOGICAL AND MODELING BACKGROUND
The modern Virchow's triad summarizes that endothelial injury is one of the risk factors leading to thromboembolism along with stasis and hypercoagulability of blood [1]. In arteries, dangerous endothelial injuries can be caused by the disruption of unstable atherosclerotic plaques. Subsequently, the exposed collagenous fiber and released adenosine diphosphate (ADP), von Willebrand factor (vWF), and tissue factor (TF) provide the environment promoting thrombus formation [17].
For the arterial mural thrombus, the head of thrombus, also known as primary thrombus, is formed by platelets through two interrelated processes called adhesion and aggregation. Primarily, platelets adherer to the subendothelial tissue, and form a single cell layer by binding platelet receptors GpIb/V/IX to vWF, and the receptors GpIa/IIa and GpVI to collagen in the subendothelial matrix of the vascular wall. After that, the activated platelets form multicellular aggregates by the binding the platelet receptors GpIIb/IIIa to fibrin and fibrinogen [18], [19]. Above the primary thrombus, white cells and red cells are captured constantly forming a mixed thrombus. Hence, the primary thrombus is the key basic part of a mix thrombus. In the above process, the vWF which is a large multimeric glycoprotein plays an important role for platelets adhesion and activation under conditions of high shear. [20]- [22].
In the past, the mechanisms and interactions of vWF and hemodynamic conditions have been studied as distinct disciplines. Lkeda et al. [23] demonstrated that platelets aggregation occurred at shear stresses above 80 dyne/cm 2 with sufficient vWF in the blood collected with hirudin as anticoagulant. Their study explanted the reason why a stenosis with high shear stress could promote thrombosis. Westein et al. [24] artificially provoked a severe thrombosis in an injured blood vessel by creating a stenosis at the injury site with a blunted needle. The experiments shown that platelet aggregation was primarily driven by the shear rate of blood flow. Furthermore, they also studied the mechanism of thrombosis in a stenosis by using a microchannel with a stenotic shape to simulate pathological changes of vessel geometry. In this experiment, they chose completely anticoagulated blood to promote that the increased activity of vWF under high local flow rate is an inducement of thrombosis.

III. PROPOSED HYBRID MODEL
In this section, an overview of the baseline hybrid model is presented at first, and then the details of the thrombus growth sub-model and the implementation of the Cartesian adaptive grid are descripted.

A. BASELINE MODEL
The proposed hybrid computational model is aimed at the process of primary thrombus formation which mainly containing platelets adhesion and aggregation. In the model, the materials taken into account are platelet, vWF, plasma, vessel and thrombus. Figure 1 is the illustration of the model in three-dimensions. The vascular wall, plasma, vWF, and thrombus are described in continuum, while the platelet is in particle. The vessel is designed as a tube with a stenosis. There is a patch of endangium located in the stenosis is assumed to be injured. The plasma flows from inlet to outlet in the vessel, and there is vWF stimulated by the endangium injury. As platelets following with plasma, the platelets near the injury are activated by the vWF and the stuck platelets adhere on the vascular wall becoming a part of thrombus. The above interaction relationships between every two materials are summarized in Figure 2.   According to the essence of hybrid model, the baseline model is divided into continuum sub-model and particles sub-model.

1) CONTINUUM SUB-MODEL
The continuum field contains vascular wall, injured endangium, plasma, vWF, and thrombus. The computational domain is denoted by containing plasma and thrombus which denoted by p and t , respectively. The interface between these two materials is thrombus surface denoted by t , and vascular wall which is the boundary of plasma is denoted by v . The continuum filed is formulated by Level Set Method that will presented in the next sub-section detailedly.
The plasma, assumed as an incompressible Newton fluid, provides a hemodynamic environment for the model. The Navier-Stokes (N-S) equation is introduced as the motion function of plasma: Equation (1) is the momentum equation of the plasma under the action of the external force from systole, and (2) is the incompressible condition of plasma. In the equations, u is plasma velocity varying with time t, ρ is plasma density, p is pressure of plasma inside the vessel, a is the acceleration from systole, and ν is plasma viscosity. Since the vascular lumen is the fluid space of plasma, as the thrombus growing, the flow status of plasma changes with the space limiting. Since vWF is in microscale relative to other continuum materials, rather than their effects to blood flow, only the interactions of vWF and platelets that controlled by high shear rate are considered. The details of how vWF effects platelets are explanted in the particles sub-model.

2) PARTICLES SUB-MODEL
In the model, the primary force acting on platelets is the mobilization force F mo from plasma: The force is the integral of the difference between the velocity of plasma u plasma and the velocity of a platelet u platelet on the particle surface S p .
Depending on the process of white thrombus growth, platelets have two states: resting and activated. Under the action of vWF, platelets are directed toward to the injured endangium, and then the resting platelets are activated after adhesion and aggregation. According to the definition of shear rate, it can be calculated by the gradient of velocity, where u is the value of plasma velocity field in the position of the platelet. The shear rate is use to decide whether platelets are affected by vWF. When the shear rate is higher than the threshold γ act , the platelet moves toward the injured endangium under the adhesion force F ad and the aggregation force F ag formulated as follows: where Dis ad (Dis ag ) is the distance from the platelet to the closest point on the surface of the injured endangium (thrombus), and the unit vector to the point is n ad (n ag ). F ag_max and F ag_max are the maximal magnitude of the adhesion force and the aggregation force, respectively. Soon after the resting platelets adhering on the injured endangium or aggregating with platelets in the thrombus, they are activated. The activation of platelets makes the adhesion and aggregation become firm. It means the activated platelets become a part of the thrombus resulting in the thrombus growing. The implementation of the thrombus growth based on LSM is described in the next sub-section.

B. THROMBUS GROWTH MODEL 1) LEVEL SET METHOD
The level set method (LSM) [25] is an effective approach for interface tracking and shape modeling that uses a flexible implicit description. The main idea behinds LSM is to describe an interface ∈ R n by the zero-contour of a higher dimensional function φ. The closed curve in 2D is ∂φ ∂t + F · ∇φ = 0 (7) One of the advantages of LSM is that it can calculate the evolving curves and surfaces on Cartesian grid without parameterization. In the applications of level set method [26], [27], the adaptive mesh is usually introduced to control the time complexity of the algorithm and ensure the resolution of the contour and the accuracy of solving the level set equation.
Since the evolution process is continuous, an improved LSM named Local Level Set Method (LLSM) [28] is proposed. In LLSM, the evolution only focuses on the Gamma through defining a narrow strip around the zero level set named narrow band as shown in Figure 4. The local level set equation is defined as The truncation function C(φ) is where T 0 = {x : |φ(x)| < δ|} being a narrow band with the half-band width δ. By this way, the LLSM can effectively reduce invalid calculations.

2) RED-BLACK TREE IN CARTESIAN ADAPTIVE MESH
The level sets formulated computational domain of model is embedded in a Cartesian grid. Cartesian adaptive mesh refinement (AMR) is a method adept in simulating multiscale and multi-physics problems. Given the platelet size is microscopic relative to other continuum material, Cartesian AMR is chosen for thrombus growth model. For the simulation using AMR, cell storage is one of the key-points that affect the simulation performance. When choosing a data structure for AMR, there are several requirements containing memory usage, data locality, and data access.
Octree is the general data structure for cell storage of Cartesian AMR in 3D deficient in memory usage and data access. A storage strategy with Z-order curve and red-black tree is proposed in Jaber's study [29] satisfying the requirements better. Red-black tree is a special self-balancing binary search tree with an extra bit for color (red or black). In their study, the Morton code representation depending on Z-order curve is used as the comparison criteria for inserting new elements in the red-black tree. Figure 3 is an example of red-black tree in 2D case. On the basis of binary search tree, the neighboring data is stored closer to each other naturally. The significant advantage of red-black tree is the high efficiency of insert and delete operations of which the time complexities are O(log N ). Meanwhile, RBT is a binary tree which needs only 2 pointers for children, the storage efficiency is higher than Octree.

3) DRLSF BASED SPEED FUNCTION
The primary thrombus grows when platelets are stuck by adhesion or aggregation. The domain of thrombus is defined by level set function φ th setting thrombus surface as φ th = 0, thrombus region as φ th < 0, and other region as φ th > 0.
The distance regularization level set function (DRLSF) is introdued to drive the thrombus growth. And the thrombus surface expansion is represented by the level set evolution processing on the red-black tree improved Cartesian adaptive mesh.
DRLSF is a re-initialization-free level set evolution method by regularizing the level set function with property |∇φ| = 1, where a double-well potential function p(s) with minimum at both s = 1 and s = 0 are suggested. To provide a more stable level set function and fewer iterations in implementation, an improved double-well potential function p(s) [31] is adopted in our model as follows:  Firstly, to expand the thrombus to cover the stuck platelets, the level set function φ th is binary initialized as follows: where D is a positive constant, is the whole continuum domain, plasma is plasma domain, and ∂ plasma Then the level set function was evolved through a forward-andbackward diffusion of the above regularization term as follows: ∂φ th ∂t = div d p (|∇φ th |) ∇φ th (12) where d p is the diffusion function of p(s). By adding curvature speed to smooth the joint between thrombus and the stuck platelets, the level set evolution is reformulated as follows: where µ and ε were constants for regularization term and curvature speed term, respectively. κ was the curvature of the expanding surface calculated from the level set function.

4) THROMBUS LEVEL SET FUNCTION EVOLUTION
To implement thrombus growth, the thrombus region expands under the proposed speed function above. As well as the surface of thrombus is traced by LLSM. To solve the level set equation efficiently, the improved Cartesian adaptive mesh is applied to refining the grid hierarchically. First of all, the minimum and maximum coordinate values of the region containing injured endangium and thrombus is used to build a cuboid. Then the refined strategy is implemented on it. The steps of an expansion are: 1. Determine the half-width of narrow band w NB , and localize the computation on the Alive Cells; 2. Configure adaptive refined mesh for Alive Cells, and mark Land Mines; 3. Solve Level Set Equation on the adaptive mesh with one timestep, and update level set; 4. Adjust the refinement according to updated level sets; 5. If zero level set hits Land Mine, loop from step 1, otherwise loop from step 3.  For adaptive refined mesh configuration, as shown in the illustration (Figure 4), the Cartesian grid is refined in 3 levels. Given the level set equation of LLSM, the half-width of narrow band is assigned the larger width δ. In the coarse mesh (level-0, L0) the cells crossed by narrow band is tagged as Alive Cells and roughly refined (level-1, L1). In this L1 refined grid, the cells in the narrow band are refined further (level-2, L2), while other cells are marked as Land Mine. Then the L2 refined cells which is crossed by zero level set are refined further (level-3, L3). The refined mesh is stored by red-black tree as described above, and the evolution of level set function is computed on the refined mesh, along with the addition and deletion of the secondary refined cells. The function of Land Mines is to prevent zero level set moving over the boundary of narrow band. Once the Land Mines are crossed by zero level set, the narrow band should be re-built, and the adaptive mesh be re-configured. The general numerical solutions for level set equation are Upwind difference method, Hamilton-Jacobi WENO method, and TVD Runge-Kutta method. Whatever method is used, the cells on the boundary of different refinement levels, called ''hanging nodes'', should be handled specially, because they have no complete set of neighbors. Therefore, to evolve level set function in a non-uniform Cartesian grids, the partial differential of hanging node need to be approximated with ghost neighboring node. As illustrated in Figure 5, the hanging node v 0 has only 3 neighboring nodes v 1 , v 2 , and

IV. NUMERICAL SIMULATION AND RESULTS
The numerical simulations which are performed to study clot development under various vascular conditions is shown in this section.

A. SIMULATION EXPERIMENT DESIGN
The simulation experiments are primarily focused on the platelet thrombus formation in stenotic vessels. Figure 6 is the design drawing of the virtual vessel from two perspectives. The grid box for mesh configuring is in blue, the vascular wall in red, and injured vascular wall in yellow. For the size of grid box and vessel, L B , W B , and L V are the box length, box width, and vascular length, respectively. To simulate the stenosis result from the press of needle, referring to 2D Gaussian Function, the height of the stenosis is determined as where d is the distance of a point on vascular wall to the press point. A and σ are the parameters of Gaussian Function that can control the height and the width of stenosis respectively. So the vascular inner dimeter (D V ) is The diameter of the vascular lumen without stenosis is the maximum diameter (D Vmax ), while the minimum diameter (D Vmin ) is located at the peak of the stenosis. Here, the maximum of H S is used to represent stenosis degree. The area where H S is larger than the size of grid cell is considered as stenosis of which the width is W S . The whole stenotic region is assumed to be injured.   Table 1 shows the parameters of vessel and grid used in the simulations.

B. RESULTS
The simulation experiments were carried out to study two contents: first, the effect of stenosis to thrombus formation; second, the effect of blood flow velocity to the thrombus formation. Figure 7 is the visualization of the thrombosis simulations in vessels with varying degrees of stenosis. In the frames, vascular wall was in red, while the injured wall and the thrombus were in yellow. The blue spheres were the platelet particles. By comparing the frames under different stenoses at the same time, it could be observed that the   growth rate of clot was related to the severity of stenosis. Figure 8 is the growth curves of the clot in vessels with different stenosis in 0.3s. The growth curves were formed by the number of particles contained in thrombus in every millisecond. It was obvious that the clot grew faster when the stenosis became more serious. This phenomenon was consistent with the experiment results in Westein's study [24].
Another phenomenon could be observed in Figure 7 was that more particles was stuck in outlet region than inlet region. Figure 9 is the growth curves of inlet region and outlet region of the thrombus when the stenosis degree was 1.2mm. It was same as in /citeref24 that the thrombus grown faster at the outlet region than inlet region. Besides, as shown in Figure 10, with the increasing of stenosis degree, the difference of particle distribution between inlet and outlet went up first and then descended.
For the stenotic vessels, the thrombus formation was also affected by the blood flow velocity as in the normal vessels. Figure 11 (a) shows the average growth rate of the thrombus forming in normal vessel and stenotic vessel for each different blood velocity. The clot grew faster with the blood flowed faster both in normal vessel and stenotic vessel. This phenomenon coincided with the most recognized results of Begent's in vivo experiments done in the venules of Golden hamsters (Figure 11 (b)) [32].
To evaluate the performance of the proposed adaptive grid, we also organized the grid cell count of the adaptive refined grid, and compared with a uniformly refined grid. In the uniformly refined grid, all the L0 cells inside the narrow band was refined into L3. The simulations in which the stenosis degree is 1.2mm, and the blood flow velocity is 0.1mm/ms were implemented with adaptive refined grid and uniformly refined grid. In Figure 12, (a) shows the growth curve of cell count in adaptive refined grid, and (b) shows the growth curve in uniformly refined grid. In 0.3s of the simulations, the cell total count (T) of the adaptive refined grid was 37 tens of millions, while that was 8.6 hundred million in the uniformly refined grid which was 10 times of adaptive refined grid. As thrombus growing, the cell count of adaptive refined grid also increased, while the cell count of uniformly refined grid no longer increased after 60ms. But in the adaptive refined grid the increase of L0-L2 cells was almost negligible, and only the number of L3 cells increased obviously which makes more sense for the computation of Level Sets.
In Figure 13 shows the cell count of adaptive refined grid for each different stenosis degree at 0.3s of simulations. The increases of grid cell count were not significant for each refined level when the stenosis became serious. For the uniformly refined grid, the grid cell count increased with the degree of stenosis increasing.

V. DISCUSSION AND CONCLUSION
This work focuses on studying primary thrombus formation in stenotic arteries though building computational model for the formation process and using the model to implement simulation experiments.
In the proposed particle-continuum computational model, we take full advantage of Level Set Method in simulating stenotic vessels and keeping thrombus morphology, as well as solving the problems of clot growing and fluid-solid coupling. Besides, a novel adaptive mesh refinement strategy on Cartesian grid for Level Set Method is proposed to prevent high computational complexity and data overload. The main innovation of the proposed strategy is the using of red-black tree as the data structure to manage the refined mesh of a moving surface. By analyzing the cell count in the simulation progress, the efficiency of the adaptive grid in decreasing simulation data is proved.
Through the simulations of primary thrombus formation, we tested variety stenotic arteries. The experiments show that serious stenosis can intensify the formation of clot. On the stenosis, the aggregation of platelets is more serious in the outlet region than that in the inlet region. There is consistence in the influence of stenosis to clot growth rate and platelets distribution in thrombus between our simulation results and the results in similar studies [24].
We also test the clot growth under different blood flow velocities in the stenotic vessels, and results show the same trend as in normal vessels that clots grow faster at the beginning of blood flow velocity increasing, and the growth rate slows down as the velocity continue increases. The effect of blood flow velocity to thrombus growth is consistent with the in vivo experiments [32].
In conclusion, the present study developed an efficient hybrid model for thrombus formation in arteries. By using Level Set Method and adaptive mesh refinement, the clot growth process in the model become more efficient. The future work will focus on extending the variety of vascular morphology and improving the consistency between simulation experiments and reality experiments.