Numerical Study on Dynamics of Blood Cell Migration and Deformation in Atherosclerotic Vessels

: A phase ﬁeld model is used to study the effect of atherosclerotic plaque on hemodynamics. The migration of cells in blood ﬂows is described by a set of multiple phase ﬁeld equations, which incorporate elastic energies and the interacting effects of cells. Several simulations are carried out to reveal the inﬂuences of initial velocities of blood cells, cellular elasticity and block rates of hemodynamic vessels. The results show that the cell deformation increases with the growth of the initial active velocity and block rate but with the decrease of the cellular elasticity. The atherosclerotic plaque not only affects the deformation and migration of cells but also can promote the variation in hemodynamic properties. The atherosclerotic plaque causes a burst in cell velocity, and the greater the block rate and cellular elasticity, the more dramatic the variation of instantaneous velocity. The present work demonstrates that the phase ﬁeld method could be extended to reveal formation atherosclerosis at the microscopic level from the perspective of hemodynamics.


Introduction
Atherosclerosis is the leading cause of cardiovascular and cerebrovascular diseases, accounting for 31% of total deaths worldwide [1][2][3][4].Thus, the investigation into factors of the effects of atherosclerotic plaque on the flow of blood is imperative for the treatment of atherosclerosis.
Previous work [5,6] on atherosclerosis has demonstrated that the formation and rupture of atherosclerotic plaque caused by variation in hemodynamic properties is a crucially important inducement for thrombosis.However, whether the atherosclerotic plaque can inversely promote the variation on hemodynamic properties, what the undergoing factors are and how the factors play their role remain barely studied.At present, not only have the hemodynamics been extensively dissected at the macroscale [7][8][9], but also research on hemodynamics at the microscale for the mechanics of atherosclerosis formation [10][11][12] has made progress.However, it is difficult to decouple and change the characteristics of plaque and blood cells quantitatively by experimental methods.Additionally, the diagrams and trend patterns of morphological features of a cell without limitations on the time and lacation are only available in numerical simulation and cannot be measured by experimental methods.With the development of computers in recent decades, it is of great significance to search for efficient and powerful methods of numerical simulation to study the mechanism of blood cell dynamics.
Several numerical attempts have been carried out to study the cell dynamics under deformation and migration in recent decades.The models developed span the sub-cellular to the supracellular scales, including lattice models, continuum models, phase-field models, active network models and particle models [13].The lattice model that is used in [14][15][16], known as the cellular Potts model, describes individual cells as domains on a lattice.Particle models, which are used in [17,18], treat each cell as one or two circular particles so that they can capture cell shape anisotropy.The phase-field model, originating in interface dynamics [19], describes each cell by a phase field ϕ.Jakob et al. [20] described the movement of many self-organized, interacting cells by presenting a phase-field model that took into account the main mechanisms of cell motility acto-myosin dynamics, as well as substratemediated and cell-cell adhesion.Akiyama et al. [21] presented a model of multicellular systems with several types of cells, based on the criteria for minimizing the surface area and retaining a certain volume.The phase-field model improved by Palmieri et al. [22] is particularly suited to describing cell dynamics.The lack of a need to track the boundary explicitly, the convenience in describing the extremely large deformations and modeling the mechanical properties and the velocities of each cell individually are the competitive advantages of this approach compared with other cell models.Based on the former, Yin et al. [23] further explored the effect of different elastic mismatches on the migration dynamics of multiple cells in a monolayer of cancer cells in more precise detail.Due to the superiority of this model, it is used in the present work, as the previous two works.Here, the model is applied concretely to simulate the influence of atherosclerosis-a specific disease.There is precedent for using numerical methods to study atherosclerosis.The first calculation of the blood flow in aortic arch geometries in the presence of plaques was approached by Assemat et al. [24].In addition, Shit et al. [25] investigated the pulsatile flow characteristics of blood when it passes through a porous overlapping stenosed artery.Much of the research on the effects of atherosclerotic plaque on blood flow is hydrodynamic-based; some other approaches include the numerical study of the macroscopic physiological(e.g., [26]) or therapeutic processes(e.g., [27]) of atherosclerotic plaque.However, very few studies have specifically sought to achieve the effect of atherosclerotic plaques on the blood cells.The work of Lázaro et al. [28,29] is worthy of reference in the numerical study of the deformation of blood cells.Here, the migration and deformation of blood cells are explored with three influencing factors of the atherosclerosis simultaneously.
The present work aims at revealing the effects of atherosclerotic plaque on hemodynamic variation within the phase-field method.Firstly, the deformation and migration of blood cells in the atherosclerotic vessel is mathematically described by using a multiple phase-field equation.After that, the phase-field model is used to investigate the effects of the cellular initial active velocity, cellular elasticity and vessel block rate on the dynamics of blood cell migration and deformation, and the results are then presented.

