A molecular signaling model of platelet phosphoinositide and calcium regulation during homeostasis and P2Y₁ activation

Platelets respond to endothelial injury or activating agonists by engaging a host of intracellular signaling events, including receptor activation, G protein signaling, second messenger generation, Ca²⁺ release, granule secretion, and cytoskeletal rearrangement.^1,2  Although the molecular details underlying individual activation events are continuously refined through focused studies, it has been difficult to develop a unified and integrated view of platelet metabolism because it involves a large number of simultaneously interacting molecular components. To this end, computational models help provide an integrated view of the many interacting components involved in cellular signaling.³  Successful applications of the “systems biology” approach have led to an improved understanding of cell surface receptors,^4,5  quantitative prediction of protease cascades in blood coagulation,^6-8  and the discovery of novel molecular interactions.⁹  When based on reliable datasets, accurate models have the potential to not only explain previously observed behaviors but also to make experimentally verifiable predictions about how cells process biologic signals.¹⁰ 

In many ways, the platelet provides an ideal test system for human systems biology. Platelets from normal donors or patients are readily available for in vitro diagnostic research and clinical monitoring. Isolated platelets, platelet-rich plasma (PRP), or whole blood are amenable to high throughput experiments to study signal transduction and associated coagulation protease cascades. In addition, numerous genetic mutations in humans¹¹  as well as knockout and transgenic mice^12,13  are known that predispose to bleeding phenotypes. From a modeling perspective, the platelet is particularly well suited for study because it lacks a nucleus, allowing one to avoid the challenges of describing whole-genome transcriptional regulation. Finally, the well-appreciated importance of platelets in mediating thrombosis and hemostasis, as well as their contribution to systemic disorders such as inflammation and cancer,¹⁴  places a high practical value on an accurate platelet model. Such a model would be useful both as a basic research tool for predicting normal activation behavior and as a strategy for the rational design of patient-specific pharmaceutical therapies.¹⁵ 

In this study, we describe a computational model of the human platelet that accurately predicts both resting and activated system behaviors. The model is based on 24 peer-reviewed studies spanning 3 decades of platelet research. Because of the inherent complexity in a model of this size, we first constructed a set of 4 distinct signaling “modules” to describe: (1) Ca²⁺ release and uptake, (2) phosphoinositide (PI) metabolism, (3) P2Y₁ G-protein signaling, and (4) protein kinase C (PKC) regulation of phospholipase C-β (PLC-β). Using fixed reaction equations and fixed kinetic parameters, each module was tuned to match a relevant set of experimental data by selecting allowable values for the resting concentrations of the module species. For the final analysis, the 4 modules were merged into a single kinetic model and fit to Ca²⁺ release data from adenosine diphosphate (ADP)–stimulated platelets. Dynamic traces of intracellular Ca²⁺ represent multiple steps in a complex signaling pathway from extracellular ligand to intracellular activation response, incorporating dynamic contributions from each module in the system. The fully integrated model was capable of reproducing population and single-cell Ca²⁺ release data as well as measurements of PI turnover and G protein activation.

Human blood was collected from healthy donors via venipuncture and anticoagulated with sodium citrate (9 parts blood to 1 part sodium citrate). PRP was incubated for 45 minutes with Fluo-4 NW (Invitrogen, Carlsbad, CA) per the manufacturer's instructions. PRP was diluted to a final concentration of 12% on a 384-well plate. A separate plate containing various concentrations of ADP was prepared on a PerkinElmer Janus (PerkinElmer Life and Analytical Sciences, Boston, MA). Intracellular Ca²⁺ concentration was quantified by measuring F(t)/F₀ using a Molecular Devices FlexStation (Sunnyvale, CA). Results are the average time-series measurements of 6 replicates. Phlebotomy was conducted in accordance with the Declaration of Helsinki and under University of Pennsylvania Institutional Review Board approval.

