Regulation of Early Steps of GPVI Signal Transduction by Phosphatases: A Systems Biology Approach

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.


S.1 Model A
In each of our models a ligand, CRP, (denoted by l) binds to a platelet's GPVI receptor (g) to form a receptor complex (G). Binding of the ligand directly leads to phosphorylation of the receptor complex (G p ) allowing the cytosolic protein tyrosine kinase Syk (s) to dock (G b 0 ) and subsequently phosphorylate (G b 1 ). We equate Syk phosphorylation to its activity state and the subsequent ability of the receptor to signal downstream. In Model A regulation of Syk activity is incorporated in the form of a simple, constituitevely active, abundant phosphatase that can dephosphorylate Syk. A network diagram representing the events incorporated in Model A is shown in Fig. 3 (upper panel) while a summary of the variables and dimensional parameters is presented in Tables A and B, respectively. Utilising mass action kinetics these reactions transform to give the following system of nonlinear ordinary differential equations (ODEs) where A v represents Avogadro's number, V e the extracellular medium per cell available under experimental conditions and V p the cytosolic volume per platelet. These parameters ease conversion of reactions that occur between proteins that are embedded in the cell surface (receptors) and those that are in the extracellular space or the cell's cytosol. To facilitate comparison to experimental data the variables representing the receptor, in all its forms, and cytosolic proteins are expressed in molecules while the ligand is expressed as a concentration.
The initial conditions prior to ligand addition are set so that i.e. the signalling cascade is inactive, all receptors are unbound and all ligands free.

S2
Conservation of the ligand, receptor and Syk is given by respectively.

S.1.1 Parametrisation
Here we describe the model's parameters and initial conditions that are either well defined in the literature or we have been able to determine experimentally (initial conditions and parameter values are summarised in Tables A and B).
The rates that the ligand, CRP, binds to and dissociates from, a GPVI receptor are given in [10], allowing us to set k 1 = 8 m 3 moles −1 s −1 and k −1 = 3.02 × 10 −2 s −1 . We obtained the initial conditions for the number of GPVI receptors and Syk molecules per platelet experimentally so that g I = 5000 molecules and s I = 2763 molecules. Determining the initial conditions for the ligand was complicated by CRP being of an indeterminate length, weight and number of binding sites for a GPVI receptor. This made it difficult to convert from the experimental units of mass concentration (measured in units of µg/mL) to the model's molar concentration (measured in units of moles/m 3 ). In all of our experiments we utilised a high saturating level of the ligand (10µg/mL) and therefore for our model estimated l I = 3 × 10 −2 moles m −3 , a level that is non-rate limiting. For every blood sample taken we measured the mean platelet volume allowing us to set V p = 7.4 × 10 −18 m 3 and in all experiments platelets were suspended in media at a concentration of 3×10 8 cells/mL [2] giving us an extracellular volume per platelet (Ve) of 3.33×10 −9 m 3 .

S.1.2 Parameter fitting
We obtained values for the remaining parameters (k 2 , k 3 , p 1 , γ 1 ) via a process of parameter fitting that is described in the main paper (see Methods). The values that these parameters could take are constrained to biologically realistic limits. These are based on estimates in the literature, for similar processes either involving different proteins and/or in different cell lines.
Estimates describing phosphorylation are reported in the literature to be in the range 0.1−100 s −1 [3,5,11]. Relaxing this range by one order of magnitude we allowed parameter values for phosphorylation to fall between a range of 10 −2 −10 3 s −1 . Binding rates for proteins are reported to be in the range of 1×10 2 −1×10 5 m 3 moles −1 s −1 [5][6][7][8]11,12] allowing us to set a range of 10 1 − 10 6 m 3 moles −1 s −1 for all protein binding. The reverse processes of dissociation and dephosphorylation are reported to be in the range 0.03 − 20 s −1 [5,7,8,11,12] allowing us to set a range of 1 × 10 −3 − 10 3 s −1 .  Simulations utilise five optimal sets of parameter values obtained from the parameter fitting process (see Table C) and demonstrate that utilising these sets results in simulations that are indistinguishable, they all describe the data equally well. Model A is unable to describe the early peak in phosphorylation displayed in the experimental observations. experimental data have on the model profiles representing Syk activity (both its steady state and the time to reach that steady state). These scores demonstrate that Model A's steady state is determined by the rate of Syk phosphorylation (p 1 ) and its reverse process (γ 1 ), the ratio of which is constant for the five 'best' fits (p 1 /γ 1 = 0.145, see Table C). The time to reach the steady state is mainly sensitive to variation in the rate that the receptor complex is phosphorylated (k 2 ). Utilising equations (S1) and (S2) we can see that at steady state, in the regime where l I >> g I > s I , all Syk molecules are bound to the receptor (s = 0) and the number of inactive and active receptors are given by

