Shear induced diffusion of platelets revisited

The transport of platelets in blood is commonly assumed to obey an advection-diffusion equation with a diffusion constant given by the so-called Zydney-Colton theory. Here we reconsider this hypothesis based on experimental observations and numerical simulations including a fully resolved suspension of red blood cells and platelets subject to a shear. We observe that the transport of platelets perpendicular to the flow can be characterized by a non-trivial distribution of velocities with and exponential decreasing bulk, followed by a power law tail. We conclude that such distribution of velocities leads to diffusion of platelets about two orders of magnitude higher than predicted by Zydney-Colton theory. We tested this distribution with a minimal stochastic model of platelets deposition to cover space and time scales similar to our experimental results, and confirm that the exponential-powerlaw distribution of velocities results in a coefficient of diffusion significantly larger than predicted by the Zydney-Colton theory.

1 Introduction

Platelets, or thrombocytes, are an essential blood constituent, from a physiological and heamodynamical point of view. Their motion is mainly a consequence of mechanical and hydrodynamic interactions with deformable red blood cells and the plasma, which makes an accurate description of their transport challenging.

Among the blood constituents, platelets are the second most numerous cells in blood, after red blood cells, with a concentration of 150–450 × 10⁹/L. They are involved in multiple physiological and pathophysiological processes such as haemostasis, thrombosis, clot retraction, vessel constriction and repair, inflammation, including promotion of atherosclerosis, host defense, and even tumor growth/metastasis (Harrison (2005)). When needed, platelets respond rapidly through activation, adhesion, aggregation, release of the materials stored in their granules, etc. Any disorder in these physiological processes results in impaired haemostasis and inappropriate thrombus formation. For example, arterial thrombi can develop within atherosclerotic lesions resulting in stroke and heart attack, two of the major causes of morbidity and mortality in the western world Harrison (2005).

From a physical point of view, platelets are small rigid suspensions interacting with other, larger, deformable suspensions, the red blood cells. It is well recognized that the physics of such systems is rich and complicated Eckstein and Belgacem (1991); Kumar and Graham (2012). Processes such as margination and segregation are observed in which platelets move towards the wall of the system, while red blood cells concentrate in the center of the vessel, producing the so-called cell free layer at the wall (Kim et al. (2009); Fedosov et al. (2010); Carboni et al. (2016)).

Adhesion and aggregation of platelets depend not only on their affinity with the endothelium or the deposition surface but also on their flow towards this surface. Platelets movement in the blood is affected by their interactions with red blood cells (D’Apolito et al. (2015)). In order to properly interpret platelet function tests and evaluate if platelets are dysfunctional in a patient, it is important to identify and separate the effect of transport properties to the wall from the intrinsic biochemical platelets properties (adhesion and aggregation). In several diseases, the shape of red blood cells is modified, affecting platelet motion Boudjeltia et al. (2020). When analyzing the way platelets deposit on a surface, their adhesion and aggregation rates can only be determined if their flow is correctly known.

In a blood flow subject to a shear rate, platelets experience an enhanced random motion in the direction perpendicular to the flow. The accepted description of this process, the so-called Zydney-Colton theory (Zydney and Colton (1988)), is that platelets are subject to a diffusion process, whose diffusion coefficient depends on the shear rate $\dot{\gamma }$ and the hematocrit H, the fraction of space occupied by red blood cells Affeld et al. (2013):

$D_{ZC}=D_{\mathit{PRP}}\left(1-H\right)+0.15d_{\mathit{RBC}}^2\left(H/4\right)\dot{\gamma }{\left(1-H\right)}^{1.8}(1)$

where D_PRP is the diffusivity of platelets in a platelet-rich plasma (i.e. without red blood cells) with a typical value of $D_{\mathit{PRP}}=\mathcal{O}(10^{-13})m^2{s}^{-1}$ and d_RBC is the typical diameter of a red blood cell. For H = 0.35 and $\dot{\gamma }=100s^{-1}$ the value of D_ZC is

$D_{ZC}=5\times 10^{-11}{m}^2{s}^{-1},(2)$

The present study was motivated by the determination of platelets adhesion and aggregation rates from in-vitro experiments with the Impact-R platelet analyzer device. It is made of a cylinder filled with whole blood and closed by two disks. The upper disk rotates to produce a prescribed shear flow and the lower disk is a deposition surface on which platelets can adhere and aggregate. In Chopard et al. (2017); Dutta et al. (2018, 2022) we show that platelet adhesion and aggregation rates can be computed from the in-vitro observations (deposition patterns in the Impact-R), by combining a mathematical model of platelets deposition and a machine learning technique. However, to explain the deposition pattern observed in our experiments, platelets must experience a large flow towards the deposition surface. More precisely, we found that with H = 0.35 and $\dot{\gamma }=100s^{-1}$, a diffusion coefficient of the order 10^–8 m²s⁻¹ is needed. This value is more than two orders of magnitude larger than the value predicted by Eq. 2.