High-resolution electron micrographs of platelets stained with glucose-6-phosphatase obtained from Ebbeling et al¹⁶  were used to quantify the relative area of the dense tubular system (DTS) with respect to cytoplasmic area. Stained areas were delineated by using a threshold value for pixel intensity. To determine the sensitivity of the technique to the chosen threshold value, we calculated the DTS area using a maximally inclusive threshold (all stained regions and some of periphery included in delineation) and a minimally inclusive threshold value (delineation entirely within stained portion of image). The range of results obtained by both threshold values did not differ by more than a factor of 2.

Ordinary differential equations (ODEs) were used to describe reaction rates corresponding to published reaction mechanisms and kinetic parameters (Table 1). To efficiently estimate the dynamic range of initial conditions (ICs), the abundance of each estimated species or compartment size was estimated over a base-10 logarithmic scale. For reactions occurring at compartmental interfaces (eg, cytosol and plasma membrane), reaction rates were adjusted by scaling the concentrations of reactants to the bulk compartment as described by Kholodenko et al.⁵  Within each of the 5 compartments, species were assumed to be well mixed. Although we do not discount the importance of spatial gradients in Ca²⁺ signaling,¹⁷  previous studies have shown that changes in [Ca²⁺]_i within platelets are effectively instantaneous.¹⁸  Simulations were performed using the SBToolbox for MATLAB (Mathworks, Natick, MA).¹⁹  Numerical integration was performed with an absolute tolerance of 10⁻⁴⁵ and a relative tolerance of 10⁻⁹.

The model topology was defined by a set of reaction equations, rate laws, and kinetic constants (Table 1). The concentrations of all enzymes and metabolites in the resting platelet at homeostasis comprise the IC space of the model. The topology of the Ca²⁺ signaling network and almost all other kinetic parameters are known from the literature (Table 1) except for those describing P2Y₁ activation, which were calibrated from measurements of P2Y₁-reconstituted proteoliposomes (“Receptor activation”). A subset of IC space for each module was established by use of homeostasis constraints and experimental data from platelets. The homeostasis constraint requires that a resting platelet remains resting such that the IC is also a steady-state solution of a given set of ODEs. We performed a dense sampling of unknown ICs over fixed reaction topologies to identify combinations of IC values that were consistent with known steady-state concentrations.

As described by Heemskerk et al,¹⁸  strong release events were identified by taking the derivative of the Ca²⁺ flux traces and setting a noise threshold above which peaks were considered to be significant (Figure S3, available on the Blood website; see the Supplemental Materials link at the top of the online article). Video-imaged platelets were found to be quite variable both in peak interval times (4-40 seconds) and peak amplitudes (20-300 nM). We also observed significant variability in the platelet model, with peak amplitudes ranging from 100 to 300 nM. We were unable to estimate a single dominant frequency. using Fourier transform analysis of the simulated Ca²⁺ traces because of the substantial noise and irregularity of spiking. Half-maximal agonist dose (1 μM) was used to generate a frequency distribution for desensitized platelets treated with antagonists prostaglandin I₂, apyrase, and aspirin before 20 μM of ADP.

The full model (Figure 1; Table 1) comprises 4 interlinked kinetic modules (Figure 2A-D). The Ca²⁺ module (Figure 2A) spans 5 compartments and functions to maintain a low intracellular Ca²⁺ concentration ([Ca²⁺]_i) by pumping ions across the PM into the extracellular space, or across the membrane of the DTS into platelet stores using a sarcoplasmic/endoplasmic reticulum Ca²⁺ ATPase (SERCA). Inositol trisphosphate receptor (IP₃R) channels release Ca²⁺ from the DTS and are regulated by inositol 1,4,5-trisphosphate (IP₃) and Ca²⁺_i. In the PI module (Figure 2B), membrane-derived signaling intermediates, such as IP₃, are continuously recycled between the plasma membrane (PM) and cytosol by a series of phosphorylation, dephosphorylation, synthesis, and hydrolysis reactions (Table 1). Among its hydrolysis products, PLC-β generates PM-bound diacylglycerol (DAG), which, along with Ca²⁺_i, modulates the activity of PKC. PKC dampens G protein-coupled receptor (GPCR)–mediated signaling by phosphorylating PLC-β, reducing its hydrolytic activity. In a module for receptor activation (Figure 2D), free ADP in the extracellular compartment binds to the GPCR P2Y₁, causing a subsequent rise in activated G protein (G_q-GTP). G_q-GTP has a basal rate of hydrolysis that is accelerated when it is bound to PLC-β. The transient complex PLC-β-G_q-GTP (PLC-β*) is the hydrolytically active form of the enzyme and regulates transmembrane signaling in the PI module.