S.2 Model B
Model B extends Model A to incorporate the hypothesis that Syk activity is regulated via a negative feedback pathway initiated through a newly introduced Syk phosphorylation site (Y323). Phosphorylation on Y323 is thought to allow the cytosolic protein c-Cbl to bind to Syk where it facilitates ubiquitination and subsequent binding of the cytosolic protein TULA-2 that is a phosphatase thought able to desphosphorylate Syk Y525. A network diagram representing the events incorporated in Model B is shown in Fig. 3 while a summary of the variables and the newly introduced dimensional parameters is presented in Tables D and E, respectively. The remaining parameters are the same as those for Model A (Table B). Incorporating these reactions into Model A leaves equations (S1a)-(S1c), and (S1g) unchanged.
Model A's variables G b 0 and G b 1 are modified to the form G k i,j which represents Syk (bound to the receptor complex) in its various states. Subscripts i and j indicate phosphorylation on sites Y323 and Y525 respectively and superscript k denotes the sequential events of Syk binding to the receptor complex  Figure C. Sensitivity of Model A to variation in the model's initial conditions. Sensitivity scores are shown for the effect of parameter variation on variable G b 1 (representing Syk activity). Parameters here are varied ± 50% and 90%.

S5
(b), c-Cbl binding (c), Syk's ubiquitination (u) and TULA-2 binding to Syk (r). In Model B the equations (S1d) and (S1e) are therefore replaced by the following ten differential equations where terms unchanged from those in Model A are shown in black and those newly introduced in blue We introduce two new equations to describte cytosolic c-Cbl and TULA-2 The initial conditions are set so that where we obtain the number of c-Cbl molecules per platelet experimentally and the number of TULA-2 molecules per platelet is given by [2], so that c I = 2581 molecules and s I = 7800 molecules. The newly introduced variables and parameters are described in Tables D and E respectively. In Model A the variable G b 1 was compared to experimental data of Syk Y525 phosphorylation.
and data describing phosphorylation of Syk on Y323 is compared to G Y 323 where  Table G and demonstrate that utilising these sets results in simulations that are indistinguishable.

S.3 Model B solutions
Model B has eight newly introduced parameters (listed in Table E). Parameters p 4 and p −4 describe ubiquitination and its reverse process. [3] provide an estimate for ubiquitination (1.67 s −1 ) that is of a similar order to estimates for phosphorylation and this leads us to set The remaining parameters (p 2 , p −2 , p 3 , p −3 , p 5 , p −5 ) describe phosphorylation, binding or their reverse processes and we therefore follow the arguments described in Section S.1.1 in setting the constraints that their values can take during the fitting process.
The number of parameters that require estimation has increased from four, in Model A, to twelve in Model B. As in Model A we fixed the parameters k 1 , k −1 , V p , V e to the values described in Table B and then infer the remaining parameters from data. The five 'best' parameter sets obtained from fitting model B to experimental data describing Syk 525 phosphorylation and then simultaneously fitting the model to data describing Syk phosphorylation on Y525 and Y323 are shown in Tables F and G respectively. Simulations utilising these values are indistinguishable from each other demonstrating that the fitting process has converged (see Fig. D for simulations that utilise the five optimal sets of parameters obtained from fitting Model B to observations of Y525 and Y323).
Model B fitted to data describing the activatory phosphorylation site (Y525) is able to describe experimental observations of Y525 but model predictions for the regulatory phosphorylation site (Y323) settle to a steady state that is approximately twice that described by data (see Fig. 4a and 6b). Model B fitted simultaneously to experimental data describing both phosphorylation sites (see Fig. 4c and 4d) is able to describe the steady states that experimental observations settle to but is unable to describe the early transient peaks displayed in the data.
Fig.E presents the sensitivity scores corresponding to the effect that variation in the parameters obtained from fitting the model to experimental data, describing phosphorylation on Y525 and Y323, have on the model profiles representing both phosphorylation sites (their steady states and the time to reach that steady state). These scores demonstrate that the time to reach the steady state in Syk Y525 (and indeed in Y323) is, like Model A, mainly sensitivie to variation in the rate that the receptor complex is phosphorylated (k 2 ). Utilising equations (S4) and (S13) we can see that at steady state, in the regime where l I >> g I > s I , all Syk molecules are bound to the receptor (s = 0) and the number of active receptors are given by where G b 1,1 (the proportion of active Syk with TULA-2 bound) reflects Syk activity participating in its own regulation. The steady states of Syk Y525 is now sensitive to variation in the rate of Syk phosphorylation and the parameters that comprise the regulatory pathway.

