# A Novel Dynamical Analysis Method for Platelet Aggregation: From Traffic Flow to Blood Flow

Authors: Min Qiu; Jiahe Ren; Zhipeng Xu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This study proposed an ordinary differential equation model to simulate platelet aggregation by considering the Stokes viscosity resistance and interaction force between each platelet. This model is similar to the car-following traffic model, and its stability condition can be derived via long-wave expansion. The phase diagram can be divided into stable, metastable, and unstable regions. The theoretical results show that the flow viscosity and platelet diameter play essential roles in thrombus formation. The mKdV equation near the critical point is derived to describe the platelet density waves’ evolution properties by applying the reductive perturbation method. Furthermore, through the simulation of the space-time evolution of the distance distributions for platelets, it is shown that a smaller platelet diameter and lower viscosity will lead to thrombus formation, and the analytical results are consistent with the simulation results. The model indicates that thrombosis is similar to traffic congestion.


---

## Body

## 1. Introduction

Platelets promote wound closure and prevent excessive blood loss when the skin or tissue is injured, but they are also the main factor affecting thrombus formation. If the receptor proteins on a platelet’s surface are activated by platelet-activating factors such as thrombin or ADP, the platelet will firmly adhere to a blood vessel and become the base of a thrombus. Subsequently, a large number of platelets will attract each other and construct a filamentous protein network through platelet-platelet interaction via GPIba-filament protein. At the same time, the platelet membrane will deform and exhibit filopodia, which raises the probability of interconnection due to increased surface area. Abnormal platelet function may lead to bleeding disorders or thrombosis, which can be life threatening. This condition affects more than 2% of the population in the United States [1].Many studies have carried out numerical simulations of platelets to supplement experimental and clinical results. These simulations typically include processes in the macroscopic scale, coagulation cascade scale, cell scale, and microscopic scale.The macroscopic scale is equal to the hydrodynamics represented by the commonly recognized blood flow. The characteristic length ranges from 10 to 500μm, including both capillaries and the aorta, where fluid equations are applied. For example, Poziridis [2] used the boundary element method to study the behavior of platelets at the boundary of the blood vessel wall. Wu et al. [3] investigated the flipping of platelets close to the blood vessel via the lattice Boltzmann method. Moreover, capillaries are a significant part of the blood-brain barrier and connect the smallest arteries to the smallest veins. Hence, nanofluid flow through capillary attracts much attention [4–11].Platelet aggregation corresponds to the coagulation cascade scale. The characteristic length of this scale is near 10μm, so this scale is capable of showing how platelets interact with filamentous proteins in the blood to construct a network structure while enveloping red blood cells and white blood cells and promoting more blood factors to participate in the reaction, finally producing a thrombus. Most studies simplify each platelet as one mass point. Yazdani et al. [12] used the mass point model to simulate thrombus formation and the interactions of platelets with red blood cells under the Morse potential. The phase-field method, percolation theory, and reaction-diffusion equations have also been used to analyze thrombus formation.The characteristic length of the cell scale is near 1μm, which is close to the size of platelets. This scale focuses on the interaction between platelets and the deformation effect. In one study, the cell scale was used to explore the interactions between platelets, red blood cells, and white blood cells, and it was found that most platelets were distributed near the blood vessel wall [13].Furthermore, the microscopic scale can be used to investigate receptor and membrane reactions. Many receptor proteins are distributed on the platelet membrane at the molecular level, and biochemical reactions involving these proteins are the fundamental causes of platelet aggregation. Regarding this scale, due to the limitation of the calculation speed of all-atom molecular dynamics, the coarse-grained molecular dynamics (CGMD) method is normally used for research. For example, Zhang [14] conducted a comparative study on platelet surface elasticity and found that the rigid-body model’s surface pressure was higher than the actual results. Prachi [15] proposed the modified Morse potential, which is calibrated with in vitro experiments, to simulate platelet aggregation. Recently, Zhang et al. [16] proposed a new model based on the results of scanning electron microscopy and established a corresponding model to simulate filopodia. Figure 1 and Table 1 show some typical methods and references that are applied to platelet aggregation, and it is easy to observe that the thrombus simulation involves biomedicine, physics, and numerical methods.Figure 1
The summary of the systematic structure of different methods, scales of platelet aggregation, and thrombus development.Table 1
Example references for each method with their characteristic length scale.
MethodScaleReferencesNavier–Stokes equationFluid dynamics[17–19]Lattice Boltzmann equationFluid dynamics[3]Molecular dynamicsFluid dynamics[14, 15]Navier–Stokes equationCoagulation cascade[19]Lattice Boltzmann equationCoagulation cascade[20]Advection-diffusion equationCoagulation cascade[21]Phase-field methodCoagulation cascade[22, 23]Molecular dynamicsCoagulation cascade[24, 25]Constitutive equationCell[26]Elasticity equationCell[27]Molecular dynamicsCell[15]Molecular dynamicsReceptor reaction[15]According to the above analysis, most numerical simulations for platelets with molecular dynamics require a large amount of calculation and complex integration of multiscale information. Hence, the objective of this study was to develop a simplified model to analyze the thrombosis phenomenon based on the traffic flow model, since traffic congestion and thrombosis have many similar features. For example, a cluster of low-speed cars will lead to traffic congestion; similarly, a thrombus will block blood flow.This paper is organized as follows. Section2 introduces Bando’s car-following model and our model considering Stokes’ law and different interaction potential functions. Section 3 presents the stability analysis for this model. Section 4 describes the reductive perturbation method used to derive the nonlinear equation. Section 5 details the numerical simulation performed to validate the theoretical analysis. Section 6 offers the conclusions of this study.