The Zydney-Colton model has been extensively validated by numerous numerical studies in which red blood cells and platelets were resolved (Zhao and Shaqfeh (2011); Zhao et al. (2012); Reasor et al. (2013); Vahidkhah et al. (2014); Mehrabadi et al. (2015, 2016); Závodszky et al. (2019)). However, divergence from Zydney-Colton has also been observed and attributed to the presence of a drift term (Eckstein and Belgacem (1991); Crowl and Fogelson (2011); Kumar and Graham (2012)). While this leads to much higher diffusion than Zydney-Colton, it remains challenging to understand why such symmetry breaking occurs, at least with the geometry of the Impact-R.

In this context, we have developed a high-fidelity numerical blood flow solver, Palabos-npFEM, described and validated in Kotsalos et al. (2019, 2020). This computational framework is based on a lattice Boltzmann fluid solver with suspensions (red blood cells and platelets), whose deformable membranes are described with the finite element method (FEM). The fluid-blood cell interaction is computed with the immersed boundary method. Although not the purpose of the present discussion, we emphasize the fact that such numerical simulations are extremely demanding in computing resources, even for relatively small system sizes (much less than a cubic millimeter) and short periods (of the order of a second). In the present paper, we revisit platelets transport and propose a new characterization based on a probability distribution of their velocity. We estimate this distribution from high-fidelity blood flow simulations resolved at typical length and time scales similar to those considered in Müller et al. (2016); Závodszky et al. (2019), and we extrapolate platelets trajectory at much greater spatiotemporal scales using a surrogate stochastic model of platelets motion. We find that the resulting random walk process leads to a diffusion coefficient whose magnitude tends towards our in-vitro observations. We also argue that the traditional ways to compute the diffusion coefficient from the particle mean square displacement is very sensitive to finite size effects, hence resulting, for small systems, in a significant underestimation of the diffusion coefficient and questioning limitations of the Zydney-Colton theory at large spatiotemporal scales.

2 Methods

2.1 Palabos-npFEM framework

Palabos-npFEM is a computational framework for the simulation of blood flow with fully resolved constituents. The software computes the movement and deformation of red blood cells and platelets, and the complex interaction between them. The tool combines the lattice Boltzmann solver Palabos for the simulation of blood plasma (fluid phase), a finite element method (FEM) solver for the resolution of blood cells (solid phase), and an immersed boundary method (IBM) for the coupling of the two phases. Palabos-npFEM provides, on top of a CPU-only version, the option to simulate the deformable bodies on Graphic Processing Units (GPUs), thus the code is tailored for the fastest supercomputers Kotsalos et al. (2021).

In more details, the framework resolves blood cells like red blood cells and platelets individually (both trajectories and deformed state), including their detailed non-linear viscoelastic behavior and the complex interaction between them.

The fluid solver is based on the lattice Boltzmann method (LBM) and solves indirectly the weakly compressible Navier-Stokes equations. The solid solver is based on the nodal projective finite elements method (npFEM) Kotsalos et al. (2019), which offers an alternative way of describing elasticity. The npFEM framework is a mass-lumped linear finite element solver that resolves both the trajectories and deformations of the blood cells with high accuracy. The solver has the capability of capturing the rich and non-linear viscoelastic behaviour of red blood cells as shown and validated in Kotsalos et al. (2019). The platelets are simulated as nearly-rigid bodies by modifying the stiffness of the material. The implicit nature of the npFEM solver renders it capable of resolving extreme deformations with unconditional stability for arbitrary time steps. The fluid-solid interaction is realized by the immersed boundary method (IBM) and more specifically by the multi-direct forcing scheme proposed in Ota et al. (2012). The IBM imposes a no-slip boundary condition, so that each point of the surface and the ambient fluid moves with the same velocity.

Collisions between blood particles, whether red blood cells or platelets, are implemented through a repulsive force acting as a spring, when the surfaces delimiting two particles are getting too close to each other. In the current study, we employ the same parameters as reported in Kotsalos et al. (2019, 2020).

2.2 Probability distribution function of platelets absolute velocities

From the fully resolved simulations we have recorded the position y_i(t), perpendicular to the flow, of the ith platelet, at every time interval δt = 10^–5 s. We define the absolute velocity as:

$v_i\left(t\right)=\frac{\left|y_i\left(t+\delta t\right)-y_i\left(t\right)\right|}{\delta t}(3)$

and computed the mean absolute velocity v_moy:

$v_{moy}=\frac{1}{N}\sum _{i=i}^N\left(\frac{1}{n}\sum _{t=0}^{(n-1)\delta t}{v}_i(t)\right)(4)$

where N is the total number of platelets and n the number of iterations.

We find that the probability distribution function of the platelets absolute velocities P(v) is well described by the relation:

