## Abstract

*Aim*. Arterial occlusion is a leading cause of cardiovascular disease. The main mechanism causing vessel occlusion is thrombus formation, which may be initiated by the activation of platelets. The focus of this study is on the mechanical aspects of platelet-mediated thrombosis which includes the motion, collision, adhesion and aggregation of activated platelets in the blood. A review of the existing continuum-based models is given. A mechanical model of platelet accumulation onto the vessel wall is developed using the dissipative particle dynamics (DPD) method in which the blood (i.e. colloidal-composed medium) is treated as a group of mesoscale particles interacting through conservative, dissipative, attractive and random forces.

*Methods*. Colloidal fluid components (plasma and platelets) are discretized by mesoscopic (micrometre-size) particles that move according to Newton's law. The size of each mesoscopic particle is small enough to allow tracking of each constituent of the colloidal fluid, but significantly larger than the size of atoms such that, in contrast to the molecular dynamics approach, detailed atomic level analysis is not required.

*Results*. To test this model, we simulated the deposition of platelets onto the wall of an expanded tube and compared our computed results with the experimental data of Karino *et al*. (*Miscrovasc. Res.* **17**, 238–269, 1977). By matching our simulations to the experimental results, the platelet aggregation/adhesion binding force (characterized by an effective spring constant) was determined and found to be within a physiologically reasonable range.

*Conclusion*. Our results suggest that the DPD method offers a promising new approach to the modelling of platelet-mediated thrombosis. The DPD model includes interaction forces between platelets both when they are in the resting state (non-activated) and when they are activated, and therefore it can be extended to the analysis of kinetics of binding and other phenomena relevant to thrombosis.

## 1. Introduction

Platelets circulate with the blood and aid in the repair of traumatic injuries of the vessel wall. Platelets are cell fragments produced from megakaryocyte during differentiation. Physiological concentration varies from 150 to 350×10^{9} l^{−1}. Owing to their small size, platelets occupy less than 0.3% of the blood's volume. Their basic shape is disc-like and they can adhere to the vessel wall as well as to each other. There are a large number of receptors of various types, which could be upregulated to the surface of the platelet membrane, participating in the interactions between vessel walls (endothelium), other platelets and different coagulation factors. The impaired function of these interactions causes excessive bleeding or pathological thrombosis within vessels (Vander *et al*. 1998).

Platelet coagulation may be initiated by the expression of tissue factor molecules and the activation of the coagulation cascade. Platelet activation can also be induced by a direct contact between platelets and collagen, enhanced by thrombin or other agonists such as ADP (Fogelson & Kuharsky 1998; Loscalzo & Schafer 2003), or when blood is exposed to a reactive surface, e.g. the synthetic surface of a cardiovascular device (Guyton 1991; Resmi *et al*. 2004). Also, other platelets in the plasma can cohere and can then adhere to the platelets already attached to the wall.

Pathological platelet aggregation is triggered when the balance between pro- and anti-coagulation factors present in the blood is not sustained. This imbalance initiates the coagulation cascade (a series of enzymatic reactions) that enhances platelet activation and leads to thrombus formation. Details about the conditions leading to thrombosis development are given elsewhere (e.g. Kuharsky & Fogelson 2001; Loscalzo & Schafer 2003). A schematic of platelet adhesion and aggregation is shown in figure 1.

It is of particular importance in the context of this study to note that platelet aggregation, adhesion and detachment are affected by factors of a mechanical nature, among which the shear rate, or the shear stress, is dominant. Numerous studies have shown that the shear stress enhances the thrombosis process and can, without any exogenous agonist, even directly cause platelet activation (e.g. Alkhamis *et al*. 1988, 1990; Strony *et al*. 1993; Bluestein *et al*. 1999; Einav & Bluestein 2004; Spieker *et al*. 2004; Yamazaki *et al*. 2005). In a flow field with high shear rate, the interaction between GPIb receptors and vWF multimers can initiate the tethering of circulating platelets to the vessel wall and to already adherent platelets (Fogelson & Guy 2004).

