Quantitative dynamics of reversible platelet aggregation: mathematical modelling and experiments

Although reversible platelet aggregation observed in response to ADP stimulation in the presence of calcium is a well-known phenomenon, its mechanisms are not entirely clear. To study them, we developed a simple kinetic mass-action-law-based mathematical model to use it in combination with experiments. Light transmission platelet aggregometry (LTA) induced by ADP was performed for platelet-rich plasma or washed platelets using both conventional light transmission and aggregate size monitoring method based on optical density fluctuations. Parameter values of the model were determined by means of parameter estimation techniques implemented in COPASI software. The mathematical model was able to describe reversible platelet aggregation LTA curves without assuming changes in platelet aggregation parameters over time, but with the assumption that platelet can enter the aggregate only once. In the model, the mean size of platelet aggregates correlated with the solution transparency. This corresponded with flow cytometry analysis and with optical density fluctuations data on aggregate size. The predicted values of model parameters correlated with ADP concentration used in experiments. These data suggest that, at the start of the aggregation, when platelet integrins switch “on”, large unstable platelet aggregates are rapidly formed, which leads to an increase in light transmission. However, upon fragmentation of these aggregates, the probability of the post-aggregate platelets’ attachment to each other decreases preventing new aggregation and resulting in the reversible aggregation phenomenon.

www.nature.com/scientificreports www.nature.com/scientificreports/ Although aggregometry is a routine straightforward technique, direct interpretation of its results is lacking. One of the unexplained phenomena is that platelet aggregation in response to the main activator ADP seems to be weaker in the presence of calcium ions in the solution (which is the correct physiological condition) than without them [10][11][12][13][14] (Fig. 1). Platelet aggregation in these conditions is always reversible (Fig. 1A,B). One of the possible explanations is the inactivation of α IIb β 3 that leads to platelet aggregate disassembly 6,9,13,15 . On the other hand, the breakdown of platelet aggregates could originate from some mechanistic reasons, such as turbulent flows in the aggregometric cuvette 1,2,16,17 . Understanding of these processes could improve the diagnostic potential of LTA.
To obtain platelet functional parameters from the complex shapes of the aggregation curves, a computational modelling approach could be utilized. The computational approach most widely used for particle aggregation tasks is Smoluchowski coagulation-fragmentation equation 18 , where some fragmentation and coagulation "kernels" are constructed and then solved by Smoluchowski equation [18][19][20][21][22][23] . Although the Smoluchowski approach could be the most "correct" one, it is quite computationally expensive when solved with deterministic methods, and still contain too many unknown parameters when solved with stochastic methods 24,25 and, therefore, computational modelling of aggregometry data with Smoluchowski equation could not be readily performed alongside with personalized parameter estimation yet 26 . Another commonly used approach is the lattice models 27,28 , which is based on a construction of a lattice with a particle in each cell of the mesh. The aggregation is modelled by the energy-profit-driven movement of particles between the cells 27 . Although lattice-based models are useful, their complexity and computational cost are too high for three-dimensional mesh, which should be applied for aggregometry modelling. The most abundant are simplistic models, based on mass-action laws of kinetics with included description of particles clustering and aggregation 29,30 . However, most of them were developed for irreversible aggregation, and therefore focus on the level of aggregation as the main parameter, and pre-define which parameter values should depend on the level of activation 31,32 .
Here we suggest that the crucial limitation of the existing models is the external introduction of the effect of activator on the laws of aggregation, and we propose to deduce parameters of platelet activation from the experimental data by means of parameter estimation techniques. In the current study, this approach allowed us to describe experimental data on the reversible platelet aggregation and conclude that the origin of the reversibility lies in the instability of large aggregates rather than in the loss of aggregatory ability as a function of time. This conclusion was confirmed by flow cytometry and optical density fluctuations analysis 33 of the aggregate size distribution in the course of platelet aggregation.