$P\left(v\right)=\left\{\begin{array}{cc}\hfill p_0\exp \left(-\lambda v\right)\hfill & \hfill \text{for}v\le v_{min}\hfill \\ \hfill p_0\exp \left(-\lambda v_{\mathit{min}}\right){\left(\frac{v}{v_{\mathit{min}}}\right)}^{-1-\alpha }\hfill & \hfill \text{for}v\ge v_{min}\hfill \end{array}\right.(5)$

where p₀, λ, v_min and α are parameters to be determined. The quantity v_min, a velocity threshold separating the exponentially decreasing bulk of the distribution to its heavy tail, and α were estimated with the python package powerlaw Alstott et al. (2014), which is the python implementation of the seminal work of Clauset et al. (2009). Note that there has to be a lower bound as power-law distributions diverge when the variable tends to zero. More generally, power laws are observed in the tail, and these distributions deviate from the power law exhibiting various behaviors below some threshold value of the measured variable Redner (1998); Ioannides and Skouras (2013). We obtain the following values:

$\alpha =3.8v_{min}=5\times 10^{-3}m{s}^{-1}{v}_{moy}=10^{-3}m{s}^{-1}(6)$

The values of p₀ and λ are obtained from the conditions

$1={\int }_0^{\infty }dvP\left(v\right)=p_0\exp \left(-\lambda v_{min}\right)\left[\frac{v_{min}}{\alpha }-\frac{1}{\lambda }+\frac{p_0}{\lambda }\right](7)$

and

$v_{moy}={\int }_0^{\infty }dvvP\left(v\right)=p_0\exp \left(-\lambda v_{min}\right)\left[\frac{v_{min}^2}{\alpha -1}-\frac{v_{min}}{\lambda }-\frac{1}{{\lambda }^2}\right](8)$

These two equations for p₀ and λ can be solved numerically, knowing α, v_min and v_moy. The solution is p₀ = 1015.24 and λ = 1017.36.

Here we assume that v can be arbitrarily large. This is of course not true in practice but adding to P(v) the numerically observed cut-off speed does not change our results because the probability of such high velocity decreases fast enough for the α value we have here. Similarly, imposing a low velocity threshold did not modify our results.

2.3 Stochastic model for platelets velocities

From the expression of P(v), one can sample the distribution and generate platelets velocities that mimic the behavior of the fully resolved blood flow simulation. We call this velocity sampling a stochastic model for platelets velocities. Concretely if r is a random number uniformly distributed in [0, 1], the value of v is sampled through the relation:

$r={\int }_0^{v}d{v}^{\prime }P\left(v^{\prime }\right)(9)$

The resulting velocity is:

$v=\left\{\begin{array}{cc}\hfill -\frac{1}{\lambda }\ln \left(1-\frac{\lambda r}{p_0}\right)\hfill & \hfill \text{if}r<\frac{p_0}{\lambda }\left(1-\exp \left(-\lambda v_{min}\right)\right)\hfill \\ \hfill v_{min}{\left[1-\frac{\alpha }{p_0\exp \left(-\lambda v_{min}\right)}\left(r-\frac{p_0}{\lambda }\left(1-\exp \left(-\lambda v_{min}\right)\right)\right)\right]}^{-1/\alpha }\hfill & \hfill \text{if}r>\frac{p_0}{\lambda }\left(1-\exp \left(-\lambda v_{min}\right)\right)\hfill \end{array}\right.(10)$

Of course, when generating a velocity from P(v) we should remember that here, v is the absolute velocity, without its sign. The data from the simulation show that both signs of the platelets velocities are equally likely. This is expected due to the symmetry of the problem with respect to the y-axis. Therefore the sign of v is chosen at random with a probability 1/2 for both directions, independently of the value v. This is obtained by generating a second random number r′.

For each stochastic simulation we used the following procedure. First, we set the initial positions y_i (0) by placing randomly N platelets on a segment of length L. Then at each iteration the velocities v_i were drawn from the distribution described previously and the platelets positions were updated accordingly 3. A platelet i with a position y_i ≤ 0 is removed from the system (deposited platelet), while the boundary y_i > L is reflective. Based on velocity auto-correlation measurements discussed below, we set δt = 0.5 ms.

3 Results and discussion

We have performed simulations of the system depicted in Figure 1A, for shear rates 100 s⁻¹ and an hematocrit of 35% as in the Impact-R experiments. We have analyzed platelets trajectories along the vertical direction, that is perpendicular to the flow. We consider these simulations for 1s of physical time and for a system of size 50μm × 50 μm horizontally, and of height L, with L ∈ {50, 100, 250, 500} μm in the vertical direction. Periodic boundary conditions are imposed horizontally along the x and z directions, while the y direction (or wall direction) is bounded by walls. The upper moving wall has a constant velocity V_mw so as to produce the desired shear flow. Although these simulations are significantly smaller in size and time than the actual Impact-R experiment, they require supercomputing capabilities.

FIGURE 1

FIGURE 1. (A) A typical configuration of our high-fidelity blood simulation, with deformable red blood cells and platelets (yellow particles) in suspension in a Newtonian flow subject to a shear rate. In this setup, fluid velocity increases linearly from 0 (at the bottom fixed wall) to V_mw (velocity of the top moving wall). The system has periodic boundary conditions along the x and z axis, while the y-direction (or wall direction) is bounded by walls. (B) Illustration of the process taking place in the Impact-R platelets analyzer. Under the action of a shear flow and their interaction with red blood cells, platelets move towards the bottom of the system where they deposit, forming aggregates displayed on the gray slab as dark clusters.

In what follows, we focus on shear rate 100s⁻¹ as in the Impact-R experiments, and on the smallest system of size L = 50 μm. The reason is two-fold. First this is the typical size that has been already considered in the literature for high-fidelity blood flow numerical simulation. Second, due to the small size, the flow reaches much faster a steady state on which reliable measurements can be considered.

A typical trajectory is shown in Figure 2A. Note that we consider the trajectory as long as the platelet stays far enough from the boundaries. Indeed, we observe in our high-fidelity simulations that the trajectory of the platelets is strongly affected by the presence of the walls (see Figure 2A). This is the so-called margination process that naturally brings and traps platelets close to walls. To obtain data not affected by the system boundaries, the region we consider goes from y = 10 μm to y = L − 10μm, which corresponds to discarding from the analysis a bottom and a top layer of about one red blood cell diameter. Note also that, due to the symmetry of the flow, the density of red blood cells is uniform along y so that the vertical platelets transport is independent of y as shown by the profiles of the velocity and its variance along the wall direction (Figure 2B). Such symmetry consideration does not hold anymore in different geometries, for instance in a tube, where the shear rate is not uniform.

FIGURE 2

FIGURE 2. (A) Sketch of the system with the region in which platelets trajectories are recorded. Trajectories y(t) of two representative platelets along the direction perpendicular to the flow. Platelets (gray) with trajectory inside the excluded region are not considered for analysis. The situation shown in this figure corresponds to a wall to wall distance of L = 50μm, hematocrit H = 0.35 and shear rate $\dot{\gamma }=100s^{-1}$. (B) Profiles of the velocity V_y (blue) and the variance of the velocity $\mathrm{\Delta }{V}_y^2$ (green) along the wall direction y. Velocity and its variance are computed in bins of 0.1 μm and bins with frequency lower than 0.1% are not represented. The total mean (resp. variance) in thick (resp. thin) dashed line is computed away from the wall (between y = 10 μm and y = 40 μm). (C) The mean square displacement MSD(t) measured within the system limits is shown, as well as the corresponding diffusion coefficient. The blue curve represents the MSD measured from the platelets that remain in the domain of interest during all the simulation time (in), while orange curve is computed from the trajectories of the platelets during their residency in the domain (in and out). The dashed straight line represents the MSD expected by Zydney-Colton theory.

From the values of y_i(t), it is traditional to compute the platelets mean square displacement MSD(t) as

$\mathrm{M}\mathrm{S}\mathrm{D}\left(t\right)={⟨{\left[y_i\left(t\right)-y_i\left(0\right)\right]}^2⟩}_i(11)$

where ⟨⋅⟩_i indicates an average over the platelets, in our case those that are still in the domain at time t. The diffusion coefficient is linked to the means square displacement, as D_MSD = MSD(t)/(2t) for 1D systems (Figure 2B).

We consider platelets that remain in the domain of interest during all the simulation time to compute a first quantity D_in. Further, we compute D_in&out from the trajectories of all the platelets until they leave the domain.

We obtain, for the present small system, D_in&out = 1.71 × 10⁻¹⁰m²s⁻¹ and D_in = 1.99 × 10⁻¹⁰m²s⁻¹ (Figure 2), which is compatible with previous numerical observations (see for instance Reasor et al. (2013)), but larger than the MSD(t) that would emerge from a diffusion constant given by D_ZC. However, as shown below, the determination of D directly from the platelets trajectories in such small spatial system is not accurate and we claim that it underestimates D. Therefore we will consider below a different procedure, relating the diffusion coefficient to the velocity probability distribution (Figure 3A and Method) and velocity auto-correlation function.

FIGURE 3

FIGURE 3. (A) Velocity probability distribution P(v) obtained from sampling the platelets trajectories at a time resolution δt = 10^–5 s (blue), and δt = 10^–3 s (orange). The straight lines represent an exponential decrease and the dashed one a powerlaw behavior, hence suggesting a velocity distribution described by Eq. 5. The distributions were computed with logarithmically spaced bins. The inset indicates how velocities are measured from trajectories samples at different values of δt. (B) Platelet velocity auto-correlation function $\mathrm{V}\mathrm{A}\mathrm{C}\mathrm{F}(\tau )={⟨v_i(t+\tau )v_i(t)⟩}_i$ as a function of the time delay τ. The fit of the data points gives a power law relation, VACF(τ) ∼ τ^−b, with b = 0.808.

Our goal with this stochastic model is to generate representative platelets trajectories, much faster than with the fully resolved blood flow simulation, and in any geometry, including scales and time span similar to the Impact-R experiment. To do so, we still need to determine the mean time between two random changes of v. This is obtained by measuring the velocity auto-correlation function ${⟨v_i(t)v_i(t+\tau )⟩}_i$ as a function of τ. The results obtained from the high-fidelity blood simulation are shown in Figure 3B. We observe that

$\mathrm{V}\mathrm{A}\mathrm{C}\mathrm{F}(\tau )=⟨v(t+\tau )v(t)⟩=a{\tau }^{-b}(12)$

with a = 0.017 × 10⁻³s (SD = 0.0015) and b = 0.808 (SD = 0.029), when τ is expressed in ms.

The fact that VACF(t) is not an exponential function indicates that platelets keep a memory of their velocity over a rather large time interval Δt. The effect of these long-range correlations is significant when sampling the trajectories at different rates. For δt = 10⁻³s, the velocity histogram captures only the exponential part of the distribution, suggesting that the tail of the distribution results from short-time correlated platelets displacements (see the inset in Figure 3A). To determine the value of Δt we compute the standard deviation σ of the distribution of values in the tail of VACF(τ) (gray box in Figure 3B). The value of τ above which measurement points match the value of σ is considered as the typical time Δt after which the velocity is randomized. Therefore we define

$\mathrm{\Delta }t=0.5\times 10^{-3}s(13)$

From the value of Δt and the distribution P(v), we can use the stochastic model to generate representative platelets trajectories at much larger spatial and temporal scales, since it requires very modest computational resources.

A typical trajectory is obtained as

$y\left(t+\mathrm{\Delta }t\right)=y\left(t\right)+v\left(t\right)\mathrm{\Delta }t(14)$

with

$v\left(t\right)\sim \pm P\left(v\right)(15)$

where the notation ∼±P(v) means that a new velocity is drawn from P(v) with an equally likely positive or negative sign.

Figure 4A compares the velocity distributions of the stochastic model, the high-fidelity model and the analytical expression. We measure a Pearson correlation coefficient 0.999 (p_val < 10^–64) between the analytical expression and the distribution of velocity of the stochastic model, and a correlation coefficient 0.995 (p_val < 10^–35) between the analytical expression and the distribution of velocity of the fully resolved model.

FIGURE 4

FIGURE 4. (A) Velocity probability distribution obtained from the fully resolved simulation, the analytical expression and the stochastic model. (B) Mean square displacement obtained with the stochastic model of trajectories, for different system sizes. The inset show of the corresponding diffusion coefficient is sensitive to the size of the system, for small systems. The red dashed line represents the diffusion coefficient estimated from the Impact-R experiments. (C) The number $\mathcal{N}$ of platelets that deposit as a function of time for a system of size L = 0.82 mm and initially 4,800 platelets uniformly distributed along L. The red dashed line represented the number of deposited platelets measured on the Impact-R device.

Figure 4B shows the mean square displacement MSD(t) corresponding to particles following this stochastic model. We investigate the finite size effects on the measure of the MSD(t) and the estimation of the diffusion coefficient D. For small system sizes, the MSD saturates, hence resulting in spurious estimation of the diffusion coefficient (inset Figure 4B). Note that the asymptotic value of D can only be found at large system sizes. Increasing only the simulation time is not sufficient as the saturation of the MSD is still present after 20s. For the largest system, 10⁻²m, the diffusion coefficient associated with MSD(t) is:

$D=0.67810^{-9}{m}^2{s}^{-1}(16)$

which is more than 65 times larger than D_zc. It is also about 4 times larger than the value obtained from a direct measurement of MSD(t) as shown in Figure 4B.

One can further use our stochastic model to simulate the number $\mathcal{N}$ of platelets that deposit within 20 s in the Impact-R. In the simulation, platelets whose trajectory reaches the deposition surface at y = 0 are counted as absorbed and removed from the system. Those reaching the upper boundary at y = L = 0.82 mm are bounced back. This process can be iterated for 20s with an initial number of 4,800 platelets uniformly distributed along L. The total number of platelets that reach the absorbing boundary y = 0 increases in time as $\sqrt{t}$, as predicted from the survival probability in a 1D diffusion-absorption process. After 20s we measure near 800 deposited platelets (see Figure 4C). This is however still significantly less than the 3,000 platelets observed to deposit in the Impact-R during this same time interval.

It is important to notice that our results are sensitive to the values α, v_min and v_moy measured from the velocity distribution given in Figure 3A. We first investigate the effect of the parameterα in Figures 5A–C. Changing the value of α only modifies the tail of the distribution, whereas the exponentially decreasing bulk of the distribution remains unchanged (see Figure 5A). We explored the set of following values α = 1.8, 2, 3, 4, and 5. This set of values includes the very peculiar Levy-flight regime that is characterized by infinite variance Viswanathan et al. (2008). We do not observe significant change, both for the number of platelets and the diffusion coefficient, for values of α above 2 (Figures 5B,C). Instead, the diffusion coefficient D increases as α decreases below 2, and we measure D = 4.4 10⁻⁹m²s⁻¹ and approximately 1300 deposited platelets. Our fully resolved simulation does not show such values for α, but we cannot exclude totally the existence of such a regime for simulated systems with size and time scales comparable to those of the Impact-R setup. Yet speculative, these results show that except for the Levy-flight regime, the parameter α does not modify significantly the transport properties of the platelets.

FIGURE 5

FIGURE 5. Velocity distribution, number of deposited platelets and finite size effect of the diffusion coefficient for different values of α (A–C) and different values of v_moy and v_min (D–F). In Figures (D–F), the orange curve represents our reference values of v_moy and v_min.

We then investigate the effect of changing v_moy and v_min. For a sake of simplicity, and also to keep the distribution of velocity smooth, we multiply concordantly v_moy and v_min by values between 0.5 and 3 (see Figure 5C). We were motivated to increase the mean velocity by the fact that the observed mean velocity in the fully resolved model is related to the actual time resolution of the model. Indeed, decreasing the time resolution results in an increase of the mean velocity. Typically, the mean velocity increases by a factor 2 when the time resolution is decreased by a factor 100. This implies that at the fluid resolution time δt = 10⁻⁷s, one could expect a platelets mean velocity twice greater than the one measured at δt = 10⁻⁵s. Simulations at lower time resolution can be envisioned for future work, however at a dramatic computational cost. According to the different velocity distributions (see Figure 5D), we observe that both the number of deposited platelets and diffusion coefficient increase with the values of v_moy and v_min. Interestingly, for the highest v_moy we tested, we find a diffusion coefficient similar to the one obtained with α < 2 (D = 6.11 10^–9 m²s⁻¹), but a number $\mathcal{N}$ of platelets twice greater. This suggests the predominant role of the platelets mean velocity to explain the number $\mathcal{N}$ observed in the Impact-R experiment.

4 Conclusion

This paper proposes a detailed analysis of the statistics of platelets velocities when subject to an imposed shear flow of red blood cells and plasma. Our motivation was to better understand the deviation between the transport properties and deposition of platelets predicted by the Zydney-Colton relation 1) and those inferred from our Impact-R experiments.

Using fully resolved blood flow numerical simulations, in which deformable red blood cells and platelets are in suspension in a shear rate flow created between two walls, we were able to reconstruct the probability distribution P(v) of the platelets velocities. Such simulations are extremely time consuming as they have to solve the flow, the deformability of red blood cells, and the interaction between the fluids and the blood cells, as well as the interaction between these blood cells. Due to this high computational requirement, our study is limited here to systems of size 50μm × 50μm × 50μm, for 1 s of physical time.

From the simulation data, we found that P(v) is made of an exponential part followed by a power law tail. We also determine the velocity auto-correlation function and its characteristic time scale Δt for velocity decorrelation. From P(v) and Δt one can propose a stochastic model to generate platelets trajectories that mimic actual trajectories. However, with this stochastic model one can consider much larger spatial and temporal scales. At such larger scales, we found that the platelets diffusion coefficient perpendicular to the flow direction is about 65 times larger than D_ZC, predicted by the Zydney-Colton theory. However, for small system sizes, the diffusion coefficient inferred from the evaluation of the MSD within the systems boundaries gives a value of the same order than predicted by the ZC relation. More in-silico experiments, yet at high computational cost, are critical to extend our findings in a broader parameter space, including higher shear rate values, and validate our description in terms of probability of velocities distribution and diffusion coefficient at larger scales.

Overall our results go in the direction of the experimental observation about the flux of platelets in the Impact-R device, which requires a much larger diffusion coefficient than D_ZC but still a diffusion constant at least one order of magnitude larger than what is found here. More analysis is still needed to clarify this discrepancy.

In view of the importance of a right characterization of platelets transport in clinical devices to correctly test platelets functionality, we hope that this study will stimulate more experimental and numerical work.

Data availability statement

The datasets generated and analyzed for this study can be found on Zenodo: https://zenodo.org/record/7028165.

Author contributions

CK, FR and BC contributed equally to this work. CK developed the code for the fully resolved blood Palabos-npFEM model and performed the simulations. FR analyzed the data and developed the code to analyze them, including the stochastic model. He contributed to the interpretation of the results, produced figures and wrote parts of the manuscript. JL contributed to the development of the Palabos-npFEM code and the interpretation of the results. RD contributed to the interpretation of the results. FD contributed to the experimental data acquisition. KZ brought the biomedical questions and knowledge, collected experimental data, interpreted the results and contributed to the design of the research and the writing of the paper. BC designed the research project proposed the interpretation of the results, contributed to the data analysis and wrote the first draft of the paper. All authors contributed to the writing and reading of the manuscript.

Funding

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 823712 (CompBioMed2 project). Furthermore, this work was supported by grants from the CHU Charleroi (Fonds pour la Chirurgie Cardiaque; Fonds de la Recherche Medicale en Hainaut) and by the Swiss PASC project “Virtual Physiological Blood: an HPC framework for blood flow simulations in vasculature and in medical devices”. Open access funding was provided by the University of Geneva.

Acknowledgments

We acknowledge support from 1) the Swiss National Supercomputing Center, project ID 323, Digital blood: a study of platelet transport in blood, 2) the National Supercomputing Center in the Netherlands (Surfsara, Cartesius supercomputer), and 3) the HPC Facilities of the University of Geneva (Baobab cluster). This work also got support from EUREKA Eurostar project 3DPlts.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Affeld, K., Goubergrits, L., Watanabe, N., and Kertzscher, U. (2013). Numerical and experimental evaluation of platelet deposition to collagen coated surface at low shear rates. J. Biomech. 46, 430–436. doi:10.1016/j.jbiomech.2012.10.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Alstott, J., Bullmore, E., and Plenz, D. (2014). Powerlaw: A python package for analysis of heavy-tailed distributions. PLoS ONE 9, e85777. doi:10.1371/journal.pone.0085777