When a platelet is activated, its shape changes and its surface becomes sticky due to the expression of surface receptors (GPIIb/IIIa; Sorensen 2002). These receptors interact with the plasma protein fibrinogen (figure 1). Fibrinogen binds to the receptors of two activated platelets and forms a molecular bond between platelets (Fogelson & Guy 2004).

The focus of this study is on the mechanical aspects of activated platelet–wall and platelet–platelet interactions. We have stated above that platelets must contact the damaged vessel tissue in order to start the process of aggregation. Fluid flow plays the main role for transport of platelets by convection and diffusion. Also, the duration of platelet activation and the proximity of the activated platelets to the damaged tissue are important. Furthermore, high-flow shear stress significantly affects the process of breaking the bonds of existing aggregated platelets. Thus, many aspects of thrombosis are intricately related to the nature of the local blood flow. Many experimental results suggest that vessel branches, recirculation zones near protruding atherosclerotic plaques and bends are zones most favourable for thrombosis (Karino & Goldsmith 1980; Karino & Motomiya 1984).

The mechanical events of platelet-mediated thrombosis include the motion, collision, adhesion and aggregation of activated platelets in the colloidal medium (blood plasma and blood cells). We assume that these mechanical events are governed by external forces (e.g. gravity, inertia) and by internal interaction forces with neighbouring platelets, plasma and other blood constituents. The fluid flow and aggregate growth are mutually coupled, since the fluid flow induces platelet aggregation and the already developed occlusion has a direct influence on the haemodynamics.

Currently, there is considerable effort to model these mechanical processes by computer simulation. The most commonly used computational models treat blood as a continuous homogeneous medium. In this approach, blood flow is usually defined by the continuity equation and the Navier–Stokes equations, and the distributions of platelets and proteins relevant to thrombosis are described by the diffusion–convection reaction equations (Fogelson 1992; Sorensen *et al*. 1999*a*,*b*). These studies use simplifications such as diffusion-limited rates of platelet–surface adhesion, simple flow fields (Poiseuille or Couette) and constant platelet–surface reactivity. These simplifications allow solutions to be obtained in an analytical closed form (Basmadjian 1990). Although these models are based on significant simplifications, they have been helpful in qualitatively describing the relative effects of convection, diffusion and surface reaction on platelet aggregation (Strong *et al*. 1987).

Numerous computational and experimental investigations have been performed to describe, qualitatively, the effects of flow characteristics, such as recirculation and separation flow region, and shear stress, on platelet deposition (Karino & Goldsmith 1979; Badimon & Badimon 1989). These studies are useful owing to the clinical evidence of high prediction of thrombosis in the areas of these flow regions (Ruggeri 2000).

Fogelson (1992) investigated platelet activation, agonist transport and bulk aggregation using a continuum model. A fictitious fluid-dynamic body force field was introduced to simulate flow disturbances due to platelet aggregation. Another discrete-platelet model of platelet–surface adhesion and bulk aggregation was proposed by the same author, but has not been experimentally validated (Fogelson 1984). Sorenson *et al*. (1999*a*,*b*) formulated a two-dimensional continuum model that includes platelet–platelet and platelet–surface adhesion, platelet activation by relevant agonists, platelet-phospholipid-dependent thrombin generation and thrombin inhibition by anti-thrombin III with heparin catalysis.

A major drawback of these types of models is that they treat ‘blood’ as a homogeneous medium and do not describe the behaviour of individual blood constituents (e.g. platelets, red blood cells (RBCs), white blood cells). Blood is a highly complex suspension of chemically and electrostatically active cells suspended in an electrolytic fluid with active proteins and organic substances (Boryczko *et al*. 2003). Especially in microvascular vessels, blood flow cannot be treated as a homogeneous fluid flow but rather the motion of RBCs lubricated by the rest of the colloidal suspension. Therefore, to get an insight into the nature of platelet activation, aggregation and adhesion, it would be instructive to track an individual platelet in the sequence of the process. Boryczko *et al*. (2003) demonstrated such a Lagrangian type computational approach; Chopard *et al*. (2006) applied the lattice Boltzmann approach to model blood flow and platelet-induced clotting; and Harrison *et al*. (2007) modelled flow-related clotting by the lattice Boltzmann method. The objective of this study is to apply a new computational dissipative particle dynamics (DPD) approach to simulate platelet-mediated thrombosis in a simple model system.