Results
Construction of the mathematical model of platelet aggregation in response to ADp. Platelet aggregation in response to ADP in the presence of calcium ions in the medium is reversible (Fig. 1) 34,35 . There was a dose-dependent increase of reversible aggregation in heparinized or hirudinated platelet rich plasma (PRP) in response to ADP, with transition reversible/irreversible aggregation at 1.25-5 µM (Fig. 1A). Washed platelets with 200 µg/ml fibrinogen reacted similarly (Fig. 1B). Calcium ions chelation by sodium citrate led to an irreversible aggregation (Fig. 1C).
Movement of individual platelets and aggregates in the cuvette of the aggregometer could not be described analytically and requires excessive computer time to model numerically 17 . However, the rate of hydrodynamic flows in an aggregometry cuvette without platelets could be calculated (Fig. S1). Our calculations together with previous reports 36 show that there is no significant dependence of the shear rate on the rpm in the working range of the aggregometer (800-1200 rpm). We also confirmed these results experimentally (Fig. S2). Surprisingly, it appears that neither stirring rate, nor the measurement height significantly affect the behaviour of the aggregation curve (Fig. S2). Accordingly, we assumed that the sole purpose of fluid flows in the cuvette is an enhancement of cell-cell interactions with each other and with aggregates. The event of aggregate formation then can be described by a probability constant. In response to ADP fibrinogen binding to a single platelet increases step-like (Fig. S3), so we can assume that the probabilities of platelet attachment to each other do not change with time. In case www.nature.com/scientificreports www.nature.com/scientificreports/ of secondary platelet activation, for example, due to thromboxane A2 synthesis 37 , the assumption of constant attachment probabilities will not hold and in such case the model will not be applicable. In the real solution, the platelets can be found in numerous states: single platelets with or without filopodia and aggregates of various size and density. However, the simplest approximation is that a significant difference exists between the states of single and aggregated platelets 16 . Then the simplest mathematical model of platelet aggregation in suspension will be: where p is the concentration of single platelets, n is the concentration of aggregates, k 2 is the probability of new aggregate formation from two single platelets, k 1 is the probability of another platelet attachment to an existing aggregate, k −1 is the probability of single platelet detachment from an aggregate, k −2 is the probability of formation of one aggregate from two existing ones and k 3 is the probability of an aggregate fragmenting into two ( Fig. 2A). All probabilities k i theoretically depend on the average size of aggregates, which can be calculated as estimation of model parameters from LtA data. We then proceeded to the estimation of the model parameters from the experimental data. We used five different parameter estimation techniques for each experimental dataset to ensure unbiased parameter estimation. In Fig. 2B, typical parameter estimation results are shown. The output parameter of the aggregometer (percentage of aggregation) could be compared with the concentration of light-scattering particles in the solution according to the following equation 30 : where parameter a is introduced to describe platelet shape change upon activation, it is considered to be linearly increasing with time in the first two seconds of the simulation.
Analysis of the experimental data shows that ADP concentration influences all parameters of the model ( Fig. 2C with legend), with the probability of aggregate formation (k 1 ) increasing with ADP concentration increase for the particular donor, and the parameters of single aggregate stability (k −2 , k 3 ) decreasing with it (note that weak response to 1.25 µM ADP slightly falls out of these trends). With the mathematical model presented here, we have tried to perform comparison between platelet aggregation in PRP and in buffer solution (Fig. S4). However, we could not reach any reliable conclusion except that the parameter of aggregate stability (k 3 ) increases for platelet-rich plasma. Analysis of the validated model revealed that the aggregate size increased from 1 to 2-10 in the first minute of experiment and then dropped to 2-1 ( Fig. 2D,E). Meanwhile the corresponding concentration of aggregates increases constantly with time and the corresponding concentration of free platelets decreases with it ( Fig. 2D,E).
Aggregate size distribution over the course of aggregation. The conclusion of the model that aggregation peak corresponds to the peak in aggregate size, while assuming that the adhering properties of single platelets do not change, can be tested by means of flow cytometry. The prediction that the aggregation peak corresponds to an increase in aggregate size, was investigated by comparison between side scattering-forward scattering dot plots of samples taken from the aggregometer (Fig. 3, S5). The distribution of fluorescence intensity (that we assume to mimic the distribution of aggregate sizes) appeared to have several characteristic ranges. At the aggregation curve maximum, about 60% of platelets remained single, 25% were in small aggregates of 2-4 cells, and 15% were in large aggregates of 4-8 cells. After that, the small aggregates began to disappear, while the large ones seemed even to increase in size (up to 10-20 cells). The large aggregates disappeared completely over the further course of disaggregation (Fig. 3), and only a small number of smaller aggregates remained. To support this conclusion, we also performed dark-field microscopy of samples from the aggregometry cuvette at various time points and confirmed that large aggregates appear in the cuvette only at the peak of aggregation (Fig. S6). In addition, we have extended the model to describe dynamics of aggregates of various sizes. It contains the same number of unknown parameters, but increased number of variables. This model demonstrates, that at the aggregational peak most platelets are in aggregates of size 2-10 and only a few in aggregates conformed of >10 cells. All these data lead us to the conclusion that size distribution of platelets is only significant in small range (1-10) and so does not dramatically affect the model, described in (1-3) (Fig. S7). simultaneous measurement of platelet aggregation and mean aggregate size. Another support for the predicted aggregate size dynamics comes from a version of the LTA technique, where the fluctuations of www.nature.com/scientificreports www.nature.com/scientificreports/ the light transmission signal are used to estimate a relative size of light-dissipating particles 38 . The LTA signal received in Biola aggregometer possesses the same features as all other LTA data ( Fig. 4A for standard platelet concentration, Fig. S8 for low platelet concentration (100 000/μl) and high activator concentrations (20-50 μM)) and could be successfully described by the model (1)-(3) (Fig. 4A, dashed curves). Predicted by the corresponding models time-courses for the aggregate size are given in Fig. 4B. Calculated by the Biola device time-course of mean size is given in Fig. 4C. It could be seen by comparison of panels B and C that although the main features of the experimental and theoretical size distribution are similar, the particulars are different. The origin of such discrepancy could be the unknown value for the platelet concentration in the solution. As the next step, we allowed it to be estimated by the model, where both aggregation data and the mean size data were described simultaneously ("dual" model in Fig. 4). It appeared that the model could successfully describe both types of data, although with less accuracy (Fig. 4, Table 1). However, only the parameters of the "dual" model highly correlate with ADP concentration (Fig. 4D,F) and with each other (Fig. 4G). This suggests that they are not independent, but rather represent a combination of a smaller number of some basic functions of platelet activation. www.nature.com/scientificreports www.nature.com/scientificreports/ To assess comparative aggregation potential of single platelets during platelet aggregation, we performed aggregation test of Fura-Red loaded washed human platelets in presence of FITC-labelled fibrinogen (Figs S9 and S10). The relative aggregation potential was assessed as the ratio of fluorescence of fibrinogen-FITC and Fura-Red for each event. The fibrinogen binding to single platelets increases upon activation and was stable over time, while fibrinogen binding to aggregated platelets is significantly higher. Together these data indicate that the process of platelet disaggregation is determined by some instability in bonds between platelets in an aggregate.