Model Description
The deformation and migration of blood cells in the atherosclerotic vessel is mathematically described by using the multiple phase field model [22].In the model, each cell is treated as a soft body and considered as a separate phase field.The time-dependent behavior of a monolayer cell is governed by the phase field equation where ϕ i (x, y, t) is the order parameter of a cell at position (x, y) and time t, v i represents the instantaneous migration velocity of the i-th cell, F(x, y, t) stands for the total free energy, and the subscript i denotes the cell number.Here, ϕ i (x, y, t) describes the inner and outer region of the cell i.When ϕ i (x, y, t) = 1, the position (x, y) is occupied by the cell i.It is outside the i-th cell at (x, y) when ϕ i (x, y, t) = 0, and ϕ i (x, y, t) ∈ (0, 1) indicates the position (x, y) is covered by the membrane boundary of the i-th cell.
In the present model, the total free energy consists of two parts: the free energy for an individual cell and the free energy F o caused by the interaction between each two cells F int , i.e., F = F o + F int .These two free energies can be respectively expressed as and where γ is the elastic coefficient of a cell, λ is the width coefficient of the cell boundary, µ represents the volumetric coefficient of a cell, R is the cell radius, and κ is the wave number, which controls the energy cost owing to the overlap of the cells i and j.In the model, a cell tends to have a circular shape with a radius R when it has no interaction with other cells.
The velocity v i is composed of two parts: i.e., v i = v i,A + v i,I .Here, v i,A is the active part of the velocity, which is triggered by the internal and external forces on each individual cell.v i,I is the inactive part, which is caused by the interaction between a cell and its neighbors.In the present model, cells are considered to be living, and thus v i,A is set as a constant in all the following simulations.The inactive velocity v i,I is computed by where ξ is the friction coefficient between every two cells.The parameters above are dimensionless, and the actual parameters need to be dimensionless before operation.Substituting Equation (4) into Equation ( 1) with the expansion of the variation δF(x, y, t)/δϕ i (x, y, t), the following dimensional equation for cell dynamics can be derived: This is a time-dependent equation that is solved by using the explicit finite difference method on a fixed grid and time stepped with a simple Euler scheme.To verify the model and code, a case with identical conditions to Yin's study [23] is carried out.Figure 1 shows the average shape change of a tagged cell over the average shape change of a normal cell.It can be seen that the points of the present data are distributed with the same trend as the reference data, within the margin of error, which proves that the model and code used here are basically correct.The data of the present work can be described as a function of