Clearly, the model does not include all of the reactions that govern platelet homeostasis and activation. Such a list is both unavailable and intractable to deploy until a core set of modules can be validated, which is the focus of this study. We have focused on the subset of reactions that govern ADP-mediated phosphoinositide signaling and mobilization of intracellular Ca²⁺ (Table 1). The P2Y₁ receptor is coupled primarily to the G_q family of G proteins,^20,21  which activate the beta-isoform of phospholipase C.^1,20  Previous work has shown that this pathway accounts for more than 90% of the ADP-mediated Ca²⁺ signaling in platelets.²²  To eliminate Ca²⁺ influx and the need to model store-operated Ca²⁺ entry, we used experimental data from studies in which extracellular Ca²⁺ was removed by ethylenediaminetetraacetic acid (EDTA) or another chelating agent. EDTA was also used in our own experiments. ADP releases Ca²⁺ only from the DTS²³ ; thus, modeling the acidic store is not required at this point.

The Ca²⁺ module addresses a fundamental question of what resting level of IP₃ is needed to balance the open probability of IP₃R channels with SERCA in a platelet of a given DTS volume and [Ca²⁺]_dts. Our strategy for modeling platelet Ca²⁺ signaling can be summarized in 2 steps: First, we considered the kinetic properties of the IP₃R and SERCA,^24-26  the resting Ca²⁺_i concentration,²  and the volume of the platelet²⁷  to be accurate as reported. These values were held fixed in the module. Second, we generated more than 10⁹ combinations of the unknown values in the module (number of IP₃R/platelet, SERCA pumps/platelet, [Ca²⁺]_dts, volume of the DTS) to find a set of “configurations” that reached an equilibrium state with the known resting [Ca²⁺]_i (“Module calibration and analysis”). Each configuration may be thought of as a unique platelet model with a characteristic compartmental structure and molecular makeup (eg, different numbers of SERCAs, IP₃Rs, etc). Note that more than one module configuration can produce the same resting [Ca²⁺]_i.

Existing kinetic, electrochemical, and physiologic data were incorporated into the module as follows: Resting platelets maintain a [Ca²⁺]_i between 40 and 100 nM^2,18  while storing concentrated Ca²⁺ ([Ca²⁺]_dts) in the DTS.²  Pumping of Ca²⁺ by SERCA isoforms (Figure 2E) was modeled according to a kinetic study of the type 3 SERCA isoform,²⁴  which is abundant in human platelets.²⁸  For the IP₃R,²⁵  each subunit of the type 2 receptor exists in one of 6 states: native, open, shut, active, or 2 inactive states (Figure 2F), where state transitions depend on [IP₃] and [Ca²⁺]_i. The channel open probability (P_o) is determined by the number of IP₃R tetramers with all subunits in either open or active conformations (Table 1). We used the Nernst equation²⁹  to relate P_o to the release of Ca²⁺ from the DTS:

Here, N is the total number of channels per platelet, γ is the single-channel conductance (10 pS) of the platelet IP₃R,²⁶  and e is the number of elementary charges (z) per second per Ampere (6.24 × 10¹⁸). The cytosolic volume of the platelet is approximately 6 fL.²⁷ 