PubMed Abstract | CrossRef Full Text | Google Scholar

Boudjeltia, K. Z., Kotsalos, C., Ribeiro, D., Rousseau, A., Lelubre, C., Sartenaer, O., et al. (2020). Spherization of red blood cells and platelet margination in COPD patients. Ann. N. Y. Acad. Sci. 1485, 71–82. doi:10.1111/nyas.14489

PubMed Abstract | CrossRef Full Text | Google Scholar

Carboni, E. J., Bognet, B. H., Bouchillon, G. M., Kadilak, A. L., Shor, L. M., Ward, M. D., et al. (2016). Direct tracking of particles and quantification of margination in blood flow. Biophys. J. 111, 1487–1495. doi:10.1016/j.bpj.2016.08.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Chopard, B., de Sousa, D. R., Lätt, J., Mountrakis, L., Dubois, F., Yourassowsky, C., et al. (2017). A physical description of the adhesion and aggregation of platelets. R. Soc. Open Sci. 4, 170219. doi:10.1098/rsos.170219

PubMed Abstract | CrossRef Full Text | Google Scholar

Clauset, A., Shalizi, C. R., and Newman, M. E. J. (2009). Power-law distributions in empirical data. SIAM Rev. Soc. Ind. Appl. Math. 51, 661–703. doi:10.1137/070710111

