Three-dimensional multi-scale model of deformable platelets adhesion to vessel wall in blood flow.

When a blood vessel ruptures or gets inflamed, the human body responds by rapidly forming a clot to restrict the loss of blood. Platelets aggregation at the injury site of the blood vessel occurring via platelet-platelet adhesion, tethering and rolling on the injured endothelium is a critical initial step in blood clot formation. A novel three-dimensional multi-scale model is introduced and used in this paper to simulate receptor-mediated adhesion of deformable platelets at the site of vascular injury under different shear rates of blood flow. The novelty of the model is based on a new approach of coupling submodels at three biological scales crucial for the early clot formation: novel hybrid cell membrane submodel to represent physiological elastic properties of a platelet, stochastic receptor-ligand binding submodel to describe cell adhesion kinetics and lattice Boltzmann submodel for simulating blood flow. The model implementation on the GPU cluster significantly improved simulation performance. Predictive model simulations revealed that platelet deformation, interactions between platelets in the vicinity of the vessel wall as well as the number of functional GPIbα platelet receptors played significant roles in platelet adhesion to the injury site. Variation of the number of functional GPIbα platelet receptors as well as changes of platelet stiffness can represent effects of specific drugs reducing or enhancing platelet activity. Therefore, predictive simulations can improve the search for new drug targets and help to make treatment of thrombosis patient-specific.

Introduction

When a blood vessel ruptures or gets inflamed, the human body responds by rapidly forming a clot to restrict the loss of blood. Blood clot (thrombi) formation is a complex biological process involving an extensive system of biochemical (coagulation) reactions, platelet hydrodynamics, platelet–platelet and platelet–blood vessel wall interactions leading to ligand–receptor adhesion bond formation and platelet activation. Platelet adhesion to the vessel wall is one of the first events associated with formation of haemostatic clots and pathological thrombi.

In this paper, a new three-dimensional multi-scale model of platelet–blood flow–vessel wall interactions combining submodels at three biological scales crucial for the early platelet aggregation is introduced and calibrated to investigate how platelet stiffness, GPIb receptor expression and platelet–platelet interaction affect platelet–wall adhesion quantified in terms of platelet pause time. We implemented a novel approach of combining a recently developed platelet hybrid membrane submodel, the subcellular element (SCE) representation of the cytoskeleton network and a continuum description of the lipid bilayer to study the very first step of blood clot formation, the rapid formation of unstable bonds which slow platelets and cause platelet-flipping and adhesion to the damaged surface. The hybrid platelet model was also coupled with the lattice Boltzmann model (LBM) of blood using the immersed boundary method (IBM) to simulate platelet motion and deformation in shear flow. The kinetics-based adhesive dynamics model was also integrated into the three-dimensional model to simulate formation and disassociation of the receptor–ligand bonds during the platelet–platelet and platelet–vessel wall interactions. Parallelized model simulations were implemented on a GPU computer cluster which speeded up simulations by a factor of 100 (table 3) in comparison with CPU implementation which allowed, for the first time, to run biologically relevant predictive simulations.

By using a novel biologically calibrated three-dimensional modelling approach, it is shown that the platelet stiffness, the number of GPIbα platelet functional receptors and mutual interaction between platelets can significantly alter the adherence of platelets at the site of vascular injury. Our results demonstrate how a comprehensive modelling approach coupling three biologically relevant scales can provide new insights into the biomedically important problem of early thrombus development. Variation of the number of functional GPIbα platelet receptors as well as changes of platelet stiffness can represent effects of specific drugs for reducing or enhancing platelet activity. This emphasizes the importance of predictive simulations as they can potentially improve the search for new drug targets and help with making treatment of thrombosis patient-specific.

Damage or alteration of a blood vessel lining can result in activation of flowing platelets and their subsequent aggregation at sites of vascular injury. The ability of platelets to tether to and translocate on injured vascular endothelium relies on the interaction between the platelet glycoprotein receptor Ibα(GPIbα) and the A1 domain of von Willebrand factor (vWF-A1) [1].

Along with biochemical activation of platelets, large shear disturbances of blood flow are one of the key factors promoting pathogenic activation of platelets and formation of thrombi. In addition, platelets flowing in a whole blood exhibit increased concentrations in the vicinity of the vessel wall, making platelet–platelet interactions more frequent near vascular surfaces. Excessive accumulation of platelets at injury sites is one of the pathological events that result in acute myocardial infarction, sudden death and ischaemic stroke. This pathological process is responsible for mortality and morbidity rates higher than for any other disease, making platelets a major target for therapeutic interventions. Thus, studying individual platelet dynamics as well as platelet–platelet interactions and platelet adhesion to a vascular or thrombus surface is of high biomedical importance and urgency.

Effects of shear flow on accumulation of platelets on various surfaces have been extensively studied in in vitro and in vivo experiments [1-7]. However, there is a limited amount of available experimental data on an individual platelet dynamics in the vicinity of the vascular surface as well as platelet–surface attachment. There is also a lack of experimental data demonstrating how platelet–surface attachment is affected by mechanical properties of a platelet as well as by platelet receptor–ligand kinetics. Better understanding of platelet aggregation requires study of the interplay among biochemical, mechanical and hydrodynamic processes occurring at different scales, including a nanometre scale (receptor–ligand kinetics), a micrometre scale (cellular level) and a millimetre scale (early platelet aggregate). Multiple characteristic scales make it difficult to experimentally discern effects of different processes involved in platelet–surface attachment and overall thrombus growth dynamics. Meanwhile, a multi-scale modelling approach can provide a useful predictive tool to aid in elucidating mechanisms of platelet–wall attachment and aggregation.