The unknown quantities in the module were the number of IP₃R and SERCA per platelet, [Ca²⁺]_dts, and the volume of the DTS. The range of values sampled for each quantity is shown in Figure 3A. A dense sampling (n = 10⁹) of these quantities generated 100 000 unique configurations of the Ca²⁺ module (with [Ca²⁺]_i = 100 ±10 nM) that were further divided into 3 groups (low-, mild-, and high-response) based on the estimated [Ca²⁺]_dts, an indicator of the configuration's response to an increase in [IP₃] (Figure 3B). To examine any molecular or structural differences among low-, mild-, and high-response groups, we compared the levels of SERCA, IP₃R, IP₃, and relative DTS volume that were estimated for each group. These values are presented as probability distributions (Figure 3C), where regions of high density reflect frequently occurring values for the estimated concentrations and compartment sizes. For example, low response configurations tended to have approximately 1000 IP₃R channels per platelet as indicated by red regions (Figure 3B). We note that, of 10⁹ sampled configurations, only 0.005% satisfied the dual constraints of steady-state Ca²⁺ homeostasis and IP₃ responsiveness (mild- or high-response configurations). This observation indicated that the kinetic properties of Ca²⁺-regulating enzymes such as IP₃R and SERCA alone place strong constraints on the physical size and molecular makeup of the platelet.

Configurations with a negligible gradient ([Ca²⁺]_dts ∼ 100 nM) lacked IP₃ responsiveness and had relatively low SERCA levels (67% had < 5000 copies) with a median estimate of 1600 IP₃R channels/platelet. By contrast, mild-response configurations ([Ca²⁺]_dts < 10 μM) responded to elevated [IP₃] with a transient rise in [Ca²⁺]_i and were characterized by high SERCA levels (∼ 10⁶/cell) and higher IP₃R abundance (median, 2000 channels/cell). High-response configurations ([Ca²⁺]_dts > 10 μM) clustered near 10⁶ SERCA pumps but showed a broad distribution of IP₃R abundance. The calculated IP₃R/SERCA ratio for low-response configurations had a loosely defined distribution, whereas the mild- and high-response configurations favored 10³ to 10⁵ more pumps than channels (1:6900 for mild-response and 1:5200 for high-response). Thus, the model predicted a very low IP₃R/SERCA ratio for functioning platelets, a hypothesis supported by noting that a single SERCA3b pump operating at a membrane with a 1000-fold Ca²⁺ gradient transports approximately 0.4 Ca²⁺ ions/second,²⁴  whereas a single type 2 IP₃R channel in the same membrane at 37° conducts approximately 3000 Ca²⁺ ions/second, assuming a steady-state P_o of 0.0005.²⁵ 

The low-response configurations harbored high resting concentrations of IP₃, inconsistent with a resting platelet.³⁰  The median IP₃ count for mild-response configurations was 750 molecules/cell, similar to the measured value (∼ 1200 molecules/platelet or ∼ 200 nM)³¹  in resting platelets. More than 80% of high-response configurations harbored less than 1000 IP₃ molecules/cell in resting platelets. The relative volume of the DTS compartment to cytosol, as predicted from the estimation procedure, was quite constrained at 0.5% to 5% (median, 2%) of the nongranular intracellular volume. To compare this result to a direct physical observation, we quantified the relative area of the platelet DTS as revealed by glucose-6-phosphatase staining.¹⁶  This calculation gave a DTS/cytosol fraction of 4.3%. Assuming the stained region marks the true DTS border, this estimate is accurate within a factor of 2 (“Platelet image analysis”). Thus, estimates obtained from independent kinetic calculations or image analysis of stained platelets gave similar values for the DTS compartment volume fraction.

PIs are continuously interconverted in platelets, even under resting conditions.^2,30,32  To capture this behavior, we used a set of synthesis and degradation reactions that continuously recycles PIs in the absence of a stimulating dose of agonist (Figure 2B; Table 1). Failure to account for PI resynthesis would have prevented prediction of steady-state concentrations. Using the same strategy for modeling platelet Ca²⁺ regulation (“Platelet Ca²⁺ balance”), we generated 10 000 unique PI module configurations that were consistent with the measured resting concentrations of PtdIns, PIP, PIP₂, DAG, and PA.^32-34  Thus, each configuration represents a potential molecular arrangement of PIs in the platelet that is consistent with published kinetic rates and resting measurements for key membrane phospholipids.