Blood Cells Deform and Gather in an Atherosclerotic Vessel
The cell migration in an atherosclerotic vessel is simulated to investigate the deformation and accumulation of cells.Initially, 80 circular cells with a uniform radius R and an initial active velocity are randomly placed in the lumen of vessel.Here, the initial active velocity is labeled as v 0 .As the vascular calcification is a hallmark of atherosclerosis, the atherosclerotic vessel is much stiffer than the normal vessel [30].Thus, the vessel is considered here as a rigid body with a constant order parameter ϕ = 0; while non-physical, this approximation has little influence on the simulation results.The configuration of the vessel is set as a straight channel of L × H in a Cartesian coordinate with an atherosclerotic plaque.The plaque is described by using the function Here, R b denotes the block rate caused by atherosclerotic plaque and is defined as the ratio of the minimum diameter to the maximum diameter (normal vascular segment) of the vessel, L is the length of the vessel, H is the height of the vessel, h v denotes the thickness of normal blood vessels shown in the simulation domain, x α determines the lateral position of the atherosclerotic plaque, and σ controls the length of vessels in atherosclerotic segment.These parameters in the function are also dimensionless and assigned as σ = 28, 800, x α = 2L/3 and h v = 10.The other parameters are set as κ = 60, λ = 7.0, ξ = 1500, µ = 40, and R = 12 (i.e., listed as Table 1).These parameters are obtained from [31] or by invoking physical approximations.In the present simulation, the channel is divided into 800 × 200 mesh grids, and periodic boundary conditions are imposed on the right and left sides of the vessels.The up and down side are considered as the solid wall boundary, as shown in Figure 2.  Figure 3 shows the evolution of the cells-vessel system with the dimensionless initial active velocity of the cell v 0 = 5, the block rate R b = 50% and the dimensionless cellular elasticity γ = 10.The cells are randomly distributed in the lumen at the beginning, and the distance between the centers of two cells is not less than 1.6R to reduce the degree of cell overlap.To analyze the migration and dynamics of an individual cell, a randomly chosen cell is tagged in dark-red color.The gray on either side is the vascular walls with protruding atherosclerosis.At the initial stage, the blood cells are almost round with little deformation and gathering, as shown in Figure 3a.Later, the tagged cell moves near to the narrowest position of the vessel, and it deforms significantly due to extrusion, as shown in Figure 3b.After that, the tagged cell almost restores its original shape, as displayed in Figure 3c.On the whole, the cells migrate simultaneously with the blood flow from the left side of the vessel to right; the atherosclerotic plaque leads to the accumulation and deformation at the narrowest position.This demonstrates that the multiple phase-field model can be used to model the dynamic behaviors of cells in a vessel with an atherosclerotic plaque.

Block Rate Effect on the Dynamics of Blood Cells
Vascular block rates would affect the migration and deformation of cells in vessels, which would imperatively impact the causes of vascular diseases.Several simulations are carried out to reveal the effects of vascular block rates on the dynamic behaviors of blood cells.In these simulations, the block rates of R b = 25%, R b = 50% and R b = 75% are selected for mild, moderate and severe stenosis, respectively.The velocity and cellular elasticity are set to be v 0 = 5 and γ = 10.Other parameters are chosen the same as those in Figure 3. Figure 4 shows that the blood cells rapidly pass through the atherosclerotic plaque with little accumulation and deformation at R b = 25%.However, with the block rate increase, more and more cells gather and deform due to the extrusion with the neighbor cells and the atherosclerotic plaque.
A dimensionless length of the cell perimeter is proposed to quantitatively characterize the cell deformation, which is expressed as It is clear that L i ∼ 1 depicts little deformation of a cell.This indicates that the shape of cell i is close to its initial shape.On the contrary, it indicates the deformation of cell i when L i > 1.Here, the L i of the tagged cell is labeled as L 0 , while the average of the dimensionless length of all cells is marked as L 1 , i.e., where N is the amount of the blood cells in the simulations.In the following simulations, L 0 and L 1 are both used to analyze the deformations of migrating cells.Figure 5 shows the evolutions of the dimensionless perimeters of the tagged cells L 0 s and the average of the dimensionless length of all cells L 1 s for the simulations of block rates R b = 25%, R b = 50% and R b = 75% with γ = 10.From Figure 5a, it can be observed that L 0 for R b = 75% rises sharply around t = 25 s, reaching the peak at t = 30 s.The maximum value of L 0 for R b = 75% exceeds 1.5.This means that the perimeter of the tagged cell increases by 50%.When R b = 25%, L 0 remains nearly constant, indicating that the tagged cell deforms barely.By combining with Figure 4, it can be found that the tagged cell just approaches the atherosclerotic plaque; meanwhile, the L 0 starts increasing at t = 10 s.Therefore, it is the atherosclerotic plaque that results in the deformation and the blockage of blood cells.While R b = 75%, L 2 tends to increase to a high value before t = 50 and then decreases.The increasing value of L 1 to its peak for R b = 75% is much larger than the values for R b = 50% and R b = 25%.The reason is that the atherosclerotic plaque causes the deformation of blood cells.The larger the block rate, the greater the cell deformation.This also demonstrates that L 1 can be considered as an indicator to describe the deformations of accumulated cells.
Furthermore, the cells are considered to have elliptic shapes with a major axis and a minor axis.The length ratio of these two axes is proposed to characterize the deformations of cells.Here, L 2 is used to represent the length ratio of these two axes of the tagged cell i. Figure 6 displays the maximum and averaged values of the dimensionless perimeters of all cells L 1,max , L 1,ave and of the ratios L 2,max and L 2,ave in the simulations of v 0 = 5, γ = 10 and various block rates R b = 25%, R b = 50% and R b = 75%.It can obviously be found from Figure 6 that all the four indicators increase simultaneously with the increase of the R b .The increased amounts of L 1 and L 2 between every two adjacent R b s are less at a lower R b than those at a higher R b .Therefore, the length ratio of the major and minor axes can also be selected as an indicator to characterize the cell deformation.Figure 7 plots the evolution of the instantaneous velocity of cell i over time for the above simulations.It shows that v i s varies more dramatically at R b = 75% than at R b = 25% and R b = 50%.It is obvious that v i almost remains constant throughout the simulation of R b = 25%.That is because the atherosclerotic plaque results in the variation of instantaneous velocity of blood cells.The more blocked the vessel, the more dramatically the cell velocity fluctuates.