In what follows, we first describe a DPD computational model of fluid flow with simple interaction forces. Besides the usual DPD interaction forces, additional platelet–wall adhesion forces are included. The modelling is implemented for a few numerical examples and verified by comparison with experimental results from the literature. These simulations involved up to 2.5 million particles. Finally, we discuss and summarize the computational results.

## 2. Methods

The aim of this study is to develop a DPD model of the platelet adhesion process. We first model fluid flow as a motion of mesoscale DPD particles, and then simulate the thrombosis development through platelet accumulation onto the wall. The basic idea in the model of the platelet adhesion process is to parametrize the bonding interaction force, involved in the process of sticking of activated platelets to the vessel wall, by a characteristic ‘effective bond stiffness’ (or spring constant), *k*_{bw}, whose magnitude is unknown. The values of the parameter *k*_{bw} are obtained by fitting the DPD model simulation results to the *in vitro* experimental data of platelet adhesion onto the wall surface in plasma shear flows.

### (a) Basic DPD equations

Flow of a colloidal fluid is viewed as a motion of the collection of DPD mesoscale particles. The motion of each DPD particle (further called ‘particle’) is described by the following Newton's law equation:(2.1)where *m*_{i} is the mass of particle ‘*i*’; is the particle acceleration as the time derivative of velocity; , and are the conservative (repulsive), dissipative and random (Brownian) interaction forces that particle ‘*j*’ exerts on particle ‘*i*’, respectively, provided particle *j* is within the radius of influence *r*_{c} of particle *i*; and is the external force exerted on particle *i*, which usually represents gradient of pressure or gravity force as a driving force for the fluid domain (Boryczko *et al*. 2003). Hence, the total interaction force (figure 2) between the two particles is(2.2)The component forces can expressed as (Español & Warren 1995)(2.3)Here, *a*_{ij} is the maximum repulsion force per unit mass; *r*_{ij} is the distance between particles *i* and *j*; is the unit vector pointing in direction from *j* to *i*; *γ* stands for the friction coefficient; and *σ* is the amplitude of the random force. Also, *w*_{D} and *w*_{R} are the weight functions for dissipative and random forces, dependent on the distance *r* from the particle *i*, and *ξ*_{ij} is a random number with zero mean and unit variance. The interaction force is equal to zero outside the domain of influence *r*_{c}, hence *F*_{ij}=0 for *r*_{ij}>*r*_{c}.

Furthermore, in order that a DPD fluid system possesses a Gibbs–Boltzmann equilibrium state, the following relation between the amplitudes of the weight functions of dissipative and random forces, *w*_{D} and *w*_{R}, must hold (Español & Warren 1995):(2.4)

Also the amplitude of the random force *σ* is related to the absolute temperature *T*,(2.5)where *k*_{B} is the Boltzmann constant. The weight functions can be expressed in a form (Groot & Warren 1997) given as(2.6)The particles used in this study represent both plasma fluid particles and platelets. The interaction forces between any two particles are defined in equation (2.2), while an additional attractive force is introduced between platelets and the wall (see §2*d*).

### (b) DPD boundary conditions

The implementation of boundary conditions at the wall in DPD is not simple and straightforward. There are several approaches for imposing wall boundary conditions. The layer of particles at a wall is assumed to be frozen, but they interact with other particles.

Free-slip and no-slip conditions at the rigid walls can be imposed by employing the specular and bounce-back approaches, respectively. The local specular reflection is expressed by the following equations:(2.7)where is the unit tensor; the minus and plus superscripts denote velocities before and after collision, respectively, with the wall; stands for the velocity of particle *j* relative to the wall velocity; and *n*_{W} is a unit vector perpendicular to the wall. Thus, only the normal component is reflected while the tangential components remain unaltered. The bounce-back reflection is expressed by the relation (Haber *et al*. 2006; Filipovic *et al*. 2008*a*,*b*)(2.8)