Several multi-scale models attempting to couple large numbers of submodels at different scales have been developed (see, among others, for reviews [8,9]). These models implemented simplified submodels in order to make simulations less computationally expensive. It is extremely difficult at this time, if not impossible, to validate predictions of multi-scale models attempting to combine submodels at all scales representing processes of blood clot formation using existing experimental data. In addition, most experimental data are currently available at the molecular level and individual platelet level. Therefore, it is important to develop detailed multi-scale models coupling two or three scales and considering only a few processes at a time. Such models when properly calibrated with available experimental data can provide useful predictive tools aiding in designing new experiments, drug design and planning new patient-specific therapeutic strategies.

Several computational models have been developed to characterize platelet and other types of blood cell motion and adhesion dynamics under hydrodynamic shear flow at cell and receptor levels (see [8,9] for a review). Analytical solutions for forces and torques exerted on a platelet treated as a rigid object in the Stokes regime in a two-dimensional case were obtained in [1,10] and compared with the data obtained using an image analysis algorithm for tracking the motion of platelets before, during and after contact with the surface. Kinetic properties of the receptor–ligand adhesion bonds, GPIbα–vWF, were quantified in [1,4] using Monte Carlo simulations and pause time analysis of transient capture/release events. This approach provided association and disassociation rate constants kon and koff, depending on the shear rate of the blood flow.

Experimental study in [11] showed that platelets have viscoelastic properties and the elastic moduli in the range of 1–50 kPa. Large deformation occurred when platelets were suspended in the shear flow [12]. To account for the elastic and viscoelastic properties of cells, a number of methods accounting for cell structural properties have been developed [13,14]. The SCE model introduced in Sandersius & Newman [13] represented each cell by a collection of elastically coupled SCEs, interacting with each other via short-range potentials. Sweet et al. [15] and Xu et al. [14] presented a three-dimensional modelling approach in which cells, modelled by SCEs, were coupled with fluid flow and substrate models by using the Langevin equation.

The fluid–structure interaction approach is an essential part of the model. Previously, the IBM introduced by Peskin [16] to investigate blood flow in the human heart has been applied to many other fluid–structure interaction problems, including platelet aggregation [17] and deformation of red blood cells [18]. Skorczewski et al. [19] developed a two-dimensional model, using a lattice Boltzmann IBM to investigate the motion of platelets near a vessel wall and close to an intravascular thrombus, in which they modelled the platelets as rigid bodies, whereas the red blood cells were represented as deformable bodies.

The results of the predictive simulations of the three-dimensional model introduced in this paper revealed that the platelet pause time strongly depends on the stiffness of the platelet as well as on the number of expressed GPIb membrane functional receptors. Additionally, we demonstrated that the platelet–platelet interaction near the surface of the vessel wall could significantly decrease the platelet paused time, and thus decrease the rate of platelet attachment to the injury site.

This paper is organized as follows. It starts with the Biological Background section. Then the Methodological Innovation is described in detail, including description of submodel at each of three space scales and of the coupling approach. This is followed by the Results section, which includes model validation and description of the predictive simulation. Biological relevance of the predictive simulations is discussed in the Discussion section. GPU implementation of the three-dimensional model is described in appendix A.

Biological background

The mechanism by which platelets bind to a damaged blood vessel wall is similar to that of leucocyte binding to activated endothelium [20], and requires two binding steps. The first step is the rapid formation of unstable catch–slip bonds which slow platelets and cause platelet-flipping along the damaged surface. (Counterintuitively, the dissociation rate first decreases with increasing force until reaching a threshold.) This is mediated by the platelet receptor component, GPIbα, forming transient bonds with the vWF exposed at the injury site. Rapid association and dissociation kinetics of the bonds result in transient tethering and subsequent flipping (or rolling) and pausing of platelets on the vessel surface [1,21]. Then, stable bonds slowly form between platelet receptors and ligands (often integrin αIIbβ3 binding with vWF or fibrinogen) bound to the damaged wall or the surface of the thrombus resulting in strong adhesion, initiating transmembrane and, subsequently, intracellular signalling. As the blood clot grows, platelet–platelet interaction becomes one of the major factors determining clot growth rate and integrity as platelets expose GPIIbIIIa receptors which permit platelet–platelet adhesion via fibrinogen. Adhesion of platelets to the injured surface is also affected by shear rates of the flow. At high shear, platelet integrin α2β1 and GPVI receptors are not sufficient to initiate binding to collagen, and binding of the GPIbα receptor to vWF immobilized on collagen becomes essential in platelet adhesion.

The stiffness of the platelet not only determines the shape and morphology of the clot, but also affects clot mechanical properties, as platelet stiffness determines cell shape when exposed to various flow conditions and contact interaction with other cells and blood vessel wall. This will affect the number of receptor–ligand pairs in platelet–platelet and platelet–substrate interactions. Platelet stiffness is also an important property reflecting platelet functioning, because it reorganizes its structure during activation or as a response to physiological or pathological conditions.