CrossRef Full Text | Google Scholar

Crowl, L., and Fogelson, A. L. (2011). Analysis of mechanisms for platelet near-wall excess under arterial blood flow conditions. J. Fluid Mech. 676, 348–375. doi:10.1017/jfm.2011.54

CrossRef Full Text | Google Scholar

D’Apolito, R., Tomaiuolo, G., Taraballi, F., Minardi, S., Kirui, D., Liu, X., et al. (2015). Red blood cells affect the margination of microparticles in synthetic microcapillaries and intravital microcirculation as a function of their size and shape. J. Control. Release 217, 263–272. doi:10.1016/j.jconrel.2015.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Dutta, R., Chopard, B., Lätt, J., Dubois, F., Zouaoui Boudjeltia, K., and Mira, A. (2018). Parameter estimation of platelets deposition: Approximate bayesian computation with high performance computing. Front. Physiol. 9, 1128. doi:10.3389/fphys.2018.01128

PubMed Abstract | CrossRef Full Text | Google Scholar

Dutta, R., Zouaoui-Boudjeltia, K., Kotsalos, C., Rousseau, A., de Sousa, D. R., Desmetand, J. M., et al. (2022). Personalized pathology test for cardio-vascular disease: Approximate bayesian computation with discriminative summary statistics learning. PLoS Comput. Biol. 18, e1009910. doi:10.1371/journal.pcbi.1009910

