Numerical analysis of in vivo platelet consumption data from ITP patients

Background Numerical methods have recently allowed quantitative interpretation of in vivo murine platelet consumption data in terms of values for the random destruction rate constant (RD), intrinsic lifespan (LS), and the standard deviation of ln LS (SD), as well as the platelet production rate (PR) and age distribution (AD). But application of these methods to data obtained in thrombocytopenic patients is problematic for two reasons. First, such data has in all cases been obtained with radiolabeled platelets, and uptake of the radio-isotope by long lived cells complicates the analysis. Second, inferred values of the platelet production rate (PR) and random destruction rate (RD) are difficult to interpret, since increased RD can occur either as a cause or a consequence of thrombocytopenia. Methods We used a numerical method to analyze in vivo platelet consumption data from a series of 41 patients with immune thrombocytopenic purpura (ITP). An additional parameter, the fraction of labeled long-lived cells (LL), was evaluated concurrently with RD, LS, and SD. To provide a basis for interpreting these values, we used an iterative interpolation process to predict their response to different pathophysiologic mechanisms. The process also generates predicted effects on the widely used immature platelet fraction (IPF). Results Optimal parameter value sets were identified in 76 % (31 of 41) of the data sets. 27 of 31 ITP patients showed no substantial homeostatic increase in platelet production, with the remaining 4 showing both augmented platelet consumption and a compensatory increase in PR. Up to 1/3 of the patients showed the degree of increased RD expected to result from reduced thrombopoiesis only. “Jacknife” resampling yielded CV values of <0.5 in over 75 % of the evaluable data sets. Predicted platelet age distributions indicate that interpretation of the IPF and absolute IPF (aIPF) is a complex function of platelet count. We found, counter-intuitively, that reduced PR can increase the IPF, and increased RD can reduce the aIPF. Conclusions Our findings support the feasibility of using numerical analysis to quantitatively interpret in vivo platelet consumption data, to identify likely etiologies of thrombocytopenias, and to assess the utility of IPF measurements in that context. Electronic supplementary material The online version of this article (doi:10.1186/s12878-015-0034-4) contains supplementary material, which is available to authorized users.


Background
In vivo platelet consumption studies have often been used to quantify the rates of both random and lifespandependent consumption processes, and to evaluate production rate, but their interpretation can be problematic. The aim of the present study is to demonstrate that a numerical analysis method can reliably quantify these rates, and thereby identify fundamental pathophysiologic features, in thrombocytopenic patients.
Methods which have been used to interpret such studies include a simple exponential decay model [1], a weighted mean method that applies an empirical mixture of separately optimized linear and exponential decay processes [2], a purely lifespan-dependent model [3], the Mills-Dornhorst equation (which includes but does not solve for a random (exponential) consumption rate constant) [4,5], the widely used multiple hit model (based on a unique consumption mechanism for which there is little experimental support) [6,7], and combined use of the latter two approaches [8,9].
None of these methods allow concurrent modeling of the random (hemostatic and phagocyte-mediated) and lifespan-dependent processes known to result in most in vivo platelet consumption. Numerical analysis models have been designed for that purpose, and their utility for the analysis of murine platelet consumption data has been demonstrated [10,11].
It is not clear, however, whether these methods can be adapted to evaluate existing clinical data, since the latter in all cases involves tracking of radiolabeled platelets and includes contributions from uptake of the radioisotope by other, longer-lived cell types. It is also unclear how to translate the resultant kinetic parameter values into useful conclusions about how individual patients became thrombocytopenic. Here we have modified a previously described numerical model [10] to successfully analyze and interpret published data from a series of patients with immune thrombocytopenic purpura (ITP) [12].

Patients
Entry criteria for the ITP patients have been described previously [12]. The study, including informed consent procedures, was performed according to the principles outlined in the Declaration of Helsinki of 1975. Briefly, 111 Indium-labeled autologous platelet consumption data from 41 consecutive adult patients with prednisone non-responsive primary ITP was reviewed. The data was obtained either (A) at the time of diagnosis, or (B) after failure to sustain a platelet count response to prednisone treatment. Diagnostic criteria for all patients included exclusion of other malignant, metabolic, or pharmacologic causes, as well as causes of "secondary" ITP such as hepatitis C virus (HCV) infection. For those in group (A), rapid consumption of autologous 111 Indium-labeled platelets (interpreted via the multiple hit model [6]) was an additional diagnostic criterion. For those in group (B), demonstration of antiplatelet antibodies on the platelet surface via indirect immunofluorescence [13] was used for this purpose.
All of the patients in this study underwent splenectomy after failing to respond adequately to prednisone treatment. Patients were deemed to have had a complete response to splenectomy if their platelet count persistently exceeded 100 × 10 9 /L thereafter, with no significant bleeding episodes. They were considered "non-responders" if their subsequent platelet count did not exceed either 30 × 10 9 /L, or twice their baseline count, or if they had persistent significant bleeding episodes.

