Parameter estimation of platelets deposition: Approximate Bayesian computation with high performance computing

by   Ritabrata Dutta, et al.
University of Geneva

A numerical model that quantitatively describes how platelets in a shear flow adhere and aggregate on a deposition surface has been recently developed in Chopard et al. (2015); Chopard et al. (2017). Five parameters specify the deposition process and are relevant for a biomedical understanding of the phenomena. Experiments give observations, at five time intervals, on the average size of the aggregation clusters, their number per mm^2, the number of platelets and the ones activated per μℓ still in suspension. By comparing in-vitro experiments with simulations, the model parameters can be manually tuned (Chopard et al. (2015); Chopard et al. (2017)). Here, we use instead approximate Bayesian computation (ABC) to calibrate the parameters in a data-driven automatic manner. ABC requires a prior distribution for the parameters, which we take to be uniform over a known range of plausible parameter values. ABC requires the generation of many pseudo-data by expensive simulation runs, we have thus developed an high performance computing (HPC) ABC framework, taking into account accuracy and scalability. The present approach can be used to build a new generation of platelets functionality tests for patients, by combining in-vitro observation, mathematical modeling, Bayesian inference and high performance computing.



There are no comments yet.


page 6


An ABC interpretation of the multiple auxiliary variable method

We show that the auxiliary variable method (Møller et al., 2006; Murray ...

Computer code validation via mixture model estimation

When computer codes are used for modeling complex physical systems, thei...

Bahadur efficiencies of the Epps–Pulley test for normality

The test for normality suggested by Epps and Pulley (1983) is a serious ...

Perceptual similarity between piano notes: Simulations with a template-based perception model

In this paper the auditory model developed by Dau et al. [J. Acoust. Soc...

Approximate Hubel-Wiesel Modules and the Data Structures of Neural Computation

This paper describes a framework for modeling the interface between perc...

Corrected Kriging update formulae for batch-sequential data assimilation

Recently, a lot of effort has been paid to the efficient computation of ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Blood platelets play a major role in the complex process of blood coagulation, involving adhesion, aggregation and spreading on the vascular wall to stop a hemorrhage while avoiding the vessel occlusion. Platelets also play a key role in the occurrence of cardio/cerebro-vascular accidents that constitute a major health issue in our societies. In 2015, Cardiovascular diseases (CVD), including disorders of the heart and blood vessels, were the first cause of mortality worldwide, causing 31% of deaths (Organization, 2015). Antiplatelet therapy generally reduces complications in patients undergoing arterial intervention (Steinhubl et al., 2002; Mehta et al., 2001). However, the individual response to dual antiplatelet therapy is not uniform and consistent studies reported that even under platelets therapy there were recurrences of atherothrombotic events (Matetzky et al., 2004; Gurbel et al., 2005; Geisler et al., 2006; Hochholzer et al., 2006; Marcucci et al., 2009; Price et al., 2008; Sibbing et al., 2009). In most cases, a standard posology is prescribed to patients, which does not take into account the inter-individual variability linked to the absorption or the effectiveness of these molecules. This was supported by a recent study (Koltai et al., 2017), reporting the high patient-dependency of the response of the antithrombotic drugs. We should also note that the evaluation of the response to a treatment by the existing tests is test-dependent.

Nowadays, platelet function testing is performed either as an attempt to monitor the efficacy of anti-platelet drugs or to determine the cause of abnormal bleeding or pro-thrombotic status. The most common method consists of using an optical aggregometer that measures the transmittance of light passing through plasma rich in platelets (PRP) or whole blood (Harrison, 2009; Born and Cross, 1963), to evaluate how platelets tend to aggregate. Other aggregometers determine the amount of aggregated platelets by electric impedance (Velik-Salchner et al., 2008) or luminescence. In specific contexts, flow cytometry (Michelson et al., 2002) is also used to assess platelet reactivity (VASP test (Bonello et al., 2009)). Determination of platelet functions using these different existing techniques in patients undergoing coronary stent implantation have been evaluated in Breet et al. (2010), which shows the correlation between the clinical biological measures and the occurrence of a cardiovascular event was null for half of the techniques and rather modest for others. This may be due to the fact that no current test allows the analysis of the different stages of platelet activation or the prediction of the in vivo behavior of those platelets (Picker, 2011; Koltai et al., 2017). It is well known that the phenomenon of platelet margination (the process of bringing platelets to the vascular wall) is dependent on the number and shape of red blood cells and their flow (Piagnerelli et al., 2007), creating different pathologies for different diseases (e.g., diabetes, End Renal Kidney Disease, hypertension, sepsis). Further, platelet margination is also known to be influenced by the aspect ratio of surrogate platelet particles (Reasor et al., 2013). Although there is a lot of data reported by recent research works (Maxwell et al., 2007) on the molecules involved in platelet interactions, these studies indicate that there is a lack of knowledge on some fundamental mechanisms that should be revealed by new experiments.

Hence, the challenge is to find parameters connecting the dynamic processes of adhesion and aggregation of platelets to the data collected from the individual patients. Recently, by combining digital holography microscopy (DHM) and mathematical modeling, Chopard et al. (2015, 2017) provided a physical description of the adhesion and aggregation of platelets in the Impact-R device. A numerical model is developed that quantitatively describes how platelets in a shear flow adhere and aggregate on a deposition surface. This is the first innovation in understanding the molecular dynamics involved in platelet interactions. Five parameters specify the deposition process and are relevant for a biomedical understanding of the phenomena. One of the main intuition is that the values of these parameters (e.g., adhesion and aggregation rates) are precisely the information needed to assess various possible pathological situations and quantify their severity regarding CVD. Further, it was shown in Chopard et al. (2017) that, by hand-tuning the parameters of the mathematical model, the deposition patterns observed for a set of healthy volunteers in the Impact-R can be reproduced.

Assuming that these parameters can determine the severity of CVD, how do we estimate the adhesion and aggregation rates of given patients by a clinical test? The determination of these adhesion and aggregation rates by hand-tuning is clearly not a solution as we need to search the high-dimensional parameter space of the mathematical model, which becomes extremely expensive and time consuming. We further notice, this has to be repeated for each patient and thus requires a powerful numerical approach. In this work, we resolve the question of estimating the parameters using Bayesian uncertainty quantification. Due to a complex stochastic nature, the numerical model for platelet deposition does not have a tractable likelihood function. We use Approximate Bayesian Computation (ABC), a likelihood-free inference scheme, with an optimal application of HPC (Dutta et al., 2017a) to provide a Bayesian way to estimate adhesion and aggregation rates given the deposition patterns observed in the Impact-R of platelets collected from a patient. Obviously, the clinical applicability of the proposed technique to provide a new platelet function test remains to be explored, but the numerical model (Chopard et al., 2017) and the proposed inference scheme here, bring the technical elements together to build a new class of medical tests.

In section 2 we introduce the necessary background knowledge about the platelet deposition model, whereas Section 3 recalls the concept of Bayesian inference and introduces the HPC framework of ABC used in this study. Then we illustrate the results of the parameter determination for platelet deposition model using ABC methodology, collectively for 7 patients in Section 4. Clearly, the same methodology can be used to determine the parameter values for each individual patients in a similar manner for a CVD clinical test. Finally, in section 5 we conclude the paper and discuss its impact from a biomedical perspective.

2 Background and Scientific Relevance

The Impact-R (Shenkman et al., 2008) is a well-known platelet function analyzer. It is a cylindrical device filled in with whole blood from a donor. Its lower end is a fixed disk, serving as a deposition surface, on which platelets adhere and aggregate. The upper end of the Impact-R cylinder is a rotating cone, creating an adjustable shear rate in the blood. Due to this shear rate, platelets move towards the deposition surface, where they adhere or aggregate. Platelets aggregate next to already deposited platelets, or on top of them, thus forming clusters whose size increase with time. This deposition process has been successfully described with a mathematical model in Chopard et al. (2015, 2017).

The numerical model (coined in what follows) requires five parameters that specify the deposition process and are relevant for a bio-medical understanding of the phenomena. In short, the blood sample in the Impact-R device contains an initial number of non-activated platelets per and a number of pre-activated platelets per

. Initially both type of platelets are supposed to be uniformly distributed within the blood. Due to the process known as shear-induced diffusion, platelets hit the deposition surface. Upon such an event, an activated platelets will adhere with a probability that depends on its adhesion rate,

, that we would like to determine. Platelets that have adhered on the surface are the seed of a cluster that can grow due to the aggregation of the other platelets reaching the deposition surface. We denote with the rate at which new platelets will deposit next to an existing cluster. We also introduce the rate at which platelets deposit on top of an existing cluster. An important observation made in Chopard et al. (2015, 2017) is that albumin, which is abundant in blood, compete with platelet for deposition. This observation is compatible with results reported in different experimental settings (Remuzzi and Boccardo, 1993; Fontaine et al., 2009; Sharma et al., 1981). As a consequence, the number of aggregation clusters and their size tends to saturate as time goes on, even though there are still a large number of platelets in suspension in the blood.

To describe this process in the model, two extra parameters, , the deposition rate of albumin, and , a factor that accounts for the decrease of platelets adhesion and aggregation on locations where albumin has already deposited, were introduced. The numerical model is described in full detail in Chopard et al. (2015, 2017). Here we simply repeat the main elements. Due to the mixing in the horizontal direction, it was assumed that the activated platelets (AP), non-activated platelets (NAP) and albumin (Al) in the bulk can be described by a 1D diffusion equation along the vertical axis


where is the density of either AP, NAP or Al, and are correspondingly the flux of particles and the shear induced diffusion. Upon reaching a boundary layer above the deposition substrate, adhesion and aggregation will take place according to


where is the number of particles in the boundary layer, a surface element on the deposition surface, and is the deposition rate, which evolves during time and varies across the substrate, according to the deposition history. For the deposition process, particles are considered as discrete entities that can attach to any position of the grid representing the deposition surface, as sketched in Fig. 1. In this figure, the gray levels illustrate the density of albumin already deposited in each cell. The picture also illustrates the adhesion, aggregation, and vertical deposition along the -axis. On the left panel, activated platelets (gray side disks) deposit first. Then in the second panel, non-activated platelets (white side disks) aggregate next to an already formed cluster. Both pre-activated and non-activated platelets can deposit on top of an existing cluster.

The deposition rules are the following. An albumin that reaches the substrate at time deposits with a probability which depends on the local density of already deposited Al. We assume that is proportional to the remaining free space in the cell,


where is a parameter and is determined by the constraint that at most 100,000 albumin particles can fit in a deposition cell of area , corresponding to the size of a deposited platelet (obtained as the smallest variation of cluster area observed with the microscope).

An activated platelet that hits a platelet-free cell deposits with a probability , where decreases as the local concentration of albumin increases. We assumed that


where and are parameters. This expression can be justified by the fact that a platelet needs more free space than an albumin to attach to the substrate, due to their size difference. In other words, the probability of having enough space for a platelet, decreases roughly exponentially with the density of albumin in the substrate. This can be validated with a simple deposition model on a grid, where small and large objects compete for deposition.

Once an activated platelet has deposited, it is the seed of a new cluster that grows further due to the aggregation of further platelets. In our model, AP and NAP can deposit next to already deposited platelets. From the above discussion, the aggregation probability is assumed to be


with another parameter.

The above deposition probabilities can also be expressed as deposition rate over the given simulation time step (see Chopard et al. (2017) for details), hence giving a way to couple the diffusion equation (1) with the 2D discrete deposition process sketched in Fig. 1. Particles that did not deposit at time are re-injected in the bulk and contribute to boundary condition of eq. (1) at .

Figure 1: Sketch of the deposition substrate, discretized in cells of area equal to the surface of a platelet.
Figure 2: The deposition surface of the Impact-R device after 300 seconds (left) and the corresponding results of the deposition in the mathematical model (right). Black dots represent the deposited platelets that are grouped in clusters.

To the best of our knowledge, except for Chopard et al. (2015, 2017) there is no model in the literature that describes quantitatively the proposed in-vitro experiment. The closest approach is that of Affeld et al. (2013), but albumin is not included, and the role of pre-activated and non-activated platelets is not differentiated. Also, we are not aware of any other study than ours that reports both the amount of platelets in suspension as a function of time and those on the deposition surface.

The validity of the proposed numerical model has been explored in detail in Chopard et al. (2017). This validation is based on the fact that the model, using hand-tuned parameters can reproduce the time-dependent experimental observations very well. We refer the readers to Chopard et al. (2017) for a complete discussion. Here we briefly recall the main elements that demonstrate the excellent agreement of the model and the simulations. We reproduce Fig. 2 from Chopard et al. (2017)

, showing the visual similarity between the actual and simulated deposition pattern. In the validation study, the evolution of the number of clusters, their average size and the numbers of pre-activated and non-activated platelets still in suspension matched quantitatively with the experimental measurements at times 20, 60, 120 and 300 s. In addition, a very good agreement between the simulated deposition pattern and the experiment was also found by comparing the distributions of the areas and volumes of the aggregates.

To be noticed, the validation reported in Chopard et al. (2017) was done using manually estimated parameters. As the main goal of this research is to propose an inference scheme to learn the parameters in a data-driven manner, a validation for the model and the inference scheme is reported in Fig. 6 below, using the inferred posterior distribution which also includes a quantification of prediction error.

For the purpose of the present study, the model is parametrized in terms of the five quantities introduced above, namely the adhesion rate , the aggregation rates and , the deposition rate of albumin , and the attenuation factor . Some additional parameters of the model, specifically, the shear-induced diffusion coefficient and the thickness of the boundary layer (Chopard et al., 2017), are assumed here to be known. Collectively, we define

If the initial values for and , as well as the concentration of albumin are known from the experiment, we can forward simulate the deposition of platelets over time using model for the given values of these parameters :


where and are correspondingly average size of the aggregation clusters, their number per , the number of non-activated and pre-activated platelets per still in suspension at time .

The Impact-R experiments have been repeated with the whole blood obtained from seven donors and the observations were made at time, 0 sec., 20 sec., 60 sec., 120 sec. and 300 sec. At these five time points, are measured. Let us call the observed dataset collected through experiment as,

By comparing the number and size of the deposition aggregates obtained from the in-vitro experiments with the computational results obtained by forward simulation from the numerical model (see Fig. 2 for an illustration), the model parameters were manually calibrated by a trial and error procedure in Chopard et al. (2017). Due to the complex nature of the model and high-dimensional parameter space, this manual determination of the parameter values are subjective and time consuming.

However, if the parameters of the model could be learned more rigorously with an automated data-driven methodology, we could immensely improve the performance of these models and bring this scheme as a new clinical test for platelet functions. To this aim, here we propose to use ABC for Bayesian inference of the parameters. As a result of Bayesian inference to this context, not only we can automatically and efficiently estimate the model parameters, but we can also perform parameter uncertainty quantification in a statistically sound manner, and determine if the provided solution is unique.

3 Bayesian Inference

We can quantify the uncertainty of the unknown parameter by a posterior distribution given the observed dataset

. A posterior distribution is obtained, by Bayes’ Theorem as,


where , and are correspondingly the prior distribution on the parameter , the likelihood function, and the marginal likelihood. The prior distribution

ensures a way to leverage the learning of parameters with prior knowledge, which is commonly known due to the availability of medical knowledge regarding cardio-vascular diseases. If the likelihood function can be evaluated, at least up to a normalizing constant, then the posterior distribution can be approximated by drawing a sample of parameter values from the posterior distribution using (Markov chain) Monte Carlo sampling schemes

(Robert and Casella, 2005). For the simulator-based models considered in Section 2, the likelihood function is difficult to compute as it requires solving a very high dimensional integral. In next Subsection 3.1, we illustrate ABC to perform Bayesian Inference for models where the analytical form of the likelihood function is not available in closed form or not feasible to compute.

3.1 Approximate Bayesian computation

ABC allows us to draw samples from the approximate posterior distribution of parameters of the simulator-based models in absence of likelihood function, hence to perform approximate statistical inference (eg., point estimation, hypothesis testing, model selection etc.) in a data-driven manner. In a fundamental Rejection ABC scheme, we simulate from the model a synthetic dataset for a parameter value and measure the closeness between and using a pre-defined discrepancy function . Based on this discrepancy measure, ABC accepts the parameter value when is less than a pre-specified threshold value .

As the Rejection ABC scheme is computationally inefficient, to explore the parameter space in an efficient manner, there exists a large group of ABC algorithms (Marin et al., 2012). As pointed in Dutta et al. (2017a), these ABC algorithms, consist of four fundamental steps:

  • (Re-)sample a set of parameters either from the prior distribution or from an already existing set of parameter samples;

  • For each of the sample from the whole set or a subset, perturb it using the perturbation kernel, accept the perturbed sample based on a decision rule governed by a threshold or repeat the whole second step;

  • For each parameter sample calculate its weight;

  • Normalize the weights, calculate a co-variance matrix and adaptively re-compute the threshold for the decision rule.

These four steps are repeated until the weighted set of parameters, interpreted as the approximate posterior distribution, is ‘sufficiently close’ to the true posterior distribution. The steps (1) and (4) are usually quite fast, compared to steps (2) and (3), which are the computationally expensive parts.

These ABC algorithms can be generally classified into two groups based on the decision rule in step (2). In the first group, we simulate

using the perturbed parameter and accept it if , an adaptively chosen threshold. Otherwise we continue until we get an accepted perturbed parameter. For the second group of algorithms, we do not have this ‘explicit acceptance’ step but rather a probabilistic one. Here we accept the perturbed parameter with a probability that depends on ; if it is not accepted, we keep the present value of the parameter. The algorithms belonging to the ‘explicit acceptance’ group are RejectionABC (Tavaré et al., 1997) and PMCABC (Beaumont, 2010), whereas the algorithms in the ‘probabilistic acceptance’ group are SMCABC (Del Moral et al., 2012), RSMCABC (Drovandi and Pettitt, 2011), APMCABC (Lenormand et al., 2013), SABC (Albert et al., 2015) and ABCsubsim (Chiachio et al., 2014). For an ‘explicit acceptance’ to occur, it may take different amounts of time for different perturbed parameters (more repeated steps are needed if the proposed parameter value is distant from the true parameter value). Hence the first group of algorithms are inherently imbalanced. We notice that an ABC algorithm with ‘probabilistic acceptance’ do not have the similar issue of imbalance as a probabilistic acceptance step takes approximately the same amount of time for each parameter.

The generation of from the model, for a given parameter value, usually takes up huge amounts of computational resources (e.g. 10 minutes for the platelets deposition model in this paper). Hence, we want to choose an algorithm with faster convergence to the posterior distribution with minimal number of required forward simulations. For this work we choose Simulated Annealing ABC (SABC) which uses a probabilistic decision rule in Step (2) and needs minimal number of forward simulation than other algorithms as shown in Albert et al. (2015). As all tasks of SABC in Step (2) can be run independently, in our recent work (Dutta et al., 2017a), we have adapted SABC for HPC environment. Our implementation is available in Python package ABCpy and shows a linear scalability.

We further note that the parallelization schemes in ABCpy were primarily meant for inferring parameters from models, for which forward simulation takes almost equal time for any values of . Due to the complex stochastic nature of the numerical model, forward simulation time for different values of , can be quite variable. To solve this imbalance in the forward simulation, additionally to the imbalance reported for ABC algorithms, we use a new dynamic allocation scheme for MPI developed in Dutta et al. (2017b).

Figure 3: Comparison of work-flow between MPI (left) and dynamic-MPI backend (right).
Figure 4: Imbalance of ABC algorithms using MPI(straight-forward) (left) and MPI(dynamic-allocation) backend (right).

3.2 Dynamic Allocation for MPI

Here we briefly discuss how a dynamic allocation strategy for map-reduce provides better balancing of ABC algorithms compared to a straightforward allocation approach.

In the straightforward approach, the allocation scheme initially distributes tasks to executors, sends the map function to each executor, which in turn applies the map function, one after the other, to its map tasks. This approach is visualized in Figure 3, where a chunk represents the set of map tasks. For example, if we want to draw samples from the posterior distribution and we have cores available, at each step of SABC we create groups of 100 parameters and each group is assigned to one individual core.

On the other hand, the dynamic allocation scheme initially distributes tasks to the executors, sends the map function to each executor, which in turn applies it to the single task available. In contrast to the straightforward allocation, the executor requests a new map task as soon as the old one is terminated. This clearly results in a better balance of the work. The dynamic allocation strategy is an implementation of the famous greedy algorithm for job-shop scheduling, which can be shown to have an overall processing time (makespan) up to twice as better than the best makespan (Graham, 1966).

This approach is illustrated in Figure 3, reproduced from Dutta et al. (2017b). The unbalanced behavior is apparent if we visualize the run time of the individual map tasks on each executor. In Figure 4, the individual map tasks processing time is shown for an ABC algorithm performing inference on a weather prediction model, reported in Dutta et al. (2017b). Each row corresponds to an executor (or rank) and each bar corresponds to the total time spent on all tasks assigned to the respective rank (row) for one map call. For the straightforward allocation strategy, one can easily verify that most of the ranks finish their map tasks in half the time of the slowest rank. This clearly leads to large inefficiencies. Conversely, using the dynamic allocation strategy, the work is more evenly distributed across the ranks. The unbalancedness is not a problem that can be overcome easily by adding resources, rather speed-up and efficiency can drop drastically compared to the dynamic allocation strategy with increasing number of executors. For a detailed description and comparison, we direct readers to Dutta et al. (2017b).

3.3 Posterior Inference

Using SABC within HPC framework implemented in ABCpy (Dutta et al., 2017a), we draw samples approximating the posterior distribution , while keeping all the tuning parameters for the SABC fixed at the default values suggested in ABCpy package, except the number of steps and the acceptance rate cutoff, which was chosen respectively as 30 and . The parallelized SABC algorithm, using HPC makes it possible to perform the computation in 5 hours (using 140 nodes with 36-core of Piz Daint Cray architecture (Intel Broadwell + NVidia TESLA P100)), which would have been impossible by a sequential algorithm. To perform SABC for the platelets deposition model, the summary statistics extracted from the dataset, discrepancy measure between the summary statistics, prior distribution of parameters and perturbation Kernel to explore the parameter space for inference are described next.

Summary statistics

Given a dataset, , we compute an array of summary statistics.

defined as following,

  • , mean over time.

  • , variance over time.

  • , auto-correlation with lag 1.

  • , correlation between different pairs of variables over time.

  • , cross-correlation with lag 1 between different pairs of variables over time.

The summary statistics, described above, are chosen to capture the mean values, variances and the intra- and inter- dependence of different variables of the time-series over time.

Discrepancy measure:

Assuming the above summary statistics contain the most essential information about the likelihood function of the simulator-based model, we compute Bhattacharya-coefficient (Bhattachayya, 1943) for each of the variables present in the time-series using their mean and variance and Euclidean distances between different inter- and intra- correlations computed over time. Finally we take a mean of these discrepancies, such that, in the final discrepancy measure discrepancy between each of the summaries are equally weighted. The discrepancy measure between two datasets, and can be specified as,

where is the Bhattacharya-coefficient (Bhattachayya, 1943) and . Further, we notice the value of the discrepancy measure is always bounded in the closed interval .


We consider independent Uniform distributions for the parameters with a pre-specified range for each of them, , , , and .

Perturbation Kernel:

To explore the parameter space of

, we consider a five-dimensional truncated multivariate Gaussian distribution as the perturbation kernel. SABC inference scheme centers the perturbation kernel at the sample it is perturbing and updates the variance-covariance matrix of the perturbation kernel based on the samples learned from the previous step.

3.4 Parameter estimation

Given experimentally collected platelet deposition dataset , our main interest is to estimate a value for . In decision theory, Bayes estimator minimizes posterior expected loss,

for an already chosen loss-function

. If we have samples from the posterior distribution , the Bayes estimator can be approximated as,


As we consider the Euclidean loss-function as the loss-function, the approximate Bayes-estimator can be shown to be .

4 Inference on experimental dataset

Figure 5: Marginal posterior distribution (black-dashed) and Bayes Estimate (back-solid) of

for collective dataset generated from of 7 patients. The smoothed marginal distribution is created by a Gaussian-kernel density estimator on 5000 i.i.d. samples drawn from the posterior distribution using SABC. The (gray-solid) line indicates the manually estimated values of the parameters in 

Chopard et al. (2017).

The performance of the inference scheme described in Section 3 is reported here, for a collective dataset created from the experimental study of platelets deposition of 7 blood-donors. The collective dataset was created by a simple average of

over 7 donors at each time-point . In Figure 5, we show the Bayes estimate (black-solid) and the marginal posterior distribution (black-dashed) of each of the five parameters computed using 5000 samples drawn from the posterior distribution using SABC. For comparison, we also plot the manually estimated values of the parameters (gray-solid) in Chopard et al. (2017)

. We notice that the Bayes estimates are in a close proximity of the manually estimated values of the parameters and also the manually estimated values observe a significantly high posterior probability. This shows that, through the means of ABC we can get an estimate or quantify uncertainty of the parameters in platelets deposition model which is as good as the manually estimated ones, if not better.

Next we do a Posterior predictive check to validate our model and inference scheme. The main goal here is to analyze the degree to which the experimental data deviate from the data generated from the inferred posterior distribution of the parameters. Hence we want to generate data from the model using parameters drawn from the posterior distribution. To do so, we first draw 100 parameter samples from the inferred approximate posterior distribution and simulate 100 data sets, each using a different parameter sample. We call this simulated dataset as the predicted dataset from our inferred posterior distribution and present the mean predicted dataset (blue-solid) compared with experimental dataset (black-solid) in Figure 6

. Note that since we are dealing with the posterior distribution, we can also quantify uncertainty in our predictions. We plot the 1/4-th quantile, 3/4-th quantile (red-dashed), minimum and maximum (gray-dashed) of the predicted dataset at each timepoints to get a sense of uncertainty in the prediction. Here we see a very good agreement between the mean predicted dataset and the experimentally observed one, while the

-th and -th quantile of the prediction being very tight. This shows a very good prediction performance of the numerical model of platelet deposition and the proposed inference scheme.

Figure 6: Posterior Prediction Check: To validate the numerical model of the platelet deposition and the inference scheme we perform a posterior prediction check by simulating 100 datasets, each using a different parameter sample drawn from the posterior distribution. Here, we plot the experimental dataset (black-solid) used for inference, mean predicted dataset (blue-solid), -th and -th quantile (red-dashed), minimum and maximum (gray-dashed) of the predicted datasets at each timepoints.

Additionally, to point the strength of having a posterior distribution for the parameters we compute and show the posterior correlation matrix between the 5 parameters in Figure 7, highlighting a strong negative correlation between , strong positive correlations between and . A detailed investigation of these correlation structure would be needed to understand them better, but generally they may point towards: a) the stochastic nature of the considered model for platelet deposition and b) the fact that the deposition process is an antagonistic or synergetic combination of the mechanisms proposed in the model.

