|Home | About | Journals | Submit | Contact Us | Français|
Patient-specific prediction of cellular response to multiple stimuli is central to evaluating clinical risk, disease progression, and response to therapy. We deployed Pairwise Agonist Scanning (PAS) to measure calcium signaling of human platelets in EDTA-treated plasma exposed to 6 different agonists (at 0.1, 1, and 10×EC50) used individually or in 135 pairwise combinations. With 154 traces, we trained a neural network (NN) model to accurately predict the entire 6-dimensional response to ADP, convulxin, U46619, SFLLRN, AYPGKF, and PGE2. The NN successfully predicted calcium responses to sequential agonist additions, all ternary combinations of [ADP]+[convulxin]+[SFLLRN] (R=0.88), and 45 different combinations of 4 to 6 agonists (R=0.88). Furthermore, PAS provided 135 pairwise synergy values that allowed a unique phenotypic scoring and differentiation of 10 donors. Training of NNs with pairs of stimuli across the dose-response regime represents a highly efficient approach to predict integration of multiple, complex signals in a patient-specific disease milieu.
Platelet activation is central to the thrombotic risks and events surrounding 1.74 million heart attacks and strokes, 1.115 million angiograms, and 0.652 million stent placements in the United States each year1. During clotting, platelets experience diverse signaling cues simultaneously. Collagen activates glycoprotein VI (GPVI)-dependent tyrosine kinase signaling, ADP is released from dense granules to activate the G-protein coupled receptors (GPCRs) P2Y1 and P2Y12, while thromboxane A2 (TxA2) is synthesized by platelet cyclooxygenase 1 (COX1) and binds TP receptors. Tissue factor at the damaged vasculature leads to the production of thrombin which cleaves the protease activated receptors PAR1 and PAR4. These activating signals occur in the context of inhibitory signals from endothelial nitric oxide and prostacyclin. While all these events occur simultaneously in vivo and platelet signaling varies spatially and temporally in growing thrombi2, few experimental or computational tools are available for building a global understanding of how a cell integrates multiple stimuli present at varying levels.
We developed a high throughput platform that measures the human platelet response to all pairwise combinations of six major agonists. Agonists tested were: convulxin (CVX; GPVI activator), ADP, the thromboxane analog U46619, PAR1 agonist peptide (SFLLRN), PAR4 agonist peptide (AYPGKF), and PGE2 (IP receptor activator). The method yields high-resolution dynamic measurements of intracellular calcium (Ca2+) which is the convergent node of platelet signaling pathways (Fig. 1a) and is critical to granule release, exposure of phosphatidylserine, actin polymerization, shape change, and integrin activation3. To determine appropriate dynamic ranges for the 6 agonists (Fig. 1a), each compound was first tested individually to determine each dose response relationship (Supplementary Fig. 1). The inhibitory response of PGE2 was studied by concomitantly stimulating the platelet with 60 μM SFLLRN. Then pairwise “inputs” were applied in all combinations to dye-loaded platelet rich plasma (Fig. 1b). The calcium traces allowed training of a 2-layer NN model to predict the dynamic response to all combinations of agonists used at any concentration (Fig. 1c).
Endogenously released ADP or TxA2 had no effect on the Ca2+ signal since 2 units/ml apyrase and/or 15μM indomethacin had no effect on individual responses (Supplementary Fig. 2 and Supplementary Tables 1 and 2). All experiments were conducted in 5 mM EDTA to chelate extracellular calcium. The removal of external calcium does not affect the ability of the receptors studied to signal, since no appreciable difference in EC50s were noted with or without external calcium (Supplementary Fig. 1a,b). Although such an experimental design does not capture the contribution of store operated calcium entry, it offers several operational advantages by (i) lowering background fluorescence without extensive platelet washing, (ii) preventing thrombin production, (iii) inhibiting granule release4, 5 as well as TxA2 formation6, and (iv) inhibiting integrin mediated signaling downstream of Ca2+ release7. Thus, the resulting traces of Ca2+ are directly dependent only on receptor mediated release from intracellular stores.
We tested all 135 pairwise combinations of low (0.1×EC50), moderate (1×EC50), and high (10×EC50) agonists concentrations (Fig. 2a). The pairwise agonist synergy score (Sij) is the scaled difference between the integrated transient (area under the curve) for the combined response and the integrated area for the individual responses (Fig. 2b) (Sij>0, synergism; Sij=0, additivity; Sij<0, antagonism). Neural networks are remarkable in learning patterns of inputs and predicting outputs by optimizing intermediate connection weights, akin to a platelet’s ability to respond to multiple thrombotic signals through coupled biochemical reactions. Motivated by the notion that a living cell is essentially a neural network whose connection weights have been selectively adjusted during evolution8, we took a “top-down” approach9 to model platelet signaling. A NN model was trained on these 154 time-course traces (135 pairwise responses, 18 single agonist responses, 1 null control response) and captured both the time-course behavior (R=0.968 for correlation between time points) and the pairwise agonist synergy (R=0.884) for correlation between Sij scores (Supplementary Fig. 3) with excellent accuracy (Fig. 2a,b).
Then, we measured all 64 ternary combinations of the agonists ADP, SFLLRN, and CVX at 0, 0.1, 1, and 10×EC50 concentrations (Fig. 3a). A CVX response requires GPVI multimerization10 and is characterized by a slow rise to a large peak signal followed by a slow decline. Gq–coupled responses (ADP or SFLLRN) produce rapid bursts that are quickly brought down to baseline. Increasing CVX for a fixed ADP level resulted in a steady increase in Ca2+ on longer timescales. In contrast, increasing ADP for a fixed CVX level bolstered early Ca2+ release. A moderate dose of both ADP and CVX (for 0 and low SFLLRN) produced a response that almost instantaneously plateaued at a steady level above baseline. The NN model, trained exclusively on the pairwise interaction data set (Fig. 2a), successfully predicted experimentally observed trends in the ternary interactions among the three agonists studied. Both the time-course behavior (R=0.844) and ternary agonist synergy scores (R=0.881) (Supplementary Fig. 4) were accurately reproduced for the 27 unique ternary conditions in this experiment that were not present in the training set.
To fully test and utilize the predictive power of the NN, we made in silico time-course and synergy predictions for the complete 6-dimensional agonist space consisting of 4077 unique agonist combinations of 2 to 6 agonists at 0.1, 1, or 10×EC50 concentrations. The full distribution of synergy predictions for all 4077 agonist combinations is shown as a vertical heat map in Fig. 3c. Based on these predictions (Supplementary Fig. 5), we selected 45 combinations of 4, 5 or 6 agonists that displayed a range of predicted synergy scores from synergy to strong antagonism and tested them experimentally in addition to no agonist and 18 single agonist controls (Fig. 3b). To prevent any bias in the selection we selected conditions that had maximal dissimilarity in the types and concentrations of agonists. We found strong agreement between both predicted and measured transient shapes (R=0.845) in Fig. 3b and Supplementary Fig. 6a, as well as between predicted and measured Sij scores (R=0.883, slope=1.08) (Fig. 3c). Artificially neglecting NN inputs typically but not always reduces predictive accuracy (Supplementary Fig. 6b), indicating that the NN does not merely rely on smaller subsets like dominant pairs. Conditions containing high levels of all agonists showed especially low synergy due to saturation of Ca2+ release. The highest synergy was observed for agonist combinations that contained high levels of the thromboxane analog U46619 with no PGE2 present (Fig. 3c, orange bar). Given that only 8 of 45 conditions had maximal U46619/PGE2 ratio, this ordering of the top 3 conditions was highly significant (p <0.004), considering there are 14,190 possible ways to order the first 3 conditions of which only 56 combinations would contain high U46619 and low PGE2. Thus, the NN model trained on pairwise data facilitated discovery of a high dimensional synergy that occurs at high U46619/PGE2 ratio (at low levels of ADP, SFLLRN and submaximal levels of AYPGKF) consistent with the known cardiovascular risks of COX2 inhibitors that prevent endothelial production of prostacyclin without affecting platelet production of thromboxane11. This points to a “high dimensional COX2 inhibition risk” of high concentrations of thromboxane, in the absence of PGI2, potentiating the effects of other agonists.
We also explored the effect of adding the agonists ADP, SFLLRN and CVX in various sequential combinations (Fig. 3d). Several interesting behaviors were accurately predicted by the NN model. Preceding a high dose of ADP with prior treatment with low dose of ADP (Fig. 3d, panel 12) desensitized the signal expected from a high dose of ADP (as in Fig. 3d, panel 1). This behavior has been observed previously12, 13 and is attributed to the internalization of P2Y1. The temporal sequence ADP-SFLLRN-CVX (Fig. 3d, panel 1) produced three distinct Ca2+ bursts, whereas the ADP response was completely abolished in the sequence SFLLRN-ADP-CVX (Fig. 3d, panel 3). This behavior points to mechanisms of cross-down regulation of ADP signaling by component(s) of the PAR1 cascade.
Prior addition of CVX abolished responsiveness to both ADP (Fig. 3d, panel 5) and SFLLRN (Fig. 3d, panel 6) again suggesting a mechanism where components of the GPVI signal are able to down regulate the Gq-coupled ADP or SFLLRN signal. Additions of any two of these agonists in combination followed by the third agonist confirm the observation that any mixture containing CVX down-regulates responsiveness to both ADP (Fig. 3d, panel 9) and SFLLRN (Fig. 3d, panel 8). CVX-mediated calcium mobilization events were unaffected by pretreatment with either ADP (Fig. 3d, panel 2) or SFLLRN (Fig. 3d, panels 4 and 3), or a binary combination of these agonists (Fig. 3d, panel 7). Activation of GPVI or thrombin receptors phosphorylates the ITIM domain of platelet PECAM14. ITIM phosphorylation inhibits response via phosphatases like SHP-215. Such mechanisms, or even agonist selective stores16, may explain the lack of ADP/SFLLRN response after prior CVX stimulation and the lack of ADP response after prior SFLLRN stimulus. Simulation traces containing CVX did not decay as observed experimentally after ~260 s (Fig. 3d, panels 2, 4–10). This limitation was expected because the NN was trained on measurements spanning only 260 s (Fig. 2a) and not the entire duration (up to 900 s). Importantly, the NN captured cross-talks of sequential additions despite being trained on purely synchronous interactions.
During the in vivo response to injury, the first adhering platelets experience strong GPVI signaling while subsequent depositing platelets contribute to the rising ADP and serotonin levels, followed by rising thromboxane levels and thrombin formation. While no single combination of agonists will replicate the dynamics of this in vivo situation, PAS-trained NN models offer the potential of using time-dependent agonist forcing functions as demonstrated in Fig. 3d (and Supplementary Fig. 7 for tests with thrombin compared to SFLLRN+AYPGKF). The ability to simulate thrombus formation under realistic hemodynamic conditions of flow17, 18 will require calculations of [Ca2+]i for platelets experiencing spatial and temporal exposures to stimuli. NN models are especially well suited for incorporation into multiscale models of patient-specific cell function in convective, reactive and dispersive flow fields19.
To investigate reproducibility and donor specificity, PAS was performed twice in a 2-week period for 10 healthy male donors (Fig. 4). The 135 conditions containing pairs of agonists in a single PAS experiment make up the synergy map for each donor experiment (Supplementary Fig. 8) and individual columns of the synergy matrix (Fig. 4). The standard errors in synergy scores across all 135 conditions were uncorrelated with the magnitude of synergy and are measures of the experimental uncertainty and day-today fluctuations in mean synergy values at these conditions. The uncertainties across all 135 conditions for a representative donor (Donor A) are shown in Supplementary Fig. 9 (mean uncertainty for this donor was ±0.0523 for Sij ranging from −1 to 1). The mean standard error in synergy scores for all 10 donors are tabulated in Supplementary Table 3 and ranged from ±0.0347 to ±0.0627. A simple hierarchical cluster tree was generated using the Euclidean distances between donor experiments. A total of 7 out of the 10 donor pair vectors (Donor pairs D, C, A, H, E, F and I) self clustered, demonstrating that in spite of intra donor variations, pronounced inter donor variations allow us to distinguish donors. This pattern of clustering was found to be highly significant (p<8×10−7) by randomizing observed donor synergies (Supplementary Fig. 10). The observed pattern of self clustering was platelet cell autonomous (and not related to donor plasma), since the PAS scans of an individual donor’s platelets with autologous or heterologous plasma self clustered (Supplementary Fig. 11). In general across all conditions and donors, highest probability of pairwise synergy was observed when moderate doses of both agonists were used. Low doses of both agonists produce additive responses, whereas high doses of both agonists skews synergy distributions towards antagonism. (Supplementary Fig. 12)
Donors separated into at least 2 major subgroups with the cluster of donor experiments D1, D2, J2, C1, C2, B1 and B2 characterized by relative lack of synergy (red) in comparison to other experiments. The cluster of experiments A1, A2, H1, H2, J1, E1, E2, F1, F2, G1, I1, I2 and G2 had marked synergy between moderate doses of SFLLRN and all doses of U46619/ADP, as well as marked synergy for moderate U46619 and high CVX. All donors showed some synergism between low and moderate doses of SFLLRN and U46619. Synergy was typically also found between AYPGKF and U46619. Moreover, synergistic/additive interactions were noted also between low and moderate doses of SFLLRN and AYPGKF. These results suggested a mechanism of synergy between thrombin and thromboxane. To test this, binary synergy maps of the physiological agonist thrombin and U46619 were constructed for donors A and E (Supplementary Fig. 13) over 7 doses spanning the active concentration ranges. To our knowledge, this is the first report of conserved synergy between thrombin and thromboxane mimetics.
Studying the combinatorial effects of pairs of agonists in low, moderate and high concentrations allowed a rapid, donor-specific phenotypic scan that was predictive of responses to 6 multiple agonists. Importantly, a single 384-well plate of data was sufficient to train a NN model (Fig. 2) capable of making accurate predictions of the global 6-dimensional agonist reaction space (Fig. 3), which is difficult to probe experimentally but fundamental to the processes of thrombosis. Synergies between platelet agonists are dependent not just on agonist pairs and doses, but also vary from donor to donor (Fig. 4). Previous studies have reported synergistic aggregation responses to combinations of multiple agonists20–22. Such unique patterns of synergisms could be used to distinguish donors and be correlated with certain risk factors.
The application of neural networks for predicting dynamic cellular signaling is beneficial because NNs are “dense” modeling structures that do not require detailed knowledge of the kinetic structure of a system. By comparison, an ordinary differential equation (ODE) model of ADP-stimulated calcium mobilization via P2Y1 required almost 80 reactions and over 100 kinetic parameters to describe just this one single pathway23. We estimate that an ODE model that describes the signaling mechanisms of the 6 agonists (Fig. 1a) in this study on a similar level of detail would require over 500 parameters, many of which are currently unavailable.
The PAS/NN approach works because unitary and binary interactions dominate and they are sampled across the full dose range of inputs. We expect the method to break down when ternary interactions in excess of summing binary interactions become strong. We show that the residual ternary synergy (Δ(ABC) = SABC-SAB-SBC-SAC) was ~ 0 in each of 27 responses of platelets to [CVX]+[ADP]+[SFLLRN] and was minimized in the NARX model training (Supplementary Fig. 14). We stimulated platelets with extracellular ligands for short times and explicitly blocked autocrinic feedbacks, which would otherwise expand the repertoire of participating factors. With these “pure inputs” (which are not broad intracellular perturbations like ionophores or phosphodiesterase inhibitors), the network architecture (Fig. 1a) necessitates rapid convergence to a common second messenger. Because of the speed of signaling and the known cell circuitry, it is difficult to envision realization of greater than second order synergies in such a situation. Thus observed ternary synergies appear to be well approximated by the linear superposition of the individual pairwise synergies. Since residual ternary interactions were undetectable, even higher ordered residuals are likely very small and potentially difficult to evolutionarily select. This exact result of unitary and binary interactions dominating with few detectable ternary interactions was reported for multiple cytokine stimulation of RAW264.7 cells24. Thus, the PAS approach is powerful for training NN models of signaling systems that accommodate multiple binary interactions where Δ(ABC) or higher ordered residuals are small.
In general, knowledge of pairwise interactions alone cannot be expected to predict response to multiple stimuli (>2) present simultaneously. However, certain characteristics of platelets and the conditions under which they were studied made such an approach feasible in this instance. These include: (i) the relative abundance of binary interactions in signaling systems with minimized ternary interactions (Supplementary Fig. 14)24, (ii) the efficient utilization of system history (Supplementary Fig. 15), (iii) the dense sampling of interactions across a full dose-response range, (iv) known intracellular wiring that rapidly converges on Ca2+, without the possibility of higher order effects from genetic regulation or other interactions on long time scales and (v) choice of well characterized extracellular ligands and careful design to avoid autocatalytic feedback.
Clinically, we anticipate PAS profiles to depend on such variables as ancestry, age, sex, pharmacology, and cardiovascular state, all of which requires further testing. Current measurements of platelet phenotype provide the coarsest measure of differences among healthy donors. For instance, platelet aggregometry for 359 subjects could classify them as “hypo and hyper” reactive to platelet agonists25; and flow cytometry of 26 individuals was able to classify them as high, medium or low responders. A combination of genetic variations may underlie differences in platelet phenotype, however linking genotype (1327 SNPs) to phenotype (flow cytometric measurement of P-Selectin exposure and fibrinogen binding) in 500 individuals26 demonstrated weak association probabilities.
Platelets are ideal “reduced” cellular systems for quantifying the effects of multiple signaling pathways because they are anucleate, easily obtained from donors, amenable to automated liquid handling, and biochemically well-characterized. In spite of the distinction from the in vivo platelet environment, the simplified in vitro conditions using EDTA were useful for the dissection of signaling pathways ‘decoupled’ from confounding autocrine effects of soluble mediators that are highly dependent on local platelet concentrations and prevailing transport processes. However, the operational advantages of using EDTA prevent direct prediction of important physiologic phenomena like granule release, integrin activation, and outside-in signaling.
Further expansion of the PAS set to include epinephrine, soluble CD40L, serotonin, and nitric oxide donors would essentially map a major portion of the entire platelet response space. The use of PAS with orthogonal pharmacological agents (indomethacin, P2Y12 inhibitors, selective PAR antagonists, quanylate cyclase or adenylate cyclase inhibitors) would allow further assessment of individual clinical risk or sensitivity/resistance to therapy. The PAS method demonstrates that sampling all dual orthogonal “axes” (every agonist pair) can successfully predict the dynamic responses and cross-talks of a higher dimensional system (6 agonists in this case).
Experimental methods and associated references are available in the online version for the paper at http://www.nature.com/naturebiotechnology/.
The authors thank Dr. Hongshe Li for suggesting the permutation test to evaluate the significance of donor clustering. This work was supported by NIH R01-HL-56621 (SLD), R33-HL-87317 (SLD, LFB), and T32-HG000046 (JEP).
Note: Supplemental information is available on the Nature Biotechnology website.
AUTHOR CONTRIBUTIONSM.S.C. designed and performed all experiments. J.E.P. constructed NN models of platelet activation. M.S.C. wrote the paper with contributions from all authors. L.F.B. advised on experimental conditions and S.L.D conceived the study.