To examine the transient behavior of the PI module under activating conditions, we set the concentration of PLC-β* in the module to 1 μM (∼ 100-fold above basal levels). Figure 4A compares the measured and simulated changes in the platelet PIs after treatment with thrombin or elevated PLC-β*, respectively (both thrombin and ADP signal through G_q-mediated stimulation of PLC-β¹ ). The simulated time course for PtdIns degradation was consistent with experiment throughout the 90-second interval after activation. Simulated levels of PIP and PIP₂ were consistent with experiments³²  through the first 25 seconds (Figure 4A). After this time, PIP and PIP₂ continued to be degraded in the simulation but were sustained, or slightly elevated, in thrombin-stimulated platelets. A recent study in murine megakaryocytes³⁵  has confirmed that resynthesis of PIP₂ is necessary to replenish basal levels after agonist stimulation, a feature that is included in the model. However, because the rate of hydrolysis exceeds the rate of PI synthesis, we observed a monotonic decrease in PIs when the level of PLC-β* was held fixed. To address this discrepancy between model and experiment, we introduced a negative-feedback module (Figure 2C) in which the activity of PLC-β* is attenuated by activated PKC (PKC*) through phosphorylation of the phospholipase.

As noted in the analysis of PI metabolism, the degradation of PIs in platelets after PLC-β activation does not decrease monotonically.^2,32  Furthermore, the sharp rise in [IP₃] and [Ca²⁺_i] after GPCR stimulation returns to near-basal levels within 30 to 60 seconds.³⁶  This negative-feedback mechanism was modeled by allowing direct inactivation of PLC-β by PKC through a phosphorylation reaction, which has been observed in several cell types.^37,38  In the PKC module (Figure 2C), activation of the kinase requires association with Ca²⁺_i and DAG and translocation to the PM.³⁹  Activated PKC (PKC*) phosphorylates PLC-β, rendering it unable to bind G_q-GTP. Clearly, this represents only one of several potential mechanisms for regulation of G_q-dependent PI hydrolysis in platelets. Other mechanisms include receptor internalization,⁴⁰  deactivation of G proteins,¹  and accelerated dephosphorylation of IP₃.⁴¹  Although these mechanisms will provide useful contributions to future iterations of the platelet model, it is interesting that the use of PLC-β regulation was sufficient to reproduce the observed attenuation in PI signaling after agonist stimulation (Figure 5C,D). Specifically, the gradual accumulation of Ca²⁺ and DAG during the first 10 to 15 seconds provides the appropriate time delay for shutting off PIP₂ hydrolysis.

P2Y₁ is a GPCR expressed on the surface of human platelets (∼ 150 copies/platelet⁴² ) that is essential for ADP-induced platelet shape change and aggregation.²⁰  Although other ADP receptors are involved in platelet Ca²⁺ signaling (eg, P2 × 1 and P2Y₁₂), P2Y₁ contributes more than 90% of the Ca²⁺ signal.²²  In our representation of P2Y₁ activation (Figure 2D,H), binding of extracellular ADP to P2Y₁ leads to activation of G_q through GDP/GTP exchange reactions. G_q-GTP is a substrate for GTPase activating proteins (GAPs), such as PLC-β and RGS4, which can accelerate G_q-GTP hydrolysis more than 1000-fold.⁴³  Bornheimer et al showed that the combination of G proteins, active GPCRs, and GAPs can form kinetic “ternary” modules with distinct signaling patterns.⁴  Thus, we used results from a published in vitro assay of P2Y₁ activity to construct a ternary model of the P2Y₁ receptor. In that assay, Waldo and Harden combined purified human P2Y₁, ADP, G_α and G_βγ subunits, PLC-β or RGS4, and radiolabeled GTP in large phospholipid vesicles and measured the agonist binding and G_q-stimulating activities of the receptor.²¹  GTPase activity was measured with and without agonist (Figure 4B). ADP dose-response was measured with and without GAP protein (Figure 4C). Using identical simulated conditions, we obtained the kinetic parameters in the P2Y₁ module from these time-course measurements. The kinetic rate constants, given in Table 1, accurately reproduced both sets of measurements (Figure 4B,C).