To date, to the best of our knowledge, cumulative effects of platelet stiffness, different levels of expression of GPIbα receptors and platelet–platelet interaction impacting the strength of platelet–substrate binding have not been systematically investigated. Our model provides a unique means for quantitatively understanding these effects, which are critical for improving our knowledge about the initial stage of the blood clot formation.

Methodological innovation of the three-dimensional modelling approach

The novelty of the three-dimensional model lies in developing a novel membrane submodel as well as in new approaches of coupling submodels of biological processes at three spatial scales (figure 1) which are crucial to early blood clot formation. At the subcellular scale (nanoscale), a kinetics-based stochastic dynamic adhesion submodel (DAM) is used to simulate vWF–GPIbα binding and GPIbα–vWF–GPIbα binding, in which individual vWF and GPIbα molecules are represented as elastic springs (table 1). This is justified by the fact that these receptor–ligand bindings are probabilistic in nature [1]. Moreover, individual filaments in the cytoskeleton network of the platelet membrane are treated as coarse-grained harmonic springs. At the cellular scale, a novel continuum description of the lipid bilayer of the cell membrane is used. We developed this new platelet membrane model to study effects of membrane stiffness on cell–substrate interaction, which was shown to strongly affect platelet–injury site adhesion. (See also §4c(i) for model prediction.) The subcellular scale and the cellular scale components are integrated by distributing GPIbα receptors at the vertices of the cytoskeleton network and by superimposing the cytoskeleton network and the lipid bilayer. At the macroscale, the dynamics of the fluid flow is represented using the LBM to facilitate parallelizing the simulation code on GPUs. The platelet model is coupled with the LBM using the IBM. (The coupling and data flow between all the submodels are demonstrated in figure 1.) We calibrate and validate this three-dimensional model by comparing simulations at different scales with either theoretical results or available experimental data at these scales. Specifically, the platelet model coupling with the LBM was validated using theoretical results and previous simulation results (see also §4a), whereas the platelet–substrate adhesion simulations were compared with experimental data to calibrate the DAM under different flow conditions (see also §4b).

Figure 1.

Flow chart of the simulation algorithm. Using F(t) hybrid membrane model, the forces F(t) acting on cell elements were calculated. The forces f acting on fluid nodes were spread from F(t) by the immersed boundary method. The velocity field v(x,t+Δt) of fluid was obtained by the lattice Boltzmann method. The velocities of cell elements dX/dt(t+Δt) were determined by immersed boundary interpolation. The cell elements positions were updated based on the velocities. Finally, the stochastic adhesion model was used to determine the force Fbond acting on receptor–ligand bonds that bind platelets to vessel walls. (Online version in colour.)

Table 1.

Biological processes and submodels at different scales.

scales	processes	submodels	coupling
<0.1 μm subcellular-level nanoscale	ligand–receptor interactions	stochastic dynamic adhesion model	1. Nano–micro scales: coupled by explicitly modelling receptors on platelet membrane nodes
approximately 1 μm cellular-level microscale	individual platelet deforming, flipping and adhering to vessel wall	hybrid membrane model	2. Micro–macro scales: coupled through immersed boundary method
>10 μm macroscale	blood flow and its Interaction with platelet	lattice Boltzmann method

At each time step of simulation, the hybrid membrane model is first used to calculate forces acting on the nodes of the Lagrangian mesh representing platelet geometry, such as bond forces resulting from stretching or compression of the cytoskeleton network, bending forces resulting from deformation of the lipid bilayer and attraction/repulsion between platelets and the environment owing to formed ligand–receptor bonds. This is followed by coupled LBM and IBM to update fluid flow and position of platelets. Finally, the MC computations of platelet adhesion to a surface expressing vWFs are performed to break the already formed bonds and to generate new bonds from unbound GPIbα and vWF.

We note that this is the first time that a detailed platelet membrane model has been developed and implemented on GPUs for studying cell–flow, cell–cell and cell–substrate interactions. Because of the speed-up gained by GPU implementation, we are able to investigate effects of these interactions and cell mechanics on platelet dynamics in a timely manner. Additionally, this model can be directly used for modelling any biological cells with membrane structures similar to these of eukaryotic cells. Here, we describe in detail individual submodels and explain how they are coupled.

Platelet membrane submodel

We simulate the motion of platelets in a three-dimensional region bounded by an infinite flat plane at z=0 (see figure 2a for example). A platelet has initial shape defined by x2/a2+y2/a2+z2/(λa)2=1, where a=1 μm is the approximate particle radius and λ=0.25 is the aspect ratio [22,23]. The Reynolds number of this system is Re=γρa2/μ=O(10−3), where γ=300 and 400 s−1 are the shear rates used in experiments [1], a=1 μm is the particle radius, ρ=1.0239 g cm−3 is the density of blood plasma and μ=1.2 cP is the viscosity of blood plasma [24].

Figure 2.

(a) Schematic diagram showing one platelet translating and rotating in shear flow near an infinite plane wall. (b) Structure of a platelet consisting of 958 SCEs. The major radius, a, and centroid height, H, are defined as is the coordinate system and flow direction. One platelet is represented by a collection of elastically linked SCEs, interacting with one another via a spring-like elastic force. The GPIbα receptors are randomly uniform distributed on the platelet membrane, and vWF ligands are distributed on the wall. (Online version in colour.)