PubMed Abstract | CrossRef Full Text | Google Scholar

Eckstein, E. C., and Belgacem, F. (1991). Model of platelet transport in flowing blood with drift and diffusion terms. Biophys. J. 60, 53–69. doi:10.1016/S0006-3495(91)82030-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Fedosov, D. A., Caswell, B., Popel, A. S., and Karniadakis, G. E. (2010). Blood flow and cell-free layer in microvessels. Microcirculation 17, 615–628. doi:10.1111/j.1549-8719.2010.00056.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, P. (2005). Platelet function analysis. Blood Rev. 19, 111–123. doi:10.1016/j.blre.2004.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ioannides, Y., and Skouras, S. (2013). Us city size distribution: Robustly pareto, but only in the tail. J. Urban Econ. 73, 18–29. doi:10.1016/j.jue.2012.06.005

CrossRef Full Text | Google Scholar

Kim, S., Ong, P. K., Yalcin, O., Intaglietta, M., and Johnson, P. C. (2009). The cell-free layer in microvascular blood flow. Biorheology 46, 181–189. doi:10.3233/BIR-2009-0530

PubMed Abstract | CrossRef Full Text | Google Scholar

Kotsalos, C., Latt, J., Beny, J., and Chopard, B. (2020). Digital blood in massively parallel cpu/gpu systems for the study of platelet transport. Interface Focus 11, 20190116. doi:10.1098/rsfs.2019.0116