The 4 signaling modules (Figure 2A-D) were merged into a single kinetic model using a bootstrap method (Figure S1). With a fixed reaction network (Figures 1,2) and fixed kinetic rate constants (Table 1), ICs were fit to dose-response [Ca²⁺]_i time-course data from ADP-stimulated platelets (Figure 5A,B). In the presence of basal ADP levels similar to human venous plasma,⁴⁴  the model maintained a resting [Ca²⁺]_i of 75 nM (Figure 5B). At higher ADP concentrations, [Ca²⁺]_i reached peak levels approximately 20 seconds after addition of agonist with a decay constant of approximately 1 minute. We found the synchronous [Ca²⁺]_i peak-response at all agonist doses to be a complex and unique feature in the data, particularly well suited for testing model validity. Achieving this characteristic shape required a transient (nonmonotonic) rise in [IP₃]. In ADP- and thrombin-stimulated platelets,^36,45  [IP₃] increases rapidly after agonist addition, peaks around 15 seconds, and then decreases to near-basal levels. The model predicted this behavior accurately (Figure 5C). We observed that PLC-β*, which comprised approximately 0.2% of the total PLC-β pool in the model, strongly controlled this rise and fall of [IP₃] (Figure 5C,D). Phosphorylation of inactive PLC-β by PKC* (Figure 5E) was sufficient to dampen G_q-mediated PI hydrolysis. The simplicity of this negative-feedback model suggests a mechanism whereby the phosphorylation of inactive PLC-β, present at much higher abundance than PLC-β*, fine-tunes the hydrolytic activity of the enzyme by sequestering it from G_q-GTP. Addition of ADP caused a sharp burst in [G_q-GTP] that saturated within 5 to 10 seconds for all doses (Figure 5F). This trend was consistent with the observed time-scale for G protein activation.⁴⁶  Maximal [Ca²⁺]_i and [IP₃] responses to increasing ADP were predicted quantitatively (Figure 5G,H).

Deterministic simulations, such as those in Figure 5B through F, resemble data obtained from suspensions of platelets (Figure 5A), in which a million or more cells may be assayed simultaneously. These conditions often produce average or “smooth” responses that are typical of cell populations. However, when platelets are monitored individually, the addition of agonist can produce a series of sharp, asynchronous Ca²⁺ spikes.^18,47  To test the single-platelet response in the model, we used the Stochastic Simulation Algorithm⁴⁸  to simulate Ca²⁺ release in a single platelet during rest or activation by ADP. Using the estimated sizes of the platelet compartments, species concentrations were converted to integer values for use in the stochastic simulation. Note that any species present at 1 nM in a 6-fL platelet constitutes less than 4 molecules, a level that can give rise to stochastic fluctuations.

Dramatic, asynchronous [Ca²⁺]_i spikes increased in frequency and amplitude after activation with ADP (Figure 6A,B). Before agonist stimulation, the resting platelet model showed occasional Ca²⁺ spikes with relatively small amplitude. Similar behavior was noted in a third of recorded platelets,⁴⁷  which spiked on average 1.6 (± 0.2) minutes without ADP stimulation. The onset of Ca²⁺ release after agonist addition was rapid in the model (∼ 5 seconds) and persisted for several minutes, also consistent with experimental results.^18,47  Key drivers in this asynchronous spiking were the channel open probability, the low copy number of IP₃R channels, the size of the cytosolic compartment (8 fL in the model), and the large Ca²⁺ gradient at the DTS membrane. To determine the frequency distribution of Ca²⁺ spiking, we followed exactly the empirical strategy used by Heemskerk et al¹⁸  to distinguish individual Ca²⁺ release events in a detailed study of video-imaged platelets (Figure S3). There was a striking resemblance between simulation and experiment for the peak interval distribution (Figure 6C), with the most common intervals occurring at 6 to 8 seconds or 11 to 13 seconds. Because platelet volume was artificially increased 100-fold in the stochastic simulation (Figure 6D), the Ca²⁺ spiking was lost and typical, smooth ADP responses (as in Figure 5B) for populations emerged. Thus, the asynchronous spiking appeared to be a result of the fundamentally stochastic nature of signal transduction in cells as small as human platelets.