The platelet membrane, which is similar to the membrane of a red blood cell, is also assumed to consists of a lipid bilayer and an attached cytoskeleton. Following ideas from [16,25], the platelet membrane surface geometry is represented by a triangular mesh consisting of a collection of N (N=958 in our simulation) points {,i∈1…N} (figure 2b). The connected edges of the mesh are used to model the cytoskeleton network of the platelet membrane, and the triangulated mesh surface represents the lipid bilayer of the cell membrane, where the cytoskeleton attaches to. The mesh points represent coarse-grained actin vertices and each edge of the mesh represents a coarse-grained filament. The Helmholtz free energy of the membrane is defined as Here, HSCE is the in-plane energy of the cytoskeleton network; Hbending is the bending energy representing the bending resistance of the lipid bilayer; Hvolume and Harea are volume and area conservation constraints, respectively, and Hwall represents the energy relating to interactions due to ligand–receptor binding (explained in detail in §3c).

We use a harmonic ‘spring’ model to simulate the elasticity of the edge connecting mesh points i and j, which mimics a coarse-grained filament. The associated potential energy functions for points i and j are where L is the rest length, =− the position vector difference for points i and j, respectively, and k=EΔx/5 the coefficient that defines the spring ‘stiffness’ [26] for elastic modulus E=25 kPa [11] and Δx=0.1 μm unit link length of the spring. The total potential energy for the cytoskeleton network is . The corresponding force vector acting on point i by point j is The area and volume conservation constraints, which account for area incompressibility of the lipid bilayer and incompressibility of the inner cytosol, respectively, are expressed as and where k, k and k are the global area, local area and volume constraint coefficients, respectively. The terms Stotal and V denote the surface area and volume for the whole platelet, whereas S0, and V 0 are the individual triangle mesh area, the total membrane area and the volume for unstressed platelet, respectively.

We adopt the energetic variational approach developed in [27] to represent the lipid bilayer of the cell membrane. Let Σ∈R3 be a smooth, closed surface representing the lipid bilayer of the platelet. The bending energy of the lipid bilayer is defined as [27] where K(x)=1/2(κ1()+κ2()) is the mean curvature, and κ1(), κ2() are the principle curvatures at the point x. We follow the finite-element method in [28] to calculate κ1() and κ2(). Briefly, let u(ξ,η) be a function defined over a triangle of the surface mesh representing the lipid bilayer and approximated as where ξ, η are the local parametric coordinates, u is the value of u at node i and N(ξ,η) are the basis functions for a quadratic six-node triangular finite element. To evaluate the membrane curvature tensor , one needs to calculate the left Cauchy–Green strain tensor, which is determined from the surface deformation gradient tensor, . For each triangular element, the surface deformation gradient tensors at the element nodes are obtained by solving the following system of equations at each node of the element, and , are its positions in the unstressed state and after deformation at time t, respectively. is the unit normal vector to the undeformed membranes. To evaluate the curvature tensor κ at a point, one needs to solve at each element node and then average over the elements sharing that node, and n is the unit normal vector to the deformed membranes. The mean curvature is given by The normal component of the elastic force-associated bending energy (3.6) is obtained by taking variational derivative and is given by bend=(ΔK+2K3).

Thus, nodal forces F are derived from the total energy as follows Computation of ∂(Hwall)/∂ is explained in §3c.

Platelet stochastic dynamic adhesion submodel

The kinetics-based stochastic DAM based on ideas of the Dembo model [29,30] is used to simulate the GPIbα of unactivated platelet binding to immobilized vWF on the vessel–wall or platelet–platelet adhesion through forming GPIbα–vWF–GPIbα bonds, in which vWF was originally in plasma. Here, we provide details of the model for GPIbα–(immobilized) vWF bond formation; modelling of GPIbα–vWF-GPIbα is treated similarly. Each platelet has approximately 10 688 GPIbα receptors distributed uniformly on its membrane surface, to achieve a surface density of approximately 1500 receptors μm−2 [31]. In our model, 5344 receptor point locations on the platelet surface are randomly distributed on the platelet membrane mesh, with each point location representing two GPIbα receptors, because there are two GPIbα receptors existing on each GPV molecule [32]. On the bottom plane of the simulation domain, z=0, immobilized vWFs are uniformly distributed, resulting in a vWF density of 25 μm−2, which is consistent with experimental conditions in Doggett et al. [1].

The following rules are used for governing the GPIbα–vWF binding [33]. (i) Two vWF molecules cannot bind to the same receptor nodes for reasons of steric blocking, and (ii) receptors from a maximum of four receptor nodes present on a platelet surface can bind a vWF molecule.

In our stochastic DAM, when an unbound GPIbα and an unbound vWF are separated by less than the length of a GPIbα–vWF bond of 128 nm [34,35], a test for forming a bond is performed. Next, the formed bonds are tested for breakage. A GPIbα–vWF bond is modelled as a linear spring.