## 2. Model

In 1998, Bando et al. proposed the optimal velocity model based on ordinary differential equations, which assumed that each driver would set their optimal speed by referring to their distance from the vehicle ahead. The equation of this model is as follows:(1)d2xntdt2=aVΔxnt−dxntdt,where xnt is the position of the nth vehicle at time t and Δxnt=xn+1t−xnt is the distance between the nth vehicle and the n+1th vehicle. In addition, the sensitivity coefficient α reflects the driver’s reaction time, and the optimal velocity VΔxnt is the driver’s expected velocity. Generally, the optimal velocity is related to the distance between the preceding vehicle and can be written as(2)VΔxnt=vmax2tanhΔxnt−hc+tanhhc.As the distance between adjacent vehicles increases, the expected speed increases monotonously with space, but it remains almost straight when the velocity approaches the maximum speedvmax. Therefore, many scholars [28–49] have extended the model to cover scenarios such as intelligent transportation, multiple lanes, and the honk effect.Here, the optimal velocity function shows the interaction behavior between two cars, which inspires us to convert the optimal velocity function into the interaction force function to analyze the movements of platelets by using Newton’s second law. Then, thrombosis can be treated as traffic congestion because the way in which a thrombus blocks a blood vessel is similar to the way in which a traffic jam blocks a road.In fluid mechanics, the force formula of an object in fluid (Stokes’ formula of viscous resistance [50]) is(3)f=6πμRv,where μ is the viscosity coefficient and v is the relative velocity. The shape of an unactivated platelet is an ellipsoid, whereas that of an activated platelet is a sphere with filopodium. R represents Stokes’ radius, which means the particle with radius R has the same dynamics parameters as an irregularly shaped platelet. Hence, platelets can be regarded as fine solids, and Stokes’ law can describe the movement of fine solids in fluid. However, the movement of platelets is different from general free particles in that there is a specific interaction force between platelets. Therefore, after considering the interaction force FΔxn, we propose the one-dimensional platelet motion model equation as follows:(4)dvntdt=6πμRmvf−dxndt−FΔxnm,where FΔx is the interaction force between platelets and m is the mass of one platelet. Yazdani et al. [12] used the Morse potential energy function to study the aggregation behavior of platelets. Therefore, we will also use the Morse force as follows:(5)FΔx=2Deβe−2βdΔx/d−1−e−βdΔx/d−1,where De is the well depth, β is the width of energy, and d is the equilibrium distance. For convenience, we set De=1,d=1,β=1 for (5), and the blue curve in Figure 2 represents the Morse force function. By expanding the Morse potential energy function near x=0, one can obtain the standard elastic potential energy, shown by the red dashed line. Meanwhile, we use the hyperbolic tangent format, shown by the black dotted line, corresponding to (2) by approaching the Morse force. In Figure 2, the upper area of each line indicates the repulsive force, and the lower region indicates the attractive force.Figure 2
The interaction force function ofFΔx.


## 3. Linear Stability Analysis