Discussion
In the present study, we developed a mathematical model of reversible platelet aggregation which incorporated a novel mechanism of disaggregation, could be used for clinical parameter estimation, and revealed a dependence of the aggregation and aggregate stability parameters on ADP concentration. The experiments performed to test these predictions added several more pieces of information, including the dynamics of aggregate size distribution. The study altogether suggested that the loss of the ability to form aggregates is maybe not the cause, but rather the consequence of disaggregation, that is supported by experiments with fibrinogen binding to single platelets (Fig. S3).
The mathematical model (eqs 1 and 2, Fig. 2A) described processes of aggregate formation and fragmenting without taking into consideration possible heterogeneity of aggregates in size and other aspects. The parameters of the model were derived for each donor from the given experimental data by means of parameter estimation techniques. Assessment of additional experimental data is to be done. Changing between PRP and washed platelets also has an impact on the parameter values. There is even a difference in platelet aggregation in PRP with hirudin or heparin (Fig. 1). There are several possibilities for the origin of the more pronounced aggregation in platelet rich plasma compared to the washed platelet suspension. To test the most obvious ones we have adjusted platelet concentrations and fibrinogen concentration in the washed platelet suspension to that of PRP levels (Fig. S11). However, platelet aggregation still did not reach the PRP levels. Other possibilities like addition of platelet poor plasma to washed platelets could be tested in the future. Never the less, we expect that, given enough data, it would be possible to establish an explicit connection between values of the model parameters and features of the sample (ADP concentration, plasma presence, device features).
ADP induced aggregation in the presence of free calcium ions is clearly distinct from the one in the presence of calcium ion chelator (citrate). Furthermore, PRP recalcification resulted in a reversible aggregation, which was observed in hirudinated of heparinated PRP (Fig. 1C). However, some differences in aggregation between hirudinated and heparinated PRP could not be explained so readily. It has been reported that heparin is capable of inducing the micro-aggregate formation in PRP, which are potentiating further platelet aggregation as well as false-positive LTA results 39 . Furthermore, it has been demonstrated that heparin is a positive mediator of platelet integrin function 40 and thus, in heparinized PRP platelet aggregation might be additionally potentiated. Based on this, it can be claimed that hirudin is a more attractive option for aggregation evaluation in PRP.
The first mathematical model prediction is that the features of the aggregating particles (platelets or aggregates) do not change with time by themselves. The reverse course of aggregation originates from the transition of platelets between single and aggregated states. These prediction is partly supported by flow cytometry experiments performed here where the level of fibrinogen binding to single platelets was shown to be unchanged during aggregation. The increase in fibrinogen binding was observed for the platelet population as a whole, and was associated with disappearance of large aggregates (Fig. 3). Similar experiments with diluted platelets also demonstrate that platelet integrins do not inactivate within minutes (Fig. S3).
The mathematical model presented here is exceptionally simple and thus is not capable to describe fine features of processes in the aggregometer. As could be deduced from data in Fig. 3, a broad range of sizes of www.nature.com/scientificreports www.nature.com/scientificreports/ aggregates are present in the mixture during aggregation, while only the mean aggregate size is calculated in the model. The non-uniformity of aggregate size distribution has been shown several decades ago experimentally as well as theoretically 16,17 , but the results of the present study suggest that this distribution is not essential for correct description of platelet aggregation. An extension of the model proposed here to assess the dynamics of aggregates of different sizes supported this conclusion (Fig. S7). Another limitation, avoided in the majority of more profound mathematical models, is the inequality of the activation of platelets situated in different parts of the cuvette due to varying shear stresses 16,17 . In future work, we propose to derive intrinsic parameters of activation from experimental data by means of parameter estimation techniques instead of introducing them externally into the model, and the shear stress will be linked to one of the parameters. Another limitation of the model is that on long scale (more than ten minutes) the features of platelet activation start to change and the model with fixed parameter values cannot describe the data any more. To avoid this limitation either long-scale experiments should not be considered, or additional parameters should be introduced into the model.