Probabilities of GPIbα–vWF bond formation and dissociation are calculated using Pf (probability of forward reaction) and Pr (probability of reverse reaction) described in [36]: where koff and kon are given in s−1 units and Δt is the simulation timestep. The reverse rate constant is calculated using the Bell model for the force-dependent dissociation rate of weak non-covalent bonds where koff(Fb) is the bond dissociation rate, k0off is the unstressed off-rate, γ is the reactive compliance, Fb is the applied force on the bond and kBT is the product of Boltzmann constant and temperature. The dependence of bond formation rate constant kon on the deviation bond length is described by [29,33] as where k0on is the intrinsic cross-linking formation rate constant, σ is the spring constant, lb is the equilibrium bond length, xb is the distance spanning the endpoints of the GPIbα receptor on the platelet surface and the vWF-A1 binding site.

The adhesion force of the GPIbα–vWF bond located at ith node of the cell membrane is calculated using a spring model as follows where σ is the spring constant, lb is the equilibrium bond length, xb is the distance spanning the endpoints of the GPIbα receptor on the platelet surface and the vWF-A1 binding site. Table 2 lists values of parameters of the DAM used in simulations.

Table 2.

Values of physical parameters used in simulations.

parameters	definition	value	references
a	platelet radius	1.0 μm	[23]
λ	platelet aspect ratio	0.25	[22]
Υ	flow shear rate	300 and 400 s−1	[1]
ρ	blood plasma density	1.0239 g cm−3	[24]
μ	blood plasma viscosity	1.2 cP	[24]
E	platelet elastic modulus	25 kPa	[11]
l0	average length of initial spring length	75 nm
ks	global area constraint coefficient	6000	[25]
kt	local area constraint coefficient	6000	[25]
kv	volume constraint coefficient	6000	[25]
k0	bending modulus	200 kBT	[27]
T	temperature	300 K
k0on	intrinsic cross-linking formation rate	10−5 s−1	[1]
k0off	unstressed disassociation rate	3.45 s−1	[1]

Lattice Boltzmann method for simulating blood flow

The LBM uses purely localized fluid particle evolution and relaxation, which in turn facilitates parallelization in computer implementation. The LBM decomposes the fluid domain into structured lattice nodes and operates on the lattice. The fluid is modelled as a group of fluid particles that are only allowed to move between lattice nodes or remain at rest. The composition of the lattice nodes depends on the chosen lattice model. In this paper, we used the three-dimensional model of a cubic lattice (16×64×16 μm with spacing h=0.2 μm) with 19 discrete velocity directions (model D3Q19, as shown in figure 3). The LBM solves the Boltzmann equation describing the dynamics of fluid from a microscopic point of view: in fluid, particles with velocities v collide with certain probability and exchange momentum. The collisions are assumed to be ideal, that is the total momentum and energy is conserved during the collisions. The Boltzmann equation describes the probability f(x,v,t) of finding a particle with velocity v at a position x and at time t evolving with time where F denotes an external body force, ∇ is the gradient in position and momentum space and Ω(f) denotes the collision operator which is chosen as a relaxation of f with a characteristic time τ to the equilibrium distribution f(eq)(v,ρ): The equilibrium distribution function depends on the local density ρ(,t) and the velocity field (,t). In the D3Q19 lattice model, 19 values f(x,t) are stored at each lattice site assigned to a lattice vector . The local density at a lattice point is obtained by summing all f, and the streaming velocity is given by where =h/Δt is the lattice speed associated with the ith direction and Δt is the time step of our simulation.

Figure 3.

Lattice Boltzmann D3Q19 (three-dimensional and 19 velocities) model. The lattice vectors c represent the velocities of the particle moving from the centre grid point to its neighbour grid point.

Using a Chapman Enskog expansion, Guo et al. [37] showed that the following lattice Boltzmann equations give a second-order accurate , the Navier–Stokes velocity in the presence of a spatially varying, time-dependent force where is a streaming velocity defined in equation (3.18), =+Δt/2ρ, and with the lattice speed of sound for the D3Q19 lattice and the lattice weights The pressure turns out to be proportional to the density and the dynamic shear viscosity is given by

To ensure convergence and stability of LBM, we follow the method in [38] to choose our parameters. Spacing h=0.2 μm was determined by our simulated fluid domain and memory size of the GPU card. Time step Δt was determined from the equation [38]: Δt=(τ−0.5)h2/(3υ), where υ=μ/ρ is kinetic viscosity, μ and ρ are fluid viscosity and density as defined in table 2. Generally speaking, a larger value of τ leads to a more stable LBM simulation, and τ must be greater than 0.5. We set τ=1.379 s in our model, such that Δt=10−8 s.

Periodic boundary conditions in x–z and y–z boundary planes (y=0, y=64 μm, x=0 and x=16 μm), are realized by propagating the f from the computational domain on the one boundary to the boundary on the opposite side of the domain. In the x–y boundary planes, we used the onsite velocity boundary conditions proposed by Hecht & Harting [39]. For instance, in the x–y boundary plane z=0, f (i=1, 2, 3, 4, 6, 7, 8, 10, 11, 12, 14, 16, 18, 19) can be obtained from the streaming step, but f (i=5, 9, 13, 15, 17) are undetermined. Following the methods of Hecht & Harting [39], we obtain Here, v, v and v are boundary velocities in x-, y- and z-directions, and the transverse momentum corrections on the z-boundary for distributions propagating in x- and y-directions, respectively and