Note finally that the posterior distribution being the joint probability distribution of the 5 parameters, we can also compute any higher-order moments, skewness etc. of the parameters for a detailed statistical investigation of the natural phenomenon.

Figure 7: Posterior correlation matrix of computed from the 5000 i.i.d. samples drawn from the posterior distribution using SABC.

5 Conclusions

Here, we have demonstrated that approximate Bayesian computation (ABC) can be used to automatically explore the parameter space of the numerical model simulating the deposition of platelets subject to a shear flow as proposed in Chopard et al. (2015, 2017)

. We also notice the good agreement between the manually tuned parameters and the Bayes estimates, while saving us from subjectivity and a tedious manual tuning. This approach can be applied patient per patient, in a systematic way, without the bias of a human operator. In addition, the approach is computationally fast enough to provide results in an acceptable time for contributing to a new medical diagnosis, by giving clinical information that no other known method can provide. The clinical relevance of this approach is still to be explored and our next step will be to apply our approach at a personalized level, with a cohort of patients with known pathologies. The possibility of designing new platelet functionality test as proposed here is the result of combining different techniques: advanced microscopic observation techniques, bottom-up numerical modeling and simulations, recent data-science development and high performance computing (HPC).

Additionally, the ABC inference scheme provides us with a posterior distribution of the parameters given observed dataset, which is much more informative about the underlying process. The posterior correlations structure shown in Fig. 7 may not have a direct biophysical interpretation, though it illustrates some sort of underlying and unexplored stochastic mechanism for further investigation. Finally we note that, although the manual estimates achieve a very high posterior probability, they are different from the Bayes estimates learned using ABC. The departure reflects a different estimation of the quality of the match between experimental observation and simulation results. As the ABC algorithms are dependent on the choice of the summary statistics and the discrepancy measures, the parameter uncertainty quantified by SABC in Section 4 or the Bayes estimates computed are dependent on the assumptions in Section 3.3 regarding their choice. Fortunately there are recent works on automatic choice of summary statistics and discrepancy measures in ABC setup (Gutmann et al., 2017), and incorporating some of these approaches in our inference scheme is a promising direction for future research in this area.