Platelet kinetics studies
These have been described in detail previously [12]. Briefly, platelet rich plasma was prepared by differential centrifugation, and platelets prepared by a subsequent high speed centrifugation were labeled with 111 Indium oxine by standard methods. Peripheral blood specimens obtained at 30 min after injection were considered "baseline" measurements (for patient 30, a 1.5 h time point was used), and all subsequent measurements were normalized to these for each patient. Equilibration with a pool of splenic platelets is thought to be complete well within this initial time frame [14]. All post-injection specimens for each patient were evaluated in a gamma counter at the same time to eliminate decay-related effects on recovery.

Numerical analysis
Data analysis was performed on desktop computers using Microsoft Excel. Baseline (initial) parameter ranges searched were RD 0-19 (resolution 1 %) %/h, LS 0-15.2 % (resolution 0.8 %), and SD 20-267 h (resolution 13 h). LL and RD ranges and resolutions were then empirically optimized to the values shown in Table 1. The equilibration metric was calculated as the net platelet count produced by the model at the midpoint of the equilibration phase divided by the net platelet count at the end of the equilibration phase, as described previously [10]. All searches achieved an equilibration metric value of > 0.997 (1000 interval equilibration phase, 0.5 h per interval). Searches of resampled data sets were performed over the same parameter ranges shown in Table 1. For cases in which all resampled data sets yielded the same parameter values as the complete data set at resolution "R", the upper limit of the standard deviation was estimated by the value obtained had one resampled data set yielded a parameter value one "R" range removed from that of the complete data set.

Modeling in vivo platelet turnover
The numerical model used here [10] posits that an in vivo platelet population can be visualized in a spreadsheet as a series of small platelet cohort concentrations. The cohorts are assumed to be produced at a constant rate (PR, K/ul/h) in short sequential time periods, and individually consumed, by both random and lifespandependent processes, at the end of each such time period. The consumption curve for individual cohorts is determined by a random destruction rate constant (RD, %/h), by the lognormally distributed cohort lifespan (LS, hr), and by the standard deviation of ln LS (SD). Population platelet consumption curves are generated by summing the cohort values at sequential time points. Optimal theoretical consumption curves generated by a large range of possible parameter values are identified via quantitative comparison to each data set (summed squared residual values, or SS).
For the 41 patient studies analyzed here, each patient's platelet consumption data was normalized to the first (baseline) measurement of circulating 111 In. Visual inspection of the data strongly suggests that many of the labeled platelet preparations contained long lived species, as others have described in similar studies [15,16]. This is evident in the plateau phase seen at late times in the consumption data (for example, patients 35 and 27, Fig. 4). To take this into account, we evaluated a fourth parameter: the fraction of the labeled cells/platelets consisting of long lived (species (LL, %). This parameter simply "shrinks" the scope of the analysis to the consumption of platelets from 100 % of the time zero value to an optimizable minimum percentage (LL). For modeling purposes, lifespans of these long lived species are assumed to be infinite.

Optimal parameter value search process
Optimal parameter value searches were performed as shown schematically in Fig. 1. For a given Fig. 2). This process is repeated over a range of 20 RD values, yielding in most cases single "volume" SS minima as shown in Fig. 3a. Finally, the entire process is repeated at a series of SD values, and the resultant volume minima are compared in order to identify a "global" minimum SS value and its associated parameter values. Distinguishable alternative volume minima showing SS values greater than those of the global minima were also seen in some data sets (see below). Searches were performed for only three SD values, as this generated a plausible range of distribution widths for the resultant lifespan-dependent consumption rates (see examples in Fig. 3b) while significantly reducing computation time. Examples of the consumption curves generated by the optimal parameter values are shown in Fig. 4.