Coupling platelet, dynamic adhesion and flow submodels

Coupling platelet and dynamic adhesion submodels

As described in §3c, GPIbα receptors are randomly distributed on the cell membrane. In each step of the simulation, forming and breaking a GPIbα–vWF bond is updated using the DAM. When a formed bond is either stretched or compressed, the bond deformation force is computed using equation (3.14). This force is exerted on the cell membrane at the place where the GPIbα receptor of the bond is located. When only the platelet membrane and vessel wall interaction is considered, the term Hwall of equation (3.1) represents the energy associated with these interactions. In particular, ∂(Hwall)/∂ corresponds to the sum of following two forces (i) adhesion forces caused by GPIbα–vWF bond and (ii) short-range repulsive forces accounting for contact of vessel wall. The short-range repulsive force is given by an empirical relationship: Frep=F0(τe−)/(1−e−), where F0=500 p Nm, τ=2000 μm−1 and ε is the separation distance between platelet membrane and vessel wall [29]. Thus, ∂Hwall/∂ term in equation (3.11) is defined as

Coupling cell and flow submodels

To couple the integrated platelet and stochastic DAM submodels with the blood flow computed by the LBM, we use the IBM [16]. In the IBM (figure 4), the Eulerian description is used for the fluid dynamics, and the Lagrangian description is used for objects immersed in the fluid. Using lowercase letters for Eulerian variables, and uppercase letters for Lagrangian variables, we have and where is time, the flow velocity, the speed of the solid object boundary, the fluid flow coordinate, the boundary coordinate, the force density on the fluid node, the force density on the solid elements and δ() the Dirac delta function.

Figure 4.

Eulerian fluid grid (black) and Lagrange solid elements (red). A Eulerian description is used for the fluid dynamics, and a Lagrangian description is used for objects immersed in the fluid. The communication between these two coordinate systems is realized by the immersed boundary method. (Online version in colour.)

Equations (3.31) and (3.32) are approximated using a regularized discrete delta function δ. The discretized forms of equations (3.31) and (3.32) using δ are as follows and where h is the fluid node spacing, =(ih,jh,kh) the coordinate of the i,j,kth Eulerian grid node, the Lagrange coordinate of the mth elements, the force density on , the force density on , the velocity of , the velocity of X. The discrete delta function δ appearing in equations (3.33) and (3.34) is a smoothed approximation to the Dirac delta function δ(). (The detailed derivation procedures in several forms have been presented in the literature [40].) We use the following common form and

To sum up, first, in each step of the simulation, equation (3.20) is solved. Then, positions of nodes of the platelet membrane are updated by equation (3.33). Finally, the MC computations are performed to break the already formed bonds and to generate new bonds from unbound GPIbα and vWF.

Results

First, the model was verified by comparing simulation results with analytical solutions and available model simulation data [24]. Next, we validated the model by comparing the simulation results with the experimental data [1] on flipping platelets flowing over a vWF-coated surface. A calibrated three-dimensional model was used to predict how the stiffness of a platelet membrane, the number of receptors on a platelet membrane and the strength of platelet–platelet adhesion affect the paused time of the platelet adhering to a vessel wall.

Validation of fluid–platelet coupling by comparing with the Jeffery orbit

Mody et al. [10] described theoretical solutions using the Jeffery orbit theory and provided predictions obtained using the analytical platelet-flipping model. This analytical solution (shown as solid red line in figure 5) did not consider the wall effect and applied only to the cases of platelet motion far from the wall (H/a>20) [24], where H is the centroid height of platelet and a is the major radius (as shown in figure 2). Mody et al. [24] modified the completed double-layer-boundary integral equation method to include a flat surface boundary that was used to compute the effects of the wall on the flow behaviour of a platelet. Platelets located as far as 2.4-fold platelet radii from the surface ‘display modified’ Jeffery orbits with periodic rotational motion in the direction of flow (green dashed line in figure 5). To verify our model, we simulated the flipping of a single platelet located at the distance of 2.4a as well as greater than 20a, from the vessel wall. Our simulations revealed that the calculated orbit of rotation (blue dashed line in figure 5) agreed well with the Mody's simulation results [24] (green dashed line) within an experimental error of 2.65% for platelet located at the distance of 2.4a, and agreed perfectly with Mody's simulation results [24] for platelets located at a distance of more than 20a. Figure 6 shows the series of snapshots from our simulations of platelet-flipping in a shear flow near the vessel wall. In our model, the platelet was modelled as an elastic cell with the elastic modulus measured by the AFM experiments [11], whereas Mody et al. [10] considered the platelet as a rigid object. By comparing our simulation results and results of Mody et al. [10], we conclude that our simulations can be successfully implemented to model the motion of individual resting platelets revealing high stiffness membrane values.

Figure 5.

Validation of the platelet dynamics model. The analytical solution for the platelet rotational trajectory (Jeffery orbit), trajectory calculated by a completed double-layer–boundary integral equation method (CDL-BIEM) and our simulation (LBM-IBM) are shown by solid and dashed colour lines (inset key). (Online version in colour.)

Figure 6.