In this study, we used the frozen layers of particles on the walls as well as bounce-back reflection for the particles that reach the walls.

Also, at the inlet boundary a parabolic velocity profile is prescribed while at the outlet boundary the periodic boundary conditions are implemented.

### (c) Integration of DPD equations

In integrating differential equations of motion (2.1) with a time step Δ*t*, we evaluate the resulting interaction force *F*_{i} for the particle *i* as(2.9)The coefficient multiplying the random force comes from the integration of the stochastic equations of motion (Español & Warren 1995). Groot & Warren (1997) gave a physical interpretation of this coefficient. They found that the spread of the random force increases if a given physical time interval is divided into more and more time steps. Thus, a step-size dependence of the random force is precisely generated by the coefficient of .

A ‘usual or common’ approach to perform the integration of differential equation (2.1) is to implement a simple Euler forward method. It was, however, shown that a so-called velocity Verlet algorithm (Groot & Warren 1997) gives more accurate results. The integration scheme for a time step *n* is as given in table 1, where the left upper indices *n* and *n*+1 correspond, respectively, to the start and end of time step. We note that the particle position at the end of time step is calculated by the Euler forward method, while the force and velocity are determined with a use of a mid-step velocity .

### (d) DPD modelling of platelet adhesion to the wall

When an activated platelet comes close to the wall and if shear rate allows, it binds to the wall. However, when adhered platelets are exposed simultaneously to other forces stronger than the binding forces, the bonds break. To incorporate platelet adhesion to vessel walls, we introduce an attractive (bonding) force, , in addition to the conservative, viscous and random interaction forces. As an approximation, we model the attractive force with a linear spring attached to the platelet's surface. The spring is attached to the vessel wall or to an already adhered platelet. The effective spring constant for platelet adhesion on the vessel wall, or to another stationary activated platelet, is denoted by *k*_{bw} (figure 3).

The additional parameter involved in the model is the size of the domain from the collagen-coated wall for which the action of attractive force needs to be considered. We take the attractive force as(2.10)where *L*_{w} is the distance of the activated platelet from the wall.

Platelets are activated when some platelet agonists reach the level that is equal to or above threshold in the concentration around the platelets. Although platelet agonists are known to act synergistically (Ware *et al*. 1987), very little quantitative information is available about the kinetics of these interactions. Fogelson (1992) introduced a single agonist in his continuum platelet aggregation model. Sorensen *et al*. (1999*a*,*b*) analysed multiple platelet agonists, rather than single, as well as the possibility of agonist inhibition via a simplistic linear rate assumption. All agonists in this model have equal weight in determining the rate of platelet activation, with the approximation of time for activation to be 1 s. We here combine finite-element (FE) convective–diffusion equation in order to simulate local agonist concentrations. The governing equation for the concentration is(2.11)where *c*_{i} denotes the concentration of agonist *i*; *v*_{x}, *v*_{y} and *v*_{z} are the velocity components in the coordinate system *x*, *y*, *z*, calculated from the Navier–Stokes equations for velocity field; and *D*_{i} is the constant diffusion coefficient of agonist *i*, while *S*_{i} is the source term for agonist *i*. Therefore, we have the concentration of agonists around the platelets moving inside the plasma as well as around the platelet already adhered to the wall (Kojic *et al*. 2008*b*).

### (e) DPD simulation protocol for modelling of platelet adhesion to the wall

Prior to the DPD calculations, we determine the agonist concentration by solving equation (2.11). We use the FE method and assume the steady-state flow conditions. The attractive force defined in equation (2.10) is scaled by a factor, between 0 and 1, proportional to the concentration at the space position of the considered platelet, where 1 corresponds to the maximum agonist concentration in the whole field.

The DPD simulation protocol was as follows.

Calculation of steady-state fluid flow by running corresponding time steps in DPD. This phase does not include platelets as particles.

After steady-state fluid flow is reached, platelet DPD particles were introduced at the inlet boundary.