Cellular Elastic Effect on Dynamics of Blood Cells
Cellular elasticity is another factor that influences the passing behaviors of blood cells crossing the atherosclerotic plaque.The following simulations are performed to study the elastic effects of cells.In the simulations, v 0 = 10 and R b = 50%.The elastic coefficients are set as γ = 5, γ = 10 and γ = 20, respectively.Other parameters are set the same as those in Figure 3.
Figure 8 displays snapshots of blood cells and the atherosclerotic vessel for the simulations of v 0 = 10 and R b = 50% at t = 10 with different cellular elasticities.It can be found that the deformation of cells in Figure 8a is more notable than that in Figure 8b,c.The deformation is caused by the extrusion force when the cells pass the atherosclerotic plaque.With the increase of the elasticity of cells, the deformation of cells becomes difficult when cells move across the atherosclerotic plaque.Figure 9 shows the time evolution of L 0 and L 1 for simulations of various elasticities.It can be observed from Figure 9a that the up and down trends of three curves of L 0 versus time t are almost the same.The increasing and decreasing rates of γ = 5 are obviously greater than those in the other two cases.The increase of L 0 is inversely related to γ.The same feature is also shown in Figure 9b, which depicts the evolution of dimensionless perimeter L 1 for the simulations.Compared with the case of higher elasticity, the cellular deformation with lower elasticity has more volatility.It can be found that all the L 1 s and L 2 s decrease with the increase of γ.The maximum and averaged values of L 1 and L 2 are negatively correlated with γ, which means that the shape change of the tagged cell decreases with the growth of the cellular elasticity.This proves that the general deformation of cell populations decreases with the increase of the cellular elasticity.It also reveals that the cell deformation decreases with increasing cellular elasticity.Therefore, cellular elasticity has a great influence on cell deformation but not obviously on cell velocity.The stiffer the blood cells, the more they deform.Figure 11 displays the evolution of instantaneous velocity.In the simulations, there are two bursts in the instantaneous velocity of the tagged cell.Both bursts are caused by the atherosclerotic plaque, which indicates the accumulation, compression and deformation of the cells.The difference in the height of the two peaks is due to the change in cell position when the tagged cell repasses the atheromatous plaque.It can also be found that the overall deformation of blood cells is increased when decreasing the cellular elasticity.Diseases such as coronary heart disease, diabetes, hypertension and hyperlipidemia may lead to the increase of elasticity of blood cells, which means a decreased cell deformation ability and that it will be difficult to pass the blocking plaques smoothly.