Data quality evaluation and optimal parameter value search results
This process is outlined in Fig. 5. Of the 41 originally reported ITP patient data sets, one was excluded due to lack of initial time point data. One patient demonstrated an initial platelet clearance rate of 46 % in the first 1.5 h of the study (>4 standard deviations faster than the mean). Because this value suggests the type of platelet activation during labeling/processing that we have on occasion seen in murine platelet clearance studies (TS, unpublished), this data was also excluded. For the remainder, quality of parameter value optimization was evaluated in terms of the ratio of SS to n (the number of data points per patient data set), where n ranged from 5 to 9 (Table 1). One case (patient 31) with an SS/n value of 143 (over four standard deviations from the mean value of 18.5) was then excluded. No other SS/n values fell beyond two sd from the mean. A single "global" minimum SS value, with its associated (optimal) parameter values, was identified in 24 of the 39 evaluable data sets. The optimal consumption curves show a large amount of inter-patient variation, as the examples in Fig. 4a demonstrate. Of those showing more than one minimum, convincing global minima were identified in four data sets on the basis of goodness-of-fit. Specifically, the global minima in these cases showed SS values which were less than 50 % of those defining the alternative (local) minima. Three data sets showed global minima for which comparison of absolute vs. squared residuals provided additional support for their significance (see Additional file 1). Seven data sets, however, showed local minima that could not be distinguished from the global minima on these bases. In sum, we were able to identify convincing global minima in 31 of the 39 evaluable data sets (79 %) shown in Table 1.
Data quality was further evaluated by performing "jackknife" resampling studies on each of the patient data sets in Table 1 [17]. Optimal parameter values were obtained for each of the n subsets for each data set via the same process  used to analyze each complete data set (at the SD value of the complete data set's global minimum). We found CV values for calculated RD and LL parameters to be under 0.5 for over 75 % of these cases ( Table 1). Quantification of platelet lifespan was more difficult, with only 52 % of our cases showing CV values for the LS parameter of under 0.5. That is expected, however, because LS value estimates showed a larger variance in cases where random destruction predominated (Fig. 4b).
Predicting the effects of reduced platelet production and increased random destruction As a guide to interpreting the parameter values in Table 1, we used the model to predict how a normal platelet population's consumption parameter values might shift in response to A) impaired production, B) increased consumption, or C) increased consumption in association with a homeostatic increase in platelet production. Our assumptions were: i) The optimal parameter values (RD 0 , LS 0 , and SD 0 ) and the associated platelet production rate (PR 0 ) obtained for the three patients in the study whose platelet counts transiently normalized in response to prednisone (patients 3, 33, and 38, Table 1) are representative of normal. ii) RD is comprised of two component processes: Hemostatic RD (HRD) and non-hemostatic RD (NHRD) (i.e. RD = HRD + NHRD). Substantial hepatic NHRD is a well characterized phenomenon [18]. iii) The absolute HRD value at a normal platelet count (aHRD 0 ) makes up a given normal fraction ("f") of absolute RD (i.e. aHRD 0 /RD 0 = f ). We do not know the normal value of f. iv) aHRD 0 is maintained, as platelet count declines, via an increase in HRD and a resultant increase in RD, as suggested by earlier studies of platelet turnover [8].  Table 1) predictive of the consumption curves shown. b Relationship between RD value and variance of LS. Values are from Table 1 v) LS is not affected by reduced platelet production. Studies of the genetic basis of platelet lifespan support this assumption [3].
The effect of reduced platelet production (PR) on RD and platelet count was modeled as shown in Fig. 6. The process begins (step A) with the optimal (baseline) parameter values for pooled data from the three patients who transiently normalized their platelet counts (Table 1), using an initial "f" value of 1.0. From this set, a "target" reduced platelet production rate (PR 1 ) is generated, corresponding to 90 % of PR 0 . Using that value, the model generates the expected (reduced) aHRD value (aHRD 1 ). We then (step B) incrementally increase HRD until the model generated value of aHRD (aHRD i ) = aHRD 0 . The associated RD i value (=HRD i + NHRD 0 ) and platelet count values are those predicted to occur at PR 1 . Finally (step C), we repeat steps A and B with a series of reduced platelet production rates (PR i ). This generates predicted HRD and platelet count values for each PR i value. We then repeated this analysis at f values of 0.5 and 0.2.
We note that for our baseline parameter values, aRD 0 (RD × platelet count) is equal to 45 % of PR 0 . We make no quantitative predictions for the effect of reducing PR Fig. 5 Data quality evaluation. Schematic showing criteria used to remove non-evaluable data sets from evaluation. The first data point obtained for the patient lacking an initial data point was obtained at 6 h after infusion of labeled platelets Fig. 6 Modeling the effect of reduced Production rate on HRD and Platelet count. Schematic of the interpolation process described in the text below aRD 0 because the assumptions underlying the current numerical analysis model may not hold in that case. This is because the number of hemostatic targets is expected to increase below platelet counts at which hemostasis begins to be impaired. That in turn would invalidate the assumption of a constant absolute HRD rate, which is one of the bases for our iterative predictive method (Fig. 5). A model incorporating a dynamic hemostatic target population will be needed to predict platelet consumption rates in these circumstances.
To model the effect of increased random platelet consumption (RD), we generated a series of incrementally reduced target platelet counts (P i )(range: 90 % to 10 % of baseline), and to achieve each we incrementally increased RD from its baseline value until the model generated value of P (P m ) was equal to P i . To model the concurrent effects of increased RD and homeostatically increased PR, we used the same series of target platelet counts (P i ), and for each we increased PR in a manner proportional to the reduction in platelet count (to a maximum of twice the baseline PR value, a conservative theoretical starting point) before, again, empirically identifying RD i .
The results of these three modeling approaches are plotted with the values obtained for the patients in Fig. 7.

Optimal patient parameter values in comparison to modeled values
Surprisingly, only four patients in the study showed a platelet production rate that is even modestly increased (>50 %) in comparison to the presumed normals (Fig. 7c). The latter showed a mean platelet production rate (2.12 K/ul/h) comparable to the 1.7 K/ul/h rate estimated for normals in a previous study [8]. The finding of predominantly low to normal production rates in the thrombocytopenic cases (Fig. 7a) is corroborated by the distribution of random destruction rates (Fig. 7b), where  Table 1 are plotted for each of the patients in the study. Projected values for thrombocytopenias due to reduced production (at three values of HRD/RD), increased consumption with no homeostatic increase in production, and increased consumption with a compensatory increase in production rate, were interpolated as described in the text. a Population turnover rate, which at equilibrium equals platelet production rate, vs. platelet count. b Random destruction rate (RD) vs. platelet count. c RD and PR values from a and b. Error bars are standard deviations arrived at via jacknife resampling (see text) rates consistent with no increase in platelet production are seen for, again, all but a handful of the patients. A surprisingly large number of cases (at least 12 of 31) fall near the RD rates predicted to result solely from impaired platelet production. The predicted rates vary significantly, however, as a function of aHRD 0 /RD 0 ('f').
Because we don't know the normal value of 'f' , our ability to predict the increase in RD at low platelet counts is limited. Future studies in patients with thrombocytopenias due to impaired platelet production could resolve that problem.

Modeling of immature platelet fraction values
An ability to take up fluorescent marker dyes such as thiazole orange (a marker of "reticulated platelets", RP) or the proprietary dyes used in Sysmex hematology analyzers (marking the "immature platelet fraction", IPF) is thought to be characteristic of those platelets which have recently been released into the bloodstream. The age threshold (T) at which "young' platelets stop taking up these marker dyes is not known. Because the numerical analysis model generates a platelet age distribution for any given set of parameter values, it can be used both to estimate T and to predict the effect of altered production and consumption rates on the fraction of platelets of age less than T (i.e. the IPF).
Specifically, the normal range for the IPF is approximately 4.5 % (each clinical laboratory typically establishes its own range; this is the value in use at the Memphis VA Medical Center). Per the age distribution predicted by the model for our normalized controls (Fig. 8a), the youngest 4.5 % of platelets corresponds to those aged less than 4 h (i.e. T = 4 h). Application of that cutoff to the age distributions generated during modeling of the effects of altered production and/or consumption (Fig. 7), generates the predicted IPF and absolute IPF (aIPF) for thrombocytopenias induced by those mechanisms, as shown in Fig. 8b.
We note that this analysis depends on the assumption that all nascent platelets below a given age (T) take up the fluorescent markers used in the RP and IPF assays. Our measurements of mass turnover for mature and reticulated murine platelets suggest that this may not be the case [19]. Fig. 8 Age distribution histograms and IPF modeling. a Histograms were generated from the optimal kinetic parameter values shown in Table 1, using bin sizes of 10 h. b Projected IPF and absolute IPF (aIPF) were generated from the parameter values for the normalized platelet count group (Table 1), as described in the text