At each time step, the displacement/velocity was calculated for each DPD particle (plasma or platelet) and those platelet particles that adhered to the wall within time step were identified.

The total number of adhered platelets was calculated along the vessel wall.

The inlet and outlet boundary conditions of the channel were updated.

Steps (iii)–(v) were repeated until the total time period of the analysis was reached.

The computed axial density distribution of adhered platelets was compared with the experimental results using the least-squares method.

Steps (iii)–(vii) were repeated by adjusting the effective spring constant (

*k*_{bw}) until the best fit with the experimental results was achieved.

## 3. Results

In this section, we first present modelling of plasma flow in a microchannel with narrowing (stenosis) and then platelet aggregation on the collagen surface inside a tubular chamber and inside a tubular expansion, using the DPD method.

### (a) Example 3-1. Plasma flow in a microchannel with narrowing

Steady blood flow between two parallel plates with narrowing is considered (figure 4*a*). We investigate applicability of DPD to the modelling of stenosis inside microvessels. For the comparison, we also present the FE solution (Kojic *et al*. 1998). Parameters of both the FE and DPD models are given in the caption of figure 4.

FE boundary conditions consist of the prescribed inlet parabolic velocity profile, which corresponds to the Poiseuille flow far from the stenosis, and zero traction at the outlet. For the DPD model, we use the prescribed inflow particle velocities of the parabolic profile, bounce-back boundary conditions for the walls and preserve a constant number of particles within the fluid domain.

Velocity profiles, in the region far and at the middle of stenosis, are shown in figure 4*b*. The solution using the DPD method compares well with the FE solution.

### (b) Example 3-2. Platelet deposition in a perfusion chamber

To test application of the DPD method to thrombosis simulation, with the assumption about the wall attractive force, we model platelet deposition in a perfusion chamber. The model corresponds to the experiment of Hubbell & McIntire (1986), who measured platelet accumulation on a collagen wall after 120 s from heparinized whole human blood flowing between parallel plates under steady conditions. The flow is characterized by wall shear rates of 500 s^{−1} and 1500 s^{−1}, while the entrance platelet concentration is 2.0×10^{8} ml^{−1}.

A constant inflow flux of plasma and the entering platelet concentration are imposed. The flow domain is represented by a simple lattice with 10 000×100 particles and periodic boundary conditions along the direction of the fluid flow. The necessary external body force in the flow direction is used to reach wall shear rates of the experiments. The conservative force parameter is taken as *a*_{ij}=25 and the dissipative force parameter as *γ*=4.5 (see §2). It is taken that be the whole chamber gap (200 μm) for the calculation of the wall attractive force (equation (2.10)). The solid walls are modelled by freezing the DPD particles at the wall surface, without possibility of breaking the bond. We use the bounce-back boundary condition for particles that reach the walls (see §2, equation (2.8)).

The experimental results and the computed results for the adhered platelet distribution after 120 s for the shear rates of 500 s^{−1} and 1500 s^{−1} are shown in figure 5. It is found that the best fit of the numerical solution and experiment is achieved for the value of spring constant *k*_{bw}=50 N m^{−1} (Filipovic *et al*. 2008*a*,*b*).

The results are also in agreement with the theoretical continuum-based solution of a diffusion-controlled transport of platelets over a very reactive surface (Kojic *et al*. 2008*b*).

### (c) Example 3-3. Adhesion of platelet to collagen inside a tubular expansion

It is well known that platelet adhesion is increased in regions of disturbed flow, such as in the vortex downstream of a sudden expansion of vessel lumen. We analyse the blood flow inside a tubular expansion using a two-dimensional DPD model. The geometry of the model is shown in figure 6*a*. This problem was studied experimentally by Karino & Goldsmith (1979). In their experiments, they used a suspension of washed platelet with 5×10^{5} cells μl^{−1} and reconstituted blood with 20% haematocrit (2.5×10^{5} cells μl^{−1}). We compare our DPD numerical results using *Re*=37.5, steady flow condition and platelet concentration of 5×10^{5} cells μl^{−1}.