Blood collection and platelet isolation. Healthy volunteers, both men and women aged between 18
and 39 years were recruited into the study. Investigations were performed in accordance with the Declaration of Helsinki and approved by CTP PCP RAS ethic committee, written informed consent was obtained from all donors. For platelet rich plasma, blood was collected into 9 ml tubes containing 3,8% sodium citrate (1:9 vol/vol) or into 9 ml tubes containing Li-heparin at 1600 g for 15 min as described previously 41 . Platelet washing procedure was described previously 42 . Briefly, blood was collected into 15 ml tubes containing 2.1 ml ACD. Platelets were purified by double centrifugation in the presence of PGI 2 , or PGI 2 and heparin in Tyrode's albumin buffer as described previously 42 . Final solution of washed platelets in Tyrode's buffer was supplemented with apyrase (0.1 U/mL) and rested for at least 30 minutes before the start of the experiments. All steps of washing procedure and all aggregometry experiments were conducted at 37 °C. PRP was diluted by PPP in order to achieve equal concentrations in all experiments. 27 healthy volunteers age from 18 to 39 (median age 25, 16 men, 11 women) were enrolled.
Aggregometry. Platelet aggregation was performed using Chrono-Log 490 and Biola LA-230 turbi-diametric aggregometers. Experiments were conducted in 250 μL (Chrono-Log) or 300 μL (Biola) aliquots of PRP or washed platelets. The platelet suspension was mixed by a stirrer with 800 rpm. ADP was added at various concentrations as the platelet activator. Washed platelets were pre-incubated with 200 µg/ml fibrinogen or without it for control sample. Before ADP activation of PRP, fibrinogen was not added. Tyrode's albumin buffer was used as a reference for washed platelets, and platelet poor plasma was used as a reference for PRP. Before measurements platelets were incubated at 37 °C. An optical signal was recorded every 0.5 s for Chrono-Log and 1 s for Biola.
Aggregate size analysis by Biola. Aggregate size analysis was performed using Biola LA-230 turbi-diametric aggregometer. The method is based on the Gabbasov'method of light transmission fluctuations 33,38 . The fluctuations arise due to a change in the number of particles in the suspension. The basic parameter of the method is calculated by the formula:  www.nature.com/scientificreports www.nature.com/scientificreports/ Number of equations = N where [j] is concentration of an aggregate of size j and other parameters are the same as described in Section 2.1. parameter estimation and model solution. Construction of computational model was based on the principles described in the Results section. The set of ordinary differential equations was integrated using the LSODA method implemented in COPASI software 43 . The method LSODA 44 is the numerical method to solve a system of ordinary differential equations with an ability to switch between explicit and implicit methods, which allows rapid integration of stiff problems (the equation dx/dt = ax presents a stiff problem). The COPASI (COmplex PAthway Simulator, copasi.org) is an open-source free software, where several methods for numerical integration and investigation of a system of ordinary differential equations are implemented. Model parameters were assesed from experimental data by means of the following teqniques implemented in COPASI software: Evolutionary Strategy (SRES) 45 , Particle Swarm 46 , Random Search, Hooke&Jeeves 47 or Levenberg-Marquardt 48 .