The configuration of the simulated platelets at different times flipping over the vessel wall for the wall shear stress of 3.0 dyn cm−2. Image sequence (a) shows projection of the platelet on the x–y plane from dimensionless time point 0 to 14. Image sequence (b) shows projection of the platelet on the y–z plane from dimensionless time point 0–14. The coordinate system is defined in figure 2a. (Online version in colour.)

Validation of the model of the platelet–substrate adhesion

To validate the kinetic submodel, we simulated flowing platelets adhering to substrate through GPIbα–vWF binding and calculated koff rates to compare with available experimental data. The model parameters used in our simulations (table 2) were obtained in biological experiments [1,11,22-25,27]. The adhesive dynamic parameters were measured in in vitro flow chamber tests [1].

Doggett et al. [1] measured the kinetics that governs platelet interactions with vWF in haemodynamic flow. In their experiment, the frequency of tethering for platelets was measured by determining the percentage of cells that paused, but did not translocate, on vWF substrates. The frequency of tethering for microspheres coated with vWF on antibody-immobilized platelet substrates was also measured. A transient tether event was defined as flowing platelets that abruptly halted forward motion for a defined period of time and subsequently released, without evidence of translocation, to resume a velocity equivalent to that of a non-interacting cell. Dissociation rate constants (koff) were determined by plotting the natural log of the number of beads that interacted as a function of pause time after the initiation of tethering (figure 7, the slope of the line is −koff).

Figure 7.

The number of tethering events as a function of the platelet paused time. The solid lines are the fitting lines of experimental data for shear stresses of 3.0 dyn cm−2 (shown in blue) and 4.0 dyn cm−2 (shown in red). The corresponding slopes of the fits (koff values) are −4.83 and −5.18. The dashed lines are the fitting lines of simulation results (shown with circles) for shear stresses of 3.0 (shown in blue) and 4.0 dyn cm−2. The corresponding slopes of the fits (koff values) are −3.31 and −3.58. (Online version in colour.)

To calculate dissociation constants, we performed simulations for various numbers of random seeds (1000–1200). The results of the simulations and experimental data are presented in figure 7 for two different flow shear rates as the natural log of the number of platelet tethering events versus the pause time. The values of dissociation rate koff were found to be 3.31 and 3.58 s−1 for flow shear stress 3.0 and 4.0 dyn cm−2, respectively. The corresponding experimental values obtained in [1] were 4.83 and 5.18 s−1.

It should be mentioned that our simulations confirmed several experimental observations. It was reported in [1] that in the range of flow rates considered, forces acting on the GPIbα–vWF bond were not sufficient to alter the rate of dissociation koff. Our simulations also demonstrated that the koff values altered only in a small range from 3.31 to 3.58 s−1. Additionally, it was reported in [1] that the forces acting on a platelet in shear flow were 14.7 and 19.6 pN for flow shear stress 3.0 and 4.0 dyn cm−2, respectively, whereas our model yielded very close force values of 12.8 and 15.6 pN, respectively.

Our simulations also confirmed that in the range of flow rates (0–4 dyn cm−2 wall shear stress), bond association and dissociation kinetics can be successfully described by the Dembo model (equations (3.12) and (3.13)).

Predictive simulations

The responses of a platelet to interactions with the environment depend, among others, on the mechanical forces that platelets experience. Here, we consider effects of platelet membrane tension, flow shear stresses and adhesion bond forces on platelet–substrate adhesion dynamics.

Effect of platelet membrane stiffness

Simulation results in §4a show that the flow dynamics of the platelet in linear shear flow can be studied by modelling platelets as rigid objects. How the stiffness of the platelets affects the platelet–substrate interaction remains to be answered. In [41], it was reported that alteration of platelet stiffness can modulate platelet aggregation. We hypothesized that softer cells lead to prolonged adhesion time and could potentially increase chances of platelets being activated after adhesion. Here, we report the simulation results indicating remarkable changes in platelet paused time as the platelet membrane stiffness changes. We varied the platelet membrane stiffness from 25 to 2.5 kPa, and performed simulations with 30 different random seeds to obtain 30 different paused times under flow shear stress of 3.0 dyn cm−2. The paused time was 6.69±0.71 s (mean±SD) for the membrane stiffness of 25 kPa, which is about twice higher than the paused time 3.15±0.69 s (mean±SD) for the membrane stiffness of 2.5 kPa (t-test, p<0.0008, figure 8). The total deviation of all the nodes in the deformed shape in figure 8a is 3.5 μm compared with the reference configuration, and in figure 8b is 0.28 μm. Thus, these simulation results indicated that softer cells have prolonged average paused time.

Figure 8.

The simulated deformations of platelet structures during their adhesion to the vessel wall for platelet stiffness of 2.5 (a) and 25 kPa (b). The effect of the platelet membrane stiffness on the platelet paused time (c). The paused time was 6.69±0.71 s (mean±SD) for the membrane stiffness of 25 kPa, which was about twice higher than the paused time of 3.15±0.69 s (mean±SD) for the membrane stiffness of 2.5 kPa. (Online version in colour.)

Effect of the number of GPIbα receptors expressed on the platelet membrane