The velocity profiles of plasma with constant platelet concentration are prescribed at the inlet and the outlet of the model in order to accelerate steady-state DPD solution for hydrodynamic macroscale behaviour. We first calculated the two-dimensional fluid flow field by the FE and DPD methods in order to verify the DPD solution for the considered tubular expansion conditions. The computational DPD model geometry is represented by a rectangular box with dimension 10 000×300 in the *x*, *y* directions, respectively, the first being the direction of the flow. At the inflow boundary, the velocity parabolic profile of plasma and platelet particles is prescribed, as shown in figure 6*a*.

At each time step, the total number of DPD particles (plasma+platelets) in the model is kept balanced. Some platelets are deposited at the bottom wall and some go out of the model. The concentration of platelets is kept constant by introducing the corresponding number of new platelets at the inlet flow at each time step. The no-slip conditions at the wall surface have been imposed using a simple bounce-back rule (Haber *et al*. 2006). The number of fluid particles used for this computation was 2 500 000, while 317 500 frozen particles were employed to simulate the walls and narrowed region (figure 6*a*). The results are obtained by subdividing the domain into 1000×400 bins. The computation was run over 300 000 time steps and the results were averaged over every 30 000 time steps. The CPU time for this calculation was 76 hours on INTEL Core 2 Duo E4500 2.2 GHz.

Karino & Goldsmith (1979) used various *Re* numbers and flow conditions (steady and pulsatile). Experimental results from Karino & Goldsmith (1979) and DPD simulation results for *Re*=37.5 and steady flow condition are shown in figure 6*c*.

The results show that immediately after tubular expansion, the number density of adhered platelets rapidly increases until approximately 4 mm and then starts to sharply decline until approximately 5.3 mm which corresponds to the reattachment point of the fluid. As a substantial flow divergence towards upstream and downstream occurs in the vicinity of the reattachment point, the number density of adhered platelets shows a local minimum. After this point, the number density of adhered platelets starts to increase again approaching an asymptotic constant value at the distal part of the tube. The spring constant of the platelet–vessel wall adhesion interaction as well as platelet–platelet interaction, which gives the best fit to the experiment, was *k*_{bw}=70 N m^{−1}. We used the least-squares criteria for the best fitting between computational and experimental results.

## 4. Concluding remarks

The most distinct feature of the DPD method, over the continuum-based methods, such as the FE method, is that the behaviour of each individual particle (plasma particle or platelet in our case) can be tracked in time sequence and thus can be analysed. This has a tremendous advantage over the traditional continuum-based methods in terms of studying its mechanical interaction, including binding, with soundings (i.e. sounding particles) in detail. The drawback of the DPD method, however, is that the simulation is currently limited within small spatial and time domains because DPD models require a very large number of particles (large number of equations). However, the DPD method may be used since platelet aggregation occurs within a local region of blood flow, where molecular and biochemical processes take place in micro time and micro length scales (e.g. Michelson 2002; Loscalzo & Schafer 2003).

DPD calculations generate a large system of equations (of the order of billions for a domain of millimetres in size), which needs to be solved successively over very small time steps (of the order of microseconds). This requires massive parallel computation and the application of the DPD method is currently limited to the simulation within small domains, such as those that exist in microvessels, and for a short period of time. An approach to overcome these limitations in the DPD thrombosis modelling is by using a multiscale modelling, where the classical FE analysis can be implemented over large blood vessel domain, while only a small clinically important zone is modelled by the DPD method (Kojic *et al*. 2008*a*).

We have demonstrated that the DPD method can be used in modelling platelet-mediated thrombosis. It is shown that not only the experimental findings of platelet accumulation to a collagen wall can be matched, but also the attractive interaction force acting on each platelet generated by the collagen surface can be quantified.

## Acknowledgments

This research was supported by the Ministry of Science of Serbia, TR6209, OI144028 and National Heart, Lung and Blood Institute grants NIH HL054885, HL070542 and HL074022.

## Footnotes

One contribution of 11 to a Theme Issue ‘The virtual physiological human: building a framework for computational biomedicine II’.

- © 2008 The Royal Society