Anomalous Platelet Transport & Fat-Tailed Distributions
Christos Kotsalos, Karim Zouaoui Boudjeltia, Ritabrata Dutta, Jonas Latt, Bastien Chopard
AAnomalous Platelet Transport & Fat-Tailed Distributions
Christos Kotsalos ∗ , Karim Zouaoui Boudjeltia , Ritabrata Dutta , Jonas Latt , and BastienChopard University of Geneva, Computer Science Department Universit´e Libre de Bruxelles and CHU de Charleroi, Laboratory of Experimental Medicine Warwick University, Department of Statistics
Abstract
The transport of platelets in blood is commonly assumed to obey an advection-diffusion equation. Here we propose adisruptive view, by showing that the random part of their velocity is governed by a fat-tailed probability distribution,usually referred to as a L´evy flight. Although for small spatio-temporal scales, it is hard to distinguish it fromthe generally accepted “red blood cell enhanced” Brownian motion, for larger systems this effect is dramatic asthe standard approach may underestimate the flux of platelets by several orders of magnitude, compromising inparticular the validity of current platelet function tests.
Keywords : platelets, anomalous transport, fat-tailed distributions, power law behaviour, cellular blood flow simu-lations, random walks, stochastic models
Platelets (PLTs) are the second most numerous cell in blood, after red blood cells (RBCs), with concentration of150-450 × /l . Their size, shape and material allow them to be optimally placed as close as possible to the vesselwall, a physical requirement for the constant inspection of the integrity of the vasculature. Platelets are polyvalententities and are involved in multiple physiological and pathophysiological processes such as haemostasis, thrombosis,clot retraction, vessel constriction and repair, inflammation including promotion of atherosclerosis, host defence,and even tumour growth/ metastasis [1]. Upon injury platelets respond rapidly (through activation, adhesion,aggregation, release reactions, etc.), and form a haemostatic plug, occluding the damaged site and preventing bloodloss. Any disorder in these physiological processes results in impaired haemostasis, and inappropriate thrombusformation. For example, arterial thrombi can develop within atherosclerotic lesions resulting in stroke and heartattack, two of the major causes of morbidity and mortality in the western world [1] † .Focal point of this analysis is the underlying mechanism of platelet transport, which has an impact on theirphysiological/pathophysiological behaviour. Platelets move within blood due to the combined effect of the plasmavelocity and the collisions with RBCs. In particular, in a shear flow, PLTs experience a random motion in thedirection perpendicular to the flow. The accepted description of this process (the so-called Zydney-Colton theory[2]) is that PLTs are subject to a diffusion process, whose diffusion coefficient is [3] D = D P RP × (1 − H ) + 0 . × ( d RBC × H/ × ˙ γ × (1 − H ) . (1)where D P RP is the diffusivity of PLTs in a platelet-rich plasma (without RBCs) and its values is typically D P RP = O (10 − ) m /s . The quantity ˙ γ is the shear rate, H is the hematocrit, and d RBC is the diameter of a RBC.Equation (1) gives that D = O (10 − ) m /s , for the situation described in Chopard et al. [4], namely the diffusionof PLTs in a shear flow with ˙ γ = 100 s − and H = 0 .
35, as created by the so-called impact-R PLT analyser. Thisdevice allows one to measure the amount of PLTs that deposit on a surface perpendicular to the flow. However, ∗ Corresponding author: [email protected] † a r X i v : . [ phy s i c s . c o m p - ph ] J un ased on such experimental evidence, Chopard et al showed that if one assumes that the concentration ρ of PLTsobeys the diffusion equation ∂ρ = D ∇ ρ (2)a value of D ≈ O (10 − ) m /s is needed to explain the number of platelets that is observed to deposit in impact-Rplatelet analyser.This result obviously raises the question of the validity of eq. (1) or the applicability of eq. (2). The Zydney-Coltonmodel results from accumulating experimental/theoretical data, and it has been extensively validated by numerousnumerical studies, in which RBCs and PLTs where resolved [5, 6, 7, 8, 9, 10, 11]. However, these studies concernspatio-temporal scales much smaller than those characterising the impact-R device. The latter considers a layer ofblood of 820 µm of thickness, rotating in a cylinder of diameter 6 . mm . The amount of platelets that disappearedfrom the bulk due to their deposition on the bottom part of the cylinder was observed after 20 s . State of the artnumerical simulations consider much smaller systems, usually less than 100 µm in size, for less than 1 s of physicaltime.A possibility to explain this 1000-fold difference between the Zydney-Colton theory and the effective observeddiffusion is to postulate the presence of a drift term in addition to the diffusion process, as suggested in [12, 5, 13], asa general mechanism for the case of transport of deformable suspensions. This drift-diffusion model, which includesa “rheological potential” (Φ), has however an ill-considered/posed origin. Furthermore, due to the symmetry ofthe problem it is hard to understand why such a symmetry breaking would appear in the case of the impact-R.As a matter of fact, simulations of PLT transport in the impact-R, including an ad hoc drift term and a Zydney-Colton diffusion, do not fit so well the in-vitro time evolution of PLTs deposition (data unpublished, B. Chopard2018).Our claim in the present paper is that it is possible to reconcile these contradictory results by assuming that PLTsdo not follow a Gaussian random walk (as is implied by the standard diffusion equation), but a random walkwith a fat-tailed distribution of velocities. We show here that a very careful analysis of fully resolved blood flowsimulations, in which deformable RBCs and platelets interact and move in a suspending fluid (the plasma) revealsthat collisions between RBCs and PLTs result in a power law probability distribution of PLT velocities, P ( v ) ∼ v − − α (3)with α around 1.5. We also show that for small systems, the value of the PLT mean square displacement (MSD),which is traditionally related to the diffusion coefficient, is compatible with the Zydney-Colton theory, explainingwhy the normal diffusion hypothesis was little questioned in the literature. But as the system size increases, weobserved that the MSD increases, a behaviour incompatible with a Gaussian random walk.It should be noted that, despite this wide consensus that Zydney-Colton theory applies, a few researchers havesuperficially pointed out events that support our current result. Vahidkhah et al. [9, 14] observed highly anisotropicRBC distribution and a “waterfall” phenomenon (cavities that act as express lanes for platelets) affecting PLTmovement (supported experimentally by Lee et al. [15]). Mehrabadi et al. [10] pointed out a possible anomalousdiffusion as PLTs get trapped in the cell free layer. Yeo et al. [16] (same approach as Gross et al. [17]) who talkedabout anomalous diffusion of wall-bounded non-colloidal suspensions, observed exponential distributions for thedensely packed spherical particles (correspondence to RBCs).It is critical to understand, why an underlying power law behaviour leads to enormous differences in transportphysics. In figure 1, we present five particles with the same average jump length (three of them exhibit power lawbehaviour, and two of them follow a Gaussian distribution). All particles start from the same point and are left toexplore the space for a thousand iterations. The particles performing L´evy flights exhibit jumps with probabilitydensity function ∼ x − − α and an average jump of size α/ ( α −
1) for α > µ = 0, σ = 1)and an average jump of size (cid:112) /π . We make sure that the average jump ( ¯ J ) is the same for every particle bymultiplying the individual jumps with a normalisation constant equal to ¯ J ∗ ( α − /α for the power laws, and equalto ¯ J ∗ (cid:112) π/
00 75 50 25 0 25 50 75 X Y Lévy Flight, Power Law =1.1, Avg Jump=1.0Lévy Flight, Power Law =1.5, Avg Jump=1.0Lévy Flight, Power Law =2.0, Avg Jump=1.0Gaussian Random Walk, Avg Jump=1.0Gaussian Random Walk, Avg Jump=1.0
Figure 1: Five particles performing random walks (three of them exhibit power law behaviour with varying exponent,and two of them follow a Gaussian distribution). The exploration of the space is largely dependent on the underlyingprobability density function.between the high fidelity blood flow simulations, which are extremely computationally demanding, and the PLTsdisplacement in the impact-R device. Our results are then presented in section 3. They consist in a careful statisticalanalysis of the random motion of the simulated PLTs, and a simulation of the platelet deposition process of theimpact-R, based on the inferred properties of the random walk. Finally, in section 4 we discuss some implicationsof our disrupting theory of PLTs motion.
Our goal is to understand whether the discrepancy found in Chopard et al [4] between the Zydney-Colton diffu-sion coefficient and the one accounting for the PLTs deposition pattern observed in the impact-R device can beexplained by an anomalous diffusion process, thus ruling out the possibility to use eq. (2). As discussed in theprevious section we made the hypothesis that PLTs may obey a fat-tailed distributed random walk instead of aGaussian one. Due to the great difficulty to measure directly the movement of platelets (individual trajectories) inwhole blood with an in-vitro approach, numerical simulations are considered. We refer these simulations as DNS(Direct Numerical Simulations) as they provide a high fidelity description of blood, integrating RBCs and PLTsas deformable suspensions in a flowing plasma. However, state of the art DNS blood solvers are still limited torather small spatio-temporal scales. For this reason we will have to perform a very careful statistical analysis of theplatelet trajectories to evidence their non Gaussian behaviour at the reachable scales.
Here we consider the Palabos-npFEM DNS blood solver [19, 20]. It offers high accuracy, flexibility and highscalability on the top fastest parallel supercomputers. This computational framework couples the lattice Boltzmannsolver Palabos [21] ∗ for the simulation of blood plasma (fluid phase), a finite element solver for the resolution of thedeformable blood cells (solid phase), and an immersed boundary method for the coupling of the two phases. Theframework resolves blood cells like RBCs and PLTs individually (both trajectories and deformed state), includingtheir detailed non-linear viscoelastic behaviour and the complex interaction between them.Collisions between blood particles, whether RBCs or PLTs, are implemented through a repulsive force acting asa spring, when the surfaces delimiting two particles are getting too close to each other. In the current study, we ∗ https://palabos.unige.ch y -axis is oriented vertically. Twohorizontal no-slip walls are positioned at locations y = 0 and y = L , with L ∈ { µm, µm, µm, µm } .A shear flow is produced in the z -direction by moving the upper wall at a proper velocity. The shear rates weconsider are ˙ γ = 100 s − and ˙ γ = 400 s − . Periodic boundary conditions are imposed along the x and z horizontalaxes. This setup is meant to approximate the impact-R geometry for which L = 820 µm and the x and z axesspan a window of size 1 mm × mm , embedded in a cylinder of diameter 6 . mm (see [4] for more details). Thehematocrit is H = 35% and the ratio RBCs/PLTs is around 20.The simulations are run for a time interval [ t , t ], which typically lasts 1 s of physical time. This interval is resolvedat the scale 10 − s , the time step of the npFEM-Palabos solver. The trajectories y i ( t ) along the y -axis of all PLTs i are recorded, based on the position of their centre of mass. From theses trajectories, one has to determine theprobability distribution of PLTs velocities. To properly sample the time series y i ( t ) of positions, one has to extractthe characteristic platelet mean free time, ∆ t , representing the average time between successive impacts/collisions.This characteristic time is important, because an under-sampling could result in missing important collision events,and thus misinterpreting the motion of platelets. On the other hand, a sampling at a too fine time scale will samplestatistically dependent velocities, before collisions could randomise the platelet movements.Multiple researchers define the sampling interval as ∼ ˙ γ − , but this formula has no robust explanation. Alter-natively, we consider the trajectories y i as seen at different time scale ∆ t , for several possible values of ∆ t ∈{ . ms, . ms, . ms, . ms } . We compute the area A formed graphically by the trajectory y i ( t ) at full DNSresolution, namely A = (cid:90) t t y i ( t ) dt Then we compute the area formed by the trajectory at scale ∆ t (i.e. y i ( t ) with t incremented by steps ∆ t ). Thesetwo areas should be identical if the platelet did not experience a collision during time ∆ t . The larger ∆ t for whichthis equality holds is chosen as the mean free time.Figure 2 presents this method. The vertical axis shows the relative deviation between these two areas, averagedover all PLTs. The inset shows the trajectory of one representative platelet at three different time scales. Wesee that ∆ t = 1 ms captures well the change in the trajectory due to collisions, while filtering out irrelevantoscillations.The set of independent platelet velocities is then obtained from the quantities (where we set t = 0 for simplic-ity) v i ( k ) = y i (( k + 1)∆ t ) − y i ( k ∆ t )∆ t (4)for all platelets i .A classical property characterising the movement of platelets is their mean square displacement (MSD), definedas MSD( t ) = (cid:10) ( y i ( t ) − y i (0)) (cid:11) (5)where the average is over all platelets i . This relation can also be written asMSD( t ) = (cid:42) t/ ∆ t (cid:88) k =1 v i ( k ) t/ ∆ t (cid:88) (cid:96) =1 v i ( (cid:96) ) (cid:43) (6)Assuming that v i ( k ) and v i ( (cid:96) ) are independent when k (cid:54) = (cid:96) , and that (cid:104) v i (cid:105) = 0 due to the symmetry along the y axis(fact confirmed numerically for the current simulation setup), we obtainMSD( t ) = t/ ∆ t (cid:88) k =1 (cid:10) v i ( k ) (cid:11) = t ∆ t (cid:104) v i (cid:105) (7)The second equality comes from the fact that in a steady state, the velocity distribution is expected to be independentof time. Note however that this hypothesis is only valid in the bulk of the sample, where RBCs are found to have a4onstant density along the y -axis. This is no longer the case in the so-called cell-free layer, near the walls at y = 0and y = L . For this reason, we separate the treatment of platelets according to their spatial location and determinethe velocity distribution for the PLTs in the homogeneous region, away from the walls.Equation (7) shows that MSD( t ) is expected to be proportional to t , but this is only the case if < v i > is finiteand well defined. We will show below that this is not the case for platelets in a shear flow. Actually the subtlepart of this observation is that, for any numerical simulation, < v i > is obviously finite as there is a finite numberof platelets in the system, for a finite number of time steps. But as the system size is increased, we observedthat (cid:104) v i (cid:105) also increases, indicating a divergence of the variance of the platelet velocities. The striking result ishowever that for small systems, namely those that are accessible to DNS blood flow simulations [5, 9, 10], thenon-converged computed value of (cid:104) v i (cid:105) is compatible with the Zydney-Colton diffusion coefficient, thus supportinga misinterpretation of a Gaussian random walk. D e v i a t i o n f r o m S i m u l a t e d P a t h ( % ) Sampling Interval (ms) P LT P o s i t i o n i n W a ll - b o und e d D i rec t i o n [ y (t)] time [t → t ] 𝐴𝑟𝑒𝑎 = න 𝑡 𝑡 𝑦 𝑡 𝑑𝑡 Figure 2: Deviation from the DNS platelet trajectories as a function of the time scale at which these trajectoriesare observed. The way the deviation is measured is illustrated in the inset, for three different sampling intervals.The black circles in the main figure show the average deviation, computed over all platelet trajectories. The “error”bars mark the minimum and maximum deviations. One observes that a sampling interval ∆ t = 10 . ms (dashedline in the inset) misses many collision events, but that ∆ t = 0 . ms is a proper choice. In this section we explain how we determined the properties of the probability distribution of PLT velocities sampledfrom the high fidelity blood flow simulations described in section 2.1.To avoid wall effects, we discard the velocities corresponding to PLTs that are less than one RBC diameter fromthe walls. In order to prove that PLT velocities are not normally distributed but have a rather high probabilityfor extreme events, we devised a set of methods originating from analyses of fat-tailed distributions for discerningand quantifying power law behaviour in empirical data (in-depth analysis by Clauset et al. [22]). A graphicalexplanation is provided in the Supplementary Material.
Family of Distributions
The term “Thick Tails” is often used to describe distributions with much higher kurtosis than the Gaussian one, and“Fat Tails” is reserved for both extreme thick tails or membership in the power law class. We avoid designations suchas “Heavy Tails” or “Long Tails” to keep ambiguity to a minimum [23]. In literature, “Heavy Tails” is commonly5inked to distributions which do not have all their moments finite/defined. The term sub-exponential distributionsis used for the ones that decay more slowly than an exponential, i.e. they are not exponentially bounded (populardistribution belonging to this class is the log-normal). To summarise, the degree of thick tailedness (ranking byseverity) is
Thick Tailed ⊃ Sub-exponential ⊃ Power Law/Fat-Tailed (Paretian) . Furthermore, the thick tailed andsub-exponential families have all their moments finite. We are particularly interested in the power law class todescribe PLT transport, because these distributions do not have all their moments finite/defined. This happensfor power law distributions whose tail decays like ∼ x − α − with α ≤ α ≤ α ≤
2, we focus on finding evidences supporting this null hypothesis.Commonly used distributions that exhibit power law behaviour (asymptotically) are the Pareto, Cauchy (half-Cauchy, for positive only data), L´evy, Dagum (or Burr Type III), Singh–Maddala (or Burr Type XII), Log-logistic(or Fisk), inverted Weibull (or Fr´echet) (see Supplementary Material: Power Law Distributions) [18]. Our approachis to check, for each DNS, the plausibility of these distributions, according to the fitting criteria presented below.Also we want to extract the asymptotic power law exponent and check how many of the distributions exhibit α ≤ Distribution Fitting
In practice, few phenomena obey power laws for all values of x (where x denotes PLT velocities). More often thepower laws apply only for values greater than some minimum value, x min . In such cases, we say that the tail ofthe distribution follows a power law [22]. Thereby, for every distribution mentioned above (power law or not), wefind this lower bound and perform the fitting at the tail. Regarding the remaining part of the data (body), thiscan be sufficiently described by the empirical distribution function (histogram - see Supplementary Material). Notehowever that this procedure does not apply when fitting the data to a Gaussian. In this case, the fit is done acrossthe whole range of velocities.Once this lower bound is known, every distribution is fitted on the tail using the Maximum Likelihood Estimate(MLE) method, and its plausibility is checked using the Kolmogorov-Smirnov (KS) test for goodness-of-fit. Given thesymmetry of the studied phenomenon (for the current simulation setup), we use the absolute value of the velocities(no discrimination between upward and downward motions of PLTs). Thus all the statistics are one-tailed. If thep-value of the KS test is greater than the significance level (10% throughout the study for goodness-of-fit tests),then the investigated distribution could be a plausible fit to the data.In addition to this goodness-of-fit test, we employ the Log-likelihood ratio (LLR) test [22] comparing the fitteddistribution with every alternative distribution proposed above (see Family of Distributions). In that case, thesmaller the p-value of the LLR test, the more confident we can be on which distribution is a good fit of the data.The significance level for the LLR test is set to 1%. Note that, to confirm or deny a distribution, a small p-value is“good” for the LLR test (shows how trustworthy is the test’s result), while it is “bad” for the KS test (shows thatthe distribution is a poor fit to the data). Thus, for the acceptance/plausibility of a distribution fitted on the tail,there exist two criteria to meet, i.e. the KS and LLR tests.The estimation of x min is an optimisation process, and is based on minimising the distance between the inves-tigated model and the empirical data [22]. The metric d , that quantifies this distance, is the widely used Kol-mogorov–Smirnov statistic, which is simply the maximum distance between the Cumulative Distribution Function(CDF) of the sampled DNS data (the empirical CDF, noted S ), and the CDF P of the fitted model. Thus d = max x ≥ x min | S ( x ) − P ( x ) | (8)Therefore, the selected x min is the one that minimises this distance. Keep in mind that for every x tested, theparameters of P are fitted using the MLE method, which requires large enough samples. Thus, the minimisationprocess described above stops when the remaining tail has less than 100 velocities.6 ridging the scale: a Random Walk description As indicated previously, we want to show that the unexpectedly high transport of PLTs observed by Chopard etal. [4] in the in-vitro impact-R device is compatible with the velocity probability distribution extracted from ourDNS blood simulation. A direct verification is not possible as the spatio-temporal scales corresponding to theimpact-R is still too hard to reach, even on the fastest supercomputers: a cubic millimetre of blood simulated for20 seconds.A solution to bridge this gap of scales is to disregard the detailed movement of RBCs at the level of the impact-R,and only consider the dynamics of PLTs in terms of a random walk, using the velocity distribution obtained fromthe statistical analysis of the fully resolved PLT trajectories in the DNS simulations. Therefore we will simulatethe PLTs deposition process taking place in the impact-R through a stochastic model implementing the determinedrandom motion of PLTs. The question is to check whether this mesoscopic transport process reproduces the numberof PLTs observed to deposit at the bottom surface of the impact-R, after 20 s of operation. More precisely, theexperimental data [4] show that about 3000 activated platelets per micro-litre of whole blood have disappearedfrom the bulk within these first 20 s .This method follows the approach of Chopard et al. [4], by replacing the 1D diffusion equations describing the bulkof the impact-R device with actual random walks of point particles. The PLTs that cross the lower boundary areconsidered as deposited. Therefore, for every DNS, we perform the statistical analysis and, for every candidate of thePLTs velocity distribution, we simulate the corresponding stochastic model and record the number of deposited PLTsit predicts after 20 s . For more details on generating data for the random walks, see Supplementary Material. Our analysis builds upon the research by Chopard et al. [4]. In more details, in the in-vitro part of their studythe researchers used the impact-R platelet function analyser to study the evolution of PLT deposition (adhesion-aggregation processes) on the substrate of the device. Impact-R is a cylindrical apparatus (820 µm height), whoselower end is a fixed disc (deposition surface, 132.7 mm ), and its upper wall a rotating disc. Due to the rotatingupper wall, the blood is subject to a pure shear flow. The imposed constant shear rate was 100 s − (inside anobservation window of 1 × mm ), the blood was extracted from seven healthy donors. The differential role ofactivated and non-activated platelets was analysed, as well as the role of albumin in the deposition process. Thein-silico counterpart of their study consisted of 1D diffusion equations describing the movement of activated platelets(AP) and non-activated platelets (NAP) in the bulk of the device, and of stochastic rules for PLT deposition on thesubstrate. The study revealed that the Zydney-Colton [2] shear induced diffusion coefficient ( D ) was significantlytoo small to explain the observed deposition rate.The fully resolved 3D cellular blood flow simulations (DNS) provide a great amount of information, that simplifiedmechanistic models or in-vivo/vitro experiments cannot provide. Numerically following individual particles andresolving the complex interactions, helped us develop an alternative theory on how platelets are transported. Wetried to reproduce numerically the in-vitro experiments performed in Chopard et al. [4], to the extent possible,given the high computational cost. To reduce the computational demand, we designed our numerical simulations inchannels with lateral dimensions of 50 µm , while resolving the wall-bounded direction at {
50, 100, 250, 500 } µm .Consider that the in-vitro counterpart consists of 1000 µm in the lateral directions, and 820 µm in the wall-boundeddirection. Numerically, we generated a constant shear rate flow regime at 100 s − , realised with the help of a movingtop wall and a fixed bottom wall. Periodic boundaries were applied in the flow and vorticity directions, and thehematocrit was 35%, as in the experiments.Regarding the platelet size and shape, we considered experiments with either activated or non-activated PLTs. TheNAP are simulated as nearly-rigid oblate ellipsoids with diameters { } µm , thicknesses { } µm , andvolumes { } f l , respectively. The AP are simulated as nearly-rigid spheres with diameters { } µm (covering the uncertainty from the complicated shape transformations), and volumes {
12, 22, 30, 60 } f l ,respectively. Upon PLT activation, negatively charged phospholipids are translocated from the inner membrane tothe external surface, leading to a more negatively charged PLT. This complex electrochemical behaviour is quantifiedby the electrophoretic mobility of platelets [25, 26, 27]. Additionally to the change of the electrophoretic mobility,there is a severe change in shape with the appearance of blebs and pseudopods, which increases the hydrodynamicvolume of AP. While we address the latter alteration through the spherical shape, the complex electrochemical7ehaviour is roughly resolved through an increase in the intensity of the collision potential between activated PLTsand RBCs, i.e. increased repulsive forces ranging from 5 to 10 times compared to the ones of NAP.As for the RBCs, the normal biconcave shape is used in the majority of the experiments, while in few of them weintroduced spherised RBCs emulating pathological conditions [28], e.g. diabetes, chronic obstructive pulmonarydisease.We performed 64 simulations on Piz Daint, the flagship system of the Swiss National Supercomputing Centre, whichis the fastest supercomputer in Europe and the 6 th worldwide ∗∗ . These simulations include for completeness, 5 casestudies in a tube of 50 µm diameter, 5 case studies at 20% hematocrit, and 6 case studies with constant shear rateat 400 s − . Forty-one simulations follow the “exact” same setup as the experiments of Chopard et al. [4], but all64 present qualitatively the same platelet transport behaviour. Out of the 41 simulations, 22 deal with AP (12 ofthem include repulsive forces), and the rest deal with NAP. The graphs and results presented below include these41 simulations. Anomalous Transport From the Velocity Distribution
A standard diffusive process is unaffected by the size of the system (at least away from the walls), i.e. the moments(e.g. variance, kurtosis) of the PLT velocities converge as the sample size increases. Figure 3 presents the diffusioncoefficient D of PLTs (AP/NAP) when extracted from the mean squared displacement, for channels of varyingsizes. Traditionally, the diffusion coefficient is linked to the MSD, as D MSD = MSD( t ) / (2 t ) for 1D systems, e.g.along the wall-bounded direction. The black square points in Figure 3 indicate the value of D MSD as a function ofthe system size, averaged over all simulations. The “error” bars denote the minimum/maximum coefficients amongthe various DNS.In Figure 3 we also evaluate the PLT diffusivity D Gauss (empty squares) when a Gaussian distribution is (poorly)fitted on the observed PLT velocities. The D Gauss is computed from the standard deviation of the correspondingnormal distribution. As expected, the latter approach returns a diffusion coefficient that is not affected by thesystem size (small variations are due to different parameters in the DNS). However, for the diffusivity emergingfrom the MSD, there is a diverging behaviour as the sample sizes increase, a signature of a fat-tailed velocity distri-bution. On top of this observation, we show in the inset of Figure 3 the mean excess kurtosis (fourth standardisedmoment) of platelet velocities, which presents a diverging behaviour as well, and values that are way higher thanthe corresponding null value of a normal distribution.Alongside with observing the moments of the data, we performed a number of normality tests to check if a Gaussiandistribution is a plausible model for the data. We deployed the Shapiro-Wilk, D’Agostino’s K and Anderson-Darling normality tests to check this hypothesis. As expected, every normality test consistently rejected thishypothesis.Summarising, the extremely high kurtosis ( (cid:29) O (10 − ) m /s . As well accepted, our simulations confirm the role of PLTs-RBCs collisions to enhance PLTs transport, but a further analysis implies that the diffusive process is anomalous.Possibly, the suggestion for extra drift terms or rheological potentials [12, 5, 13] comes from a misinterpretation ofthis anomalous diffusive process. Power Law Emergence
From the statistical analysis of the sampled DNS output, we tried to find evidences that PLT velocities follow fat-tailed distributions, and more specifically power laws with exponent α ≤
2. Out of the 41 numerical experimentsonly 2 of them did not show evidence of power laws, i.e. no valid fitting of power laws on the data was obtained.For the rest, the tails of PLT velocities can be described with distributions which asymptotically behave as powerlaws (see Family of Distributions in section 2, and Supplementary Material). ∗∗ .0E-111.2E-102.3E-103.4E-104.5E-105.6E-106.7E-107.8E-10 D i ff u s i o n C o e ff i c i e n t (f r o m M S D ) [ m / s ] System Size along the Wall-bounded Direction [µm]
PLT Diffusion Coefficient (MSD DNS) PLT Diffusion Coefficient (Gaussian Random Walk)0100200300400 E x ce ss K u r t o s i s PLT DNS velocities Kurtosis
Figure 3: PLT diffusion coefficient D (black squares) estimated from the means square displacement of PLTs,computed from the DNS. The values of D are given as a function of the size L along the y axis obtained by anaverage over all the simulations. The “error” bars denote the minimum/maximum coefficients among the variousDNS. The white square indicates the diffusion constant that would correspond to a Gaussian random walk byforcing a normal distribution fit on the measured PLT velocities. The inset presents the mean excess kurtosis ofPLT velocities coming from the DNS, as a function of L . The “error” bars correspond to the min/max values of thekurtosis. The diverging moments clearly indicate fat-tailed distributions with infinite/undefined moments (powerlaw behaviour with exponent α ≤ α averaged over all the fitted and accepted powerlaws with exponent ≤
2. The min/max indicated in table 1 denote the minimum and maximum values of theexponent from the fitting of the various distributions (Pareto, half-Cauchy, L´evy, Dagum, Singh–Maddala, Log-logistic, and inverted Weibull) for each PLT type.It is interesting to observe that the shape and the electrophoretic properties of platelets are reflected on the exponentof the power law. The lower the power law exponent, the higher the mobility through “extreme” tail events. Theexponents observed for AP-rep are consistently smaller. In the Supplementary Material, we provide a table thatsummarises the majority of the performed experiments with key quantities per DNS for the completeness of thestudy, e.g. one can find per DNS the power law exponents (considering all fitted distributions with exponent α ≤ NAP AP - no rep AP - rep
Mean Power Law α (min/max) 1.56 (1.35/1.75) 1.48 (1.23/1.75) 1.40 (1.12/1.65)Mean Deposited PLTs (power laws) 366 321 416Max Deposited PLTs (power laws) 1759 1512 2270Mean Deposited PLTs Normal Dist. 246 174 287Mean Deposited PLTs Exponential Dist. 250 175 290In a second step, we used the fitted distributions as generative mechanisms for simulating random walks and the9ransport of PLTs in the impact-R device. Table 1 presents the amount of deposited PLTs (the ones that cross thebottom wall) after running the stochastic simulations for 20 s of physical time, with L = 820 µm . Chopard et al.[4] report that 3125 AP have disappeared from the bulk during these 20 s , with an initial concentration of 4808 APper micro-litre.Given that the phenomenon is symmetric along the y -axis (cross-checked with the DNS), and due to the developedcell free layer close to the walls (trapping the crossing PLTs), we expect at least 1500 platelets out of 4808 (1500 ∼ /
2) to cross the bottom wall of the system.Table 1 shows that this expected number of deposited platelets is compatible with the proposed fat-tailed velocitydistribution, without the need to invoke special drift terms or a rheological potential. The table also presents theamount of deposited PLTs when forcing a Gaussian or exponential fit on the velocity distributions. Clearly, thedeposition values are way too small to describe the deposition rates observed in the in-vitro experiments. Therefore,the use of standard diffusive models heavily underestimates the deposition rate, compared to the L´evy flights thatproduce a 2 to 10 times higher amount of deposited PLTs.As shown in Fig. 4, it is also interesting to note that as opposed to the MSD which diverges with the system size,the exponent α keeps a consistent value when varying L from 50 µm to 500 µm (scale-invariance of power laws, i.e.the phenomena are expected to occur without a characteristic size or scale). P o w er L a w α System Size along the Wall-bounded Direction [µm]
AP - rep AP - no rep NAP
Figure 4: Mean power law exponents per cellular blood flow simulation. The power laws are independent of thesystem size (no diverging behaviour as with the moments of the velocity distribution). Any variation is due to thedifferent parameters per simulation (shapes, sizes, repulsive forces). The inset shows a few simulation snapshotsfrom systems of size 50 & 100 µm along the wall-bounded direction y . Combining the power law fitting on the tails with the normality tests that fail in all the experiments, the divergingmoments of the velocity distribution, the very high kurtosis, and the good agreement with the impact-R depositiondata, we give convergent evidences that PLTs do not follow a Gaussian random walk, but rather L´evy Flights. Thestandard diffusion equation does not apply to describe PLTs transport. Fractional differential equations might beneeded to account for such an anomalous diffusion processes at the macroscopic scale.As mentioned by Kumar and Graham [13], no clear and systematic mechanistic explanation was yet proposed for thesegregation and margination phenomena of platelets. In particular, no simplified mathematical description (such10s a set of transport equations or a simple stochastic process model) has emerged that captures the phenomena. Inthis study, we prove that PLT velocities, more specifically the tail of their distribution, can be described by powerlaws ( P ( v ) ∼ v − α − ) with exponent α ≤
2. We found no evidence of normally distributed PLT velocities, and thusPLTs cannot possibly exhibit standard diffusion, which is the norm when describing PLT transport.The new stochastic process model that we introduce does not need additional terms to describe margination, as isthe case in other studies [12, 5, 10]. A striking observation is that while our results are compatible with the stan-dard models (diffusion extracted from MSD/variance/second central moment), a further investigation (statisticalanalysis) reveals the anomalous behaviour of PLT transport. This can be explained from the fact that the morefat-tailed a distribution, the more statistical information lies in the tail, and the moments (on which standard dif-fusive models are built upon) become uninformative and unreliable. The majority of the numerical research in thisfield is limited to case studies that deal with too few blood cells in the computational domain, and this is due to thehigh computational cost of the simulations. Our highly scalable numerical framework allows us to investigate caseswith dimensions approaching the ones of the in-vitro experiment. We would like to emphasise that the ability tostudy sizes of clinical relevance (at least resolving the direction of interest), allowed us to observe the idiosyncraticbehaviour of PLT transportation, and to capture the anomalous characteristics that manifest in setups of largersizes.In addition to proposing a disruptive view of PLT transport physics in blood, our results have a concrete impacton the design of new and efficient platelet function tests. Those can be a vital part for the detection of car-dio/cerebrovascular diseases in clinical practice. Nowadays, there exist numerous platelet function tests for thediagnosis of disorders or the monitoring of anti-platelet therapies, but with limited prognostic capacity in clinicalpractice, due to contradicting results. As several studies show [29, 30, 31], there is a problem interpreting theresults and mapping them to patient risk and disease. We strongly believe that the emerging stochastic modelsfrom the present analysis will offer a paradigm shift for developing the next generation tests, as the understandingof PLT transport is inextricably associated with the success of platelet function testing. A first step in this directionwas proposed by Dutta et al. [32], a study in which clinically important properties, such as platelet adhesion andaggregation rates, can only be correctly inferred provided that PLT transport can be correctly described.
Source Code
The source code for the statistical analysis is available in GitHub: https://github.com/kotsaloscv/PLTs-FatTails . Thetools for performing the cellular blood flow simulations, as presented in Kotsalos et al. [19, 20], are available as part ofPalabos open-source library (npFEM specialised module): https://palabos.unige.ch . Acknowledgements
We acknowledge support from the Swiss National Supercomputing Centre (CSCS, Piz-Daint supercomputer), the NationalSupercomputing Centre in the Netherlands (Surfsara, Cartesius supercomputer), and the HPC Facilities of the University ofGeneva (Baobab cluster).
Funding
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grantagreement No 823712 (CompBioMed2 project). upplementary Material: Statistical Analysis (Graphical Explanation) For every cellular blood flow simulation, we produce from the sampled output the histogram of platelet velocities(velocities along the wall-bounded direction for channel flow, and radial velocities for tubular flow). Our goal isto fit a model on this histogram and use it for simulating random walks. Below some lower bound ( x min ), theoriginal data can be sufficiently described by the empirical distribution function (histogram), but above x min , dueto the fractionate information that we get from the tail, we need to apply more sophisticated statistical models.Since we are particularly interested on investigating the power law hypothesis, we fit a family of distributions thatasymptotically exhibit power law behaviour [18], e.g. Pareto, half-Cauchy, L´evy, Dagum, Singh–Maddala, Log-logistic, and inverted Weibull. Every distribution is fitted separately on the data (just at the tail), and thereforethe lower bound and power law exponent α vary per distribution (see figure below). The distributions with α ≤ x min : optimal & relaxed (sub-optimal but larger sample size) Empirical Distribution Function (EDF) - histogram Tail distributions (power laws): ➢ Pareto ➢ Half-Cauchy ➢ Lévy ➢ Burr ➢ Fisk ➢ Fréchet
Different pair (x min , α) per distribution
Try all x and choose the one that gives the best fit at the tail (KS-test): per distribution x: PLT velocity f ̴ x -α-1 , x → ∞ Asymptotic behavior
Statistical analysis per direct numerical simulation. Per investigated distribution, we define x min with the optimisation process as described in Methods, and α through the MLE method. To generate data for the random walks, the DNS sampled output (per simulation) is split into the body and thetail (see figure above). Regarding the body, we use the empirical distribution function (ECDF) for generating newvelocities. For the generated velocities that are above the x min (varies per fitted distribution), we use the fitteddistributions to re-generate velocities belonging to the tail. For the new data to be as close as possible to theoriginal (from DNS), the fitted distributions should pass the goodness-of-fit and LLR tests. Assuming a plausibledistribution for the tail, it may be the case that the tail sample size is small ( n ≈ x min threshold, i.e. reduce the lower bound and thus increase the sample size on which to fit. Forthis relaxed version, we need to perform an additional goodness-of-fit test [22]. For this test, we use as metric theKuiper’s statistic, which is a variation of the Kolmogorov-Smirnov, but more sensitive on capturing differences atthe tails. In short, we generate data from the ECDF, then synthetic data with the optimal x min , and we comparethe two data sets using the Kuiper’s test, from which we get the statistic. Following, we generate another setof synthetic data with the relaxed x min , and we compare it with the data from the ECDF, extracting anotherstatistic. If the last statistic is smaller than the one of the first comparison, then the relaxed model is a validalternative. By repeating this process more than a thousand times [22], we get a p-value coming from the fractionof valid/accepted relaxed models to the overall repetitions. Keep in mind, that the relaxed tail gives an additionalfitting for consideration, thus for each distribution we have the optimal and relaxed fittings . Therefore, per DNSand per fitted distribution, we construct a generative mechanism that “feeds” the random walk models.12 upplementary Material: Power Law Distributions The probability densities presented here are in the “standardised” form, i.e. the location parameter is zero and thescale parameter is one. However, the analysis is not affected by this “standardised” form, since these parametersare only for shifting and scaling, respectively, the investigated distributions. The power law exponent α is givenby the asymptotic analysis of the distributions (probability density function-pdf), which decay like f ∼ x − α − .The presentation of the pdf with − F = P r ( X > x ),complementary cumulative distribution function) decays like ¯ F ∼ x − α . Pareto
The probability density function for Pareto is: f ( x, b ) = bx b +1 for x ≥ b >
0. For x → ∞ , f ( x, b ) ∼ x − b − . Thevariance is ∞ for b ≤ ∞ for b ≤
1. The power law exponent is b . Half-Cauchy
The probability density function for half-Cauchy is: f ( x ) = π (1+ x ) for x ≥
0. For x → ∞ , f ( x ) ∼ x − . The meanand variance are undefined. The power law exponent is fixed to 1. L´evy
The probability density function for L´evy is: f ( x ) = √ πx exp (cid:0) − x (cid:1) for x ≥
0. For x → ∞ , f ( x ) ∼ x − . . Boththe mean and the variance are ∞ . The power law exponent is fixed to 0 . Dagum (or Burr Type III)
The probability density function for Dagum is: f ( x, c, d ) = cdx − c − / (1 + x − c ) d +1 for x ≥ c, d >
0. For x → ∞ , f ( x, c, d ) ∼ x − c − . The mean is undefined for c ≤
1, and the variance is undefined for c ≤
2. The powerlaw exponent is c . Singh–Maddala (or Burr Type XII)
The probability density function for Singh–Maddala is: f ( x, c, d ) = cdx c − / (1 + x c ) d +1 for x ≥ c, d >
0. For x → ∞ , f ( x, c, d ) ∼ x c − /x cd + c ∼ x − cd − . The mean is undefined for cd ≤
1, and the variance is undefined for cd ≤
2. The power law exponent is cd . Log-logistic (or Fisk)
The probability density function for Log-logistic is: f ( x, c ) = cx − c − (1 + x − c ) − for x ≥ c >
0. For x → ∞ , f ( x, c ) ∼ x − c − . The mean is undefined for c ≤
1, and the variance is undefined for c ≤
2. The power law exponentis c . Inverted Weibull (or Fr´echet)
The probability density function for inverted Weibull is: f ( x, c ) = cx − c − exp ( − x − c ) for x > c >
0. For x → ∞ , f ( x, c ) ∼ x − c − . The mean is undefined for c ≤
1, and the variance is undefined for c ≤
2. The power law exponentis c . 13 upplementary Material: Summary of DNS and Statistical Analysis PLT : Activated/Non-Activated platelets, rep : no repulsive forces at 20 (weight of collision potential - DNS related), exp : ei-ther scaled impact-R (box) or tube, dir : system size along the wall-bounded direction/diameter [ µm ], Ht : hematocrit [%], SR : ShearRate [ s − ], num RBCs/PLTs : number of RBCs/PLTs in DNS, D (DNS) : Diffusion Coefficient - MSD of DNS [ m /s ], D (Gaussian) :Diffusion Coefficient - Random Walk (Gaussian fitted on DNS) [ m /s ], fat-tails : number of accepted power laws (considering only theones with α ≤ ) for both optimal and relaxed tail fitting, min sample size : minimum sample size of one of the accepted power laws, avg sample size : average sample size considering all accepted power laws, avg exp : average power law exponent considering all acceptedpower laws, max depo : deposited PLTs (after 20s physical time [4]) using the random walks (Impact-R at full scale) corresponding tothe min exp . Simulations that seem to repeat are with/without the Particle In Kernel (PIK) technique [20] (validating its consistency). PLT rep exp dir Ht SR numRBCs numPLTs D(DNS) D(Gaussian) excesskurtosis fat-tails minsamplesize avgsamplesize avgexp minexp maxexp maxdepo
AP 200 box 500 35 100 4765 953 3.00E-10 1.12E-10 23.3 2 640 822 1.522 1.35 1.70 367AP 200 box 500 35 100 4765 953 3.11E-10 1.01E-10 31.4 3 362 619 1.302 1.00 1.63 384AP 200 box 50 35 100 476 95 1.40E-10 1.27E-10 17.5 4 181 772 1.319 1.00 1.59 640AP 200 box 50 20 100 272 54 1.54E-10 4.87E-11 5.4 3 104 540 1.314 1.00 1.49 413AP 200 box 100 35 100 953 190 2.87E-10 1.77E-10 26.8 2 508 3081 1.325 1.30 1.35 879AP 200 box 50 35 100 476 95 1.71E-10 1.28E-10 39.2 8 102 1194 1.359 1.00 1.62 2270AP 200 box 50 20 100 272 54 1.11E-10 5.65E-11 4.6 3 101 880 1.735 1.50 1.86 375AP 200 tube 52 35 200 374 74 1.06E-10 1.07E-10 20.0 6 136 402 1.418 1.00 1.96 4659AP 200 tube 52 35 200 374 74 7.97E-11 1.67E-10 21.8 6 144 530 1.345 1.00 1.71 2509AP 100 box 100 35 100 953 190 1.27E-10 5.33E-11 7.3 3 271 464 1.338 1.00 1.76 325AP 100 box 500 35 100 4765 953 1.09E-10 3.92E-11 41.0 6 178 513 1.490 1.00 1.97 743AP 100 box 50 35 100 476 95 7.32E-11 4.87E-11 5.1 3 108 919 1.368 1.00 1.61 270AP 100 box 50 20 100 272 54 8.94E-11 2.96E-11 2.6 6 102 505 1.489 1.00 1.91 1315AP 100 box 100 35 100 953 190 1.89E-10 6.35E-11 46.0 5 301 1633 1.282 1.00 1.59 565AP 100 box 100 35 100 953 190 1.13E-10 5.77E-11 14.2 5 283 740 1.546 1.00 1.95 731AP 100 box 250 35 100 2382 476 4.12E-10 1.60E-10 84.5 1 184 184 1.454 1.45 1.45 398AP 100 box 50 35 100 476 95 1.08E-10 5.56E-11 9.0 2 108 159 1.473 1.39 1.56 237AP 20 box 50 35 100 476 47 6.06E-11 3.01E-10 125.8 4 2155 2580 0.765 0.72 0.86 4782AP 20 box 50 35 100 476 95 5.81E-11 1.86E-11 1.6 2 141 1155 1.615 1.59 1.64 210AP 20 box 50 35 100 476 95 5.81E-11 2.57E-11 3.8 5 130 838 1.338 1.00 1.59 215AP 20 box 50 35 100 476 95 6.17E-11 2.61E-11 1.2 4 120 1176 1.454 1.00 1.64 204AP 20 box 50 35 100 476 95 7.20E-11 5.14E-11 44.7 7 457 1860 1.336 1.00 1.91 1512AP 20 box 50 35 100 476 95 7.91E-11 4.00E-11 15.2 4 107 4085 1.646 1.40 1.92 365AP 20 box 50 20 100 272 54 8.78E-11 1.35E-11 1.8 0AP 20 box 100 35 100 953 190 8.99E-11 6.47E-11 155.9 5 109 4567 1.588 1.30 1.94 573AP 20 box 250 35 100 2382 476 1.29E-10 6.15E-11 111.5 3 166 8197 1.580 1.51 1.70 345AP 20 box 500 35 100 4765 953 5.92E-11 1.91E-11 107.3 3 229 10127 1.527 1.28 1.75 223AP 20 box 50 35 100 476 95 8.17E-11 3.34E-11 39.1 5 387 1197 1.254 1.00 1.62 1312AP 20 box 50 20 100 272 54 4.11E-11 1.69E-11 1.5 0AP 20 tube 52 35 200 374 74 6.11E-11 2.74E-11 4.5 4 102 544 1.358 1.00 1.73 1093NAP 20 box 50 35 100 476 95 1.18E-10 8.74E-11 42.3 2 7090 7933 1.736 1.72 1.75 363NAP 20 box 100 35 100 953 190 9.82E-11 6.22E-11 70.8 4 756 4044 1.472 1.23 1.77 558NAP 20 box 250 35 100 2382 476 1.37E-10 6.36E-11 295.1 3 1731 8573 1.568 1.29 1.84 375NAP 20 box 500 35 100 4765 953 1.09E-10 5.05E-11 862.8 1 3867 3867 1.177 1.18 1.18 464NAP 20 box 50 35 100 476 95 5.70E-11 2.72E-11 4.7 4 110 1216 1.480 1.00 1.76 279NAP 20 box 50 35 100 476 95 1.12E-10 4.67E-11 16.1 2 1030 2223 1.708 1.59 1.83 256NAP 20 box 50 35 100 370 74 3.06E-10 1.46E-10 21.0 0NAP 20 box 50 35 400 476 95 1.72E-10 1.50E-10 28.8 8 454 1118 1.490 1.00 1.98 2180NAP 20 box 50 35 400 476 95 1.61E-10 2.07E-10 15.2 3 115 2879 1.552 1.30 1.85 1342NAP 20 box 50 35 400 370 74 2.73E-10 2.89E-10 10.5 2 1180 1567 1.753 1.66 1.85 839NAP 20 box 50 35 100 476 95 7.93E-11 7.36E-11 28.3 8 102 3119 1.467 1.00 1.94 1759NAP 20 box 50 35 100 476 95 1.02E-10 4.89E-11 21.9 4 498 2012 1.562 1.28 1.88 432NAP 20 box 50 35 100 370 74 1.46E-10 7.13E-11 14.0 2 869 1901 1.719 1.68 1.76 393NAP 20 box 50 35 400 476 95 1.62E-10 1.73E-10 9.7 0NAP 20 box 50 35 400 476 95 1.44E-10 1.09E-10 2.7 3 189 1145 1.333 1.00 1.52 536NAP 20 box 50 35 400 370 74 2.06E-10 2.47E-10 10.3 2 776 1308 1.757 1.69 1.83 619NAP 20 box 100 35 100 953 190 1.91E-10 9.95E-11 36.1 2 9172 9529 1.516 1.50 1.53 471NAP 20 box 500 35 100 4765 953 1.44E-10 5.95E-11 65.6 2 110 9762 1.820 1.72 1.92 301NAP 20 box 50 35 100 476 95 1.25E-10 4.84E-11 14.7 6 105 1707 1.597 1.37 1.88 513NAP 20 tube 52 35 200 374 74 5.31E-11 3.09E-11 5.3 3 108 474 1.770 1.55 1.91 268NAP 20 tube 52 35 200 374 74 4.39E-11 4.48E-11 7.1 7 123 298 1.423 1.00 1.91 4797NAP 20 box 50 35 100 476 95 1.36E-10 1.11E-10 34.1 5 1044 2692 1.546 1.22 1.88 718NAP 20 box 50 35 100 476 95 1.06E-10 5.32E-11 3.1 4 113 1053 1.438 1.00 1.73 296NAP 20 box 50 35 100 476 95 6.32E-11 4.07E-11 12.9 4 147 738 1.313 1.00 1.47 360NAP 20 box 50 35 100 476 95 7.55E-11 3.65E-11 20.0 5 107 1049 1.325 1.00 1.73 1386NAP 20 box 50 35 100 476 95 1.26E-10 1.17E-10 38.2 2 11308 12913 1.879 1.84 1.92 421NAP 20 box 50 35 100 476 95 1.18E-10 9.06E-11 42.3 2 7090 7933 1.736 1.72 1.75 378 eferences [1] P. Harrison, Platelet function analysis, Blood Reviews 19 (2) (2005) 111–123. doi:10.1016/j.blre.2004.05.002 .[2] A. L. Zydney, C. K. Colton, Augmented Solute Transport in the Shear Flow of a Concentrated Suspension.,PCH. Physicochemical hydrodynamics 10 (1) (1988) 77–96.[3] K. Affeld, L. Goubergrits, N. Watanabe, U. Kertzscher, Numerical and experimental evaluation of plateletdeposition to collagen coated surface at low shear rates, Journal of Biomechanics 46 (2) (2013) 430–436. doi:10.1016/j.jbiomech.2012.10.030 .[4] B. Chopard, D. R. de Sousa, J. L¨att, L. Mountrakis, F. Dubois, C. Yourassowsky, P. Van Antwerpen, O. Eker,L. Vanhamme, D. Perez-Morga, G. Courbebaisse, E. Lorenz, A. G. Hoekstra, K. Z. Boudjeltia, A physicaldescription of the adhesion and aggregation of platelets, Royal Society Open Science 4 (4) (2017) 170219.[5] L. Crowl, A. L. Fogelson, Analysis of mechanisms for platelet near-wall excess under arterial blood flow con-ditions, Journal of Fluid Mechanics 676 (2011) 348–375. doi:10.1017/jfm.2011.54 .[6] H. Zhao, E. S. Shaqfeh, Shear-induced platelet margination in a microchannel, Physical Review E - Statistical,Nonlinear, and Soft Matter Physics 83 (6) (2011) 061924. doi:10.1103/PhysRevE.83.061924 .[7] H. Zhao, E. S. Shaqfeh, V. Narsimhan, Shear-induced particle migration and margination in a cellular suspen-sion, Physics of Fluids 24 (1) (2012) 011902.[8] D. A. Reasor, M. Mehrabadi, D. N. Ku, C. K. Aidun, Determination of critical parameters in platelet margina-tion, Annals of Biomedical Engineering 41 (2) (2013) 238–249. doi:10.1007/s10439-012-0648-7 .[9] K. Vahidkhah, S. L. Diamond, P. Bagchi, Platelet dynamics in three-dimensional simulation of whole blood,Biophysical Journal 106 (11) (2014) 2529–2540. doi:10.1016/j.bpj.2014.04.028 .[10] M. Mehrabadi, D. N. Ku, C. K. Aidun, A Continuum Model for Platelet Transport in Flowing Blood Basedon Direct Numerical Simulations of Cellular Blood Flow, Annals of Biomedical Engineering 43 (6) (2015)1410–1421. doi:10.1007/s10439-014-1168-4 .[11] M. Mehrabadi, D. N. Ku, C. K. Aidun, Effects of shear rate, confinement, and particle parameters on margina-tion in blood flow, Physical Review E 93 (2) (2016) 023109. doi:10.1103/PhysRevE.93.023109 .[12] E. C. Eckstein, F. Belgacem, Model of platelet transport in flowing blood with drift and diffusion terms,Biophysical Journal 60 (1) (1991) 53–69. doi:10.1016/S0006-3495(91)82030-6 .[13] A. Kumar, M. D. Graham, Margination and segregation in confined flows of blood and other multicomponentsuspensions, Soft Matter 8 (41) (2012) 10536–10548. doi:10.1039/c2sm25943e .[14] K. Vahidkhah, P. Bagchi, Microparticle shape effects on margination, near-wall dynamics and adhesion ina three-dimensional simulation of red blood cell suspension, Soft Matter 11 (11) (2015) 2097–2109. doi:10.1039/c4sm02686a .[15] T. R. Lee, M. Choi, A. M. Kopacz, S. H. Yun, W. K. Liu, P. Decuzzi, On the near-wall accumulation ofinjectable particles in the microcirculation: Smaller is not better, Scientific Reports 3 (1) (2013) 1–8. doi:10.1038/srep02079 .[16] K. Yeo, M. R. Maxey, Anomalous diffusion of wall-bounded non-colloidal suspensions in a steady shear flow,EPL (Europhysics Letters) 92 (2) (2010) 24008. doi:10.1209/0295-5075/92/24008 .[17] M. Gross, T. Kr¨uger, F. Varnik, Fluctuations and diffusion in sheared athermal suspensions of deformableparticles, EPL (Europhysics Letters) 108 (6) (2014) 68006. doi:10.1209/0295-5075/108/68006 .[18] C. Kleiber, S. Kotz, Statistical Size Distributions in Economics and Actuarial Sciences, Wiley Series in Prob-ability and Statistics, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2003.[19] C. Kotsalos, J. Latt, B. Chopard, Bridging the computational gap between mesoscopic and continuum modelingof red blood cells for fully resolved blood flow, Journal of Computational Physics 398 (2019) 108905. doi:10.1016/j.jcp.2019.108905 . 1520] C. Kotsalos, J. Latt, J. Beny, B. Chopard, Digital Blood in Massively Parallel CPU/GPU Systems for theStudy of Platelet Transport, Interface Focus (accepted for publication) (Computational Biomedicine) (2020).[21] J. Latt, O. Malaspinas, D. Kontaxakis, A. Parmigiani, D. Lagrava, F. Brogi, M. B. Belgacem, Y. Thorimbert,S. Leclaire, S. Li, F. Marson, J. Lemus, C. Kotsalos, R. Conradin, C. Coreixas, R. Petkantchin, F. Ray-naud, J. Beny, B. Chopard, Palabos: Parallel Lattice Boltzmann Solver, Computers and Mathematics withApplications (apr 2020). doi:10.1016/j.camwa.2020.03.022 .[22] A. Clauset, C. R. Shalizi, M. E. Newman, Power-law distributions in empirical data (jun 2009). doi:10.1137/070710111 .[23] N. N. Taleb, Statistical Consequences of Fat Tails: Real World Preasymptotics, Epistemology, and Applications(Technical Incerto), 1st Edition, Vol. 1, STEM Academic Press, 2020.[24] J. Alstott, E. Bullmore, D. Plenz, Powerlaw: A python package for analysis of heavy-tailed distributions, PLoSONE 9 (1) (2014) e85777. doi:10.1371/journal.pone.0085777 .[25] W. Jy, L. L. Horstman, D. Homolak, Y. S. Ahn, Original article: Electrophoretic properties of platelets fromnormal, thrombotic and ITP patients by doppler electrophoretic light scattering analysis, Platelets 6 (6) (1995)354–358. doi:10.3109/09537109509078471 .[26] J. J. Betts, J. P. Betts, J. T. Nicholson, Significance of ADP, plasma and platelet concentration in plateletelectrophoretic studies (1968). doi:10.1038/2191280b0 .[27] J. R. Hampton, J. R. Mitchell, Modification of the electrokinetic response of blood platelets to aggregatingagents, Nature 210 (5040) (1966) 1000–1002. doi:10.1038/2101000a0 .[28] K. Z. Boudjeltia, C. Kostalos, D. Ribeiro, A. Rousseau, C. Lelubre, O. Sartenaer, M. Piagnerelli, J. Dohet-Eraly, F. Dubois, N. Tasiaux, B. Chopard, A. V. Meerhaeghe, Spherization of red blood cells and plateletsmargination in COPD patients., medRxiv (2020) doi:10.1101/2020.04.03.20051748 .[29] N. J. Breet, J. W. Van Werkum, H. J. Bouman, J. C. Kelder, H. J. Ruven, E. T. Bal, V. H. Deneer, A. M.Harmsze, J. A. Van Der Heyden, B. J. Rensing, M. J. Suttorp, C. M. Hackeng, J. M. Ten Berg, Comparisonof platelet function tests in predicting clinical outcome in patients undergoing coronary stent implantation,JAMA - Journal of the American Medical Association 303 (8) (2010) 754–762. doi:10.1001/jama.2010.181 .[30] S. M. Picker, In-vitro assessment of platelet function, Transfusion and Apheresis Science 44 (3) (2011) 305–319. doi:10.1016/j.transci.2011.03.006doi:10.1016/j.transci.2011.03.006