PubMed Abstract | CrossRef Full Text | Google Scholar

Kotsalos, C., Latt, J., and Chopard, B. (2019). Bridging the computational gap between mesoscopic and continuum modeling of red blood cells for fully resolved blood flow. J. Comput. Phys. 398, 108905. doi:10.1016/j.jcp.2019.108905

CrossRef Full Text | Google Scholar

Kotsalos, C., Latt, J., and Chopard, B. (2021). Palabos-npfem: Software for the simulation of cellular blood flow (digital blood). J. Open Res. Softw. 9, 16–12. doi:10.5334/JORS.343

CrossRef Full Text | Google Scholar

Kumar, A., and Graham, M. D. (2012). Margination and segregation in confined flows of blood and other multicomponent suspensions. Soft Matter 8, 10536–10548. doi:10.1039/c2sm25943e

CrossRef Full Text | Google Scholar

Mehrabadi, M., Ku, D. N., and Aidun, C. K. (2015). A continuum model for platelet transport in flowing blood based on direct numerical simulations of cellular blood flow. Ann. Biomed. Eng. 43, 1410–1421. doi:10.1007/s10439-014-1168-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehrabadi, M., Ku, D. N., and Aidun, C. K. (2016). Effects of shear rate, confinement, and particle parameters on margination in blood flow. Phys. Rev. E 93, 023109. doi:10.1103/PhysRevE.93.023109

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, K., Fedosov, D. A., and Gompper, G. (2016). Understanding particle margination in blood flow - a step toward optimized drug delivery systems. Med. Eng. Phys. 38, 2–10. doi:10.1016/j.medengphy.2015.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ota, K., Suzuki, K., and Inamuro, T. (2012). Lift generation by a two-dimensional symmetric flapping wing: Immersed boundary-lattice Boltzmann simulations. Fluid Dyn. Res. 44, 045504. doi:10.1088/0169-5983/44/4/045504