Author Contribution

Design of the research: RD, BC, AM; Performed research: RD; Experimental data collection: KZB, FD; writing of the paper: RD, BC; Contribution to the writing: AM, KZB, JL, FD, AM; design and coding of the numerical forward simulation model: BC, JL.


This study conforms with the Declaration of Helsinki and its protocol was approved by the Ethics Committee of ‘CHU de Charleroi’ (comité d’éthique OM008). All volunteers gave their written informed consent.


Ritabrata Dutta and Antonietta Mira are supported by Swiss National Science Foundation Grant No. 105218_163196 (Statistical Inference on Large-Scale Mechanistic Network Models). We thanks CADMOS for providing computing resources at the Swiss Super Computing Center. We acknowledge partial funding from the European Union Horizon 2020 research and innovation programme for the CompBioMed project ( under grant agreement 675451.


We thank Dr. Marcel Schoengens, CSCS, ETH Zurich for helps regarding. We thank CHU Charleroi for supporting the experimental work used in this study.


  • Affeld et al. [2013] Klaus Affeld, Leonid Goubergrits, Nobuo Watanabe, and Ulrich Kertzscher. Numerical and experimental evaluation of platelet deposition to collagen coated surface at low shear rates. Journal of biomechanics, 46(2):430–436, 2013.
  • Albert et al. [2015] Carlo Albert, R. Künsch Hans, and Andreas Scheidegger. A simulated annealing approach to approximate Bayesian computations. Statistics and Computing, 25:1217–1232, 2015.
  • Beaumont [2010] Mark A Beaumont. Approximate bayesian computation in evolution and ecology. Annual review of ecology, evolution, and systematics, 41:379–406, 2010.
  • Bhattachayya [1943] A Bhattachayya. On a measure of divergence between two statistical population defined by their population distributions. Bulletin Calcutta Mathematical Society, 35(99-109):28, 1943.
  • Bonello et al. [2009] Laurent Bonello, Laurence Camoin-Jau, Sébastien Armero, Olivier Com, Stéphane Arques, Caroline Burignat-Bonello, Marie-Paule Giacomoni, Roland Bonello, Frédéric Collet, Philippe Rossi, et al. Tailored clopidogrel loading dose according to platelet reactivity monitoring to prevent acute and subacute stent thrombosis. American Journal of Cardiology, 103(1):5–10, 2009.
  • Born and Cross [1963] GVR Born and MoJ Cross. The aggregation of blood platelets. The Journal of physiology, 168(1):178–195, 1963.
  • Breet et al. [2010] Nicoline J Breet, Jochem W van Werkum, Heleen J Bouman, Johannes C Kelder, Henk JT Ruven, Egbert T Bal, Vera H Deneer, Ankie M Harmsze, Jan AS van der Heyden, Benno JWM Rensing, et al. Comparison of platelet function tests in predicting clinical outcome in patients undergoing coronary stent implantation. Jama, 303(8):754–762, 2010.
  • Chiachio et al. [2014] Manuel Chiachio, James L Beck, Juan Chiachio, and Guillermo Rus. Approximate bayesian computation by subset simulation. SIAM Journal on Scientific Computing, 36(3):A1339–A1358, 2014.
  • Chopard et al. [2015] B. Chopard, D. Ribeiro de Sousa, J. Latt, F. Dubois, C. Yourassowsky, P. Van Antwerpen, O. Eker, L. Vanhamme, D. Perez-Morga, G. Courbebaisse, and K. Zouaoui Boudjeltia. A physical description of the adhesion and aggregation of platelets. ArXiv e-prints, 2015.
  • Chopard et al. [2017] Bastien Chopard, Daniel Ribeiro de Sousa, Jonas Lätt, Lampros Mountrakis, Frank Dubois, Catherine Yourassowsky, Pierre Van Antwerpen, Omer Eker, Luc Vanhamme, David Perez-Morga, et al. A physical description of the adhesion and aggregation of platelets. Royal Society Open Science, 4(4):170219, 2017.
  • Del Moral et al. [2012] Pierre Del Moral, Arnaud Doucet, and Ajay Jasra. An adaptive sequential monte carlo method for approximate bayesian computation. Statistics and Computing, 22(5):1009–1020, 2012.
  • Drovandi and Pettitt [2011] Christopher C Drovandi and Anthony N Pettitt. Estimation of parameters for macroparasite population evolution using approximate bayesian computation. Biometrics, 67(1):225–233, 2011.
  • Dutta et al. [2017a] R Dutta, M Schoengens, J.P. Onnela, and Antonietta Mira. Abcpy: A user-friendly, extensible, and parallel library for approximate bayesian computation. In Proceedings of the Platform for Advanced Scientific Computing Conference. ACM, June 2017a.
  • Dutta et al. [2017b] Ritabrata Dutta, Marcel Schoengens, Avinash Ummadisingu, Jukka-Pekka Onnela, and Antonietta Mira. Abcpy: A high-performance computing perspective to approximate bayesian computation. arXiv preprint arXiv:1711.04694, 2017b.
  • Fontaine et al. [2009] Eustace Fontaine, Richard Warwick, Priya Sastry, and Michael Poullis. Effect of foreign surface pacification with albumin, aprotinin, propofol, and high-density lipoprotein. The Journal of extra-corporeal technology, 41(1):3, 2009.
  • Geisler et al. [2006] Tobias Geisler, Harald Langer, Magdalena Wydymus, Katrin Göhring, Christine Zürn, Boris Bigalke, Konstantinos Stellos, Andreas E May, and Meinrad Gawaz. Low response to clopidogrel is associated with cardiovascular outcome after coronary stent implantation. European heart journal, 27(20):2420–2425, 2006.
  • Graham [1966] Ronald L Graham. Bounds for certain multiprocessing anomalies. Bell Labs Technical Journal, 45(9):1563–1581, 1966.
  • Gurbel et al. [2005] Paul A Gurbel, Kevin P Bliden, Kirk Guyer, Peter W Cho, Kazi A Zaman, Rolf P Kreutz, Ashwani K Bassi, and Udaya S Tantry. Platelet reactivity in patients and recurrent events post-stenting: results of the prepare post-stenting study. Journal of the American College of Cardiology, 46(10):1820–1826, 2005.
  • Gutmann et al. [2017] Michael U Gutmann, Ritabrata Dutta, Samuel Kaski, and Jukka Corander. Likelihood-free inference via classification. Statistics and Computing, pages 1–15, 2017.
  • Harrison [2009] P Harrison. Assessment of platelet function in the laboratory. Hämostaseologie, 29(01):25–31, 2009.
  • Hochholzer et al. [2006] Willibald Hochholzer, Dietmar Trenk, Hans-Peter Bestehorn, Benjamin Fischer, Christian M Valina, Miroslaw Ferenc, Michael Gick, Angelika Caputo, Heinz Joachim Büttner, and Franz-Josef Neumann. Impact of the degree of peri-interventional platelet inhibition after loading with clopidogrel on early clinical outcome of elective coronary stent placement. Journal of the American College of Cardiology, 48(9):1742–1750, 2006.
  • Koltai et al. [2017] Katalin Koltai, Gabor Kesmarky, Gergely Feher, Antal Tibold, and Kalman Toth. Platelet aggregometry testing: Molecular mechanisms, techniques and clinical implications. International journal of molecular sciences, 18(8):1803, 2017.
  • Lenormand et al. [2013] Maxime Lenormand, Franck Jabot, and Guillaume Deffuant. Adaptive approximate bayesian computation for complex models. Computational Statistics, 28(6):2777–2796, 2013.
  • Marcucci et al. [2009] Rossella Marcucci, Anna Maria Gori, Rita Paniccia, Betti Giusti, Serafina Valente, Cristina Giglioli, Piergiovanni Buonamici, David Antoniucci, Rosanna Abbate, and Gian Franco Gensini. Cardiovascular death and nonfatal myocardial infarction in acute coronary syndrome patients receiving coronary stenting are predicted by residual platelet reactivity to adp detected by a point-of-care assay: a 12-month follow-up. Circulation, 119(2):237–242, 2009.
  • Marin et al. [2012] Jean-Michel Marin, Pierre Pudlo, ChristianP. Robert, and RobinJ. Ryder. Approximate Bayesian computational methods. Statistics and Computing, 22(6):1167–1180, 2012. ISSN 0960-3174. doi: 10.1007/s11222-011-9288-2. URL
  • Matetzky et al. [2004] Shlomi Matetzky, Boris Shenkman, Victor Guetta, Michael Shechter, Roy Beinart, Ilan Goldenberg, Ilya Novikov, Hanna Pres, Naphtali Savion, David Varon, et al. Clopidogrel resistance is associated with increased risk of recurrent atherothrombotic events in patients with acute myocardial infarction. Circulation, 109(25):3171–3175, 2004.
  • Maxwell et al. [2007] Mhairi J Maxwell, Erik Westein, Warwick S Nesbitt, Simon Giuliano, Sacha M Dopheide, and Shaun P Jackson. Identification of a 2-stage platelet aggregation process mediating shear-dependent thrombus formation. Blood, 109(2):566–576, 2007.
  • Mehta et al. [2001] Shamir R Mehta, Salim Yusuf, Ron JG Peters, Michel E Bertrand, Basil S Lewis, Madhu K Natarajan, Klas Malmberg, Hans-Jürgen Rupprecht, Feng Zhao, Susan Chrolavicius, et al. Effects of pretreatment with clopidogrel and aspirin followed by long-term therapy in patients undergoing percutaneous coronary intervention: the pci-cure study. The Lancet, 358(9281):527–533, 2001.
  • Michelson et al. [2002] Alan D Michelson, Marc R Barnard, Lori A Krueger, AL Frelinger III, and Mark I Furman. Flow cytometry. In Michelson AD, editor, Platelets, pages 297–315. Academic Press, 2002.
  • Organization [2015] World Health Organization., 2015.
  • Piagnerelli et al. [2007] Michael Piagnerelli, K Zouaoui Boudjeltia, Dany Brohée, Anne Vereerstraeten, Pietrina Piro, Jean-Louis Vincent, and Michel Vanhaeverbeek. Assessment of erythrocyte shape by flow cytometry techniques. Journal of clinical pathology, 60(5):549–554, 2007.
  • Picker [2011] S.M. Picker. In-vitro assessment of platelet function. Transfus Apher Sci, 44:305–19, 2011.
  • Price et al. [2008] Matthew J Price, Sarah Endemann, Raghava R Gollapudi, Rafael Valencia, Curtiss T Stinis, Justin P Levisay, Alissa Ernst, Neil S Sawhney, Richard A Schatz, and Paul S Teirstein. Prognostic significance of post-clopidogrel platelet reactivity assessed by a point-of-care assay on thrombotic events after drug-eluting stent implantation. European Heart Journal, 29(8):992–1000, 2008.
  • Reasor et al. [2013] Daniel A Reasor, Marmar Mehrabadi, David N Ku, and Cyrus K Aidun. Determination of critical parameters in platelet margination. Annals of biomedical engineering, 41(2):238–249, 2013.
  • Remuzzi and Boccardo [1993] A Remuzzi and P Boccardo. Albumin treatment reduces in vitro platelet deposition to pmma dialysis membrane. The International journal of artificial organs, 16(3):128–131, 1993.
  • Robert and Casella [2005] Christian P. Robert and George Casella. Monte Carlo Statistical Methods (Springer Texts in Statistics). Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2005. ISBN 0387212396.
  • Sharma et al. [1981] NC Sharma, SF Mohammad, HY Chuang, and RG Mason. Albumin-igg complexes in human serum and plasma that inhibit blood platelet adhesion. Proceedings of the National Academy of Sciences, 78(12):7750–7753, 1981.
  • Shenkman et al. [2008] B Shenkman, Y Einav, O Salomon, D Varon, and N Savion. Testing agonist-induced platelet aggregation by the impact-r [cone and plate (let) analyzer (cpa)]. Platelets, 19(6):440–446, 2008.
  • Sibbing et al. [2009] Dirk Sibbing, Siegmund Braun, Tanja Morath, Julinda Mehilli, Wolfgang Vogt, Albert Schömig, Adnan Kastrati, and Nicolas von Beckerath. Platelet reactivity after clopidogrel treatment assessed with point-of-care analysis and early drug-eluting stent thrombosis. Journal of the American College of Cardiology, 53(10):849–856, 2009.
  • Steinhubl et al. [2002] Steven R Steinhubl, Peter B Berger, J Tift Mann III, Edward TA Fry, Augustin DeLago, Charles Wilmer, Eric J Topol, Credo Investigators, et al. Early and sustained dual oral antiplatelet therapy following percutaneous coronary intervention: a randomized controlled trial. Jama, 288(19):2411–2420, 2002.
  • Tavaré et al. [1997] Simon Tavaré, David J Balding, Robert C Griffiths, and Peter Donnelly. Inferring coalescence times from dna sequence data. Genetics, 145(2):505–518, 1997.
  • Velik-Salchner et al. [2008] Corinna Velik-Salchner, Stephan Maier, Petra Innerhofer, Werner Streif, Anton Klingler, Christian Kolbitsch, and Dietmar Fries. Point-of-care whole blood impedance aggregometry versus classical light transmission aggregometry for detecting aspirin and clopidogrel: the results of a pilot study. Anesthesia & Analgesia, 107(6):1798–1806, 2008.