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.

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 Effects of shear flow on accumulation of platelets on various surfaces have been extensively studied in in vitro and in vivo experiments [1][2][3][4][5][6][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 k on and k off , 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 plateletflipping 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 plateletsubstrate 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 cellsubstrate 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). 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.

(a) 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 x 2 /a 2 + y 2 /a 2 + z 2 /(λ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 = γρa 2 /μ = 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].
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 {X i , 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, H SCE is the in-plane energy of the cytoskeleton network; H bending is the bending energy representing the bending resistance of the lipid bilayer; H volume and H area are volume and area conservation constraints, respectively, and H wall 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 where L ij is the rest length, R ij = X j − X i 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 H SCE = i,j for all edges U e ij . 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 s , k t and k v are the global area, local area and volume constraint coefficients, respectively. The terms S total and V denote the surface area and volume for the whole platelet, whereas S 0 , 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 Σ ∈ R 3 be a smooth, closed surface representing the lipid bilayer of the platelet. The bending energy of the lipid bilayer is defined as [27] H bending = k 0 where K(x) = 1/2(κ 1 (x) + κ 2 (x)) is the mean curvature, and κ 1 (x), κ 2 (x) are the principle curvatures at the point x. We follow the finite-element method in [28] to calculate κ 1 (x) and κ 2 (x). 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 i is the value of u at node i and N i (ξ , η) 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, A. 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 X,X are its positions in the unstressed state and after deformation at time t, respectively.n 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 Thus, nodal forces F i are derived from the total energy as follows

(b) 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 plateletplatelet 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]. 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 P f (probability of forward reaction) and P r (probability of reverse reaction) described in [36]: where k off and k on 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 k off (F b ) is the bond dissociation rate, k 0 off is the unstressed off-rate, γ is the reactive compliance, F b is the applied force on the bond and k B T is the product of Boltzmann constant and temperature. The dependence of bond formation rate constant k on on the deviation bond length is described by [29,33] as where k 0 on is the intrinsic cross-linking formation rate constant, σ is the spring constant, l b is the equilibrium bond length, x b 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, l b is the equilibrium bond length, x b 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.

(c) 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-        where F denotes an external body force, ∇ x,p 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 ρ(x, t) and the velocity field v(x, t). In the D3Q19 lattice model, 19 values f i (x, t) are stored at each lattice site assigned to a lattice vector c i . The local density at a lattice point is obtained by summing all f i , (3.17) and the streaming velocity is given by where c i = h/ t is the lattice speed associated with the ith direction and t is the time step of our simulation. Using a Chapman Enskog expansion, Guo et al. [37] showed that the following lattice Boltzmann equations give a second-order accurate v, the Navier-Stokes velocity in the presence of a spatially varying, time-dependent force where u is a streaming velocity defined in equation (3.18), v = u + F t/2ρ, and  The pressure p = c 2 s ρ 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)h 2 /(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.  [39]. For instance, in the x-y boundary plane z = 0, f i (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 (i = 5, 9,13,15,17) are undetermined. Following the methods of Hecht & Harting [39], we obtain (3.27) and Here, v x , v y and v z are boundary velocities in x-, y-and z-directions, N z x and N z y the transverse momentum corrections on the z-boundary for distributions propagating in x-and y-directions, respectively and (d) Coupling platelet, dynamic adhesion and flow submodels (i) 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 H wall of equation (3.1) represents the energy associated with these interactions. In particular, ∂(H wall )/∂x i 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: where F 0 = 500 p Nm, τ = 2000 µm −1 and ε is the separation distance between platelet membrane and vessel wall [29]. Thus, ∂H wall /∂x i term in equation (3.11) is defined as

(ii) 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 (3.31) and and where h is the fluid node spacing, x ijk = (ih, jh, kh) the coordinate of the i, j, kth Eulerian grid node, X m the Lagrange coordinate of the mth elements, f ijk the force density on x ijk , F m the force density on X m , u ijk the velocity of x ijk , U k the velocity of X m . The discrete delta function δ h appearing in equations (3.33) and (3.34) is a smoothed approximation to the Dirac delta function δ(r). (The detailed derivation procedures in several forms have been presented in the literature [40].) We use the following common form and

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 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 doublelayer-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.

(b) 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 k off rates to compare with available experimental data. The model parameters used in our simulations (table 2) were obtained in biological experiments [1,11,[22][23][24][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 (k off ) 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 −k off ).
To calculate dissociation constants, we performed simulations for various numbers of random seeds (1000-1200) 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 k off . Our simulations also demonstrated that the k off 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)).

(c) 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.

(i) 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.

(ii) 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  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].