In the ordinary differential equation dynamical system, linear stability analysis has been widely applied to investigate the perturbation response. In other words, stability refers to free flow, which means each platelet moves individually without the cluster. Instability means there exists a thrombus blocking blood flow. In this study, we use the same linear stability analysis as in the classical car-following model to investigate the impact of the evolution of thrombus on blood flow. To do so, we consider blood flow during the stable state of platelet movement. In uniform flow, all the platelets are moving through a blood vessel with the same velocity. The steady-state solution for this equation is as follows:(6)xn0t=bn+vf−FΔxn6πμR,b=LN,where L is the length of the blood vessel and N is the number of platelets. In the steady state, the interactive force between each platelet is zero, and all platelets move with the same velocity.By introducing a small deviation into the steady-state solution, the uniform solution can be written as(7)xnt=xn0t+ynt.Then, we can obtain the linearized equation by substituting (6) and (7) into (4):(8)d2yntdt2+6πμRmF′bΔynt6πμR+dyntdt=0,where Δynt=yn+1t−ynt and F′b=dFΔxn|dΔxn|Δxn=b. For convenience, some coefficients in (8) can be replaced with c1=6πμR|m and c2=1/6πμR. By setting ynt=expikn+zt for the Fourier models, we can obtain(9)z2+c1c2F′b+c1zeik−1=0,and then we can obtain the first-order and second-order terms of ik by substituting z=z1ik+z2ik2+⋯ into (9) as follows:(10)z1=−c2F′b,z2=−c1c2F′b+2z122c1.Whenz2>0, the evolution of uniform steady-state flow still remains stable for long-wavelength modes. In the case of z2<0, the flow will become unstable. Finally, we can obtain the neutral stability condition for this simplified dynamical model of platelets as follows:(11)c1=−2F′bc2,and(12)36π2μ2Rm=−2F′b.To maintain that there exists no thrombosis in a blood vessel, the parameters should satisfy(13)36π2μ2Rm>−2F′b.Here, it can be seen clearly that decreasing viscosity leads to thrombosis. Increasing density with smallerR/m also promotes platelet adhesion, which means thrombosis may occur. These findings are consistent with the findings in [51].Figure3 shows the stable region of α for different interaction force functions. The lower part below the solid line corresponds to the unstable region, and the upper part over the solid line is the stable region. The region below the red dot-dash line is the infeasible part because α<0 has no physical meaning. The space between the coexistence curve and the neutral stability curve is the metastable region mixed with free particles and thrombi.Figure 3
Phase diagram in theΔx−α space.


## 4. Nonlinear Analysis

The nonlinear behavior near the critical point of the system is analyzed, especially for the functionFΔx=−2tanhΔx−1. By introducing a small scaling parameter, the gap between each platelet is set as follows:(14)X=εn+btandT=ε3t,0<ε≤1,Δxnt=hc+εRX,T,where hc is the mean distance between each platelet. For convenience, with c1=6πμR/m and c2=1/6πμR, (4) can also be rewritten as(15)d2△xntdt2=c1c2FΔxn−c2FΔxn+1−d△xntdt.By expanding each term in equation (16) to the fifth order of ε, we obtain(16)d△xndt=ε2b∂XR+ε4∂TR,(17)d2△xndt2=ε3b2∂X2R+2ε5b∂X∂TR,(18)Δxn+1t=hc+εR+ε2∂XR+12ε3∂X2R+16ε4∂X3R+124ε5∂X4R,(19)FΔxn=Fhc+F′Δxn−hc+16F‴Δxn−hc3,(20)FΔxn+1=Fhc+εF′R+ε2F′∂XR+ε3F′2∂X2R+F′′′6R3+ε4F′6∂X3R+F′′′6∂XR3+ε5F′24∂X4R+F′′′12∂X2R3,where(21)F′=F′hc=dFΔxndΔxn|Δxn=hc,F‴=F‴hc=d3FΔxndΔxn3|Δxn=hc.By substituting (17)–(21) into (15) and expanding to the fifth order, the partial differential equation can be obtained as follows:(22)ε2c1b+c1c2F′∂XR+ε3b2+c1c2F′2∂X2R+ε4c1∂TR+c1c2F′6∂X3R+c1c2F‴6∂XR3+ε52b∂X∂TR+c1c2V′24∂X4R+c1c2F‴12∂X2R3=0.At the critical pointΔxc,ac, let ac=a1+ε2 and take b=−c2F′. Then, the second-order and third-order terms of ε in (22) can be eliminated and written as(23)ε4∂TR−g1∂X3R+g2∂XR3+ε5g3∂X2R+g4∂X4R+g5∂X2R3=0,where(24)g1=−c2F′6,g2=c2F‴6,g3=−c2F′2,g4=−c2F′23ac+c2F′24,g5=c22F′F‴3ac+c2F‴12.To obtain the standard mKdV equation, we introduce the following transformation:(25)T=T′g1,RX,T=g1g2R′X,T′,and then we obtain the standard mKdV equation as(26)∂T′R′−∂X3R′+∂XR′3+εMR′=0,where(27)R0′X,T′=ctanhc2X−cT′.To determine the wave velocity in (27), we use the following solvability condition:(28)R0′,MR0′≡∫−∞+∞dXR0′MR0′=0.By integrating (28), we obtain the wave velocity c:(29)c=5g2g32g2g4−3g1g5.Therefore, the standard kink-antikink wave solution of the mKdV equation is(30)RX,T=g1cg2tanhc2Xcg1T.Finally, the kink-antikink wave solution of the platelet spacingΔxn can be obtained as(31)Δxn=hc+g1cg2aca−1tanhc2aca−1×n+1−cg1aca−1t,where the amplitude of the solution A can be written as(32)A=g1cg2aca−1.The kink-antikink solution of the mKdV equation means that there are coexisting phases in the system; in the context of our study, it means a small-scale thrombus and free platelets appear simultaneously. Figure3 shows the coexistence curve when FΔx=−2tanhΔx−1.