We deployed available kinetic data, electrochemical calculations, measurements of cell size, and modular organization of signaling function to estimate the physical structure and dynamical properties of ADP-stimulated Ca²⁺ release in human platelets. For the first time, a single model has accounted for several fundamental attributes of known platelet behavior: (1) In the presence of basal ADP levels, the model remained at steady state and accurately predicted the resting concentrations of Ca²⁺_i, IP₃, DAG, PA, PI, PIP, and PIP₂. Under these conditions, the model underwent slow turnover of signaling molecules (eg, IP₃, PLC*, and G_q-GTP), yet the rates of synthesis and degradation of these metabolites were exactly balanced. At the same time, Ca²⁺ was continuously leaked from the DTS and actively pumped back into stores, maintaining a resting platelet [Ca²⁺]_i near 100 nM. This “active” homeostatic state is representative of the situation in vivo because inhibition of SERCA pumps causes a gradual depletion of Ca²⁺ stores.²  (2) Increased [ADP] caused a transient rise in [Ca²⁺]_i and [IP₃] with the correct dose-response and time-course profiles. (3) Stochastic simulation of the platelet model generated noisy Ca²⁺ release and reuptake behavior with a nearly instantaneous onset. (4) Calculation of peak-to-peak interval times revealed a distribution of intervals favoring 6- to 8- and 11- to 13-second gaps. These results were strikingly similar to what was observed in video-imaged platelets and demonstrate that both averaged and stochastic behaviors may be accurately captured by the same molecular model. Importantly, the model fulfilled the dual requirement of maintaining homeostasis under resting conditions while remaining capable of generating an activation response in the presence of sufficiently high agonist concentration. This is a critical quality to consider when modeling platelets, which must be sensitive enough to respond to injury without activating prematurely.^1,2 

Our analysis also produced several novel predictions: (1) The number of SERCA pumps was predicted to outnumber IP₃R channels (< 5000 IP₃R copies were predicted in most configurations). This was largely a consequence of the kinetic properties of the enzymes, which transport Ca²⁺ at significantly different rates. (2) The relative volume of the platelet DTS was predicted to be 2% to 4% of the nongranular intracellular volume, both by kinetic modeling and image analysis of a DTS-stained platelet. (3) Lower resting [IP₃] was associated with greater [Ca²⁺]_dts and more acute Ca²⁺ release. (4) Restoration of basal PI levels after PLC-β* hydrolysis was predicted to occur in part through a negative-feedback mechanism in which PKC* phosphorylates PLC-β and thereby diminishes its hydrolytic activity. (5) By artificially increasing the volume of the simulated platelet, we were able to quench the fluctuation of Ca²⁺ spikes, achieving a response more similar to that of megakaryocytes, in which regular oscillations in [Ca²⁺]_i are observed with fixed amplitude and frequency.⁴⁹  Because mature megakaryocytes are more than 1000 times larger than platelets,⁵⁰  these observations suggest that the asynchronous and stochastic Ca²⁺ response in platelets may simply be the consequence of a small intracellular volume.