(iii) 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 plateletwall adherence.

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 fluidcell 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.

Appendix A. GPU implementation
The CPU used for our simulation is the Intel Xeon CPU L5520 with clock rate 2.27 GHz. Our NVIDIA graphics card is the GeForce GTX 480, Clock rate: 1.45 GHz, CUDA driver version: 5.50, CUDA runtime version: 5.50, CUDA capability version: 2.0. GPUs are separate devices with their own processors, and memory devices which do not have direct access to the CPUs or CPU memory units. The communication pathway for transferring data between CPU memory and GPU memory has a relatively slow bandwidth capability compared with direct access to memory devices. Thus, it is necessary to minimize communication as much as possible. The typical GPU code is composed of three main parts: (i) initialization, (ii) execution, and (iii) cleanup. During initialization, model data are firstly allocated and initialized in the CPU memory. The CPU code then initializes connection to GPU device and allocates GPU memory for the model simulation data. The model data are copied from CPU memory to GPU memory units. During execution, GPU kernel functions are called and occasionally copy data between CPU and GPU memory devices. When a simulation is finished, both CPU and GPU memory units are freed and connection to GPU device is shut down for the clean-up. Figure 1 shows the flow chart of our simulation algorithm. For a single step of the simulation, its execution on the GPU starts with the hybrid membrane model. Assuming that a platelet consists of N triangle mesh elements and P nodes (each node represents an SCE; figure 2b), The GPU kernel function to calculate forces in equation (  where blocksPerGrid and threadsPerBlock are determined according to the block and thread distribution on the GPU card, and P = blocksPerGrid*threadsPerBlock. If this is implemented on a single CPU in serial configuration, there will be a total of P iterations for one step of simulation. In our GPU implementation, this simulation is performed simultaneously on P GPU threads. Hence, it reduces the complexity of execution time from O(P) to O(1) for one step of simulation. Similarly, we use the following GPU functions to calculate forces owing to bending, area constraint and volume constraint sem_Bending_kernel<<< blocksPerGrid, threadsPerBlock >>>(grids->devImage){ //Calculation of forces due to bending . . . . . . The forces acting on solid nodes spread to the neighbour fluid nodes using the IBM. The GPU kernel function has the form fluid3d_force_distribute_kernel<<<blocksPerGrid3D,threadsPerBlock3D>>>(grids->devImage){ //Implementation of IBM . . . . . . } where both blocksPerGrid3D and threadsPerBlock3D have three-dimensional structure similar to a three-dimensional space coordinate. Let blocks per grid 3D = (x b , y b , z b ) and threads per block 3D = (x t , y t , z t ), and the fluid lattice has the size of X × Y × Z, then X = x b × x t , Y = y b × y t and Z = z b × z t . It reduces the complexity of execution time from O(XYZ) to O(1) comparing with CPU code. Similar strategy is applied to lattice Boltzmann method implementation.
There are M receptors in each of the N triangle mesh elements of a platelet. The GPU kernel function for dynamic adhesion model has the form: simulation on both CPU and GPU for three different fluid grid sizes. As table 3 shows, the execution time of GPU code is only about 1/100 of the CPU version for the fluid grid size we used in this paper. While the CPU execution time increased linearly with fluid grid size, the GPU execution time increased much slower than CPU.