S.4 Model C
In formulating Model C we explored the following biological hypothesis: H1 : TULA-2 (once bound to the receptor complex) is able to dephosphorylate not only the Syk molecule to which it is associated but any nearby Syk molecule (at rate γ 2 ), H2 : Syk activity (phosphorylation on Y525) increases the rate of phosphorylation of Y323 (denoted by q 2 ), H3 : H1 and H2 incorporated simultaneously.
A network diagram representing Model C is shown in Fig. 3 (lower panel). The newly introduced parameters (q 2 , γ 2 ) are presented in Table H. All other parameters and all variables are the same as those for Model B (see Tables B, E and D). Incorporating the hypotheses into the equations for Model B leaves equations (S1a)-(S1c), (S1g), (S4e) and (S4j-S4l) unchanged. Incorporating H1, the ability of TULA-2 to dephosphorylate nearby receptors (rate γ 2 ), requires additional terms in equations (S4a-S4d, S4f-S4i) and incorporating H2 (the increase in the rate that Y323 is phosphorylated (by q 2 ) results in the modifications to equations (S4f) and (S4g). These equations (with the modifications for H1 and S9 H2 depicted in blue and green respectively) now read The initial conditions and conservation laws for Model C are as described for Model B (equations (S5)-(S13)).

S.5 Model C solutions
Model C has two newly introduced parameters (listed in Table H). The parameter γ 2 describes dephosphorylation of Y525 but, unlike the parameter (γ 1 ) introduced in Model B, is second order. As there are no estimates for the second order form (units molecules −1 s −1 ) in literature we set a wide range that the parameter is allowed to take in the parameter fitting process (1.0 × 10 −3 to 1 × 10 4 molecules s −1 ). The parameter q 2 is dimensionless, it increases the rate that of Syk phosphorylation p 2 and we therefore allow it to take the range 1.0 to 1 × 10 4 .
We incorporate each hypothesis and fit the model to experimental data. The model is fitted twice, initially to data describing phosphorylation on Y525 and then to data describing phosphorylation on Y525 and Y323. The five 'best' parameter sets obtained from fitting Model C (with H1, H2 and then H3 incorporated) are shown in Tables I and J. Simulations utilising the optimal parameter sets from fitting Model C simultaneously to data describing both phosphorylation sites are shown in Fig 5c and 5d. These demonstrate that Model C with hypothesis H3 is best able to accurately describe Syk phosphorylation on Y525 and Y323. Model C with H1 or H2 incorporated fail to capture the full dynamics displayed in the data. If experimental data is restricted to one phosphorylation site (Y 525) then Model C, like Model B, is able to describe Syk Y525 phosphorylation but predictions for Y323 phosphorylation are inaccurate.
When fitting Model C, H3 a sample size of N=10000 was used, this larger sample size being required to fit the more complicated model. Fig. F demonstrates that simulations utilising the five best parameter sets (Table J, Table J (H3, 1 − 5) and demonstrate that utilising these sets results in simulations that are indistinguishable.
sets that describe the data equally well. and Y323 (both their steady states and the time to reach that steady state). These scores demonstrate that the time to reach both steady states is predominantly influenced by the rate that the receptor complex is phosphorylated (k 2 ) with the time to reach the steady state of Y323 also being strongly influenced by the initial conditions of GPVI and Syk. Of the paremeters that were held fixed during the parameter fitting process the volume has the most influence on the model outputs. The steady state in Syk Y525 phosphorylation is influenced by all parameters that comprise the regulatory pathway and the levels of the regulatory proteins (c-Cbl and TULA-2). Utilising Model C equations (S15) and (S13) we can see that at steady state, in the regime where l I >> g I > s I , all Syk molecules are bound to the receptor (s = 0) and the number of active receptors are given by where (G b 1,0 + G b 1,1 ) represents TULA-2 bound to the receptor complex and (G b 0,1 + G b 1,1 + G c 1,1 + G u 1,1 ) represents the portion of Syk phosphorylated on Y525 but without TULA-2 bound.