Active Velocity Effect on Dynamics of Blood Cells
The active velocity is another important factor that can influence the dynamics of blood cells during their migration crossing the atherosclerotic plaque.Several simulations are performed at various initial active velocities: v 0 = 5, v 0 = 10 and v 0 = 20.Here, the cellular elasticity and block rate are set to be γ = 10 and R b = 50%.Other parameters are chosen the same as those in Section 3.1.
Figure 12 shows the simulated multiple blood cells gathering in the atherosclerotic vessel at various initial active velocities.Despite the different velocities, cell accumulation appears when passing through the atherosclerotic plaque.At a lower velocity, the accumulation is less significant than that at a higher velocity.That is because the accumulation can be attributed to both the obstruction of atherosclerotic plaque and the deformation of cells.A large deformation indicates large extrusion force between two neighboring cells.For example, the tagged cell i migrates with a slow velocity, as shown in Figure 12a.Less formation occurs in such a condition.In the condition of a high velocity, as shown in Figure 12c, a large deformation occurs, and stronger interactions give rise to significant accumulation.Therefore, the initial active velocity can change the accumulation of a single cell and the group behavior of cells.Figure 13 displays the evolution of L 0 s and L 1 s for the simulations of various initial active velocities.As Figure 13a shows, despite the different velocities, peak values of L 0 s could be found in all the simulations.The peak value of L 0 at v 0 = 20 appears earliest, and that at v 0 = 5 comes latest.The peak value of L 0 S at v 0 = 20 is largest in the three cases, and there are several steep climbs in the condition of the largest initial active velocity.This implies that a large migrating velocity can result in a large deformation for a cell.When increasing the initial active velocity, the deformation and accumulation would be enhanced.Figures 13b shows that the L 1 goes up and down with the migration of cells.When increasing the initial active velocity, L 1 can be enhanced.Figure 14 shows the maximum and averaged values of the dimensionless perimeters of L 1 and of the length ratios of L 2 for the simulations with R b = 50%, γ = 10 and various initial active velocities.Here, the velocities are set as v 0 = 5, v 0 =7.5, v 0 = 10, v 0 = 12.5, v 0 = 15, v 0 = 17.5 and v 0 = 20.Generally, the average and maximum of L 2 is enhanced when increasing v 0 .The increasing trends of L 1,ave are similar to those of L 2,ave and L 2,max .However, when increasing the initial active velocity, L 1,max can be increased more obviously, especially in the condition of a high value of v 0 .The difference in the trends between L 1 and L 2 may be caused by the behaviors of a cell and its surrounding cells.The shape change of the cell can be generally increased with the increase of the initial active velocity, but such an effect is less obviously observed than that observed in the condition of a group of cells.The interactions of cells and the group behavior can affect and enlarge the deformations of cells.Figure 15 shows the time evolution of migrating velocities for cases of v 0 = 5, v 0 = 10 and v 0 = 20.It can be found that the migrating velocity burst more than three times during the simulation in the case of v 0 = 20, whereas the velocity varies gently and with low frequency in contrast for v 0 = 5.This implies that the deformation of a cell is increased with the increase of the initial active velocity.The higher the value of v 0 , the more obvious the effect is.According to Bernoulli's principle, where the pressure is higher, the blood flows faster, meaning that the blood pressure and pulsatility will be reflected in the blood flow rate, with the initial active velocity of the blood cells at a microscopic level.

Conclusions
The blood cell migration in atherosclerotic vessels was numerical studied by using an extended phase field model.The multiple phase-field equations were used in the model, which include free energies and interacting energies of cells.Simulations were performed with various cell initial velocities, cell elasticities and block rates of vessels, respectively.A dimensionless length and axis length ratios are proposed to characterize quantitatively the deformation of cells.The results reveal that atherosclerotic plaque leads to cell accumulation and deformation, and the initial active velocity of blood cells, the cell elasticity and the block rate could significantly change the dynamical behaviors of blood cells.The deformation and accumulation grow with the increase of the block rate, and meanwhile, the instantaneous velocity varies more drastically.The cell deformation is positively related to the initial active velocity and block rate but negatively related to the elasticity.In addition, the greater the block rate and cellular elasticity, the more dramatic the variation of instantaneous velocity.By using the explicit finite difference method and local parallelism calculation, the program operates efficiently.The new proposed dimensionless length and axis length ratios are useful to describe the behavior of migrating and deforming cells.The present work studied the function of atherosclerotic plaque on individual cell deformation and cell population migration and accumulation and explored the effects of three factors on these effects, providing a microcosmic perspective of the causes and consequences of atherosclerosis and a numerical hemodynamic approach for the research of atherosclerosis and its complications.This will be extended to further studies with the consideration of extracellular fluid flows.