## 5. Numerical Simulation

To verify the theoretical analysis results and simulate thrombus formation in a one-dimensional system for lengthL=200 and platelet number N=200, we use the fourth-order Runge–Kutta method to solve the equation. Here, we use the periodic boundary condition and set Δt=0.1, and we set the initial conditions as follows:(33)Δxn0=1,Δxn1=Δxn0=4n≠100,101,Δx100=1+0.01,Δx101=1−0.01,N=200.Figures4–6 show the time-space evolution diagrams corresponding to different functions. Although the types of forces are different, the system becomes increasingly stable with the increase of α. According to stability analysis, the phase transition occurs at α>4 when Δx=1. Therefore, the system becomes unstable with the decrease of α, so some platelets aggregate, which indicates the formation of a thrombus. For FΔx=−2tanhΔx−1, Figure 6 exhibits the kink-antikink wave structure, a typical solution of the mKdV equation. For the force of FΔx=2e−2Δx−1−e−Δx−1, although it shows characteristics of noise, no thrombus exists. Hence, thrombus formation is not related to the viscosity coefficient, and the shape of the platelet is also affected by aspects of the external environment, such as average distribution density.Figure 4
Space-time evolution of the headway for different values ofα when FΔx=2e−2Δx−1−e−Δx−1. (a) α = 3.5. (b) α = 3.7. (c) α = 3.9. (d) α = 4.1.
(a)(b)(c)(d)Figure 5
Space-time evolution of the headway for different values ofα when FΔx=−2Δx−1. (a) α = 3.5. (b) α = 3.7. (c) α = 3.9. (d) α = 4.1.
(a)(b)(c)(d)Figure 6
Space-time evolution of the headway for different values ofα when FΔx=−2tanhΔx−1. (a) α = 3.5. (b) α = 3.7. (c) α = 3.9. (d) α = 4.1.
(a)(b)(c)(d)Figures7–9 show the numerical results for the phase space Δx,α, which verify the stable analysis result. In these figures, we use color to represent maxΔxn−minΔxn when t=20000. Here, the step size of each parameter of Δx,α is 0.05. The dashed lines in Figures 7–9 are consistent with the results of the theoretical analysis, as shown in Figure 3. However, we also find that there exists an unstable area when Δx>1.6 in Figure 4 because the corresponding stability curve is negative at this region. This result shows that even if the platelet density is low, thrombus formation may still be caused by nonlinear effects, which explains why thrombosis still occurs under the low platelet of vaccine for COVID-19 theoretically, which was reported in [52].Figure 7
Heatmap ofmaxΔxn−minΔxn for t=20000 in the phase space of Δx,α, where FΔx=2e−2Δx−1−e−Δx−1.Figure 8
Heatmap ofmaxΔxn−minΔxn for t=20000 in the phase space of Δx,α, where FΔx=−2Δx−1.Figure 9
Heatmap ofmaxΔxn−minΔxn for t=20000 in the phase space of Δx,α, where FΔx=−2tanhΔx−1.


## 6. Conclusions

In this paper, we introduced additional interaction forces into the Stokes formula of viscous force and analyzed the dynamic mechanism of thrombosis. Theoretical analysis for the car-following model is applied because thrombosis is similar to traffic congestion. In addition, the theoretical analysis results are consistent with the experimental results, indicating that platelet density impacts thrombus formation. Nonlinear analysis shows that in the case of force function similar to optimal velocity function, the model can be reduced to the mKdV equation, and it exhibits the characteristics of kink-antikink wave solutions. Numerical results show that thrombus formation is related to fluid and platelet parameters, and phase-space simulation indicates thrombus formation can still occur when platelet density is low. Some of the conclusions are consistent with experimental data and phenomena. In summary, the one-dimensional model we proposed can effectively simulate thrombosis without requiring excessive calculation.



---

*Source: 3150451-2022-04-15.xml*