We present a data-driven mathematical model of a key initiating step in platelet activation, a central process in the prevention of bleeding following Injury. In vascular disease, this process is activated inappropriately and causes thrombosis, heart attacks and stroke. The collagen receptor GPVI is the primary trigger for platelet activation at sites of injury. Understanding the complex molecular mechanisms initiated by this receptor is important for development of more effective antithrombotic medicines. In this work we developed a series of nonlinear ordinary differential equation models that are direct representations of biological hypotheses surrounding the initial steps in GPVI-stimulated signal transduction. At each stage model simulations were compared to our own quantitative, high-temporal experimental data that guides further experimental design, data collection and model refinement. Much is known about the linear forward reactions within platelet signalling pathways but knowledge of the roles of putative reverse reactions are poorly understood. An initial model, that includes a simple constitutively active phosphatase, was unable to explain experimental data. Model revisions, incorporating a complex pathway of interactions (and specifically the phosphatase TULA-2), provided a good description of the experimental data both based on observations of phosphorylation in samples from one donor and in those of a wider population. Our model was used to investigate the levels of proteins involved in regulating the pathway and the effect of low GPVI levels that have been associated with disease. Results indicate a clear separation in healthy and GPVI deficient states in respect of the signalling cascade dynamics associated with Syk tyrosine phosphorylation and activation. Our approach reveals the central importance of this negative feedback pathway that results in the temporal regulation of a specific class of protein tyrosine phosphatases in controlling the rate, and therefore extent, of GPVI-stimulated platelet activation.
Platelets, small anuclear cells, are quiescent in undamaged blood vessels. They respond to injury by activating, triggering blood to clot. Whilst platelet activation is essential to prevent excessive bleeding at sites of injury, their inappropriate activation, for example as a consequence of vascular disease, can lead to the formation of clots within the circulation, or thrombosis which triggers heart attacks and strokes [1]. Damage to blood vessels results in the exposure of extracellular matrix proteins, particularly collagens, which form structural components within the vessel wall. Collagens provide an initiating signal for platelet activation. Indirect interactions between platelets and collagen, in the arterial circulation, are mediated by von Willebrand factor (vWF) which binds to collagen and to GPIb on the platelet' s surface. This causes platelets to slow down and roll along the site of vessel damage allowing direct binding to platelet collagen receptors, including integrin α2β1, which largely supports adhesion, and glycoprotein (GP) VI which stimulates cell signalling and activation [2]. Platelet activation is marked by a dramatic change in platelet shape, the secretion of various prothrombotic factors and conformational change in the integrin αIIbβ3. These secreted factors initiate a second wave of signalling. Integrin αIIbβ3 binding to fibrinogen, which in itself stimulates signalling, supports platelet aggregate formation allowing the assembly of a platelet thrombus (or haemostatic plug) to stem the loss of blood [2,3]. GPVI is a member of the immunoglobulin family of receptors that shares aspects of its mechanisms of activation with immunoreceptors, including B and T cell antigen receptors [4]. These receptors possess or associate with transmembrane proteins that have cytoplasmic domains that contain immunoreceptor tyrosine-based activation motifs (ITAMs). The binding of a ligand leads to associated Src Family Kinases (SFKs) phosphorylating conserved tyrosine residues on the ITAM. In the case of GPVI, Syk (spleen tyrosine kinase), or the related protein ZAP70 in T cells, is recruited to the receptor complex through direct binding to the phosphorylated ITAM. Recruitment of Syk in turn provokes its phosphorylation [5], pivotal to downstream signalling events: in platelets this leads to activation, shape change and aggregation (Fig 1) [6]. The development of new drugs to suppress platelet function and thereby prevent thrombosis has been shown to be an effective strategy [7–9]; however, current anti-platelet therapies are ineffective in many patients and are associated with side effects. A more detailed understanding of the signalling pathways that lead to platelet activation is needed to develop more sophisticated, effective and safer anti-platelet therapies [10]. GPVI is a candidate anti-platelet drug target [11,12] and the central role of Syk in immune cell signalling also makes it an attractive therapeutic option, with implications for the treatment of allergy, auto-immune diseases and haematological malignancies [13]. In other cell types mathematical and computational models have aided in the understanding of signalling pathways and assisted further experimental work [14–16]. The formation of blood clots and thrombi has long been of mathematical interest but, probably due to the ready availability of data, mathematical modelling has focussed on the coagulation cascade [17–23], i. e. the network of interactions between blood proteins that results in the formation of fibrin polymers that stabilise a clot. These models show promise in translating an individuals' coagulant profile into a measure of thrombotic risk [24,25]. Even though platelets are an ideal cell type to model (they lack a nucleus thereby removing the complexities of gene transcription and translation) it is only recently that the platelet has become of mathematical interest. The main focus of recent developments have been the interactions between a platelet and the coagulation cascade [26,27] and the role of platelets in a thrombus shaped by blood flow [28–31], there being several excellent reviews [32–34]. Where platelets have been studied in more detail investigations focus on the functional response of platelets to agonists [35] and the mechanics of signalling pathways [36–39]. Dandekar and co-workers [36,37] focused on platelet signalling mechanisms and the effects of prostacyclin, a molecule secreted from endothelial cells that contributes to platelets remaining in a quiescent state in the undamaged circulation. Their studies were based on platelet-specific data and simulated and tested different pharmacological conditions and the contributions of various platelet receptors. A detailed mathematical model [38] that captures current knowledge of platelet signal transduction through the ADP receptor is able accurately to predict calcium flux and ADP dose-response. This work was used as a basis for a model of PAR-1 mediated platelet activation [39] that includes the critical step of linking inside-out to outside-in integrin signalling. These latter two studies focus on the events downstream of G-protein coupled receptors through to intracellular calcium mobilisation, a critical step in platelet activation that is common to all activation pathways. While GPVI activation is critical in health and disease no mathematical study has focussed on the platelet' s response to collagen. The pathway of reactions downstream of the GPVI receptor is considerably more complex than those stimulated by G-protein coupled receptors. We therefore focus our investigations on the early reactions, proximal to the GPVI receptor, that lead to Syk activity that is thought to be a fundamental switch that controls GPVI mediated signal transduction. Studies to dissect the component parts and associated processes within the GPVI signalling pathway have provided a picture of the series of steps that occur [40] but these studies often lack detailed kinetic analysis and molecular detail. For example, they may measure the level of tyrosine phosphorylation of a specific protein, rather than assessing the phosphorylation of specific amino acid residues, despite these individual phosphorylated residues being linked to specific functions. Furthermore, previous attempts to model signalling pathways in this and other cell types have often been based on qualitative data, which may make it difficult to detect differences between experimental data and model predictions. Here, these limitations in available data were resolved by the collection of high density quantitative experimental observations that capture the temporal increase and decrease in tyrosine phosphorylation of key proteins at specific amino acid residues in the early steps of activation of the GPVI signalling pathway. Existing knowledge of platelet signalling pathways tends to focus on the sequential forward reactions that lead to platelet activation. While tight negative regulation of the pathways appears crucial for appropriate and rapid platelet activation, little effort has focussed on understanding what the key regulators are that modulate these reactions and how they integrate to precisely control platelet activation. Stimulation of GPVI results in the initiation of a signalling pathway that is highly dependent on the activities of a range of protein kinases. A protein tyrosine kinase transfers a phosphate group to a protein. Such phosphorylation may alter the function of a given protein, for example through inducing a conformational change, creating docking sites, causing intracellular relocation, and modulation of enzymatic activity. Protein tyrosine phosphatases reverse these protein modifications since they are able to remove phosphate groups returning proteins to their original state. Historically protein tyrosine kinases have been the focus of more research than protein tyrosine phosphatases. Recent studies indicate that human platelets possess at least 18 different protein tyrosine phosphatases, numbering more than 52,000 copies in total per platelet [41]. While traditionally phosphatases have been thought to be promiscuous in activity, constitutively active and able to remove phosphate moieties from proteins, it is now evident that this view is too simplistic and that phosphatases like kinases are actively recruited and regulated within specific pathways [10,42,43]. The identities and roles of phosphatases within the GPVI signalling pathway are not well understood, although it is clear that some play a fundamental role in the control of platelet function following stimulation with collagen [44,45,46]. Biological knowledge of components and interactions of signalling pathways is often expressed by cartoons or network diagrams. These descriptive static models provide little information about how the many components evolve over time, making it difficult to predict system behaviour. Here we reinterpret network diagrams, representing current biological knowledge and hypotheses surrounding the early regulation of the GPVI signalling pathway, into a series of dynamic mathematical models that include various possible modes of regulation by phosphatases. This allows direct comparison between knowledge and hypotheses and high density quantitative, temporal experimental data generated specifically to inform the development of the models. Through a process of parameter estimation, model validation and comparison with experimental data we gain understanding of the temporal dynamics that control the initiation of the GPVI signalling pathway allowing us to explore how perturbations of different kinds (e. g. ligand availability, negative regulation processes, and expression of key regulatory proteins) effect the ability of GPVI to signal downstream.
In all models we treated GPVI, its associated Fc receptor γ-chain and SFKs as a single entity (see Fig 2). The ligand (CRP and indeed collagen) is a complex one being of indeterminant length, weight and number of binding sites for the GPVI receptor. While ligand mediated clustering is thought to be an important attribute the difficulty with obtaining experimental data covering these very early timepoints leads us to neglect these complications and assume an abundant ligand source. A ligand binds to this receptor bundle, which is subsequently phosphorylated, allowing for recruitment and activation of the cytosolic protein-tyrosine kinase Syk. Binding to the receptor leads to auto-phosphorylation of Syk, on tyrosine 525 (Y525), that allows the receptor to signal downstream [48–50]. The early reactions of ITAM phosphorylation and Syk binding are unlikely to be rapidly reversible. The binding affinity of Syk to the ITAM is high, occurring through two tyrosine residues [5] that protect the ITAM from dephosphorylation. Experimental data shows quantitative differences between protein copy numbers and equilibrium levels of Syk phosphorylation, pointing to regulatory activity. Syk phosphorylation is an obvious candidate for regulation, most likely through a phosphatase. Phosphatases have traditionally been viewed as promiscuous housekeeping molecules that act quickly on exposed phosphoproteins [43]. Our initial model (A) incorporated the continuous presence of a simple abundant generic phosphatase that is able to dephosphorylate Syk on Y525, thereby restricting the ability of Syk (and therefore of the GPVI receptor) to signal downstream. Full details of the model and the parameter values obtained are given in S1 Text. Model solutions describing Syk phosphorylation monotonically rise and settle to the steady state seen in experimental data (see Fig 4, top row). The steady state is solely determined by Syk' s rate of phosphorylation and its regulation, such that Bound, inactive Syk =γ1p1+γ1sI, Active Syk =p1γ1+p1sI, (1) Where sI denotes Syk copy numbers, p1 the rate that Syk is phosphorylated and γ1 the rate at which it is dephosphorylated. Local sensitivity analysis (see Methods) reveals that the model' s steady state is insensitive to variation in the rate of ITAM phosphorylation and the rate that Syk binds to the phosphorylated ITAM. While Model A is able to describe the equilibrium to which Syk activity settles, it is unable to capture its early transient peak in phosphorylation. This early peak in Syk phosphorylation is reminiscent of negative regulatory behaviour affecting Syk activity and points to the constitutively activated phosphatase being insufficient to explain this behaviour. Therefore the model was adapted, replacing the' housekeeping' phosphatase with a more complex pathway that involves a specific phosphatase, TULA-2. While it is known that TULA-2 can interact with Syk, in a ubiquitin-dependent manner [51], the purpose of the interaction had not previously been established. This interaction becomes the template for model adaption. The mathematical model was extended to incorporate a negative feedback, which leads from Syk phosphorylation to its own regulation, with the intention of explaining the early peak in Syk activity seen in experimental data. This newly introduced regulatory pathway centres on the proteins c-Cbl and TULA-2. c-Cbl is a member of the Cbl (Casitas B-lineage lymphoma) family of adaptor proteins and is known to be found in platelets [50,52]. It has been shown to play a role in regulating signals by similar receptors in other cell types [53,54] and it has been implicated in GPVI regulation in platelets [52]. c-Cbl is a ubiquitin ligase that is able to associate with phosphorylated tyrosine kinases, such as Syk, whereby it can promote their ubiquitination [55]. Syk is known to be ubiquitinated rapidly upon activation by CRP and collagen [50]. c-Cbl has no phosphatase domain but has been suggested as a necessary intermediate scaffold between Syk and its phosphatase [50]. The TULA (T-cell Ubiquitin Ligand) family of proteins are histidine tyrosine phosphatases that have been shown to be important negative regulators in immune cells [56]. TULA proteins are able to bind ubiquitinylated proteins and they play a role in regulating T-cells, decreasing Syk related protein (ZAP-70) phosphorylation [57,58]. TULA-2 exists in platelets [41,51] where it has been shown to associate with, and subsequently dephosphorylate, Syk [56]. Model B incorporates an additional Syk phosphorylation site (Y323) that is a known binding site for c-Cbl [59]. In the model, binding of c-Cbl to Y323 leads to ubiquitination of Syk, which allows TULA-2 to associate with Syk and dephosphorylate its activatory phosphorylation site (Y525), returning the receptor complex to an inactive state. All reactions are assumed to be reversible. A schematic of these newly introduced reactions is given in Fig 2, a network diagram in Fig 3 and the corresponding equations and parameters obtained from the parameter fitting process are given in S1 Text. Experimental data specifically describing the newly introduced phosphorylation site was collected and Model B was fitted to both sets of experimental data simultaneously (Fig 5a and 5b). Model profiles describe accurately the steady state representing Syk phosphorylation on Y525 and the newly introduced regulatory site (Y323). They failed, however, to capture the early transient peak seen in experimental observations of both phosphorylation sites. If Model B is fitted to data restricted to one phosphorylation site (Y525) then while model solutions capture the full dynamics seen in Syk Y525 phosphorylation predictions for phosphorylation on Y323 are inaccurate, settling to a level that is three fold higher than that seen experimentally (Fig 5c and 5d). The model steady state representing Syk phosphorylation on Y525 is given by Active Syk =sI−γ1p1G11r, (2) where G11r denotes the phosphatase TULA-2 bound to Syk, and hence able to dephosphorylate it. Local sensitivity analysis (see Methods) revealed that the time to reach equilibrium was still predominantly influenced by the rate that the ITAM is phosphorylated; the level of steady state (for the activatory and regulatory phosphorylation site) is most strongly influenced by the strength of the parameters that comprise the feedback loop. The level of Syk Y525 at steady state is also sensitive to the rate of Syk phosphorylation on Y525 and its reversal. This is reflected in the change in parameter values when fitted to both sets of data; the rate of all reactions within the feedback pathway and the rate that Y525 are increased to allow the model to fit the data. Model B accurately describes the steady states seen in the experimental observations but, though showing promise, fails to describe fully their kinetics. This led us to investigate biologically plausible refinements that would enable the model to recapture the early transitory behaviour we had seen when the model was fitted to data solely describing phosphorylation in Syk Y525. In formulating Model C two biologically plausible modifications to our model were explored. The first modification (H1) allows the TULA-2, when bound to the receptor complex, to dephosphorylate other Syk molecules. The second modification (H2) assumes that the increase in Syk activity (phosphorylation on Y525) increases the rate at which Y323 is phosphorylated—increases in Syk activity following Y525 phosphorylation have previously been reported [47]. The simultaneous inclusion of both modifications is denoted by Model C, H3. These reactions are depicted in Figs 2 and 3 and the corresponding equations are given in S1 Text. Model C, H1 and H2 profiles (Fig 6a and 6b broken and dotted lines respectively) fail to describe transitory behaviour in both sets of experimental data. Model C, H3 (Fig 6a and 6b, solid line) incorporates H1 and H2 simultaneously and is able to describe accurately both sets of experimental observations. These optimal fits were dependent on the models being simultaneously fitted to both sets of experimental data. Like Model B, if fitted to one set of observations (i. e. Y525), model predictions for Y323 were inaccurate (Fig 7c and 7d). The steady state that model profiles representing Syk phosphorylation on Y525 is determined by Active Syk =sI−γ1p1G11r− (TULA-2 bound) (GY525−G11r), (3) where GY525 denotes the portion of Syk molecules that are bound to the receptor complex and phosphorylated on Y525. Local sensitivity analysis (see Methods) revealed (see Fig G in S1 Text) that the steady state is influenced by parameters that comprise the whole regulatory pathway reflected in the above by the proportion of TULA-2 that is bound to the receptor complex. In agreement with Thomas et al [51] the above expression (3) and sensitivity analysis (see Fig J in S1 Text) confirm that a decrease in the levels of TULA-2 lead to hyperphosphorylation of Syk. The time to reach the above steady state rate is shown to be predominantly influenced by the rate of ITAM phosphorylation while the time to reach the steady state of the regulatory phosphorylation site (Y323) is influenced by the rate of ITAM phosphorylation and by protein copy numbers of GPVI and Syk. Of the parameters that were held fixed during the parameter fitting process (platelet volume, extracellular volume and rates of ligand binding and dissociation) only the platelet volume influences Syk activity. In summary, the ability of the model to describe the full range of dynamics seen in the experimental observations requires that Syk activity increases the rate of phosphorylation on Y323 and that TULA-2 is able to dephosphorylate not only the receptor to which it is bound but also any nearby receptor complex. To understand if the extension of each model is justified through a better fit to data we use the Akaike information criterion (AIC) and a modified form that takes into account small experimental sample sizes (AICc) —a description of the metric and how we apply it is given in Methods and the results are shown in Table 3. The model with the lowest AICc provides the best balance between its replication of data and the added complexity introduced to achieve it. Model C (H3) has the lowest AICc; it provides the most faithful replication of both sets of data (a total of forty-four experimental observations) and is not overly compromised by its additional complexity when compared to Model B. If data are restricted to that describing one phosphorylation site (Y525, twenty-two time-points) then Model A has the lowest AICc; Model B and C are too complicated (too many parameters) to be inferred from such a set of observations. Overall, these results confirm Model C, H3 as our preferred model and this was used for all further simulations.
A major success over the last 20 years has been the identification of the key components of the major signal transduction pathways in cells [2–7]. What is less well understood is how these components come together in space and time in a regulated manor to bring about appropriate cellular responses. Whilst the activation of protein kinases has been studied in some detail, much less work has been done on the regulation of phosphatases in part due to the lack of specific inhibitors. A search of the literature shows that there are 20 times the number of papers on protein kinases in platelets as there are on protein phosphatases. However as noted earlier most cells, including platelets, contain many different types of protein-tyrosine phosphatases (there are roughly the same number of different protein-tyrosine phosphatases as there are protein tyrosine kinases in platelets [40]), that show distinct subcellular compartmentalisation and regulation. Some appear to be simple abundant cytosolic phosphatases without obvious regulatory domains which may be constitutively active, some are transmembrane receptor type phosphatases such as CD148, whilst others are highly regulated containing domains that allow them to be temporally recruited to specific signalling complexes either via SH2 domain phosphotyrosine interactions (eg. , SHP1 and SHP2) or through ubiquitin-dependent binding (TULA-1 and TULA-2) [40–45]. Thus the challenge now is to determine which of these phosphatases are required and sufficient to explain the regulation of the key determining steps in individual signalling cascades. The objective was to develop a mathematical model of GPVI proximal signalling that allows the comparison of current biological knowledge and hypotheses to high-density temporal experimental data explicitly collected to inform the model. The experimental observations describe in particular phosphorylation of Syk, a key protein that is recruited to the GPVI receptor and which initiates downstream signalling. It was anticipated that this approach would lead to increased knowledge surrounding the interaction of the components of this key step in platelet activation and, importantly, how this is controlled. The model was developed incrementally. Starting with a simple model (A) that captured current biological knowledge of the components and interactions that occur just downstream of the GPVI receptor more complicated interactions were incorporated until a good description of experimental data was achieved. While the components that comprise the forward steps of platelet signalling have been well studied, regulatory processes have not. The experimental data pointed to the pathway being regulated and, in consequence, a simple housekeeping phosphatase that limits the ability of Syk (and therefore the receptor) to signal downstream was incorporated into Model A. The model was unable to explain the full dynamics seen in experimental observations and led us to conclude that the inclusion of a simple phosphatase, despite their abundance in platelets [40], is insufficient to explain data. A potential regulatory pathway was selected, supported by current literature, which focuses on the Syk protein and its ability to participate in its own regulation through proteins c-Cbl and TULA-2 [49–55]. Our initial attempts to explain the experimental data while promising did not capture the full range of dynamics displayed in observations of the central protein Syk on a phosphorylation site that we equate to activity and one that initiates the ability of Syk to regulate itself [46–47]. With further model refinements (Model C) we found that a model with c-Cbl/TULA-2 incorporated could provide not only a good quantitative fit to experimental data but could capture its full dynamics. As demonstrated in sensitivity analysis the inclusion of the negative feedback loop protects Syk activity from variation in the number of Syk molecules. Interestingly, if the experimental data is limited to that describing a single Syk phosphorylation site (Y525) then the inclusion of the more complex regulatory mechanism (centring on c-Cbl and TULA-2) cannot be justified. It is only when data describing a second phosphorylation site was collected that Model B and Model C could be successfully differentiated. This emphasises the need for high density data to allow the development of such complex models that incorporate more subtle processes (and many parameters). The model predicts that, on average, very few molecules of c-Cbl and TULA-2 need to be bound to the receptor complex at any point in time. This in turn suggests a significant stochastic component to the behaviour. This, and the prediction that regulation of activated Syk molecules relies on the ability of the bound phosphatase, TULA-2, being able to dephosphorylate many activated Syk molecules, point to receptors being clustered at the plasma membrane. This fits well with the localisation of GPVI to specific membrane domains and the known ability of the multivalent ligand, collagen, to drive localised clustering of GPVI at the platelet surface [65] but also highlights the importance of further spatial mathematical investigations. To be of wider use it is essential that the model can accurately predict experimental signalling outputs under varying conditions. The model was tailored to an environment where the ligand is abundant but is able to predict Syk activity in response to decreasing levels of ligand availability. A decrease in the availability of the ligand in the model results in a reduction of the early transient peak and a delay to reach the steady state which shows good agreement with independently generated experimental data. Platelet signalling is likely to demonstrate a wide degree of variation across a healthy population but little is known about the extent and at what point this leads to potential pathologies. This normal variation is seen in the levels of the individual proteins involved in our model, both in our data and the recent quantitative human platelet proteome [40] and also in the signalling responses of individual donors studied experimentally (Fig 7). In the first instance the model was fitted to data based on triplicate samples from one healthy donor. By comparison to data describing Syk phosphorylation, on both the activatory (Y525) and regulatory (Y323) sites, in nine additional donors we found the model was able to replicate phosphorylation profiles observed in a wider healthy population. The model was also utilised to investigate signalling responses in a population that has reduced levels of GPVI receptors [1,2] with normal population variability of expression levels of GPVI and other signalling components involved (Fig 10). This suggested that the time to peak of Syk activation in normal individuals seems to be tightly controlled whilst the actual peak maxima of Syk Y525 phosphorylation is more variable, as was also seen experimentally in our normal panel of nine healthy donors who all had very similar times to peak but quite variable absolute levels of Syk Y525 phosphorylation (Fig 7). The simulated population with reduced GPVI levels however had a longer time to reach this peak and the timing seems to be less tightly controlled within the simulated population. This may suggest that temporal regulation of this signalling pathway may be more important to a functional response than the absolute magnitude of the response. In summary, this study aims to bridge the gap between biological knowledge and hypotheses and temporal experimental data in aiding understanding of how the components of a key initiating step in platelet activation interact to control precisely a platelet' s response to its environment. Platelet activation has traditionally been thought of as a series of forward steps resulting in the platelet switching on. While the fine details of even this initiating step will continue to be elucidated, the principal of negative pathways playing a central role is clear and highlights the complexity of these signalling pathways. Indeed there is evidence of a number of inhibitory receptors and signalling mechnisms including those stimulated by nitric oxide and prostacyclin that serve to prevent un-required activation [69]. Some mechanisms such as those stimulated by the adhesion receptor PECAM-1 are known to modulate signalling downstream of the processes explored in this study, although we cannot rule the potential impact of these mechanisms on the initiation of GPVI signalling. Given the outcomes of this study, this will be the focus for future work. While our approach here has focused on elucidating a key regulatory step in a platelet signalling pathway we believe that this approach of integrating biological knowledge with time-rich kinetic data sets and mathematical models is well suited to the analysis of the regulation of any signalling cascade for which suitable reagents are available.
Here we briefly describe recently established methods utilised for the measurement and quantification of protein copy numbers and their post translational modifications. Please see [70] for further details. The quantification method is based on the relationship between a known amount of a specific protein to a quantitative signal emitted by fluorescently conjugated antibodies. Individual phosphorylated proteins were first isolated by immunopreciptation (IP) using suitable site-specific antibodies. Identical serial dilutions of the IP were loaded onto two SDS-PAGE gels. One gel was also loaded with a serial dilution of a known concentration of a corresponding non-phosphorylated recombinant protein. The resulting immunoblot was probed with an antibody recognising the total protein (i. e. non phosphorylated). The second gel loaded with the IP serial dilutions was blotted then probed with the same phospho-specific antibody used for the IP. Both gels also contained serial dilutions of known amount of IgG necessary to calibrate with subsequent experiments. The different immunoblots were then treated with appropriate fluorescent conjugated secondary antibodies and scanned with a fluorescence imaging system. The resulting values for the recombinant protein dilution were then used to construct a standard curve and determine the concentration of molecules immunoprecipitated. Combined with the results from the second immunoblot allows the direct relationship between the amount of protein actually present on the immunoblot and fluorescent levels obtained when probing with a phosphospecific antibody to be determined. Further experimental samples were loaded on SDS-PAGE gels together with identical IgG serial dilutions for normalisation with the reference datasets and quantification of the amount of phosphorylated proteins. The amount of total protein was established by comparing the normalised fluorescence intensity of the samples and the recombinant protein standard curve when probed with the antibody recognising the total protein. Blood samples were obtained from healthy volunteers that had given consent, using procedures approved by University of Reading Research Ethics Committee. Washed platelets were prepared by differential centrifugation, as described previously [70], and resuspended in Tyrodes buffer containing 0. 4U/ml Apyrase, 1mM EGTA and 100 μM Indomethacin; to suppress secondary signalling and secretion. Samples were stimulated with CRP-XL (provided by Prof Richard Farndale, University of Cambridge, UK) at a final concentration of 10 μg/ml, then lysed, denatured and loaded onto 10% SDS-PAGE gels. CRP-XL, unlike collagen, is a GPVI selective agonist. Reference datasets were constructed using platelets treated with pervanadate [71] (10 μg/ml) CRP for 1 min. Immunoprecipitation was carried out using PureProteome Protein A magnetic beads (Millipore). Immunoprecipitates and experimental samples were loaded alongside the corresponding recombinant protein (Syk, GST-tagged, Abnova; c-CBL proprietary tag, Abcam) and a series of dilutions of IgG (Rabbit and Murine Isotype controls, US Biological). Westernblot transfer to Immobilon-FL membrane (Millipore) was performed using a semi-dry blotter (Bio-Rad). The membranes were blocked with 5% (w/v) BSA and probed as instructed by the manufacturer using the Millipore SNAP i. d. protein detection system and the appropriate antibodies (Anti-Syk N-19, Santa-Cruz Biotechnology Inc; anti-c-Cbl, BD Biosciences; anti-Syk phospho Y525+Y526, anti-c-Syk Y323). Immunoblots were treated with either a fluorescent Cy5 dye-labeled goat anti-Rabbit or a Alexa-Fluor 647 dye-labeled donkey anti-mouse (Life Technologies). Digital scan values of the fluorescence emission were obtained using a Typhoon Trio scanner (GE Healthcare). Quantification was performed using the ImageQuant TL software and data analysed using R [72]. We estimated all unknown parameter values for our models by utilising a constrained local optimization routine (MATLAB' s fmincon) that varies all unknown parameters to minimize the differences between the model and experimental data via a cost function Sum Squared due to Error (SSE) = (yi (θ) −Datai) 2 (4) where yi (θ) is the model’s prediction for the relevant model variable (which depends on the parameters θ) and Datai represents the experimental observations at time points i. Given two sets of parameters, the one with the smaller cost function is the one that provides the better description of experimental observables. The values that parameters can take are restricted within limits that are guided by biological knowledge and literature (S1 Text includes a summary of the limits used and how they were set). The optimization algorithm takes steps that successively decrease the value of the cost function, beginning from an initial guess of the parameter values. To avoid parameter sets being selected that fit local, rather than global, minima we utilize a multi-start approach where repeated optimization runs are carried out with multiple initial guesses. We sampled initial guesses from the complete allowable parameter range utilising Latin Hypercube Sampling (LHS). The range that parameter values can take is discretised into N bins. Where parameter limits extend several orders of magnitude, a logscale is used, which leads to a better representation of parameter space. Each interval in the parameter range is sampled once (without replacement), so that the entire range for each parameter is explored. The resulting N samples are used as initial guesses. Parameter estimation is performed utilising samples N = 1000. The five parameter sets with the lowest SSE N = 1000 are reported in each case to provide confidence that the best fit is converging, for the final model (C, H3) samples of N = 10000 are reported, and needed to converge. We performed a local sensitivity analysis, where each parameter was varied by fifty and ninety percent above and below their initially estimated value, and computed the normalised local sensitivity of either the final steady state or the time to reach this steady state according to Sensitivity Score = (Oa−Oi) Oa (5) where Oi and Oa represent the model output (e. g. the steady state and time to reach the peak in Syk activity) in respect of the initial parameter set (obtained from parameter fitting) and the adapted parameter respectively. To aid model comparisons we used Akaike’s Information Criterion (AIC) [73]. This is an information-theoretic criterion for model comparison, which incorporates not only the cost function value from the parameter fitting process (SSE) but also a penalty based on the number of parameters (K) in the model. We used a modified criterion [68] that takes into account the experimental sample size (n) by increasing the relative penalty for model complexity with small datasets defined by AICc=AIC+2K (K+1) (n−K−1) (6) where AIC=n (ln (SSE/n) ) +2K. (7) AICc converges to AIC as n becomes large with respect to K. The value of AICc has no meaning in isolation, its relevance only becoming apparent when it is used to compare (and rank) models fitted to the same experimental data. In general a better fit to experimental data is achieved when the complexity of a model is increased, simply because the number of free parameters increases. The AICc (and indeed the AIC) aims to protect against this by a penalty based on the number of parameters in the model with the lowest AICc being the better model.
Platelets are blood cells that, upon injury, trigger the blood to clot. Following blood vessel damage platelets encounter the extracellular matrix protein collagen to which they respond. They become activated, aggregating to form a major component of blood clots. The platelet collagen receptor GPVI stimulates platelet activation through a complex signalling pathway, and while many of the molecules involved in the activation of this pathway have been identified, their specific roles in determining the rate and extent of the exceptionally rapid platelet response have not been determined. Furthermore, while signalling proteins responsible for forward reactions are known, reverse or negative feedback elements are not well understood. Platelets also trigger thrombosis in diseased arteries, causing heart attacks and strokes, and therefore platelets, and particularly the GPVI signalling pathway, are therapeutic targets. To begin to understand the components in the GPVI signalling pathways that may represent tractable therapeutic targets we have developed a mathematical model of the key initiating events that occur upon stimulation of GPVI. In so doing, we have established the importance of a specific phosphatase-controlled negative feedback in determining the rate of initiation of platelet activation.
Regulation of Early Steps of GPVI Signal Transduction by Phosphatases: A Systems Biology Approach