Figure 1 .
Figure 1.Plot of the average shape change of a tagged cell over the average change of a normal cell for different values of γ tag /γ normal .The present work uses |v a | = 0.3, compared with Yin's work of three different active velocity |v a | values.The fit curve of the present work is a power law function.

Figure 2 .
Figure 2. Diagram of the simulation, where the blood cells migrate to the right in the blood vessel with atherosclerotic plaque, and the periodic boundary conditions and solid wall conditions are imposed at different boundaries.

Figure 3 .
Figure 3. Snapshots of multiple blood cells and atherosclerotic vessel over time for v 0 = 5, R b = 50% and γ = 10 at (a) t = 0.5s, (b) t = 20s, (c) t = 40s, respectively.The cell filled in red is the tagged one, and the gray areas represent blood vessel walls and the atherosclerotic plaque.

Figure 5 .
Figure 5. Evolutions of (a) dimensionless perimeters of the tagged cell L 0 s and (b) the average of the dimensionless length of all cells L 1 s for the simulations of block rates R b = 25%, R b = 50%, R b = 75% with v 0 = 5 and γ = 10.

Figure
Figure 5b displays time histories of L 1 s in the simulations of v 0 = 5, γ = 10 and block rates of R b = 25%, R b = 50%, R b = 75%.It can be found that L 1 is around its initial value

Figure 6 .
Figure 6.Variation of (a) maximum and averaged dimensionless perimeters of all cells and the (b) maximum and averaged length ratios of the major and minor axes of the tagged cell against block rate R b in the conditions of v 0 = 5, γ = 10.

Figure 7 .
Figure 7. Evolution of instantaneous velocity v i of the tagged cell for the simulations with v 0 = 5, γ = 10 and different block rates R b = 25%, R b = 50%, R b = 75%.

Figure 9 .
Figure 9. Evaluations of (a) the dimensionless perimeters of the tagged cell L 0 s and (b) the average of the dimensionless length of all cells L 1 s for the simulations of cellular elasticity γ = 5, γ = 10 and γ = 20 with v 0 = 10 and R b = 50%.

Figure 10
Figure10depicts the maximum and averaged values of the dimensionless perimeters L 1,max and L 1,ave and the maximum and averaged values of the ratios L 2,max and L 2,ave in the simulations of v 0 = 5 and R b = 50% and cellular elasticities γ = 5, γ = 10 and γ = 20.It can be found that all the L 1 s and L 2 s decrease with the increase of γ.The maximum

Figure 10 .
Figure 10.Evolution of the average and maximum of (a) the dimensionless lengths of the perimeter of all cells (L 1,ave , L 1,max ) and (b) the length ratios of the tagged cell i (L 2,ave , L 2,max ) as functions of γ in the condition of v 0 = 10 and R b = 50%.

Figure 13 .
Figure 13.Evolution of (a) the dimensionless perimeters of the tagged cell L 0 s and (b) the average of dimensionless length of all cells L 1 s for the simulations of various initial active velocities v 0 = 5, v 0 = 10, v 0 = 20 with R b = 50%, γ = 10.

Figure 14 .
Figure 14.Evolution of the average and maximum of (a) the dimensionless perimeters of all cells (L 1,ave , L 1,ave ) and (b) the length ratios of tagged cell (L 2,ave , L 2,ave ) as functions of the initial active velocity v 0 with γ = 10 and R b = 50%.

Table 1 .
Numerical values used for the dimensionless parameters in the model.