The model also provided some insight into why thrombin is a more potent agonist than ADP. Both agonists act through G_q-coupled GPCRs,¹  yet there are about 10 times as many thrombin receptors as ADP receptors.^1,42  To determine whether the increased potency of thrombin may be the result of greater receptor surface expression, we increased the number of P2Y₁ receptors in the model and examined the ADP response. A 10-fold increase in [P2Y₁] caused a 2-fold higher peak response and more sustained elevation of [Ca²⁺]_i (Figure S4), similar to the observed difference between thrombin- and ADP-stimulated platelets. Although a more thorough analysis will involve examining the coupling efficiency to G_q and the duration of the receptor active state, these calculations imply that the low copy number of P2Y₁ may be a limiting factor for ADP signaling.

Our estimation procedure generated 10 individual platelet models (unique sets of ICs representing the resting concentrations for the 70 species in the model) that each satisfied homeostasis constraints and fit all available data (Figure S5). These models may be thought of as individual cells in a platelet population that vary in their exact molecular makeup but are each capable of generating the same activation response. For certain species such as Ca²⁺_i, DAG, and PIP₂, the similarity among the models can be easily explained because these quantities were held fixed (with a small variability term) to the values reported in the literature. Other similarities in the concentrations, as for GTP and PLC-β*, are more interesting because these values were not explicitly restricted by the estimation procedure. Rather, these uniformities reflect implicit constraints in the platelet model that emerged because of other species concentrations or signaling reactions that strongly control the production or consumption of these species. As such, these quantities are predicted to be rigid, or inflexible, “nodes” in the model⁸  and are subject to further investigation as a new hypothesis revealed by the analysis. By contrast, we observed several instances of variability in concentrations (node flexibility) among the 10 models. For example, the concentrations of PKC and G_q-GTP both varied over 3 orders of magnitude. Similarly, the values for P2Y₁, SERCA, PMCA, PLC-β, Ins, IP₂, and CDP-DAG varied appreciably. The large variability observed for many of the estimated phospholipids (eg, Ins, IP₂, and CDP-DAG) was the result of allowable node flexibility and/or insufficient experimental data to constrain the estimation of the resting concentrations. Thus, whereas the essential steady-state and activation properties were accurately captured by the module as a whole (Figure 4A), additional experimental data will be necessary to generate reliable estimates for individual phosphoinositides. Finally, the estimated number of IP₃R channels was generally low (< 2500 channels/cell) as expected from the analysis of the Ca²⁺ module. Higher [Ca²⁺]_dts appeared to be compensated by higher SERCA levels.

Through a computational analysis of the human platelet, we have demonstrated the successful integration of heterogeneous datasets and legacy knowledge to construct a predictive model of platelet homeostasis and activation. This assembly of molecular detail represents, for the first time, a substantial description of the platelet signaling apparatus and was capable of reproducing diverse experimental observations as well as providing specific, testable hypotheses regarding platelet enzymatic regulation, the physical structure of the cell, and the dynamics of intracellular calcium release. Specifically, this is the first work to provide a quantitative molecular explanation of the asynchronous calcium spiking observed in stimulated human platelets. By incorporating additional experimental data, future versions of the model can include mechanisms for store-operated Ca²⁺ influx mechanisms, the acidic Ca²⁺ store, integrin activation, granule secretion, phosphatidylserine exposure, as well as additional receptors for prostacyclin, nitric oxide, thromboxane, epinephrine, thrombin, fibrinogen, and collagen. As experimental investigation of platelet function continues, new data will be used to challenge, confirm, and expand the platelet signaling model presented here.

The authors thank Ravi Radhakrishnan for help with model reduction and estimation methods.

This work was supported by the National Institutes of Health (grants R01-HL-56 621 and R33-HL-87 317; S.L.D., L.F.B.).

National Institutes of Health

Contribution: S.L.D., L.F.B., and J.E.P. designed the research and prepared the manuscript; M.S.C. measured ADP-induced Ca²⁺ release; and J.E.P. performed model constructions and simulations.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Scott L. Diamond, 1020 Vagelos Research Laboratories, 3340 Smith Walk, Philadelphia, PA 19104; e-mail: sld@seas.upenn.edu.