CrossRef Full Text | Google Scholar

Reasor, D. A., Mehrabadi, M., Ku, D. N., and Aidun, C. K. (2013). Determination of critical parameters in platelet margination. Ann. Biomed. Eng. 41, 238–249. doi:10.1007/s10439-012-0648-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Redner, S. (1998). How popular is your paper? An empirical study of the citation distribution. Eur. Phys. J. B 4, 131–134. doi:10.1007/s100510050359

CrossRef Full Text | Google Scholar

Vahidkhah, K., Diamond, S. L., and Bagchi, P. (2014). Platelet dynamics in three-dimensional simulation of whole blood. Biophys. J. 106, 2529–2540. doi:10.1016/j.bpj.2014.04.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Viswanathan, G., Raposo, E., and da Luz, M. (2008). Lévy flights and superdiffusion in the context of biological encounters and random searches. Phys. Life Rev. 5, 133–150. doi:10.1016/j.plrev.2008.03.002

CrossRef Full Text | Google Scholar

Závodszky, G., Rooij, B. V., Czaja, B., Azizi, V., Kanter, D. D., and Hoekstra, A. G. (2019). Red blood cell and platelet diffusivity and margination in the presence of cross-stream gradients in blood flows. Phys. Fluids 31, 031903. doi:10.1063/1.5085881

CrossRef Full Text | Google Scholar

Zhao, H., Shaqfeh, E. S., and Narsimhan, V. (2012). Shear-induced particle migration and margination in a cellular suspension. Phys. Fluids. 24, 011902. doi:10.1063/1.3677935

CrossRef Full Text | Google Scholar

Zhao, H., and Shaqfeh, E. S. (2011). Shear-induced platelet margination in a microchannel. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 83, 061924. doi:10.1103/PhysRevE.83.061924

PubMed Abstract | CrossRef Full Text | Google Scholar

Zydney, A. L., and Colton, C. K. (1988). Augmented solute transport in the shear flow of a concentrated suspension. PCH. Physicochem. Hydrodyn. 10, 77–96.

Google Scholar