The interaction between platelet glycoprotein (GP) Ib–IX–V complex and vWF is the first step of the haemostatic response to vessel injury. As resting platelets interact with vWF, binding of vWF to GPIbα initiates platelet activation [42]. Meanwhile, in platelet-type von Willebrand disease, mutations of GPIb functional receptors can compromise haemostasis by increasing the affinity for vWF [43,44]. Some studies demonstrated that abnormalities in the concentrations of GPIb membrane proteins are present in patients with myeloproliferative disorders. In particular, decreased GPIb concentrations were found in patients with thrombocythaemia and leukaemia [45,46]. How the platelet–substrate adhesion dynamics and subsequent platelet activation are affected by the number of GPIb is not clear. The objective of our simulations performed in this section was to gain insights into this problem. We varied the platelet receptor number from 10 688 (normal) to 5344 (insufficient), and performed simulations with 30 different random seeds to obtain 30 different paused times under flow shear stress of 3.0 dyn cm−2. The results of our simulations revealed that the paused time in the case of decreased receptor number was 2.07±0.41 s (mean±SD), which was significantly lower than 3.15±0.69 s (mean±SD) for the normal receptor number group (t-test, p<0.02, figure 9). Our simulations predicted that as the number of GPIb on the platelet membrane decreased, the paused time of platelet adhesion to vessel wall also decreased. Thus, the results of our model suggest that the number of functional GPIb is an important factor determining platelet adhesion and subsequent activation. This has important biological consequences, as controlling the number of functional GPIb receptors can provide means for development of novel anti-thrombotic drugs. The mechanism of these drugs is based on inhibiting/promoting the function of platelet GPIb receptors to decrease/increase adhesion of platelets to vWF to control blood clot growth [47].

Figure 9.

The effect of the number of platelet receptors on the platelet–vessel wall paused time. The platelet paused time for a decreased number of GPIb functional receptors was 2.07±0.41 s (mean±SD), which was significantly lower than the paused time of platelets having the normal number of receptors (3.15± 0.69 s, mean±SD). (Online version in colour.)

Effect of the platelet–platelet adhesion

To study how platelet–platelet interaction affects platelet adhesion to the blood vessel wall, we modelled dynamics of two platelets near the surface of the vessel (figure 10a). In the model, the two platelets interacted with each other and one of them adhered to the vessel wall. Our simulations revealed that the platelet paused time was 1.61±0.46 s (mean±SD) in the case of two adhesive platelets, which was significantly lower than the pause time of 3.15±0.69 s (mean±SD) calculated for a single platelet interacted with the wall (t-test, p<0.02, figure 10b). These results indicate an important mechanism by which a single platelet adhesion can be affected owing to interaction with neighbouring cells. These findings have direct biological consequences and help to explain how the increased platelet concentration in blood can affect platelet–wall adherence.

Figure 10.

(a) Initial configuration of two platelets used in simulations studying the effect of mutual interaction of platelets on platelet–wall adhesion. (b) Platelet–vessel wall paused time as a function of the number of interacting platelets. The platelet paused time was 1.61±0.46 s (mean±SD) for two adhesive platelets, which was significantly lower than the stopping time of platelets having the normal number of receptors (3.15±0.69 s, mean±SD). (Online version in colour.)

Discussion

This paper described a novel three-dimensional model coupling processes at three biologically important spatial scales critical for early blood clot development and uses models to provide predictive simulations. First, our model provides a comprehensive representation of mechanical properties of a platelet based on the implementation of a hybrid membrane submodel to describe mechanical behaviour of the cytoskeleton network and the lipid bilayer of the platelet. In previous studies, platelets were modelled as rigid bodies [24,33,48]. However, it has been experimentally shown [11] that platelets exhibited both elastic and viscoelastic behaviour and that they undergo large deformation in shear flow [12].

Experimental studies demonstrated [1,49,50] that flow shear stress could increase both bond formation and dissociation rates during platelet adhesion to the vessel wall. Additionally, estimates for the forces acting on platelet–substrate bonds were provided in [10]. However, Xu et al. [8] did not describe a detailed computational model to simulate the binding dynamics under various flow conditions. By combining three-dimensional multi-scale models with microfluidic experiments, we provided a methodology to quantify in detail single platelet-flipping in blood flow and platelet tethering to the injured vessel wall. It results in a two-way coupled fluid–cell interaction submodel combined with a stochastic submodel of formation/breakage of individual receptor–ligand bonds. This approach provided a biologically justified description of platelet dynamics, which can be used to simulate dynamics of platelets under more complex flow conditions.

By incorporating physiological parameter values characterizing cellular membrane mechanics, our method provides explicit representation for the structure of the cytoskeleton and simulation of cellular dynamics. Thus, our model allows one to examine how the mobility of cells is affected by their membrane structural and mechanical properties and, hence, aids in providing prognostic assessment in blood cell disorders outcome. The model developed in this paper can be also used for simulating important biomedical problems which involve description of dynamics and deformation of cells in fluid flow, including (patho)physiological inflammation involving leucocyte and platelet tethering to the vessel wall. Other important applications of the model include studying cell aggregate formation in blood, metastasis of tumour cells as well as stem cell attachment to the target tissues.

Table 3.

Execution time of 10 000 simulation steps on CPU and GPU for different fluid grid sizes.

fluid grid size	CPU (s)	GPU (s)
80×300×20	36 251	268
40×150×10	2393	231
20×75×5	232	51